在不平衡面板上使用 plm 对杂波样本进行随机效应估计
Random Effects Esimation on cluter sample using plm on a unbalanced panel
我通常在广义最小二乘框架内工作,估计 Wooldridge 的 Introductory (2013) 调用的随机效应和固定效应模型对由个人和时间索引的纵向数据的影响维度.
我一直在使用plm包中plm()
中的可行GLS估计来估计随机效应模型——一些统计数据文学术语混合模型。 plm()
函数采用 index
参数,我在其中指示个人和时间索引。但是,我现在面对一些数据,其中每个人在每个时间点都有多个度量,即多好的分组结构。
我发现可以使用 lme4 包中的 lmer()
来建立这样的模型,但是我对行话的差异以及可能性感到有点困惑框架,我想知道是否正确指定了模型。由于我不熟悉框架和这个术语,我担心我可能会忽略更多内容。
我可以使用 lmer()
复制我常用的 plm()
模型,但我不确定如何添加分组。我试图在下面说明我的问题。
我找了一些和我的数据有点像的数据来说明我的情况。首先是一些需要的包,
install.packages(c("mlmRev", "plm", "lme4", "stargazer"), dependencies = TRUE)
然后是数据
data(egsingle, package = "mlmRev")
egsingle
是一个不平衡的小组,由 1721 名学童组成,分为 60 所学校,跨越五个时间点。这些数据最初与 HLM 软件包(Bryk、Raudenbush 和 Congdon,1996)一起分发,但可以在 mlmrev 软件包中找到,详情请参阅 ? 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
以下是我如何使用 plm()
、
设置 随机效应模型 而没有 schoolid
library(plm)
reg.re.plm <- plm(math~Female+size+year, dta, index = c("childid", "year"), model="random")
# summary(reg.re.plm)
我可以像这样重现这些结果lme4
require(lme4)
dta$year <- as.factor(dta$year)
reg.re.lmer <- lmer(math~Female+size+year+(1|childid), dta)
# summary(reg.re.lmer)
现在,从阅读 Bates (2010) 的第 2 章“lme4:混合效应建模
使用 R” 我相信这就是我指定模型的方式,包括集群级别,schoolid
,
reg.re.lmer.in.school <- lmer(math~Female+size+year+(1|childid)+(1|schoolid), dta)
# summary(reg.re.lmer.in.school)
但是,当我查看结果时,我不太相信我实际上已经正确指定了它(见下文)。
在我的实际数据中,重复的措施是在个人内部,但我认为我可以以这个数据为例。对于如何进行的任何建议,我将不胜感激。也许是对 notation/terminology 的工作示例的引用,与 Wooldridge (2013) 中使用的示例相差不远。而且,我如何逆向工作并编写 reg.re.lmer.in.school
模型的规范?
# library(stargazer)
stargazer::stargazer(reg.re.plm, reg.re.lmer, reg.re.lmer.in.school, type="text")
#> =====================================================================
#> Dependent variable:
#> -------------------------------------------------
#> math
#> panel linear
#> linear mixed-effects
#> (1) (2) (3)
#> ---------------------------------------------------------------------
#> Female -0.025 -0.025 0.008
#> (0.046) (0.047) (0.042)
#>
#> size -0.0004*** -0.0004*** -0.0003
#> (0.0001) (0.0001) (0.0002)
#>
#> year-1.5 0.878*** 0.876*** 0.866***
#> (0.059) (0.059) (0.059)
#>
#> year-0.5 1.882*** 1.880*** 1.870***
#> (0.059) (0.058) (0.058)
#>
#> year0.5 2.575*** 2.574*** 2.562***
#> (0.059) (0.059) (0.059)
#>
#> year1.5 3.149*** 3.147*** 3.133***
#> (0.060) (0.059) (0.059)
#>
#> year2.5 3.956*** 3.954*** 3.939***
#> (0.060) (0.060) (0.060)
#>
#> Constant -2.671*** -2.669*** -2.693***
#> (0.085) (0.086) (0.152)
#>
#> ---------------------------------------------------------------------
#> Observations 7,230 7,230 7,230
#> R2 0.735
#> Adjusted R2 0.735
#> Log Likelihood -8,417.815 -8,284.357
#> Akaike Inf. Crit. 16,855.630 16,590.720
#> Bayesian Inf. Crit. 16,924.490 16,666.460
#> F Statistic 2,865.391*** (df = 7; 7222)
#> =====================================================================
#> Note: *p<0.1; **p<0.05; ***p<0.01
经过研究 Robert Long's great answer on stats.stackexchange 我发现模型的正确规格是嵌套设计,即 (1| schoolid /childid)
。然而,由于数据的编码方式(unniqe childid
在 schoolid
中),交叉设计或规范,即我上面使用的 (1|childid)+(1|schoolid)
,产生相同的结果。
这是使用与上述相同数据的插图,
data(egsingle, package = "mlmRev")
dta <- egsingle
dta$Female <- with(dta, ifelse(female == 'Female', 1, 0))
require(lme4)
dta$year <- as.factor(dta$year)
重新运行交叉的 design-model, , reg.re.lmer.in.school
, 以进行比较
reg.re.lmer.in.school <- lmer(math~Female+size+year+(1|childid)+(1|schoolid), dta)
这里是嵌套结构
reg.re.lmer.nested <- lmer(math~Female+size+year+(1| schoolid /childid), dta)
最后使用惊人的 texreg 包比较两个模型,
# install.packages(c("texreg"), dependencies = TRUE)
# require(texreg)
texreg::screenreg(list(reg.re.lmer.in.school, reg.re.lmer.nested), digits = 3)
#> ===============================================================
#> Model 1 Model 2
#> ---------------------------------------------------------------
#> (Intercept) -2.693 *** -2.693 ***
#> (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 16590.715
#> BIC 16666.461 16666.461
#> Log Likelihood -8284.357 -8284.357
#> Num. obs. 7230 7230
#> Num. groups: childid 1721
#> Num. groups: schoolid 60 60
#> Var: childid (Intercept) 0.672
#> Var: schoolid (Intercept) 0.180 0.180
#> Var: Residual 0.334 0.334
#> Num. groups: childid:schoolid 1721
#> Var: childid:schoolid (Intercept) 0.672
#> ===============================================================
#> *** p < 0.001, ** p < 0.01, * p < 0.05
我通常在广义最小二乘框架内工作,估计 Wooldridge 的 Introductory (2013) 调用的随机效应和固定效应模型对由个人和时间索引的纵向数据的影响维度.
我一直在使用plm包中plm()
中的可行GLS估计来估计随机效应模型——一些统计数据文学术语混合模型。 plm()
函数采用 index
参数,我在其中指示个人和时间索引。但是,我现在面对一些数据,其中每个人在每个时间点都有多个度量,即多好的分组结构。
我发现可以使用 lme4 包中的 lmer()
来建立这样的模型,但是我对行话的差异以及可能性感到有点困惑框架,我想知道是否正确指定了模型。由于我不熟悉框架和这个术语,我担心我可能会忽略更多内容。
我可以使用 lmer()
复制我常用的 plm()
模型,但我不确定如何添加分组。我试图在下面说明我的问题。
我找了一些和我的数据有点像的数据来说明我的情况。首先是一些需要的包,
install.packages(c("mlmRev", "plm", "lme4", "stargazer"), dependencies = TRUE)
然后是数据
data(egsingle, package = "mlmRev")
egsingle
是一个不平衡的小组,由 1721 名学童组成,分为 60 所学校,跨越五个时间点。这些数据最初与 HLM 软件包(Bryk、Raudenbush 和 Congdon,1996)一起分发,但可以在 mlmrev 软件包中找到,详情请参阅 ? 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
以下是我如何使用 plm()
、
schoolid
library(plm)
reg.re.plm <- plm(math~Female+size+year, dta, index = c("childid", "year"), model="random")
# summary(reg.re.plm)
我可以像这样重现这些结果lme4
require(lme4)
dta$year <- as.factor(dta$year)
reg.re.lmer <- lmer(math~Female+size+year+(1|childid), dta)
# summary(reg.re.lmer)
现在,从阅读 Bates (2010) 的第 2 章“lme4:混合效应建模
使用 R” 我相信这就是我指定模型的方式,包括集群级别,schoolid
,
reg.re.lmer.in.school <- lmer(math~Female+size+year+(1|childid)+(1|schoolid), dta)
# summary(reg.re.lmer.in.school)
但是,当我查看结果时,我不太相信我实际上已经正确指定了它(见下文)。
在我的实际数据中,重复的措施是在个人内部,但我认为我可以以这个数据为例。对于如何进行的任何建议,我将不胜感激。也许是对 notation/terminology 的工作示例的引用,与 Wooldridge (2013) 中使用的示例相差不远。而且,我如何逆向工作并编写 reg.re.lmer.in.school
模型的规范?
# library(stargazer)
stargazer::stargazer(reg.re.plm, reg.re.lmer, reg.re.lmer.in.school, type="text")
#> =====================================================================
#> Dependent variable:
#> -------------------------------------------------
#> math
#> panel linear
#> linear mixed-effects
#> (1) (2) (3)
#> ---------------------------------------------------------------------
#> Female -0.025 -0.025 0.008
#> (0.046) (0.047) (0.042)
#>
#> size -0.0004*** -0.0004*** -0.0003
#> (0.0001) (0.0001) (0.0002)
#>
#> year-1.5 0.878*** 0.876*** 0.866***
#> (0.059) (0.059) (0.059)
#>
#> year-0.5 1.882*** 1.880*** 1.870***
#> (0.059) (0.058) (0.058)
#>
#> year0.5 2.575*** 2.574*** 2.562***
#> (0.059) (0.059) (0.059)
#>
#> year1.5 3.149*** 3.147*** 3.133***
#> (0.060) (0.059) (0.059)
#>
#> year2.5 3.956*** 3.954*** 3.939***
#> (0.060) (0.060) (0.060)
#>
#> Constant -2.671*** -2.669*** -2.693***
#> (0.085) (0.086) (0.152)
#>
#> ---------------------------------------------------------------------
#> Observations 7,230 7,230 7,230
#> R2 0.735
#> Adjusted R2 0.735
#> Log Likelihood -8,417.815 -8,284.357
#> Akaike Inf. Crit. 16,855.630 16,590.720
#> Bayesian Inf. Crit. 16,924.490 16,666.460
#> F Statistic 2,865.391*** (df = 7; 7222)
#> =====================================================================
#> Note: *p<0.1; **p<0.05; ***p<0.01
经过研究 Robert Long's great answer on stats.stackexchange 我发现模型的正确规格是嵌套设计,即 (1| schoolid /childid)
。然而,由于数据的编码方式(unniqe childid
在 schoolid
中),交叉设计或规范,即我上面使用的 (1|childid)+(1|schoolid)
,产生相同的结果。
这是使用与上述相同数据的插图,
data(egsingle, package = "mlmRev")
dta <- egsingle
dta$Female <- with(dta, ifelse(female == 'Female', 1, 0))
require(lme4)
dta$year <- as.factor(dta$year)
重新运行交叉的 design-model, , reg.re.lmer.in.school
, 以进行比较
reg.re.lmer.in.school <- lmer(math~Female+size+year+(1|childid)+(1|schoolid), dta)
这里是嵌套结构
reg.re.lmer.nested <- lmer(math~Female+size+year+(1| schoolid /childid), dta)
最后使用惊人的 texreg 包比较两个模型,
# install.packages(c("texreg"), dependencies = TRUE)
# require(texreg)
texreg::screenreg(list(reg.re.lmer.in.school, reg.re.lmer.nested), digits = 3)
#> ===============================================================
#> Model 1 Model 2
#> ---------------------------------------------------------------
#> (Intercept) -2.693 *** -2.693 ***
#> (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 16590.715
#> BIC 16666.461 16666.461
#> Log Likelihood -8284.357 -8284.357
#> Num. obs. 7230 7230
#> Num. groups: childid 1721
#> Num. groups: schoolid 60 60
#> Var: childid (Intercept) 0.672
#> Var: schoolid (Intercept) 0.180 0.180
#> Var: Residual 0.334 0.334
#> Num. groups: childid:schoolid 1721
#> Var: childid:schoolid (Intercept) 0.672
#> ===============================================================
#> *** p < 0.001, ** p < 0.01, * p < 0.05