如何在R的nlme::lme中指定固定的单位方差-协方差矩阵？

核心思路 #

要得到固定为单位矩阵的残差VCV，需要同时满足三个约束：

• 组内残差的相关矩阵为单位矩阵（非对角元素全为0）
• 每个测量点的残差方差固定为1
• 禁用模型默认的总体残差尺度参数估计，避免方差被额外缩放

完整实现代码 #

第一步：生成测试数据（符合描述的结构） #

library(nlme)
set.seed(123)
# 构造20个受试者、3个重复测量水平的长表
df <- expand.grid(S = factor(1:20), X = c("a", "b", "c"))
df$Y <- rnorm(nrow(df))
# 新增全1列，用于指定固定方差
df$ones <- 1

第二步：拟合固定单位VCV的模型 #

model <- lme(
  fixed = Y ~ X,
  random = ~1|S,
  data = df,
  # 1. 指定相关结构为单位矩阵，所有非对角相关为0
  correlation = corIdent(form = ~1|S),
  # 2. 固定残差方差为全1
  weights = varFixed(~ones),
  # 3. 强制残差尺度参数为1，不参与优化
  control = lmeControl(sigma = 1, fixedSigma = TRUE)
)

第三步：查看结果 #

getVarCov(model, type = "conditional")

输出结果

Conditional variance covariance matrix
     1    2    3
1 1.0000 0.0000 0.0000
2 0.0000 1.0000 0.0000
3 0.0000 0.0000 1.0000

关键参数说明 #

• corIdent：属于nlme内置的相关结构类，直接指定组内相关矩阵为单位矩阵，无待估参数，不需要手动设置fixed约束。
• varFixed(~ones)：指定残差方差与输入变量的取值成正比，这里输入变量为全1列，因此残差方差的基值固定为1。
• lmeControl(sigma = 1, fixedSigma = TRUE)：nlme默认会估计一个全局残差缩放参数σ，设置该参数后会强制σ为1且不参与模型优化，避免方差被额外缩放。

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