如何在nlme中拟合组间方差不同的个体水平随机效应模型？

数据老陈

2026-8-13

问题描述

我有一份个体测量数据集，每个个体隶属于A或B组，测量类型分为1和2两类，部分个体同时具有两种类型的测量值。数据中不同测量类型的方差存在差异，且个体水平随机效应在两组中的方差差异显著。我已使用nlme拟合了如下混合效应模型：

init <- c(-1.2, 120, 2, 100)

model1 <- nlme(y ~ a,
               data = dat,
               fixed = list(a ~ group : measure + 0),
               random = a ~ 1,
               groups = ~ subject,
               start = init,
               weights = varIdent(form = ~ 1 | measure))

但我希望让随机效应在不同组中具有不同的方差，推测可通过相关结构实现但未成功。由于实际模型为非线性且更为复杂，无法使用lmer的交叉随机效应解决，请问该如何正确拟合这类模型？

示例数据集如下：

dat <- structure(list(subject = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 
                       12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 
                       28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 
                       44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 
                       60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 
                       76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 
                       92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 101, 102, 102, 103, 
                       103, 104, 104, 105, 105, 106, 106, 107, 107, 108, 108, 109, 109, 
                       110, 110, 111, 111, 112, 112, 113, 113, 114, 114, 115, 115, 116, 
                       116, 117, 117, 118, 118, 119, 119, 120, 120, 121, 121, 122, 122, 
                       123, 123, 124, 124, 125, 125, 126, 126, 127, 127, 128, 128, 129, 
                       129, 130, 130), 
           group = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 
                               2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
                               2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
                               2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
                               2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
                               2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
                               2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 
                               2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 
                               2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 
                               2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 
                               2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L), .Label = c("A", "B"), class = "factor"), 
           measure = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
                                 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
                                 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
                                 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
                                 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
                                 2L), .Label = c("1", "2"), class = "factor"), 
           y = c(-1.71, 
                 121.74, -1.57, 109.96, -0.64, 101.67, -0.13, 120.64, 1.47, 
                 101.99, -4.51, 133.18, -2.9, 117.95, -0.97, 126.94, -1.44, 
                 105.1, -1.52, 122.2, -2.29, 130.17, -0.35, 133.14, -0.94, 
                 112.68, -0.89, 105.37, -2.49, 126.75, -2.61, 139.25, -2.13, 
                 113.43, 0.61, 140.76, -0.75, 129.17, 1.94, 139.4, -0.49, 
                 119.03, -2.09, 89.97, -2.76, 107.85, 1.61, 136.31, -0.55, 
                 128.6, 0.41, 86.66, 0.54, 100.03, 2.46, 115.37, 6.94, 109.34, 
                 3.78, 102.34, -4.46, 104.06, 1.48, 105.06, 3.98, 85.21, 1.31, 
                 103.17, -3.35, 110.83, 2.75, 98.38, -2.43, 101.57, 2.2, 120.45, 
                 -4.06, 101.25, 3.85, 99.38, 2.17, 108, 9.27, 100.76, 3.27, 
                 110.3, 1.22, 98.91, 1.62, 105.65, 4.64, 113.07, 8.14, 108.75, 
                 6.84, 73.08, 1.42, 99.41, -0.5, 95.25, 1.42, 3.76, 102.95, 
                 85.45, -2.71, -0.48, 137.34, 114.61, -0.42, 1.71, 98.82, 
                 83.06, -3.51, -0.32, 109.66, 91.99, -0.46, -1.35, 113.88, 
                 97.32, -0.93, 1.17, 111.26, 103.9, -4.11, 6.78, 106.36, 88.22, 
                 -0.85, -6.56, 137.39, 112.19, -0.91, 3.26, 122.53, 105.18, 
                 -0.61, 4.25, 111.01, 95.85, -2.68, 3.1, 142.26, 114.44, -0.31, 
                 3.76, 127.61, 102.26, -1.82, 4.01, 116.61, 97.1, -3.61, 0.9, 
                 107.73, 90.6, -0.13, 3.78, 108.73, 95.12)), 
      row.names = c(NA, -160L), class = "data.frame")

解决方案

在nlme中，可以通过给随机效应添加varIdent方差结构，实现个体水平随机效应的方差随组（group）变化。具体做法是在random参数中指定varIdent(form = ~1 | group)，让模型估计不同组的随机效应方差比值，同时保留原有的残差异方差结构（随measure变化）。

注意需要调整初始值，新增一个参数对应varIdent的比例系数（初始设为1，代表两组随机效应方差初始相等）。

拟合代码

init <- c(-1.2, 120, 2, 100, 1)  # 新增第5个初始值给组间随机效应方差比例

model2 <- nlme(y ~ a,
               data = dat,
               fixed = list(a ~ group : measure + 0),
               # 指定随机效应方差随group变化
               random = a ~ varIdent(form = ~1 | group) | subject,
               weights = varIdent(form = ~1 | measure),
               start = init,
               # 可选：优化器设置，提升收敛稳定性
               control = nlmeControl(opt = "nlminb"))

模型解释

random = a ~ varIdent(form = ~1 | group) | subject：指定每个个体的随机截距方差由其所属的组决定，模型会估计参考组（默认是group的第一水平，即A组）的随机效应方差，以及另一组（B组）相对于参考组的方差比例。
weights = varIdent(form = ~1 | measure)：保留原模型中残差方差随测量类型变化的设置。
新增的初始值1表示初始时B组随机效应方差与A组相同，模型会根据数据调整这个比例。

验证结果

拟合完成后，可以通过summary(model2)查看随机效应的方差估计值，其中varIdent的参数会显示B组相对于A组的方差比例，结合参考组的方差即可得到两组各自的随机效应方差。

内容的提问来源于stack exchange，提问作者Lars Lau Raket

如何在nlme中拟合组间方差不同的个体水平随机效应模型？

拟合代码 #

模型解释 #

验证结果 #

拟合代码

模型解释

验证结果