Início / coding / Perguntas / 79598915
Accepted
eee ppp
Asked: 2025-04-30 01:36:27 +0800 CST

R predFit(): Problema com intervalo de predição para nls em duas variáveis

  • 772

Estou ajustando uma regressão não linear em um conjunto de dados com duas variáveis ​​explicativas Daye Treatmentmodelando uma resposta Amountusando {nls}. Quero obter os valores previstos para , Amountjuntamente com os intervalos de previsão de 95%, para um novo conjunto de dados. No entanto, o problema também ocorre se eu estiver apenas reutilizando o mesmo conjunto de dados. Portanto, estou incluindo o exemplo sem definir o newdataargumento -.

DS <- data.frame(Day=rep(1:5,t=2),
                 Amount=c(65,17,11,3.5,1.2,85,23,15,5,1.7),
                 Treatment=rep(c(0,1),e=5))

Model <- nls(Amount ~ (a+b*Treatment) * exp((c+d*Treatment)*Day) + e,
  data=DS, start = list(a=10,b=10,c=-0.1,d=-0.1,e=0.1))

Embora eu possa prever os Amountvalores para cada ponto de dados com predict()e também com investr::predFit(),

> predict(Model)
> [1] 64.814965 18.614433  7.150980  4.306624  3.600871 84.626103 25.804353  9.562934  5.078476  3.840261
> predFit(Model)
> [1] 64.814965 18.614433  7.150980  4.306624  3.600871 84.626103 25.804353  9.562934  5.078476  3.840261

Recebo a seguinte mensagem de erro quando tento incluir intervalos de previsão.

predFit(Model,interval="prediction",level=0.95)

Erro em eval(form[[3]]): objeto 'Dia' não encontrado. Além disso: Mensagem de aviso: Em assign(xname, newdata[, xname]): apenas o primeiro elemento é usado como nome da variável.

No entanto, se eu reduzir a complexidade do modelo usando apenas Day(no Treatment) na fórmula em {nls}, tudo funciona bem e predFit()encontra Day. O que estou fazendo errado?

1 respostas

  1. Best Answer

    O culpado está assign(xname, newdata[, xname])aqui .

    Quando há duas ou mais variáveis ​​independentes, o código acima se torna assign(c('variable1','variable2'), newdata[, c('variable1','variable2')]). Isso não faz sentido e a maneira pretendida é atribuí-las separadamente.

    Tudo o que precisamos fazer é substituir o culpado por

    for (i in 1:length(xname)) { 
      assign(xname[i], newdata[, xname[i]])  
    }

    O código está aqui:

    predFit.nls <- function(object, newdata, se.fit = FALSE,
                        interval = c("none", "confidence", "prediction"), 
                        level = 0.95, 
                        adjust = c("none", "Bonferroni", "Scheffe"), k, 
                        ...) {
  
  # Match arguments
  interval <- match.arg(interval)
  adjust <- match.arg(adjust)
  
  # Make sure se.fit is set to TRUE if intervals are requested
  compute.se.fit <- if (se.fit || (interval != "none")) {
    TRUE
  } else {
    FALSE
  }
  
  # No support for the Golub-Pereyra algorithm for partially linear 
  # least-squares models
  if (interval != "none") {
    if (!is.null(object$call$algorithm) && object$call$algorithm == "plinear") {
      stop(paste0("The Golub-Pereyra algorithm for partially linear least-", 
                  "squares models is currently not supported."), call. = FALSE)
    }
  }
  
  # Prediction data
  newdata <- if (missing(newdata)) {
    eval(stats::getCall(object)$data, envir = parent.frame()) 
  } else {
    as.data.frame(newdata) 
  }
  if (is.null(newdata)) {
    stop("No data available for predictions.", call. = FALSE)
  }
  
  # Name of independent variable
  xname <- intersect(all.vars(stats::formula(object)[[3]]), colnames(newdata)) 
  
  # Predicted values
  pred <- object$m$predict(newdata)
  
  # Compute standard error
  if (compute.se.fit) {
    
    # Assign values to parameter names in current environment
    param.names <- names(stats::coef(object))  
    for (i in 1:length(param.names)) { 
      assign(param.names[i], stats::coef(object)[i])  
    }
    
    # Assign values to independent variable name
    
    
    for (i in 1:length(xname)) { 
      assign(xname[i], newdata[, xname[i]])  
    }
    
    # Calculate gradient (numerically)
    form <- object$m$formula()
    rhs <- eval(form[[3]])
    if (is.null(attr(rhs, "gradient"))) {
      f0 <- attr(stats::numericDeriv(form[[3]], param.names), "gradient")
    } else {  # self start models should have gradient attribute
      f0 <- attr(rhs, "gradient")
    }
    
    # Calculate standard error
    R1 <- object$m$Rmat()
    # v0 <- diag(f0 %*% solve(t(R1) %*% R1) %*% t(f0))
    v0 <- diag(f0 %*% tcrossprod(solve(crossprod(R1)), f0))  # slightly faster
    se_fit <- sqrt(stats::sigma(object)^2 * v0)
    
  }
  
  # Compute results
  if (interval == "none") {
    
    # Vector of fitted/predicted values
    res <- pred    
    
  } else { 
    
    # Adjustment for simultaneous inference
    crit <- if (adjust == "Bonferroni") {  # Bonferroni adjustment 
      
      stats::qt((level + 2*k - 1) / (2*k), stats::df.residual(object))
      
    } else if (adjust == "Scheffe") {  # Scheffe adjustment
      
      if (interval == "confidence") {
        p <- length(stats::coef(object))  # number of regression parameters
        # sqrt(p * stats::qf((level + 1) / 2, p, stats::df.residual(object))) 
        sqrt(p * stats::qf(level, p, stats::df.residual(object))) 
      } else {
        # sqrt(k * stats::qf((level + 1) / 2, k, stats::df.residual(object))) 
        sqrt(k * stats::qf(level, k, stats::df.residual(object))) 
      }     
      
    } else {  # no adjustment   
      
      stats::qt((level + 1) / 2, stats::df.residual(object))   
      
    }
    
    # Interval calculations
    if (interval == "confidence") {  # confidence limits for mean response
      lwr <- pred - crit * se_fit  # lower limits
      upr <- pred + crit * se_fit  # upper limits
    } else {  # prediction limits for individual response
      lwr <- pred - crit * sqrt(stats::sigma(object)^2 + se_fit^2)  # lower limits
      upr <- pred + crit * sqrt(stats::sigma(object)^2 + se_fit^2)  # upper limits
    }
    
    # Store results in a matrix
    res <- cbind("fit" = pred, "lwr" = lwr, "upr" = upr)
    
  }
  
  # If standard errors of fitted values are requested, convert results to a list
  # and store additional information
  if (se.fit) {
    res <- list("fit" = if (interval != "none") res else pred, #res,
                "se.fit" = se_fit,
                "df" = stats::df.residual(object),
                "residual.scale" = stats::sigma(object))
  }
  
  # Return results
  return(res)
  
}
    • 2

relate perguntas