我对估计具有两个随机分量的混合效应模型感兴趣(我很抱歉符号有些不精确。我对这类模型有些陌生)。最后，我还想要随机分量方差的标准误差。这就是为什么我对使用包lme有些不屑一顾的原因。原因是我发现了关于如何计算这些标准误差的描述，而且有趣的是，这些方差函数的标准误差链接。
我相信我知道如何使用包lmer。我终于对model2感兴趣了.对于model1，两个命令产生相同的估计值。但是，model2lme与model2lme4包中的lmer产生不同的结果。你能帮我解决如何为lme设置随机组件吗？这将不胜感激。谢谢。请随函附上我的MWE。
最好
丹尼尔
#### load all packages #####
loadpackage <- function(x){
for( i in x ){
# require returns TRUE invisibly if it was able to load package
if( ! require( i , character.only = TRUE ) ){
# If package was not able to be loaded then re-install
install.packages( i , dependencies = TRUE )
}
# Load package (after installing)
library( i , character.only = TRUE )
}
}
# Then try/install packages...
loadpackage( c("nlme", "msm", "lmeInfo", "lme4"))
alcohol1 <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/alcohol1_pp.txt", header=T, sep=",")
attach(alcohol1)
id <- as.factor(id)
age <- as.factor(age)
model1.lmer <-lmer(alcuse ~ 1 + peer + (1|id))
summary(model1.lmer)
model2.lmer <-lmer(alcuse ~ 1 + peer + (1|id) + (1|age))
summary(model2.lmer)
model1.lme <- lme(alcuse ~ 1+ peer, data = alcohol1, random = ~ 1 |id, method ="REML")
summary(model1.lme)
model2.lme <- lme(alcuse ~ 1+ peer, data = alcohol1, random = ~ 1 |id + 1|age, method ="REML")
编辑(15/09/2021)：
按如下方式估计模型，然后通过nlme::VarCorr返回估计值，会给我不同的结果。虽然估计值似乎在球场上，但这是因为它们在组件之间切换。
model2a.lme <- lme(alcuse ~ 1+ peer, data = alcohol1, random = ~ 1 |id/age, method ="REML")
summary(model2a.lme)
nlme::VarCorr(model2a.lme)
Variance StdDev
id = pdLogChol(1)
(Intercept) 0.38390274 0.6195989
age = pdLogChol(1)
(Intercept) 0.47892113 0.6920413
Residual 0.08282585 0.2877948
编辑(16/09/2021)：
由于 Bob 促使我更多地考虑我的模型，我想提供一些额外的信息。请注意，我在 MWE 中使用的数据与我的真实数据不符。我只是将其用于说明目的，因为我无法上传 myy 真实数据。我有一个家庭面板，里面有收入、人口统计信息和家长指标。
我对代际流动感兴趣。永久收入的兄弟姐妹相关性是一个行业标准。至少，同时期的观察是非常糟糕的永久收入代理。由于暂时性冲击，即经典测量误差，这些估计肯定会衰减。出于这个原因，我们利用了数据的纵向维度。
对于兄弟姐妹相关性，这相当于假设收入过程如下：
$$Y_{ijt} = \beta X_{ijt} + \epsilon_{ijt}.$$
Y是个人i在t年来自家庭j的收入，X包括年龄和调查年份指标，以考虑调查年份的生命周期效应和宏观条件。Epsilon 是一个包含随机个体和家庭分量以及瞬态分量(测量误差和短寿命冲击)的复合项。它看起来如下：
$$\epsilon_{ijt} = \alpha_i + \gamma_j + \eta_{ijt}.$$
收入差异为：
$$\sigma^2_\epsilon = \sigma^2_\alpha + \sigma^2\gamma + \sigma^2\eta.$$
我们感兴趣的数量是
$$\rho = \frac(\sigma^2\gamma}{\sigma^2_\alpha + \sigma^2\gamma}，$$
这反映了兄弟姐妹之间共享家庭(和其他特征)在永久收入变化中所占的比例。
B.t.w.：挣扎只是因为我想对所有估计和\rho有一个标准误差。
这是交叉与嵌套随机效果的示例。(请注意，您引用的示例拟合了一种不同类型的模型，一个随机斜率模型，而不是具有两个不同分组变量的模型......
如果你尝试with(alcohol1, table(age,id))你可以看到每个id都与每个可能的年龄(14,15,16)相关联。或subset(alcohol1, id==1)例如：
id age coa male age_14 alcuse peer cpeer ccoa
1 1 14 1 0 0 1.732051 1.264911 0.2469111 0.549
2 1 15 1 0 1 2.000000 1.264911 0.2469111 0.549
3 1 16 1 0 2 2.000000 1.264911 0.2469111 0.549
对于具有age(按i索引)和id(按j索引)随机效应的模型，有三种可能的模型可以拟合
- 交叉(
(1|age) + (1|id))：Y_{ij} = beta0 + beta1*peer + eps1_i + eps2_j +epsr_{ij}; 酒精使用因人而异，并且独立地因年龄而异(此模型不会很好地工作，因为数据集中只有三个不同的年龄，通常需要更多级别) - 嵌套在年龄内的ID(
(1|age/id)=(1|age) + (1|age:id))：Y_{ij} = beta0 + beta1*peer + eps1_i + eps2_{ij} + epsr_{ij};酒精使用因年龄而异，并且在年龄内因人而异(请参阅上面关于级别数的说明)。 - 嵌套在ID中的年龄(
(1|id/age)=(1|id) + (1|age:id))：Y_{ij} = beta0 + beta1*peer + eps1_j + eps2_{ij} + epsr_{ij};酒精使用因人而异，并且因个人内部的年龄而异
这里eps1_i、eps2_{ij}和epsr_{ij}是正常的恶魔;epsr是残差项。
后两个模型在这种情况下实际上没有意义;因为每个年龄/ID组合只有一个观测值，嵌套方差(eps2)与残差方差(epsr)完全混淆。lme没有注意到这一点;如果您尝试在lmer中适合其中一个嵌套模型，它将给出一个错误
每个分组因子的水平数必须<观测值数(问题：ID：AGE)>
(尽管如果您尝试根据model1.lme计算置信区间，则会收到错误"无法获取 var-cov 成分的置信区间：非正确定近似方差-协方差"，这表明存在问题。
你可以把这个问题重述为残差变异，以及个体内年龄之间的变异，是共同无法识别的(在统计上不能彼此分离)。
此处更新的答案显示了如何从lmer模型中获取方差分量的标准误差，因此您不应该被lme所困扰(但您应该仔细考虑您真正尝试拟合的模型......
GLMM常见问题解答也可能有用。
更一般地说，标准误差
rho = (V_gamma)/(V_alpha + V_gamma)
将很难准确计算，因为这是模型参数的非线性函数。 您可以应用delta 方法，但最可靠的方法是使用参数化自举：如果您有一个拟合的模型m，那么这样的东西应该可以工作：
var_ratio <- function(m) {
v <- as.data.frame(sapply(VarCorr(m), as.numeric))
return(v$family/(v$family + v$id))
}
confint(m, method="boot", FUN =var_ratio)
您应该使用/而不是+来指定lme中的随机效果
通过lmer
model2.lmer <-lmer(alcuse ~ 1 + peer + (1|id) + (1|age), data = alcohol1)
summary(model2.lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: alcuse ~ 1 + peer + (1 | id) + (1 | age)
Data: alcohol1
REML criterion at convergence: 651.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.0228 -0.5310 -0.1329 0.5854 3.1545
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.08078 0.2842
age (Intercept) 0.30313 0.5506
Residual 0.56175 0.7495
Number of obs: 246, groups: id, 82; age, 82
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.3039 0.1438 2.113
peer 0.6074 0.1151 5.276
Correlation of Fixed Effects:
(Intr)
peer -0.814
按lme
model2.lme <- lme(alcuse ~ 1+ peer, data = alcohol1, random = ~ 1 |id/age, method ="REML")
summary(model2.lme)
Linear mixed-effects model fit by REML
Data: alcohol1
AIC BIC logLik
661.3109 678.7967 -325.6554
Random effects:
Formula: ~1 | id
(Intercept)
StdDev: 0.4381206
Formula: ~1 | age %in% id
(Intercept) Residual
StdDev: 0.4381203 0.7494988
Fixed effects: alcuse ~ 1 + peer
Value Std.Error DF t-value p-value
(Intercept) 0.3038946 0.1438333 164 2.112825 0.0361
peer 0.6073948 0.1151228 80 5.276060 0.0000
Correlation:
(Intr)
peer -0.814
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.0227793 -0.5309669 -0.1329302 0.5853768 3.1544873
Number of Observations: 246
Number of Groups:
id age %in% id
82 82
好的，最后。只是为了勾勒我的机密数据：我有一个个人小组。数据包括通过mnr识别的兄弟姐妹。收入是收入，波浪调查年份，年龄因素。女性是性别的一个因素，PID是识别个体的因素。
m1 <- lmer(income ~ age + wavey + female + (1|pid) + (1 | mnr),
data = panel)
vv <- vcov(m1, full = TRUE)
covvar <- vv[58:60, 58:60]
covvar
3 x 3 Matrix of class "dgeMatrix"
cov_pid.(Intercept) cov_mnr.(Intercept) residual
[1,] 2.6528679 -1.4624588 -0.4077576
[2,] -1.4624588 3.1015001 -0.0597926
[3,] -0.4077576 -0.0597926 1.1634680
mean <- as.data.frame(VarCorr(m1))$vcov
mean
[1] 17.92341 16.86084 56.77185
deltamethod(~ x2/(x1+x2), mean, covvar, ses =TRUE)
[1] 0.04242089
最后一个标量应该是我解释为永久收入兄弟姐妹的共同背景。
感谢@Ben Bolker为我指明了这个方向。