Полосы доверия регрессии от моделей PLM

Question

В настоящее время я использую модель с фиксированными эффектами с использованием функции plm из одноименного пакета. У меня есть кластеризация в данных, поэтому я использую функцию vcovCR из пакета clubSandwich для получения стандартных ошибок оценочных коэффициентов регрессии.

Я хочу сделать следующее: Создать график, который показывает прогнозируемое значение Y с учетом X1 (когда другие переменные поддерживаются постоянными при их средних значениях), а также добавляет 95% доверительные полосы к графику.

Что я пробовал:

1. Использование функции ggpredict() из пакета ggeffects. Это не сработало из-за того, что пакет plm явно проблематичен c для ggeffects.
2. Вручную с использованием «группового» подхода оценки фиксированных эффектов, а затем запустив обычную функцию lm() по центру группы данных. Выполнение ggpredict() на этих данных дает мне хорошие доверительные интервалы, но тогда значения X и Y центрируются по группам (удивление, удивление!), И я хочу, чтобы предсказания и значения X были в исходном масштабе, а не в средняя центрированная шкала.

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

Answer 1

Я пытался что-то подобное, однако я не мог найти прямого решения. Я оцениваю исследование событий с лагами и лидами на модели fe. Для меня сработало следующее:

ES1<-plm(total_thefts ~ lag6 + lag5 + lag4 + lag3 + lag2 + E + lead1 + lead2 + lead3 + lead4 + lead5 + lead6 + factor(month) + factor(year), data = panel, model="within")

Поскольку у plm нет прямого способа оценить надежные стандартные ошибки, я делаю это вручную и заменяю исходные:

G <- length(unique(panel$id))
c <- G/(G - 1)
robust_se <- coeftest(ES1, c * vcovHC(ES1, type = 'HC1', cluster = 'group'))
ES1<-summary(ES1)
ES1$coefficients[,2:4] <- robust_se[,2:4]

После этого я разделяю свои коэффициенты интереса и оцениваю полосы CI (95%), следуя этому post :

coefs <- tidy(ES1[["coefficients"]])
coefs$conf.low <- coefs$Estimate+c(-1)*coefs$Std..Error*qt(0.975,42)
coefs$conf.high <- coefs$Estimate+c(1)*coefs$Std..Error*qt(0.975,42)

Наконец, после этого я сохраняю свои параметры интереса и график исследование события:

interest <- c("lag6","lag5","lag4","lag3","lag2","E","lead1","lead2","lead3","lead4","lead5","lead6")
coefs <- subset(coefs,coefs$.rownames%in%interest)
coefs$time <- c(seq(-6,-2,1),seq(0,6,1))

ggplot(coefs, aes(time, Estimate))+
           geom_line() +
           geom_point()+
           geom_pointrange(aes(ymin = conf.low, ymax = conf.high))+
           geom_hline(yintercept = 0, linetype = 2) +
           ggtitle(paste0("N: ",nobs(ES1),", Sample: ",pdim(ES1)$nT$n,", Buff: 1000m"))) +
           theme(plot.title = element_text(hjust = 0.5)))