如何在R的lmer函数中为男女分组设置独立残差方差？

我正在针对以下数据集运行线性混合效应模型：

maleID  femaleID  dyadID  session  wai  male  female 
55      100       55_100  1        4.75 1     0
55      100       55_100  2        5    1     0
55      100       55_100  1        3.25 0     1
55      100       55_100  2        4.25 0     1
55      101       55_101  1        4    1     0
55      101       55_101  2        6    1     0
55      101       55_101  1        5    0     1           
55      101       55_101  2        4    0     1
59      101       59_101  1        5.25 1     0
59      101       59_101  2        3.5  1     0
59      101       59_101  1        4.5  0     1
59      101       59_101  2        6.5  0     1
59      102       59_102  1        5    1     0
59      102       59_102  2        6    1     0
59      102       59_102  1        5.75 0     1
59      102       59_102  2        5.25 0     1

其中wai是因变量，male和female是虚拟变量，maleID、femaleID和dyadID为唯一标识。

我用以下R代码拟合混合模型：

lmer(wai ~ 0 + male + female + session + 
       (0 + male + female|maleID) + 
       (0 + male + female|femaleID),
     data=data, na.action=na.exclude, REML=TRUE)

male和female作为虚拟变量，已经让模型为男性评级和女性评级生成两个独立的随机效应，但模型仅输出单个残差项。我想为模型的误差/残差项做同样处理，即分别为男性评级和女性评级设置独立的残差方差。

在SPSS中，通过以下语法的最后一条REPEATED语句可实现该功能：

MIXED wai WITH male female BY maleID femaleID dyadID
    /FIXED = male female session | NOINT
    /METHOD=REML
    /PRINT=SOLUTION TESTCOV
    /RANDOM=male female | SUBJECT (dyadID) COVTYPE (csh)
    /RANDOM=male female | SUBJECT (maleID) COVTYPE (csh)
    /RANDOM=male female | SUBJECT (femaleID) COVTYPE (csh)
    /REPEATED=male female | SUBJECT (dyadID*maleID*femaleID*session) COVTYPE (vc)

注：通过使用“NOINT”命令和前置的male、female虚拟变量，我抑制了全局截距，从而获取男女分组的独立方差。请问在R中如何实现这一功能？

在R中实现男女分组的独立残差方差，有两种常用方法，分别基于nlme和glmmTMB包，具体实现如下：

方法1：使用nlme包的lme函数 #

nlme支持通过weights参数指定异方差结构，用varIdent为不同性别组设置独立残差方差：

1. 先创建性别分组变量：

data$gender <- ifelse(data$male == 1, "male", "female")

2. 拟合模型：

library(nlme)

model_lme <- lme(wai ~ 0 + male + female + session,
                 random = list(maleID = pdDiag(~0 + male + female),
                               femaleID = pdDiag(~0 + male + female)),
                 weights = varIdent(form = ~1 | gender),
                 data = data,
                 na.action = na.exclude,
                 method = "REML")

• pdDiag(~0 + male + female)指定随机效应为对角矩阵，对应原模型中独立的随机效应结构
• varIdent(form = ~1 | gender)指定残差方差随gender分组变化，即男性和女性观测拥有独立的残差方差

方法2：使用glmmTMB包 #

glmmTMB支持通过分散公式直接指定残差方差的分组结构，无需额外创建分组变量：

library(glmmTMB)

model_glmmtmb <- glmmTMB(wai ~ 0 + male + female + session +
                           (0 + male + female|maleID) +
                           (0 + male + female|femaleID),
                         dispformula = ~0 + male + female,
                         data = data,
                         REML = TRUE,
                         na.action = na.exclude)

• dispformula = ~0 + male + female指定分散参数（残差方差）由male和female虚拟变量预测且无截距，会为男性和女性分别估计独立的残差方差

结果验证 #

拟合完成后，可通过以下命令查看残差方差估计值：

• 对于nlme模型：运行summary(model_lme)，在Variance function部分会展示不同性别组的残差方差乘数
• 对于glmmTMB模型：运行summary(model_glmmtmb)，在Dispersion model部分会看到male和female对应的分散参数估计值，转换后即可得到分组残差方差

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