主页 / coding / 问题 / 79598915
Accepted
eee ppp
Asked: 2025-04-30 01:36:27 +0800 CST

R predFit()：两个变量中 nls 的预测区间问题

  • 772

我正在对一个包含两个解释变量的数据集进行非线性回归拟合Day，并使用Treatment建模响应。我想获取 的预测值以及新数据集的 95% 预测区间。但是，如果我重复使用同一个数据集，也会出现同样的问题。因此，我添加了示例，但没有设置- 参数。Amount{nls}Amountnewdata

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))

尽管我可以用以及 来预测Amount每个数据点的值，predict()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

当我尝试包含预测区间时，收到以下错误消息。

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

eval(form[[3]]) 中的错误：未找到对象“Day”此外：警告消息：在assign(xname, newdata[, xname])中：仅第一个元素用作变量名

但是，如果我在公式中仅使用Day（无）来降低模型的复杂度，则一切正常，并且会得到。我做错了什么？Treatment{nls}predFit()Day

1 个回答

  1. Best Answer

    罪魁祸首就assign(xname, newdata[, xname])在这里。

    当有两个或多个独立变量时，上述代码将变为assign(c('variable1','variable2'), newdata[, c('variable1','variable2')])。这没有任何意义，预期的方式是分别对它们进行赋值。

    我们需要做的就是用

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

    代码如下：

    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

相关问题