使用 lm() 和 predict() 进行滚动回归和预测

Question

我需要将lm() 应用于我的数据框dat 的一个扩大子集，同时对下一次观察进行预测。比如我在做：

fit model      predict
----------     -------
dat[1:3, ]     dat[4, ]
dat[1:4, ]     dat[5, ]
    .             .
    .             .
dat[-1, ]      dat[nrow(dat), ]

我知道我应该为特定子集做什么（与此问题相关：predict() and newdata - How does this work?）。例如要预测最后一行，我会这样做

dat1 = dat[1:(nrow(dat)-1), ]
dat2 = dat[nrow(dat), ]

fit = lm(log(clicks) ~ log(v1) + log(v12), data=dat1)
predict.fit = predict(fit, newdata=dat2, se.fit=TRUE)

如何为所有子集自动执行此操作，并可能将我想要的内容提取到表中？

• 来自fit，我需要summary(fit)$adj.r.squared；
• 来自predict.fit 我需要predict.fit$fit 值。

谢谢。

Answer 1

（高效）解决方案

这是你可以做的：

p <- 3  ## number of parameters in lm()
n <- nrow(dat) - 1

## a function to return what you desire for subset dat[1:x, ]
bundle <- function(x) {
  fit <- lm(log(clicks) ~ log(v1) + log(v12), data = dat, subset = 1:x, model = FALSE)
  pred <- predict(fit, newdata = dat[x+1, ], se.fit = TRUE)
  c(summary(fit)$adj.r.squared, pred$fit, pred$se.fit)
  }

## rolling regression / prediction
result <- t(sapply(p:n, bundle))
colnames(result) <- c("adj.r2", "prediction", "se")

注意我在bundle 函数中做了几件事：

• 我已使用 subset 参数来选择要适合的子集
• 我使用model = FALSE 不保存模型框架，因此我们保存了工作空间

总体来说没有明显的循环，但是使用了sapply。

• 拟合从p 开始，这是拟合具有p 系数的模型所需的最小数据数；
• 拟合终止于nrow(dat) - 1，因为我们至少需要最后一列进行预测。

---

测试

示例数据（包含 30 个“观察结果”）

dat <- data.frame(clicks = runif(30, 1, 100), v1 = runif(30, 1, 100),
                  v12 = runif(30, 1, 100))

应用上面的代码得到results（总共27行，5行截断输出）

            adj.r2 prediction        se
 [1,]          NaN   3.881068       NaN
 [2,]  0.106592619   3.676821 0.7517040
 [3,]  0.545993989   3.892931 0.2758347
 [4,]  0.622612495   3.766101 0.1508270
 [5,]  0.180462206   3.996344 0.2059014

第一列是拟合模型的调整后的 R.squared 值，而第二列是预测值。 adj.r2 的第一个值是 NaN，因为我们拟合的第一个模型有 3 个数据点的 3 个系数，因此没有可用的合理统计数据。 se 也有同样的情况，因为拟合线没有 0 残差，所以预测没有不确定性。

Answer 2

我只是编造了一些随机数据用于这个例子。我将对象称为data，因为在我编写此解决方案时，它在问题中被称为（随便你怎么称呼它）。

（高效）解决方案

data <- data.frame(v1=rnorm(100),v2=rnorm(100),clicks=rnorm(100))

data1 = data[1:(nrow(data)-1), ]
data2 = data[nrow(data), ]

for(i in 3:nrow(data)){
  nam  <- paste("predict", i, sep = "")
  nam1 <- paste("fit", i, sep = "")
  nam2 <- paste("summary_fit", i, sep = "")

  fit = lm(clicks ~ v1 + v2, data=data[1:i,])
  tmp  <- predict(fit, newdata=data2, se.fit=TRUE)
  tmp1 <- fit
  tmp2 <- summary(fit)
  assign(nam, tmp)
  assign(nam1, tmp1)
  assign(nam2, tmp2)
}

您想要的所有结果都将存储在由此创建的数据对象中。

例如：

> summary_fit10$r.squared
[1] 0.3087432

您在 cmets 中提到您想要一个结果表。您可以通过编程方式从 3 种类型的输出文件创建结果表，如下所示：

rm(data,data1,data2,i,nam,nam1,nam2,fit,tmp,tmp1,tmp2)
frames <- ls()

frames.fit     <- frames[1:98] #change index or use pattern matching as needed
frames.predict <- frames[99:196]
frames.sum     <- frames[197:294]

fit.table <- data.frame(intercept=NA,v1=NA,v2=NA,sourcedf=NA)
for(i in 1:length(frames.fit)){
  tmp <- get(frames.fit[i])
  fit.table              <- rbind(fit.table,c(tmp$coefficients[[1]],tmp$coefficients[[2]],tmp$coefficients[[3]],frames.fit[i]))
}

fit.table

> fit.table
             intercept                   v1                   v2 sourcedf
2  -0.0647017971121678     1.34929652763687   -0.300502017324518    fit10
3  -0.0401617893034109   -0.034750571912636  -0.0843076273486442   fit100
4   0.0132968863522573     1.31283604433593   -0.388846211083564    fit11
5   0.0315113918953643     1.31099122173898   -0.371130010135382    fit12
6    0.149582794027583    0.958692838785998   -0.299479715938493    fit13
7  0.00759688947362175    0.703525856001948   -0.297223988673322    fit14
8    0.219756240025917    0.631961979610744   -0.347851129205841    fit15
9     0.13389223748979    0.560583832333355   -0.276076134872669    fit16
10   0.147258022154645    0.581865844000838   -0.278212722024832    fit17
11  0.0592160359650468    0.469842498721747   -0.163187274356457    fit18
12   0.120640756525163    0.430051839741539   -0.201725012088506    fit19
13   0.101443924785995     0.34966728554219   -0.231560038360121    fit20
14  0.0416637001406594    0.472156988919337   -0.247684504074867    fit21
15 -0.0158319749710781    0.451944113682333   -0.171367482879835    fit22
16 -0.0337969739950376    0.423851304105399   -0.157905431162024    fit23
17  -0.109460218252207     0.32206642419212   -0.055331391802687    fit24
18  -0.100560410735971    0.335862465403716  -0.0609509815266072    fit25
19  -0.138175283219818    0.390418411384468  -0.0873106257144312    fit26
20  -0.106984355317733    0.391270279253722  -0.0560299858019556    fit27
21 -0.0740684978271464    0.385267011513678  -0.0548056844433894    fit28