Как я запускаю множественную регрессию в R, каждый раз добавляя 100 дополнительных строк

Вот минимальный пример использования набора данных 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")

^{Created on 2019-04-07 by the reprex package (v0.2.1)}

Рассмотрите возможность разделения набора данных по фирмам, а затем итеративного запуска lm с использованием последовательности seq(1, 4000, by=100) для подмножества списка разделенных фреймов данных:

# 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

The second task, is I want to save the beta coefficient of each regression (i.e. each regression after I add extra 100 firms), and then plot beta on Y and number of firms on x (from 100 to 4000) to observe how beta changes when adding firms.

Вы можете использовать мой пакет rollRegres. Это почти идентично примеру в эта виньетка:

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.