在 R 中使用 lapply() 时，predFit() 生成大小不合理的预测文件？

Question

遵循用户 tpetzdoldt 在对 () 的回答中给出的示例结构，

library(tidyverse)
 library(investr)
 data <- tibble(date = 1:7,
                cases = c(0, 0, 1, 4, 7, 8.5, 8.5))

    model <- nls(cases ~ SSlogis(log(date), Asym, xmid, scal), data= data )
    new.data <- data.frame(date=seq(1, 10, by = 0.1))
    interval <- as_tibble(predFit(model, newdata = new.data, interval = "confidence", level= 0.9)) %>% 
      mutate(date = new.data$date)

然后我尝试将这些相同的概念应用于我自己的数据（此处生成的可复制版本）：

#Trying to create a reproducible example:

string_temp <- c(5, 12, 43, 12, 0.5, 11, 16, 15, 10, 8)
string_resp <- c(22, 15, 106, 18, 9, 14, 32, 11, 1, 4)
string_id <- c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", 
               "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V")

temp <- rep(string_temp, 220)
resp <- rep(string_resp, 220)
id <- rep(string_id, 100)

data_model <- data.frame(temp, resp, id)

#Data for predictions:

predictions <- runif(122735)

predictions <- data.frame(predictions)
predictions <- predictions %>% rename(temp = predictions)

#Split by identity:


data_model_split <- data_model %>% split(data_model$id)


#model:


model <- lapply(data_model_split, function(d) nls(resp ~ a * exp(b * temp), 
                                                  start = list(a = 0.8, b = 0.1), 
                                                  data = d))

#results:

results <- lapply(1:2, function(i) { 
  predFit(model[[i]], newdata = predictions, interval = "confidence", level = 0.9)})

我收到以下错误：

 Error: cannot allocate vector of size 112.2 Gb

这些调整会生成那种大小的数据框，这似乎很奇怪。上面示例中生成的数据框只有 4 列宽。我正在为“模型”生成 122,000 行的 22 个模型提供数据，但我仍然感到震惊，并且确信它生成的假设的 4 列 x 3,000,000 数据帧不应该接近 1 Gb。在这种情况下，我对 lapply() 的应用有问题吗？由于数据集非常大，我为我个人示例中缺乏可重现性而道歉，但我希望问题可能出在我的代码中而不是我的数据集中。如果有帮助，我可以尝试为我的数据生成一个可重现的代理。

Answer 1

错误是由predFun代码中的以下行引起的：

v0 <- diag(f0 %*% tcrossprod(solve(crossprod(R1)), f0))

它试图构建一个 122735 x 122735 矩阵（以你的例子为例），并取它的对角线。在基数 R 中构建这样大小的矩阵可能需要很多 space。但是，注意前面的函数等价于：

library(magrittr)

v0 <- lapply(1:nrow(f0), function(rw){
  f0[rw,,drop=F] %*% tcrossprod(solve(crossprod(R1)), f0[rw,,drop=F])
}) %>% do.call(c,.)

即，如果我们可以按行求和，则立即 we don't need the whole matrix。

备注：

我确定有一种 easier/quicker 方法可以实现与原始代码尝试相同的效果。对于您正在寻找的任务，可能有一些替代库可以更有效地实现同样的目标，但这些不是本答案的重点。
如果您打算专门使用 predFun，则可以更正代码并覆盖原始功能。

对于覆盖，我不是专家，必须存在一种 cleaner/more 优雅的方法来做到这一点。但是，一个示例可能是通过提取原始代码并解决手头的问题：

'predFit.nls_custom' <- function (object, newdata, se.fit = FALSE, interval = c("none", 
                                                        "confidence", "prediction"), level = 0.95, adjust = c("none", 
                                                                                                              "Bonferroni", "Scheffe"), k, ...) 
{
  require(magrittr)
  adjust <- match.arg(adjust)
  compute.se.fit <- if (se.fit || (interval != "none")) 
    TRUE
  else FALSE
  if (object$call$algorithm == "plinear") {
    stop(paste("The Golub-Pereyra algorithm for partially linear least-squares \n               models is currently not supported."), 
         call. = FALSE)
  }
  newdata <- if (missing(newdata)) {
    eval(getCall(object)$data, envir = parent.frame())
  }
  else {
    as.data.frame(newdata)
  }
  if (is.null(newdata)) {
    stop("No data available for predictions.", call. = FALSE)
  }
  xname <- intersect(all.vars(formula(object)[[3]]), colnames(newdata))
  pred <- object$m$predict(newdata)
  if (compute.se.fit) {
    param.names <- names(coef(object))
    for (i in 1:length(param.names)) {
      assign(param.names[i], coef(object)[i])
    }
    assign(xname, newdata[, xname])
    form <- object$m$formula()
    rhs <- eval(form[[3]])
    if (is.null(attr(rhs, "gradient"))) {
      f0 <- attr(numericDeriv(form[[3]], param.names), 
                 "gradient")
    }
    else {
      f0 <- attr(rhs, "gradient")
    }
    R1 <- object$m$Rmat()

    # Applied fix below: 
    v0 <- lapply(1:nrow(f0), function(rw){
      f0[rw,,drop=F] %*% tcrossprod(solve(crossprod(R1)), f0[rw,,drop=F])
    }) %>% do.call(c,.) 
    # --- End of fix
    
    se_fit <- sqrt(Sigma(object)^2 * v0)
  }
  interval <- match.arg(interval)
  if (interval == "none") {
    res <- pred
  }
  else {
    crit <- if (adjust == "Bonferroni") {
      qt((level + 2 * k - 1)/(2 * k), df.residual(object))
    }
    else if (adjust == "Scheffe") {
      if (interval == "confidence") {
        p <- length(coef(object))
        sqrt(p * qf(level, p, df.residual(object)))
      }
      else {
        sqrt(k * qf(level, k, df.residual(object)))
      }
    }
    else {
      qt((level + 1)/2, df.residual(object))
    }
    if (interval == "confidence") {
      lwr <- pred - crit * se_fit
      upr <- pred + crit * se_fit
    }
    else {
      lwr <- pred - crit * sqrt(Sigma(object)^2 + se_fit^2)
      upr <- pred + crit * sqrt(Sigma(object)^2 + se_fit^2)
    }
    res <- cbind(fit = pred, lwr = lwr, upr = upr)
  }
  if (se.fit) {
    res <- list(fit = res, se.fit = se_fit, df = df.residual(object), 
                residual.scale = Sigma(object))
  }
  return(res)
}

下一步，一种方法是包含以下代码，用我们的自定义变体 predFit.nls_custom 覆盖 predFit.nls 方法。 (See here for other ways to override).

assignInNamespace("predFit.nls",predFit.nls_custom,ns="investr")
Sigma <- investr:::Sigma
Sigma.nls <- investr:::Sigma.nls

并重新运行原代码：

results <- lapply(1:2, function(i) { 
  predFit(model[[i]], newdata = predictions, interval = "confidence", level = 0.9)}
  )

现在应该可以正常工作了。如果没有，则应用覆盖可能会出现问题。

在 R 中使用 lapply() 时，predFit() 生成大小不合理的预测文件？

predFit() generates an unreasonably sized prediction file when using lapply() in R?

r

nls

predict

dplyr