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

Question

我需要将lm()应用于我的数据框dat的一个不断增大的子集，同时对下一个观测值进行预测。例如，我正在执行以下操作：

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

我知道如何为特定的子集（与此问题相关： predict()和newdata - 这是如何工作的？）提供帮助。举个例子，如果要预测最后一行，我会做如下操作：

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的值。

谢谢。

- JohnnyDeer

---

我添加了一个示例，展示如何从输出中创建结果表格，请查看一下。我只为我的解决方案生成的3种输出对象中的1种做了这个示例，但是相同的方法适用于所有3种对象。 - Hack-R

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))

应用上述代码会得到结果（总共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平方值，而第二列是预测值。第一个adj.r2的值为NaN，因为我们拟合的第一个模型有3个系数用于3个数据点，因此没有可靠的统计数据可用。同样的情况也发生在se上，因为适配线没有0残差，所以预测是没有不确定性的。

- Zheyuan Li

---

你能回答我最后一个问题吗？ - JohnnyDeer

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

您在评论中提到您需要一个结果表格。您可以通过编程从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

- Hack-R

---

看起来还不错，但是我如何查看预测的数值或表格呢？通过一些研究，rollapply可能也是一个选择，但遗憾的是我对此了解甚少。 - JohnnyDeer

@JohnnyDeer 嗯，你可以使用一些滚动回归函数，但实际上你并不需要它们来完成这个任务。你想要的结果包含在这些数据集中，并且你可以通过 str(predict3) 来查看它们，或者直接访问整个对象或其元素。如果你想要一个结果表格，那也没问题，但请单独提出一个问题。使用 str 函数可以描述结构，以便你知道如何以编程方式访问所需的元素。 - Hack-R

@JohnnyDeer 我正准备更新，以便在输出数据集中提供更多信息。现在看一下。 - Hack-R

@JohnnyDeer 再次更新，为您提供结果表格示例。 - Hack-R