在R中拟合分段线性混合效应回归模型如何避免奇异性错误？
问题场景与原模型问题
我用线性混合效应模型拟合两臂临床试验的纵向血红蛋白（Hb）计数长格式数据，数据子集如下：
pid site arm pday Hb pday2 I 1 MOB2004_001 Boane SP 0 8.5 0 0 2 MOB2004_001 Boane SP 3 8.0 9 0 3 MOB2004_001 Boane SP 7 8.9 49 1 4 MOB2004_001 Boane SP 14 8.9 196 1 5 MOB2004_002 Boane SP 0 11.5 0 0 6 MOB2004_002 Boane SP 3 9.1 9 0
观察均值响应曲线后，我希望建模Hb与时间的二次关系，且在第3天前后使用不同参数，因此创建指示变量I（pday<3时为0，否则为1）。同时需要评估研究组间差异，必须包含研究组-时间交互项。
原模型设定如下：
lme(data=dat.mod2, fixed=Hb~ arm*pday + pday*I + pday2 + pday2*I , random=~1|site/pid, na.action = na.omit)
对应的模型公式：
$$ Hb_{ik} = \beta_{00} + \beta_{01}*I + \beta_{10}*day + \beta_{11}*dayI + \beta_{20}day^2 + \beta_{21}day^2I + \beta_3arm_i + \beta_4 *arm_iday + u_{ik} + v_{k} $$
拟合时出现共线性错误：
Error in MEEM(object, conLin, control$niterEM) :
Singularity in backsolve at level 0, block 1
单独包含pday*I或pday2*I时模型可正常拟合，但同时包含就触发错误；用lme4的lmer拟合也出现秩亏缺问题。我需要估计节点前后的二次项，不希望使用平滑处理。
问题根源
核心错误是指示变量I与原时间项的交互引发多重共线性：
- 当
I=0时，pday*I和pday2*I均为0，时间项仅由pday和pday2构成； - 当
I=1时，时间项变为pday + pday*I + pday2 + pday2*I = pday*(1+I) + pday2*(1+I)，本质是原时间项的线性缩放，导致变量间存在完全线性依赖，模型矩阵秩亏缺。
正确建模思路
要实现节点前后的分段二次模型，需重构时间变量，避免共线性：
方法1：构建分段时间变量（推荐）
创建以第3天为节点的分段时间变量，将节点前后的时间拆分为独立变量：
- 生成新变量：
t1：pday<3时取pday值，否则为0（节点前的时间）t2：pday>=3时取pday-3（节点后从0开始计数的时间），否则为0- 同时生成平方项
t1_sq = t1^2、t2_sq = t2^2
- 模型设定：包含分组、分段时间的线性/二次项，以及分组与各时间项的交互，确保组间差异可评估。
R代码示例：
# 生成分段时间变量 dat.mod2 <- dat.mod2 %>% mutate( t1 = ifelse(pday < 3, pday, 0), t2 = ifelse(pday >= 3, pday - 3, 0), t1_sq = t1^2, t2_sq = t2^2 ) # 拟合lme模型 fit_lme <- lme( data = dat.mod2, fixed = Hb ~ arm + t1 + t1_sq + t2 + t2_sq + arm:t1 + arm:t1_sq + arm:t2 + arm:t2_sq, random = ~1|site/pid, na.action = na.omit )
该设定下，t1和t2不会同时非零，彻底消除共线性，且能分别估计节点前后的二次趋势及组间差异。
方法2：中心化时间变量后结合指示变量
若坚持使用原指示变量I，需先将时间变量中心化到节点（第3天），减少共线性风险：
- 生成中心化时间变量
cent_pday = pday - 3及平方项cent_pday_sq = cent_pday^2 - 模型中纳入分组与中心化时间的交互、指示变量与中心化时间的交互
R代码示例：
# 生成中心化时间变量 dat.mod2 <- dat.mod2 %>% mutate( cent_pday = pday - 3, cent_pday_sq = cent_pday^2 ) # 拟合模型 fit_lme <- lme( data = dat.mod2, fixed = Hb ~ arm*cent_pday + arm*cent_pday_sq + I*cent_pday + I*cent_pday_sq, random = ~1|site/pid, na.action = na.omit )
需注意：拟合前建议用cor()检查变量相关性，或用car::vif()计算方差膨胀因子（VIF<5为可接受范围），确保无强共线性。
验证步骤
- 拟合前检查固定效应变量间的相关性，相关系数>0.7提示强共线性；
- 用方差膨胀因子（VIF）量化共线性程度，VIF>10需调整模型；
- 逐步添加变量，定位触发秩亏缺的变量组合，验证模型设定合理性。
内容的提问来源于stack exchange，提问作者modellingmania
