将 lm() 和 predict() 应用于数据框中的多个列

Question

我在下面有一个示例数据集。

train<-data.frame(x1 = c(4,5,6,4,3,5), x2 = c(4,2,4,0,5,4), x3 = c(1,1,1,0,0,1),
                  x4 = c(1,0,1,1,0,0), x5 = c(0,0,0,1,1,1))

假设我想根据列x1 和x2 为列x3、x4、x5 创建单独的模型。例如

lm1 <- lm(x3 ~ x1 + x2)
lm2 <- lm(x4 ~ x1 + x2)
lm3 <- lm(x5 ~ x1 + x2) 

然后我想将这些模型应用到使用 predict 的测试集，然后创建一个矩阵，将每个模型结果作为一列。

test <- data.frame(x1 = c(4,3,2,1,5,6), x2 = c(4,2,1,6,8,5))
p1 <- predict(lm1, newdata = test)
p2 <- predict(lm2, newdata = test)
p3 <- predict(lm3, newdata = test)
final <- cbind(p1, p2, p3)

这是一个简化版，你可以一步一步做，实际数据太大了。有没有办法创建一个函数或使用 for 语句将其组合成一个或两个步骤？

Answer 1

我倾向于将您的问题作为与Fitting a linear model with multiple LHS 重复的问题来结束，但遗憾的是，预测问题在那里没有得到解决。另一方面，Prediction of 'mlm' linear model object from lm() 谈论预测，但与您的情况有点远，因为您使用的是公式接口而不是矩阵接口。

我没能在"mlm" tag 中找到完美的重复目标。所以我认为为这个标签提供另一个答案是个好主意。正如我在链接问题中所说，predict.mlm 不支持se.fit，目前，这也是“mlm”标签中缺少的问题。所以我会借此机会填补这个空白。

---

这是一个获得预测标准误差的函数：

f <- function (mlmObject, newdata) {
  ## model formula
  form <- formula(mlmObject)
  ## drop response (LHS)
  form[[2]] <- NULL
  ## prediction matrix
  X <- model.matrix(form, newdata)
  Q <- forwardsolve(t(qr.R(mlmObject$qr)), t(X))
  ## unscaled prediction standard error
  unscaled.se <- sqrt(colSums(Q ^ 2))
  ## residual standard error
  sigma <- sqrt(colSums(residuals(mlmObject) ^ 2) / mlmObject$df.residual)
  ## scaled prediction standard error
  tcrossprod(unscaled.se, sigma)
  }

对于您给定的示例，您可以这样做

## fit an `mlm`
fit <- lm(cbind(x3, x4, x5) ~ x1 + x2, data = train)

## prediction (mean only)
pred <- predict(fit, newdata = test)

#            x3          x4         x5
#1  0.555956679  0.38628159 0.60649819
#2  0.003610108  0.47653430 0.95848375
#3 -0.458483755  0.48014440 1.27256318
#4 -0.379061372 -0.03610108 1.35920578
#5  1.288808664  0.12274368 0.17870036
#6  1.389891697  0.46570397 0.01624549

## prediction error
pred.se <- f(fit, newdata = test)

#          [,1]      [,2]      [,3]
#[1,] 0.1974039 0.3321300 0.2976205
#[2,] 0.3254108 0.5475000 0.4906129
#[3,] 0.5071956 0.8533510 0.7646849
#[4,] 0.6583707 1.1077014 0.9926075
#[5,] 0.5049637 0.8495959 0.7613200
#[6,] 0.3552794 0.5977537 0.5356451

我们可以验证f是否正确：

## `lm1`, `lm2` and `lm3` are defined in your question
predict(lm1, test, se.fit = TRUE)$se.fit
#        1         2         3         4         5         6 
#0.1974039 0.3254108 0.5071956 0.6583707 0.5049637 0.3552794 

predict(lm2, test, se.fit = TRUE)$se.fit
#        1         2         3         4         5         6 
#0.3321300 0.5475000 0.8533510 1.1077014 0.8495959 0.5977537 

predict(lm3, test, se.fit = TRUE)$se.fit
#        1         2         3         4         5         6 
#0.2976205 0.4906129 0.7646849 0.9926075 0.7613200 0.5356451