Создайте объект lm-класса с загрузочной загрузкой кластера (R) или: загрузочную загрузку кластера для sensemakr()

Обновлено: Мне бы хотелось помочь закрыть мой вопрос. Этот вопрос изначально был опубликован на Cross Validated, и должен был остаться там, потому что ответ - это скорее вопрос статистики, а не программирования, но там он был закрыт, и его попросили переопубликовать здесь. Я так и сделал, но ответ на самом деле вопрос статистики, поэтому я опубликую ответ здесь. Если кто-то захочет закрыть это, я буду признателен.


Я пытался выполнить кластерную загрузку с помощью lm, которую затем можно было бы использовать в пакете sensemakr(). Если вы с ним не знакомы, то это реализация некоторых идей Cinelli & Hazlett (2020, в Журнале Королевского статистического общества: Статистическая методология, Серия B).

Принцип работы sensemakr() заключается в том, что вы запускаете модель регрессии, которая создает объект класса lm или feols, который затем можно использовать в sensemakr().

library(palmerpenguins)
library(sensemakr)

penguin_dat<-penguins

model_lm<-lm(flipper_length_mm ~ sex + body_mass_g, data=penguin_dat)
summary(model_lm)

sensitivity1<-sensemakr(model=model_lm, treatment = "body_mass_g", benchmark_covariates = "sexmale")
summary(sensitivity1)

Мне нужно найти способ передать в sensemakr модели с более сложной структурой ошибок, в частности модель с кластерной начальной загрузкой. В настоящее время он использует модели lm, но я читал, что вы можете попробовать передать модели feols из fixest в него, если хотите использовать версию для разработки (см.: Используйте sensemakr с моделью fixest feols (R)). Однако сам я еще не пробовал.

Вот варианты, которые я рассмотрел:

  1. Используйте мой предыдущий объект класса lm (model_lm в предыдущем примере) для lm.boot() в simpleboot, затем используйте этот объект в sensemakr. Проблема: я не могу найти какую-либо функциональность для кластерной загрузки в простой загрузке. Боюсь, я что-то упускаю из виду. Надеюсь, я просто упустил из виду эту функцию, поскольку это, вероятно, самое простое решение.
  2. Используйте загрузочную загрузку дикого кластера из fwildclusterboot. Я не уверен, какой объект он создает, и он был удален из CRAN из-за проблем с зависимостями, поэтому я не могу его установить (или посмотреть справочное руководство).
  3. Используйте пакет начальной загрузки кластера, например clusbootglm или lmeresampler, а затем попытайтесь привести результат к объекту типа lm. Это несколько выходит за рамки моих R-возможностей, но если пункт 1 не сработает, думаю, это единственный вариант. В качестве примера того, как будет выглядеть #3, вот пример использования coeftest() из lmtest и vcovBS() из сэндвича, но я не могу понять, как можно вернуть это обратно в объект lm или даже если бы это была хорошая идея:
library(palmerpenguins)
library(sandwich)
library(sensemakr)
library(lmtest)

penguin_dat<-penguins

model_lm<-lm(flipper_length_mm ~ sex + body_mass_g, data=penguin_dat)

clustered_bootstrap_model_lm<-coeftest(model_lm, vcov = vcovBS, cluster = penguin_dat$species, R = 1000)

Ирония в том, что авторы sensemakr() также реализовали это в Stata, и очень просто выполнить кластерную загрузку, а затем использовать стандартные ошибки в другой команде.

sysuse auto2.dta
regress price mpg weight rep78, vce(bootstrap, reps(100)) cluster(foreign)
est sto m1

Но реализация Stata не позволяет вам передавать модель, вместо этого вам приходится запускать регрессию внутри самого sensemakr... и когда вы это делаете, она запрещает более сложные структуры ошибок.

Есть ли у кого-нибудь идеи о том, как это сделать? Я очень ценю любой совет, который вы можете дать.

Цитаты: Синелли, К., и Хэзлетт, К. (2020). Понимание чувствительности: расширение смещения пропущенной переменной. Журнал Королевского статистического общества, серия B: Статистическая методология, 82 (1), 39–67. https://doi.org/10.1111/rssb.12348

Синелли, К., Ферверда, Дж., и Хэзлетт, К. (2020). sensemakr: инструменты анализа чувствительности для OLS в R и Stata. Доступно на SSRN: https://ssrn.com/abstract=3588978 или http://dx.doi.org/10.2139/ssrn.3588978


54
2

Ответы:

Довольно просто «взломать» sensemakr и включить произвольные стандартные ошибки в анализ чувствительности, но неясно, является ли это статистически достоверным. Многие выводы в Cinelli & Hazlett (2020) основаны на стандартных результатах линейной регрессии с обычными стандартными ошибками. Пакет выдает предупреждение при использовании fixest с любыми стандартными ошибками, кроме обычных стандартных ошибок iid. Степени свободы, используемые в расчетах, предполагают, что единицы измерения независимы.

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


Еще раз подчеркиваю: не используйте этот код в своем анализе; результаты недействительны.

sensemakr() извлекает статистику из lm объектов, используя summary(). В частности, он извлекает стандартные ошибки и t-статистику из выходных данных summary(). Я смог узнать об этом, просматривая исходный код пакета и наблюдая, где используются стандартные ошибки и как они извлекаются.

Что можно сделать, так это дать объекту lm новый класс и метод summary() из этого класса, который использует предоставленную вами матрицу отклонений вместо вычисления стандартных ошибок с использованием предположений МНК.

Для этого мы напишем новый summary() метод для класса lm_custom. Этот метод извлекает ковариационную матрицу из компонента vcov предоставленного объекта.

summary.lm_custom <- function(object, ...) {
  s <- NextMethod("summary", object)
  
  if (!is.null(object[["vcov"]])) {
    rdf <- s$df[2]
    
    s$coefficients[, "Std. Error"] <- sqrt(diag(object[["vcov"]]))
    s$coefficients[, "t value"] <- s$coefficients[, "Estimate"] / s$coefficients[, "Std. Error"]
    s$coefficients[, "Pr(>|t|)"] <-  2 * pt(abs(s$coefficients[, "t value"]), rdf, lower.tail = FALSE)
  }
  
  s
}

Далее вам нужно вычислить ковариационную матрицу, присвоить ее компоненту vcov объекта lm и присвоить объекту новый класс.

model_lm2 <- model_lm
model_lm2$vcov <- sandwich::vcovBS(model_lm, cluster = ~species)
class(model_lm2) <- c("lm_custom", class(model_lm2))

Запуск summary() для этого объекта дает стандартную ошибку, t-статистику и значения p с использованием включенной ковариационной матрицы.

summary(model_lm2)

Если вы вызовете sensemakr() для нового объекта, он будет использовать стандартные ошибки, вычисленные на основе предоставленной ковариационной матрицы. Повторяю, этот анализ недействителен; Доказательством этого является то, что выходные данные, которые теоретически не должны зависеть от используемых стандартных ошибок (например, результаты, рассчитанные с использованием показателей R2), изменяются при использовании новых стандартных ошибок. Это связано с тем, что пакет полагается на связи между t-статистикой и R2, которые действительны только при использовании обычных стандартных ошибок OLS. НЕ ИСПОЛЬЗУЙТЕ, НЕ СООБЩАЙТЕ ИЛИ ИНТЕРПРЕТИРУЙТЕ ЭТОТ АНАЛИЗ.


Решено

Я разговаривал с основным автором программы, и предыдущий ответ на этот вопрос технически верен, но не то, что мне нужно.

Настоящий ответ — это вопрос статистики, а не программирования, и поэтому его не следует закрывать на Cross Validated. Вы можете изменить объект, но на самом деле в этом нет никакого смысла, поскольку sensemakr() почти во всех своих вычислениях опирается на точечные оценки, а не на стандартные ошибки. Большими исключениями являются (1) значения устойчивости для проверки нулевой гипотезы и (2) статистика стандартной ошибки/доверительного интервала/теста для скорректированной оценки.

Для этих двух правильных ответов будет загрузить сам sensemakr и извлечь в вектор следующие оценки: rv_q (нижняя граница которого даст вам самозагрузочный эквивалент rv_qa) и скорректированную_эстимату.

Большую часть этого кода мне предоставил Карлос Синелли (ведущий автор пакета sensemakr), и я очень ему за это благодарен.

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

adjusted_estimate_boot <- rep(NA, B) # vector to store results

Затем повторно дискретизируйте набор данных много раз, каждый раз извлекая скорректированные оценки и добавляя их в вектор.

B <- 5e2; # number of bootstrap samples
<- nrow(darfur) # number of observations in the full data

# boostrap loop
for(i in 1:B){
  cat(i,"out of", B, "\n")
  
  # resample data
  idx_boot <- sample(1:n, size = n, replace = T)
  dat_boot <- darfur[idx_boot, ]
  
  # fit model
  my.ols_boot <- lm(peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
                      pastvoted + hhsize_darfur + female + village, data = dat_boot)
  
  # sensemakr
  sense.out_boot <- sensemakr(my.ols_boot, treatment = "directlyharmed",
                              benchmark_covariates = "female",
                              kd = 1)
  
  # save the estimate you want
  adjusted_estimate_boot[i] <- sense.out_boot$bounds[1,"adjusted_estimate"]
}

А затем вычислите новый доверительный интервал:

sig.level <- 0.05
quantile(adjusted_estimate_boot, c(sig.level/2, 1-sig.level/2))

Вы можете сделать то же самое, чтобы извлечь rv_q, и нижняя граница rv_q в этой последней части эквивалентна загруженному rv_qa.