基于nlme的三层混合模型：Level3随机斜率设置问题

问题背景 #

数据为三层嵌套结构：选举年（Level3，5个类别）→选区（Level2，299个选区，共1495个年度-选区组合）→候选人（Level1，共5800个案例）。核心需求：

• 预测变量x1仅在Level3设置随机斜率，Level2不设置随机斜率
• 保留Level3的随机截距
• 后续需引入x1与Level3的交互项

此前尝试的问题 #

1. 使用~x1-1|Level3时，会移除Level3的随机截距，不符合需求
2. 将-1替换为+1（即~1+x1|Level3）时，触发报错：

   Error in logLik.lmeStructInt(lmeSt, lmePars) : 
     NA/NaN/Inf in foreign function call (arg 3)
   

   伴随大量警告，本质原因是Level3仅5个组，同时估计随机截距、x1的随机斜率及两者的协方差，组数量过少导致方差协方差矩阵估计奇异。

解决方案 #

方案1：改用lme4包的lmer函数 #

lme4对小样本组的随机效应估计更鲁棒，无需手动指定协方差结构即可处理：

library(lme4)
# 基础模型（保留Level3随机截距+随机斜率，Level2仅随机截距）
m <- lmer(y ~ x1 + x2 + xn + (1 + x1 | Level3) + (1 | Level2), 
          data = data, subset = x2 == 0, REML = TRUE, na.action = na.omit)
summary(m)

# 引入x1与Level3交互项的模型（Level3为分类变量时，交互项为固定效应）
m_interact <- lmer(y ~ x1 + x2 + xn + x1:Level3 + (1 + x1 | Level3) + (1 | Level2), 
                   data = data, subset = x2 == 0, REML = TRUE, na.action = na.omit)
summary(m_interact)

方案2：在nlme中指定对角协方差结构 #

若坚持使用nlme，可通过pdDiag()指定Level3的随机效应协方差为对角矩阵（强制截距与x1随机斜率的协方差为0），减少待估计参数，避免奇异问题：

library(nlme)
# 基础模型
m <- lme(y ~ x1 + x2 + xn,
         random = list(Level3 = pdDiag(~1 + x1), Level2 = ~1),
         data = data, subset = x2 == 0, method = "REML", na.action = na.omit)
summary(m)

# 引入交互项的模型
m_interact <- lme(y ~ x1 + x2 + xn + x1:Level3,
                  random = list(Level3 = pdDiag(~1 + x1), Level2 = ~1),
                  data = data, subset = x2 == 0, method = "REML", na.action = na.omit)
summary(m_interact)

补充说明 #

若引入x1:Level3交互项（固定效应），需注意：该交互项已直接建模不同选举年x1的固定差异，此时Level3的x1随机斜率捕捉的是年度间x1效应的额外随机变异。若理论上认为年度间x1效应的差异是固定的，可考虑移除Level3的x1随机斜率，仅保留随机截距。

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