我如何在 R 中运行多元回归，同时每次添加 100 行答案

【问题标题】：How I run multiple regression in R while adding 100 additional rows every time我如何在 R 中运行多元回归，同时每次添加 100 行
【发布时间】：2019-04-08 00:06:10
【问题描述】：

我希望在R中的以下问题上得到帮助。

我有 4 个变量，firm ID、sales、size、date，用于近 4,000 家公司。

我想运行这个回归：

lm(size~sales)，同时从 4000 家一次添加 100 家。

因此，第一次回归将有 100 家公司，第二次将有 200 家，第三次将有 300 家……直到最后一次回归包括所有公司（4000 家）。

第二个任务，我想保存每个回归的 beta 系数（即我添加额外 100 家公司后的每个回归），然后在 Y 上绘制 beta 和在 x 上绘制公司数量（从 100 到 4000）以观察添加公司时 beta 如何变化。

我是否需要某种用于回归的循环、用于保存 beta 的循环和用于绘图的循环？感谢阅读

【问题讨论】：

标签： r loops regression

【解决方案1】：

考虑按公司拆分数据集，然后使用序列 seq(1, 4000, by=100) 迭代运行 lm 以子集拆分数据框列表：

# BUILD A LIST OF DATA FRAMES (SIZE = 4,000)
firms_df_list <- split(df, df$firm_id)

# FUNCTION TO CALL lm() AND EXTRACT RESULTS
lm_results <- function(n, df) {

  model <- lm(sales ~ size, data = df)
  res <- summary(model)

  p <- res$fstatistic
  c(num_of_firms = n,
    sales = res$coefficients[2,1],
    std_err = res$coefficients[2,2],
    t_stat = res$coefficients[2,3],
    t_pvalue = res$coefficients[2,4],
    r_sq = res$r.squared,
    adj_r_sq = res$adj.r.squared,
    f_stat = p[['value']],
    f_pvalue = unname(pf(p[1], p[2], p[3], lower.tail=FALSE))
  )
}

# BUILD MATRIX RESULTS WHERE ROWS ARE MODEL RUNS AND COLS ARE RESULT ESTIMATES
mat_results <- t(sapply(seq(1, 4000, by=100), function(i) {
     # COMBINE FIRM SUBSETS BY RANGE
     curr_df <- do.call(rbind, firms_df_list[1:i])

     # CALL MODEL AND RETRIEVE RESULTS
     lm_results(i, curr_df)
}))

# PLOT ALL SALES BETAS AND NUMBER OF FIRMS
plot(mat_results[,"num_of_firms"], mat_results[,"sales"], type="b", 
     col="blue", lwd=1, pch=16, xlab="Number of Firms", ylab="Sales Estimate")

考虑到年份和月份的细分，考虑将by（类似于split + lapply）按年份和月份与内部split（类似于上述过程）进行子集，其中每次迭代运行所需的模型.然后，在每个月和年级别绑定矩阵以获得最终的大矩阵。注意：lm_results 现在接收两个用于指标月份和年份矩阵列的参数。

# FUNCTION TO CALL lm() AND EXTRACT RESULTS
lm_results <- function(n, df, yy, mm) {

  model <- lm(sales ~ size, data = df)
  res <- summary(model)

  p <- res$fstatistic
  c(year = yy,
    month = mm,
    num_of_firms = n,
    sales = res$coefficients[2,1],
    std_err = res$coefficients[2,2],
    t_stat = res$coefficients[2,3],
    t_pvalue = res$coefficients[2,4],
    r_sq = res$r.squared,
    adj_r_sq = res$adj.r.squared,
    f_stat = p[['value']],
    f_pvalue = unname(pf(p[1], p[2], p[3], lower.tail=FALSE))
  )
}    

# BUILD A LIST OF MONTHLY MATRICES BY YEAR
firms_mat_list <- by(df, df$yy, function(sub_year){

  # BUILD A LIST OF FIRM MATRICES BY MONTH
  month_mat_list <- by(sub_year, sub_year$mm, function(sub_month){

    firms_df_list <- split(sub_month, sub_month$firm)

    # BUILD MATRIX RESULTS WHERE ROWS ARE MODEL RUNS AND COLS ARE RESULT ESTIMATES
    mat_results <- t(sapply(seq(1, 4000, by=100), function(i) {
      # COMBINE FIRM SUBSETS BY RANGE
      curr_df <- do.call(rbind, firms_df_list[1:i])

      # CALL MODEL AND RETRIEVE RESULTS
      lm_results(i, curr_df, curr_df$yy[1], curr_df$mm[1])
    }))

  })

  do.call(rbind, month_mat_list)
})

firms_matrix <- do.call(rbind, firms_mat_list)

firms_matrix

【讨论】：

谢谢，它奏效了，但是，每个公司都是按月观察的（我在最初的问题中没有提到），所以我有一个专栏多年（有些是从 2000 年到 2016 年观察到的，有些时间跨度较小）和对应于每年的月份列。考虑到年份和月份，我认为它会影响回归输出，您能否建议在运行model <- lm(sales ~ size, data = df) 时如何考虑每个公司的年份和月份？我需要几个月甚至几年的循环吗？
哇！我以为我把你丢到互联网的另一边了！有时我们会得到消失的 OP。但是，几个月和几年的帐户到底是什么意思。您想进一步按年/月划分公司规模还是将它们包含在回归中？
谢谢！！我的意思是，从公司 1 到公司 4000 的每家公司多年来每月观察一次，mm 变量从 1 到 12 和 yy 从 2005 年到 2016 年。我的变量是：firm mm yy size sale 基本上，我想运行回归以包括所有每年的月份，然后一次添加 100 家公司，即涵盖所有年份的所有月份和所有 4000 家公司。我认为您提供的代码会立即获取所有观察结果......可能我们必须在数月和数年内复制回归以考虑时间？再次感谢您，期待。
在使用split 运行建模之前，请使用by 查看按年和月划分的更新。

【解决方案2】：

这是一个使用mtcars 数据集的最小示例。我建立了一个回归，一次添加一行。我将结果向量预分配到右侧，然后遍历行并存储系数结果。

results <- vector(length = nrow(mtcars))
for (j in 1:nrow(mtcars)){
  results[j] <- coef(lm(mpg ~ hp, data = mtcars[1:j, ]))[2]
}

plot(x = 1:nrow(mtcars), y = results, type = "p")

^{由reprex package (v0.2.1) 于 2019-04-07 创建}

【讨论】：

【解决方案3】：

第二个任务，我想保存每个回归的 beta 系数（即我添加额外 100 家公司后的每个回归），然后在 Y 上绘制 beta 和在 x 上绘制公司数量（从 100 到 4000）以观察添加公司时 beta 如何变化。

您可以使用我的rollRegres 包。这与this vignette 中的示例几乎相同：

set.seed(65731482)
ngrp <- 40L
n_per_g <- 100L
# create group variable
grp <- c(sapply(1:ngrp, rep, times = n_per_g))
n <- n_per_g * ngrp
p <- 1L
X <- matrix(rnorm(p * n), n, p)
y <- drop(X %*% 1.5) + rnorm(n)

library(rollRegres)
out <- roll_regres(y ~ X, do_downdates = FALSE, width = 100L)
beta <- out$coefs

# check result
tail(out$coefs, 2)
#R      (Intercept)    X
#R 3999    -0.00552 1.51
#R 4000    -0.00571 1.51
coef(lm(y ~ X))
#R (Intercept)           X 
#R    -0.00571     1.51405 

# plot 
plot(out$coefs[, 2], xlab = "Time", ylab = "slope", type = "l")

它为您提供了所有 40000 - 99 的值，但速度很快，因此您可能不会关心额外的计算

microbenchmark::microbenchmark(
  roll_regres(y ~ X, do_downdates = FALSE, width = 100L))
#R Unit: microseconds
#R                                                   expr min  lq mean median  uq  max neval
#R roll_regres(y ~ X, do_downdates = FALSE, width = 100L) 740 750  771    763 772 1090   100

然后你可以子集beta。

【讨论】：