如何用R的nlme包拟合含多随机效应的重复测量时间序列混合模型

需求说明 #

• 用nlme包拟合带异质方差结构的时间序列混合模型，因变量为Yield
• 固定效应：Genotype + Table + Genotype×Week + Table×Week + Akima（Akima为协变量）
• 随机效应：Col + Row + Col×Week + Row×Week
• 重复测量要求：
  1. 以Plant为个体，按Week重复测量（需指定无结构方差协方差矩阵）
  2. 以Col×Row×Block为个体，按Week重复测量，固定和随机效应与前者一致

现有成功模型（仅固定+随机效应） #

已完成不含重复测量结构的模型拟合，代码如下：

data$Plant <- as.factor(data$Plant)
data$Week <- as.factor(data$Week)
data$Table <- as.factor(data$Table)
data$Block <- as.factor(data$Block)
data$Row <- as.factor(data$Row)
data$Col <- as.factor(data$Col)
data$Genotype <- as.factor(data$Genotype)

data$br <- with(data, Row:Week)
data$bc <- with(data, Col:Week)

# 仅含固定与随机效应的模型
model1 <- lme(Yield ~ Genotype + Table + Genotype:Week  + Table:Week + Akima,
              random = list(Col=~ 1, Row=~ 1, br=~ 1, bc=~ 1),
              data = data)

报错的重复测量模型及错误信息 #

尝试添加重复测量结构时出现错误，错误代码与提示如下：

model2 <- lme(Yield ~ Genotype + Table + Genotype:Week  + Table:Week + Akima,
              random = list(Col=~ 1, Row=~ 1, br=~ 1, bc=~ 1, Plant = pdDiag(~Week)),
              data = data)

错误提示：

Error in lme.formula(Yield ~ Genotype + Table + Genotype:Week +  : 
  fewer observations than random effects in all level 5 groups

问题分析与解决方案 #

错误原因 #

报错核心是部分Plant个体的Week观测数量少于无结构协方差矩阵所需估计的参数数。无结构矩阵的参数数量为k(k+1)/2（k为Week的水平数），若某个Plant的观测Week数不足，模型无法估计足够参数，进而报错。另外，之前的代码错误地将重复测量结构放在random参数中，pdDiag(~Week)是定义随机效应的方差结构，而非重复测量的协方差结构。

修正步骤 #

1. 筛选完整观测数据
   先检查每个Plant的Week覆盖情况，筛选出拥有全部Week观测的个体：

# 查看每个Plant的Week观测分布
table(data$Plant, data$Week)
# 筛选出覆盖所有Week的Plant
total_weeks <- length(unique(data$Week))
complete_plants <- names(which(table(data$Plant) == total_weeks))
data_complete <- subset(data, Plant %in% complete_plants)

2. 正确指定重复测量的无结构协方差结构
   在lme中，重复测量的协方差结构需通过correlation参数指定，使用corSymm()对应无结构矩阵，同时保留原有随机效应：

model2 <- lme(Yield ~ Genotype + Table + Genotype:Week + Table:Week + Akima,
              random = list(Col=~1, Row=~1, br=~1, bc=~1),
              # 指定Plant为个体，Week间采用无结构协方差
              correlation = corSymm(form = ~ Week | Plant),
              data = data_complete)

以Col×Row×Block为个体的重复测量模型 #

先创建复合个体变量，再按同样逻辑拟合模型：

# 创建Col×Row×Block的复合个体变量
data$CRB <- interaction(data$Col, data$Row, data$Block, drop = TRUE)
# 筛选完整观测的CRB个体
complete_crb <- names(which(table(data$CRB) == total_weeks))
data_crb_complete <- subset(data, CRB %in% complete_crb)

# 拟合模型
model_crb <- lme(Yield ~ Genotype + Table + Genotype:Week + Table:Week + Akima,
                 random = list(Col=~1, Row=~1, br=~1, bc=~1),
                 correlation = corSymm(form = ~ Week | CRB),
                 data = data_crb_complete)

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