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

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

从 fit 开始，我需要 summary(fit)$adj.r.squared;
来自 predict.fit 我需要 predict.fit$fit 值。

谢谢。

Answer 1

我只是编造了一些随机数据用于此示例。我将对象称为 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

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

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

Answer 2

(高效)解决方案

这是你可以做的：

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 "observations"）

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 残差，所以预测是在没有不确定性的情况下完成的。

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

Rolling regression and prediction with lm() and predict()

regression

r

linear-regression

predict

lm