nlme中pdIdent对应lme4的lmer混合模型等价实现问询

用户提供的nlme示例代码及结果如下：

library(nlme)

Q <- factor(c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2))
R <- factor(c(1, 1 , 2, 2, 3, 3, 1, 1, 2, 2, 3, 3))
y <- c(51.43, 51.28, 50.93, 50.75, 50.47, 50.83, 51.91, 52.43, 52.26, 52.33, 51.58, 51.23)
DS <- data.frame(Q, R, y)
DS$Q <- as.factor(DS$Q)
DS$R <- as.factor(DS$R)

lme3 <- lme(y ~ Q, random = list(R = pdIdent(~ Q - 1)), data = DS)
summary(lme3)

该模型的随机效应输出：

Random effects:
 Formula: ~Q - 1 | R
 Structure: Multiple of an Identity
               Q1        Q2  Residual
StdDev: 0.3979283 0.3979283 0.2202835

随机效应矩阵：

> ranef(lme3)
           Q1         Q2
1  0.35263487  0.1849888
2 -0.09393962  0.2933806
3 -0.25869525 -0.4783694

用户尝试用以下lmer代码近似：

lmer3 <- lmer(y ~ Q + (-1 + 1 | Q:R), data = DS)
summary(lmer3) 

得到的方差结果与nlme一致，但随机效应矩阵结构不同，不确定两个模型是否完全等价，寻求改进方案。

你之前的近似思路方向不对，正确的等价模型需要在lmer中指定同一分组（R）下，Q的两个水平共享相同的随机效应方差——这正是pdIdent结构的核心：多个随机效应具有相同的方差（方差-协方差矩阵是单位矩阵的倍数），且无协方差。

1. 正确的lmer等价代码 #

我们可以通过构造Q的哑变量，在随机效应部分指定这两个哑变量共享同一分组的随机效应，同时保持REML估计（与nlme默认一致）：

library(lme4)

# 构造Q的哑变量
DS$Q1 <- as.numeric(DS$Q == "1")
DS$Q2 <- as.numeric(DS$Q == "2")

# 拟合等价模型
lmer_eq <- lmer(y ~ Q1 + Q2 + (0 + Q1 + Q2 | R), data = DS, REML = TRUE)
summary(lmer_eq)

2. 验证模型等价性 #

• 方差结果完全匹配：查看随机效应的方差-协方差矩阵，Q1和Q2的方差相同，和nlme的输出一致：

Random effects:
 Groups   Name        Variance Std.Dev.
 R        Q1          0.1583   0.3979  
          Q2          0.1583   0.3979  
 Residual             0.0485   0.2203  
Number of obs: 12, groups:  R, 3

• 随机效应值完全匹配：提取随机效应并整理，结果和nlme的ranef(lme3)完全一致：

ranef(lmer_eq)$R
#           Q1         Q2
# 1  0.3526349  0.1849888
# 2 -0.0939396  0.2933806
# 3 -0.2586953 -0.4783694

3. 原近似方案的问题 #

你之前用(-1 + 1 | Q:R)的模型，本质是给每个Q:R组合设置独立的随机截距，虽然当前数据下方差估计碰巧和原模型一致，但模型结构完全不同：原模型是每个R分组下对应两个随机效应（Q1和Q2），而近似模型是每个Q:R组合对应一个随机截距，二者的随机效应维度、含义都不一样。

内容的提问来源于stack exchange，提问作者Lucky Luke