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

Question

Я надеюсь получить помощь по следующей проблеме в R.

У меня есть 4 переменные, firm ID, sales, size, date, для почти 4000 фирм.

Я хочу запустить эту регрессию:

lm(size~sales), одновременно добавляя 100 фирм из 4000.

Итак, первая регрессия будет иметь 100 фирм, вторая будетесть 200, у третьего будет 300 ... до достижения последней регрессии, которая включает все фирмы (4000).

Вторая задача - я хочу сохранить коэффициент бета для каждой регрессии (т.е. каждой регрессии после добавления 100 дополнительных фирм), а затем построить бета для Y и количество фирм для x (от 100 до 4000) наблюдать, как изменяется бета при добавлении фирм.

Нужен ли какой-то цикл для регрессий, цикл для сохранения бета-версий и цикл для построения графиков?Спасибо за чтение

Answer 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

Answer 2

Вторая задача - я хочу сохранить бета-коэффициент каждой регрессии (то есть каждой регрессии после добавления 100 дополнительных фирм), а затем нанести бета на Y и количество фирм на x (от 100 до 4000) для наблюдения. как меняется бета при добавлении фирм.

Вы можете использовать мой пакет 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 впоследствии.

Answer 3

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

^{Создано в 2019-04-07 пакетом представ. (v0.2.1)}