在R的glmmTMB/sommer中构建含固定方差成分的混合模型

一、glmmTMB实现 #

关键调整细节 #

1. 观测方差的正确传递：glmmTMB里高斯族的离散参数对应的是方差的对数，所以得把给定的标准差SD先平方转成方差，再取对数作为初始值——你之前用的是log(SD)，这是错的，得改成log(SD²)。
2. 随机效应的简洁写法：A + A:B + A:C可以简化成A*(B + C)，功能完全一致，代码更清爽。
3. 收敛失败的核心原因：你的示例数据集只有20个观测，但这个随机效应结构要估计一个8×8的协方差矩阵（光这就36个参数），参数数量远超样本量，肯定收敛不了——不是模型指定错了，是数据撑不起这么复杂的结构。

修正后的代码 #

# 构造示例数据
data <- structure(list(A = c("H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J"),
                       B = c("i", "i", "k", "k", "k", "i", "i", "i", "k", "i", "k", "i", "i", "k", "k", "i", "i", "k", "k", "k"),
                       C = c(5L, 5L,7L, 10L, 2L, 6L, 5L, 4L, 2L, 3L, 10L, 4L, 3L, 9L, 8L, 5L,1L, 10L, 3L, 9L),
                       X = c("x1", "x2", "x3", "x4", "x5", "x1","x2", "x3", "x6", "x7", "x3", "x3", "x1", "x8", "x9", "x9","x10", "x10", "x11", "x12"),
                       Y = c(7.7, 2.9, 24, 3.6, 7.8,7.1, 73, 5.7, 18, 5.2, 43, 4.7, 3.4, 12, 5.8, 88, 6.9, 9.4, 1.1, 31),
                       SD = c(2.566666667, 0.966666667, 8, 1.2, 1.95, 2.366666667, 18.25, 2.85, 9, 1.3, 10.75, 1.566666667, 1.133333333,4, 1.45, 44, 1.725, 2.35, 0.275, 15.5)), 
                  row.names = c(NA,-20L), class ="data.frame")

# 计算观测方差V=SD²，生成观测ID
data$V <- data$SD^2
data$obs <- factor(seq(nrow(data)))

# 拟合模型
library(glmmTMB)
model_glmmTMB <- glmmTMB(
  log(Y) ~ A + A:B + A:C + (A*(B + C) | X),  # 和原随机效应结构等价，写法更简洁
  family = gaussian,
  dispformula = ~ 0 + obs,  # 给每个观测分配独立的离散参数（即固定方差）
  start = list(betad = log(data$V)),  # 初始值设为方差的对数
  map = list(betad = factor(rep(NA, nrow(data)))),  # 固定这些离散参数，不做估计
  data = data
)

# 查看结果（示例数据因参数过多大概率还是不收敛，实际数据集够大的话没问题）
summary(model_glmmTMB)

二、sommer实现 #

sommer可以用加权最小二乘的方式直接固定观测方差，权重设为1/V就行，同时能方便地指定随机效应结构。

代码示例 #

library(sommer)

# 计算权重w=1/V，权重对应1/观测方差，实现固定残差方差
data$w <- 1/data$V

# 拟合模型
model_sommer <- mmer(
  log(Y) ~ A + A:B + A:C,
  random = ~ vs(A*(B + C), grp = X),  # 指定试验X和处理的交互随机效应
  weights = w,
  data = data
)

# 查看结果
summary(model_sommer)

注意点 #

• sommer的vs()函数用来定义方差结构，vs(A*(B + C), grp=X)和你要的A + A:B + A:C随机效应完全一致，每个试验X的水平对应一组随机参数。
• 用weights=1/V本质就是让每个观测的残差方差固定为V，刚好符合元分析的需求。

三、解决收敛问题的实用建议 #

1. 简化随机效应结构：如果实际数据集样本量不算大，别搞这么复杂的随机效应，比如先保留(A | X)和(A:C | X)，或者用双竖线(A + A:B + A:C || X)让随机效应间独立，这样要估计的协方差参数会少很多。
2. 凑更多样本：元分析嘛，尽量多找些试验/观测数据，样本量上去了，复杂模型才跑得动。
3. 调整初始值：可以给随机效应的方差参数设个合理的初始值（比如用start参数指定theta的初始值），帮模型找收敛方向。

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