在数据框中应用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) 

我希望能够使用“预测”将这些模型应用于测试集，并创建一个矩阵，其中每个模型的结果作为一列。请保留HTML标记。

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语句将其合并为一两个步骤？

- SRB

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 

- Zheyuan Li

---

谢谢，这很有帮助。另外，如果我想使用另一个模型，比如glmnet，我该用什么作为'y'值呢？我尝试了上面的方法，但那种形式不被接受。 - undefined