Оценить модель случайных эффектов повторных измерений с вложенной структурой, используя `plm()`
Можно ли оценить модель случайных эффектов повторных измерений с вложенной структурой, используя plm() из пакета plm?
Я знаю, что это возможно с lmer() из пакета lme4. Тем не мение, lmer() полагаться на рамки вероятности, и мне любопытно сделать это с plm(),
Вот мой минимальный рабочий пример, вдохновленный этим вопросом. Сначала некоторые необходимые пакеты и данные,
# install.packages(c("plm", "lme4", "texreg", "mlmRev"), dependencies = TRUE)
data(egsingle, package = "mlmRev")
набор данных egsingle это несбалансированная группа, состоящая из 1721 школьника, сгруппированных в 60 школах, по пяти временным точкам. Подробнее см. ?mlmRev::egsingle
Некоторое легкое управление данными
dta <- egsingle
dta$Female <- with(dta, ifelse(female == 'Female', 1, 0))
Кроме того, фрагмент соответствующих данных
dta[118:127,c('schoolid','childid','math','year','size','Female')]
#> schoolid childid math year size Female
#> 118 2040 289970511 -1.830 -1.5 502 1
#> 119 2040 289970511 -1.185 -0.5 502 1
#> 120 2040 289970511 0.852 0.5 502 1
#> 121 2040 289970511 0.573 1.5 502 1
#> 122 2040 289970511 1.736 2.5 502 1
#> 123 2040 292772811 -3.144 -1.5 502 0
#> 124 2040 292772811 -2.097 -0.5 502 0
#> 125 2040 292772811 -0.316 0.5 502 0
#> 126 2040 293550291 -2.097 -1.5 502 0
#> 127 2040 293550291 -1.314 -0.5 502 0
Теперь, в значительной степени полагаясь на ответ Роберта Лонга, я оцениваю модель случайных эффектов с повторными измерениями с помощью вложенной структуры, используя lmer() из пакета lme4,
dta$year <- as.factor(dta$year)
require(lme4)
Model.1 <- lmer(math ~ Female + size + year + (1 | schoolid /childid), dta)
# summary(Model.1)
Я заглянул в справочную страницу для plm() и у него есть команда индексации, index, но это займет всего один индекс и время, т. е. index = c("childid", "year"), игнорируя schoolid модель будет выглядеть так,
dta$year <- as.numeric(dta$year)
library(plm)
Model.2 <- plm(math~Female+size+year, dta, index = c("childid", "year"), model="random")
# summary(Model.2)
Подвести итог вопроса
Как я могу, или вообще возможно, указать модель случайных эффектов с повторными измерениями с вложенной структурой, например
Model.1, с помощьюplm()из пакета plm?
Ниже приведены фактические результаты оценки из двух моделей,
# require(texreg)
names(Model.2$coefficients) <- names(coefficients(Model.1)$schoolid) #ugly!
texreg::screenreg(list(Model.1, Model.2), digits = 3) # pretty!
#> ==============================================================
#> Model 1 Model 2
#> --------------------------------------------------------------
#> (Intercept) -2.693 *** -2.671 ***
#> (0.152) (0.085)
#> Female 0.008 -0.025
#> (0.042) (0.046)
#> size -0.000 -0.000 ***
#> (0.000) (0.000)
#> year-1.5 0.866 *** 0.878 ***
#> (0.059) (0.059)
#> year-0.5 1.870 *** 1.882 ***
#> (0.058) (0.059)
#> year0.5 2.562 *** 2.575 ***
#> (0.059) (0.059)
#> year1.5 3.133 *** 3.149 ***
#> (0.059) (0.060)
#> year2.5 3.939 *** 3.956 ***
#> (0.060) (0.060)
#> --------------------------------------------------------------
#> AIC 16590.715
#> BIC 16666.461
#> Log Likelihood -8284.357
#> Num. obs. 7230 7230
#> Num. groups: childid:schoolid 1721
#> Num. groups: schoolid 60
#> Var: childid:schoolid (Intercept) 0.672
#> Var: schoolid (Intercept) 0.180
#> Var: Residual 0.334
#> R^2 0.004
#> Adj. R^2 0.003
#> ==============================================================
#> *** p < 0.001, ** p < 0.01, * p < 0.05
1 ответ
На основании комментария Helix123 я написал следующую спецификацию модели для модели случайных эффектов с повторными измерениями с вложенной структурой, в plm() из пакета plm, используя метод Уоллеса и Хуссейна ( 1969), т.е. random.method = "walhus" для оценки компонентов дисперсии,
p_dta <- pdata.frame(dta, index = c("childid", "year", "schoolid"))
Model.3 <- plm(math ~ Female + size + year, data = p_dta, model = "random",
effect = "nested", random.method = "walhus")
Результаты, увиденные в Model.3 ниже, как близко к идентичным, к оценкам в Model.1, как я мог ожидать. Только перехват немного отличается (см. Вывод ниже).
Я написал выше на основе примера из Baltagi, Song and Jung ( 2001), представленного в
?plm, В примере Baltagi, Song and Jung ( 2001) компоненты дисперсии оцениваются сначала с использованием Swamy и Arora ( 1972), т.е.random.method = "swar"и второй с использованием Уоллеса и Хуссейна ( 1969). Только преобразование Nerlove ( 1971) не сходится с использованием данных Song и Jung ( 2001). Принимая во внимание, что только метод Уоллеса и Хуссейна ( 1969) мог сходиться, используяegsingleнабор данных.Любые авторитетные ссылки на это будут оценены. Я буду продолжать работать над этим.
names(Model.3$coefficients) <- names(coefficients(Model.1)$schoolid)
texreg::screenreg(list(Model.1, Model.3), digits = 3,
custom.model.names = c('Model 1', 'Model 3'))
#> ==============================================================
#> Model 1 Model 3
#> --------------------------------------------------------------
#> (Intercept) -2.693 *** -2.697 ***
#> (0.152) (0.152)
#> Female 0.008 0.008
#> (0.042) (0.042)
#> size -0.000 -0.000
#> (0.000) (0.000)
#> year-1.5 0.866 *** 0.866 ***
#> (0.059) (0.059)
#> year-0.5 1.870 *** 1.870 ***
#> (0.058) (0.058)
#> year0.5 2.562 *** 2.562 ***
#> (0.059) (0.059)
#> year1.5 3.133 *** 3.133 ***
#> (0.059) (0.059)
#> year2.5 3.939 *** 3.939 ***
#> (0.060) (0.060)
#> --------------------------------------------------------------
#> AIC 16590.715
#> BIC 16666.461
#> Log Likelihood -8284.357
#> Num. obs. 7230 7230
#> Num. groups: childid:schoolid 1721
#> Num. groups: schoolid 60
#> Var: childid:schoolid (Intercept) 0.672
#> Var: schoolid (Intercept) 0.180
#> Var: Residual 0.334
#> R^2 0.000
#> Adj. R^2 -0.001
#> ==============================================================
#> *** p < 0.001, ** p < 0.01, * p < 0.05#>