Применение «функций кластеризации» к ряду линейных моделей - PullRequest
Новая
       36

Применение «функций кластеризации» к ряду линейных моделей

1 голос
/ 03 мая 2019

Я хочу перебрать список линейных моделей и применить «кластерные» стандартные ошибки к каждой модели, используя функцию vcovCL.Моя цель - сделать это как можно более эффективно (я использую линейную модель для множества столбцов данных).Моя проблема заключается в попытке указать дополнительные аргументы внутри анонимной функции.Ниже я симулирую некоторые поддельные данные.Участки представляют мое поперечное сечение;месяцы представляют собой мое измерение времени (5 единиц наблюдалось за 4 месяца).Переменная int является фиктивной, когда происходит вмешательство. 

df <- data.frame(
  precinct = c( rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4), rep(5, 4) ),
  month = rep(1:4, 5),
  crime = rnorm(20, 10, 5),
  int = c(c(0, 1, 1, 0), rep(0, 4), rep(0, 4), c(1, 1, 1, 0), rep(0, 4))
  )

df[1:10, ]

outcome <- df[3]
est <- lapply(outcome, FUN = function(x) { lm(x ~ as.factor(precinct) + as.factor(month) + int, data = df) })

se <- lapply(est, function(x) { sqrt(diag(vcovCL(x, cluster = ~ precinct + month))) })

При добавлении аргумента cluster внутри функции vcovCL я получаю следующее сообщение об ошибке. 

Error in eval(expr, envir, enclos) : object 'x' not found

Единственный способ обойти это, по моим оценкам, - проиндексировать фрейм данных, то есть df$, а затем указать переменные кластеризации.Можно ли этого достичь, указав дополнительный аргумент для df внутри вызова функции? Эффективен ли этот код?

Возможно, формальное определение уравнения модели - лучший способ, я полагаю.

Любые мысли / комментарии всегда полезны:)

1 Ответ

0 голосов
/ 03 мая 2019

Вот один из подходов, который позволяет получать кластерные стандартные ошибки для нескольких моделей: 

library(sandwich)

# I am going to use the same model three times to get the "sequence" of linear models. 
mod <- lm(crime ~ as.factor(precinct) + as.factor(month) + int, data = df)

# define function to retrieve standard errors:
robust_se <- function(mod) {sqrt(diag(vcovCL(mod, cluster = list(df$precinct, df$month))))}

# apply function to all models:
se <- lapply(list(mod, mod, mod), robust_se)

Если вы хотите настроить весь вывод, может быть полезно следующее: 

library(lmtest)
adj_stats <- function(mod) {coeftest(mod, vcovCL(mod, cluster = list(df$precinct, df$month)))}

adjusted_models <- lapply(list(mod, mod, mod), adj_stats)

Для решения проблемы с несколькими столбцами:

В случае, если вам приходится работать с линейными моделями над несколькими столбцами, может быть полезно следующее.Все вышеперечисленное останется прежним, за исключением того, что вы передадите свой список моделей в lapply.

Сначала давайте воспользуемся этим фреймом данных: 

df <- data.frame(
  precinct = c( rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4), rep(5, 4) ),
  month = rep(1:4, 5),
  crime = rnorm(20, 10, 5),
  crime2 = rnorm(20, 10, 5),
  crime3 = rnorm(20, 10, 5),
  int = c(c(0, 1, 1, 0), rep(0, 4), rep(0, 4), c(1, 1, 1, 0), rep(0, 4))
)

Давайте определим столбцы результатов: 

outcome_columns <- c("crime", "crime2", "crime3")

Теперь давайте запустим регрессию с каждым результатом: 

models <- lapply(outcome_columns, 
         function(outcome) lm( eval(parse(text = paste0(outcome, " ~ as.factor(precinct) + as.factor(month) + int"))), data = df) )

И тогда вы просто позвоните 

adjusted_models <- lapply(models, adj_stats)

Относительно эффективности:

Приведенный выше код эффективен тем, что его легко настроить и быстро записать.В большинстве случаев это будет прекрасно.Для эффективности вычислений обратите внимание, что ваша матрица проектирования одинакова во всех случаях, то есть, предварительно вычислив общие элементы (например, inv(X'X)*X'), вы можете сохранить некоторые вычисления.Однако вы бы упустили удобство многих встроенных функций.

...