Two-level 在 R 中使用 lme 建模
Two-level modelling with lme in R
我对估计具有两个随机分量的混合效应模型很感兴趣（对于有些不精确的符号，我很抱歉。我对这类模型有些陌生）。最后，我还想要随机分量方差的标准误差。这就是为什么我有点愿意使用包 lme。原因是我发现了关于如何计算这些标准误差的描述，而且也很有趣，这些方差函数的标准误差 link.
我相信我知道如何使用包 lmer。我终于对model2感兴趣了。对于 model1，这两个命令都会产生相同的估计值。但是 model2 和 lme 产生的结果与 model2 和 lme4 包中的 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 中使用的数据与我的真实数据不符。我只是将它用于说明目的，因为我无法上传我的真实数据。我有一个包含收入、人口统计信息和 parent 指标的家庭面板。
我对代际流动很感兴趣。永久收入的兄弟姐妹相关性是一个行业标准。至少，同时期的观察是永久收入的非常糟糕的代表。由于短暂的冲击，即经典测量误差，这些估计肯定会减弱。出于这个原因，我们利用数据的纵向维度。
对于兄弟姐妹相关性，这相当于假设收入过程如下：
$$Y_{ijt} = \beta X_{ijt} + \epsilon_{ijt}.$$
其中 Y 是第 t 年来自家庭 j 的个人 i 的收入。 X 包括年龄和调查年份指标，以说明调查年份的 life-cycle 影响和宏观经济状况。 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.
有一个标准误差
您应该使用 / 而不是 +
在 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
这是 交叉 与 嵌套 随机效应的示例。 （请注意，您提到的示例适用于不同类型的模型，随机斜率模型而不是具有两个不同分组变量的模型...）
如果您尝试 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 嵌套在 age (
(1|age/id) = (1|age) + (1|age:id)):
Y_{ij} = beta0 + beta1*peer + eps1_i + eps2_{ij} + epsr_{ij};酒精使用因年龄而异，并且因年龄段内的个体而异（请参阅上面关于级别数的注释）。
- age 嵌套在 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 是残差项。
后两种模型在这种情况下实际上没有意义；因为每个 age/id 组合只有一个观察值，所以嵌套方差 (eps2) 与残差 (epsr) 完全混淆。 lme 没有注意到这一点；如果您尝试拟合 lmer 中的嵌套模型之一，它将给出一个错误
number of levels of each grouping factor must be < number of observations (problems: id:age)
（尽管如果您尝试根据 model1.lme 计算置信区间，您会收到错误“无法获得 var-cov 分量的置信区间：非正定近似方差-协方差”，这是提示有问题。）
您可以将此问题重述为残差和个体年龄间的差异共同无法识别（统计上不能彼此分开） .
更新答案展示了如何从lmer模型中获取方差分量的标准误差，所以你不应该被卡住使用 lme（但你 应该 仔细考虑你真正想要适应的模型......）
GLMM FAQ 也可能有用。
更一般地说，
的标准误差
rho = (V_gamma)/(V_alpha + V_gamma)
将很难准确计算，因为这是模型参数的非线性函数。您可以应用 delta 方法，但最可靠的方法是使用 parametric bootstrapping：如果您有拟合模型 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)
好吧，终于。只是勾勒出我的机密数据：我有一个小组。数据包括通过 mnr 识别的兄弟姐妹。 income 是收入，wavey 调查年份，年龄年龄因素。 female 是性别的一个因素，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 为我指出了这个方向。
我对估计具有两个随机分量的混合效应模型很感兴趣（对于有些不精确的符号，我很抱歉。我对这类模型有些陌生）。最后，我还想要随机分量方差的标准误差。这就是为什么我有点愿意使用包 lme。原因是我发现了关于如何计算这些标准误差的描述，而且也很有趣，这些方差函数的标准误差 link.
我相信我知道如何使用包 lmer。我终于对model2感兴趣了。对于 model1，这两个命令都会产生相同的估计值。但是 model2 和 lme 产生的结果与 model2 和 lme4 包中的 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 中使用的数据与我的真实数据不符。我只是将它用于说明目的，因为我无法上传我的真实数据。我有一个包含收入、人口统计信息和 parent 指标的家庭面板。
我对代际流动很感兴趣。永久收入的兄弟姐妹相关性是一个行业标准。至少，同时期的观察是永久收入的非常糟糕的代表。由于短暂的冲击，即经典测量误差，这些估计肯定会减弱。出于这个原因，我们利用数据的纵向维度。
对于兄弟姐妹相关性，这相当于假设收入过程如下：
$$Y_{ijt} = \beta X_{ijt} + \epsilon_{ijt}.$$
其中 Y 是第 t 年来自家庭 j 的个人 i 的收入。 X 包括年龄和调查年份指标，以说明调查年份的 life-cycle 影响和宏观经济状况。 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.
有一个标准误差您应该使用 / 而不是 +
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
这是 交叉 与 嵌套 随机效应的示例。 （请注意，您提到的示例适用于不同类型的模型，随机斜率模型而不是具有两个不同分组变量的模型...）
如果您尝试 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 嵌套在 age (
(1|age/id)=(1|age) + (1|age:id)):Y_{ij} = beta0 + beta1*peer + eps1_i + eps2_{ij} + epsr_{ij};酒精使用因年龄而异，并且因年龄段内的个体而异（请参阅上面关于级别数的注释）。 - age 嵌套在 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 是残差项。
后两种模型在这种情况下实际上没有意义；因为每个 age/id 组合只有一个观察值，所以嵌套方差 (eps2) 与残差 (epsr) 完全混淆。 lme 没有注意到这一点；如果您尝试拟合 lmer 中的嵌套模型之一，它将给出一个错误
number of levels of each grouping factor must be < number of observations (problems: id:age)
（尽管如果您尝试根据 model1.lme 计算置信区间，您会收到错误“无法获得 var-cov 分量的置信区间：非正定近似方差-协方差”，这是提示有问题。）
您可以将此问题重述为残差和个体年龄间的差异共同无法识别（统计上不能彼此分开） .
更新答案lmer模型中获取方差分量的标准误差，所以你不应该被卡住使用 lme（但你 应该 仔细考虑你真正想要适应的模型......）
GLMM FAQ 也可能有用。
更一般地说，
的标准误差rho = (V_gamma)/(V_alpha + V_gamma)
将很难准确计算，因为这是模型参数的非线性函数。您可以应用 delta 方法，但最可靠的方法是使用 parametric bootstrapping：如果您有拟合模型 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)
好吧，终于。只是勾勒出我的机密数据：我有一个小组。数据包括通过 mnr 识别的兄弟姐妹。 income 是收入，wavey 调查年份，年龄年龄因素。 female 是性别的一个因素，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 为我指出了这个方向。