伙计们,
我正在尝试使用 {nlme} 包中的 lme() 在 R 中拟合一个模型,该模型具有许多固定效果和附加随机效果(必须是 {nlme},因为我还想稍后包含一个 AR(1) 相关矩阵) . 但是,该模型非常慢(由于大量的固定效应)。
现在,我们知道运行带有假人的模型在数值上但不是在计算上等同于在去均值数据上运行模型(在估计范围内)。因此,以下两个模型会产生相同的结果,但第一个模型运行得更快:
# Load library
library(nlme) # For lme()
library(lfe) # For demeanlist()
# Create reproducible dummy dataset
set.seed(1)
dat <- data.frame(y=rnorm(1000), x=rnorm(1000), cluster1=sample(LETTERS, 1000, replace=TRUE), cluster2=sample(LETTERS, 1000, replace=TRUE))
# Create de-meaned dataset
dat_demean <- demeanlist(dat[,c("y", "x")], list(dat$cluster1))
dat_demean$cluster2 <- dat$cluster2 # Copy cluster2 column over to de-meaned data
# Model 1: Fixed effect on cluster1
orig <- lm( y ~ x + cluster1, data=dat)
deme <- lm( y ~ x - 1, data=dat_demean)
# This works as expected. Coefficients are exactly identical (SE are wrong and need to be fixed).
all.equal( coef(orig)["x"], coef(deme)["x"] )
R> [1] TRUE
当我们想在第二个分组变量上额外包含随机效应时,可以做类似的事情吗?以下是我的尝试,它不起作用。
# Model 2: Fixed effect on cluster1 and additional random effect for cluster2
orig <- lme( y ~ x + cluster1, random=~1|cluster2, data=dat)
deme <- lme( y ~ x -1, random=~1|cluster2, data=dat_demean)
# This does not work - coefficients are different
fixef(orig)["x"]
R> x
R> -0.005885881
fixef(deme)["x"]
R> x
R> -0.006169549
all.equal(fixef(orig)["x"], fixef(deme)["x"])
R> [1] "Mean relative difference: 0.04819478"
我错过了什么或我做错了什么?或者,将随机效应模型拟合到贬义数据上不是那么简单吗?
干杯