линейная (log-log) модель с помощью 'lm': как получить прогнозную дисперсию суммы прогнозируемых значений

Question

Я подгоняю модель мощности к набору данных, применяя простую линейную модель с функцией R lm после преобразования log-log, как в примере ниже (вместо непосредственного подбора модели мощности, например, применяя функция nls). Я мог бы использовать функцию predict.lm, чтобы применить модель к новым данным и вычислить интервалы прогнозирования.

data(stackloss); dat <- stackloss[c(2, 4)]; colnames(dat) <- c("x","y")
dat.lm <- lm(log(y) ~ log(x), data = dat)

new <- data.frame(x = seq(0, 30, 1))
pred <- predict.lm(dat.lm, new, interval = "prediction", level = 0.95)
matplot(new$x, exp(pred), type = "l", col = 1, lty = c(1, 2, 2)); points(dat$x, dat$y)

Теперь мне нужно сложить n прогнозируемых значений (что довольно просто после применения функции exp), а также рассчитать агрегированную дисперсию и интервалы прогнозирования. Последнее было описано для простой линейной модели в следующих вопросах и ответах: линейная модель с `lm`: как получить прогнозную дисперсию суммы прогнозируемых значений . В этом интересном ответе для простой линейной модели были введены следующие функции lm_predict (что позволяет вычислить полную дисперсионно-ковариационную матрицу прогнозируемых значений) и agg_pred.

lm_predict <- function (lmObject, newdata, diag = TRUE) {
  ## input checking
  if (!inherits(lmObject, "lm")) stop("'lmObject' is not a valid 'lm' object!")
  ## extract "terms" object from the fitted model, but delete response variable
  tm <- delete.response(terms(lmObject))      
  ## linear predictor matrix
  Xp <- model.matrix(tm, newdata)
  ## predicted values by direct matrix-vector multiplication
  pred <- c(Xp %*% coef(lmObject))
  ## efficiently form the complete variance-covariance matrix
  QR <- lmObject$qr   ## qr object of fitted model
  piv <- QR$pivot     ## pivoting index
  r <- QR$rank        ## model rank / numeric rank
  if (is.unsorted(piv)) {
    ## pivoting has been done
    B <- forwardsolve(t(QR$qr), t(Xp[, piv]), r)
    } else {
    ## no pivoting is done
    B <- forwardsolve(t(QR$qr), t(Xp), r)
    }
  ## residual variance
  sig2 <- c(crossprod(residuals(lmObject))) / df.residual(lmObject)
  if (diag) {
    ## return point-wise prediction variance
    VCOV <- colSums(B ^ 2) * sig2
    } else {
    ## return full variance-covariance matrix of predicted values
    VCOV <- crossprod(B) * sig2
    }
  list(fit = pred, var.fit = VCOV, df = lmObject$df.residual, residual.var = sig2)
  }

agg_pred <- function (w, predObject, alpha = 0.95) {
  ## input checing
  if (length(w) != length(predObject$fit)) stop("'w' has wrong length!")
  if (!is.matrix(predObject$var.fit)) stop("'predObject' has no variance-covariance matrix!")
  ## mean of the aggregation
  agg_mean <- c(crossprod(predObject$fit, w))
  ## variance of the aggregation
  agg_variance <- c(crossprod(w, predObject$var.fit %*% w))
  ## adjusted variance-covariance matrix
  VCOV_adj <- with(predObject, var.fit + diag(residual.var, nrow(var.fit)))
  ## adjusted variance of the aggregation
  agg_variance_adj <- c(crossprod(w, VCOV_adj %*% w))
  ## t-distribution quantiles
  Qt <- c(-1, 1) * qt((1 - alpha) / 2, predObject$df, lower.tail = FALSE)
  ## names of CI and PI
  NAME <- c("lower", "upper")
  ## CI
  CI <- setNames(agg_mean + Qt * sqrt(agg_variance), NAME)
  ## PI
  PI <- setNames(agg_mean + Qt * sqrt(agg_variance_adj), NAME)
  ## return
  list(mean = agg_mean, var = agg_variance, CI = CI, PI = PI)
  }

Однако их нельзя применять напрямую для надлежащего агрегирования дисперсии в случае регрессии log-log. Возможно, мне следует преобразовать дисперсию в вывод lm_predict, но я не мог понять, как поступить. Заранее благодарю за любую помощь.