针对58名受试者4个时间点的lmer模型构建咨询

数据概况 #

• 样本：58名受试者，每人在4个时间点重复测量progesterone（孕酮）和factor1两个指标
• 已知趋势：所有受试者的progesterone在时间点1-3呈上升趋势，时间点4出现下降

初始模型评估 #

你设定的初始模型为：

mod <- lmer(Factor1 ~ timepoint + progesterone + (1 | ID ), data = df)

该模型已通过(1 | ID)纳入个体随机截距，能处理重复测量的非独立性，但结合已知的progesterone时间趋势，还可以进一步优化。

模型优化建议 #

1. 加入时间点与孕酮的交互项
   由于progesterone的时间趋势并非线性（1-3上升、4下降），建议添加timepoint和progesterone的交互项，检验不同时间点下孕酮对factor1的预测效应是否存在差异：

   mod_interaction <- lmer(Factor1 ~ timepoint * progesterone + (1 | ID ), data = df)
   

   之后可用似然比检验（anova(mod, mod_interaction)）对比两个模型的拟合优度，判断交互项是否有必要保留。

2. 尝试更灵活的随机效应结构
   如果个体间的孕酮效应存在差异，可以加入随机斜率，允许每个个体的孕酮对factor1的影响程度不同：

   mod_random_slope <- lmer(Factor1 ~ timepoint * progesterone + (1 + progesterone | ID ), data = df)
   

   若模型无法收敛，可简化随机结构，或通过lmerControl()调整优化器。

3. 调整时间点的编码方式
   考虑到时间点的非线性趋势，除了用连续型timepoint，还可以：

   • 将时间点设为分类变量，直接检验每个时间点与基线（如时间点1）的差异
   • 引入二次项捕捉非线性趋势
     示例代码：

   # 分类变量编码
   df$timepoint_cat <- factor(df$timepoint)
   mod_cat <- lmer(Factor1 ~ timepoint_cat * progesterone + (1 | ID ), data = df)
   
   # 二次项编码
   df$timepoint_sq <- df$timepoint^2
   mod_quad <- lmer(Factor1 ~ timepoint + timepoint_sq + progesterone + timepoint:progesterone + timepoint_sq:progesterone + (1 | ID ), data = df)
   

4. 模型诊断
   拟合模型后，需检查残差是否符合正态性、方差齐性假设：

   # 绘制残差图
   plot(resid(mod_interaction))
   qqnorm(resid(mod_interaction))
   qqline(resid(mod_interaction))
   

   同时检查随机效应的分布是否合理，避免出现奇异拟合。

模拟数据代码（中文注释） #

set.seed(123)

# 设置受试者数量和时间点数量
num_part <- 58
num_timepoints <- 4

# 生成孕酮数据
progesterone <- rnorm(num_part * num_timepoints, mean = 0, sd = 1)
progesterone[is.na(progesterone)] <- -3.8  # 将缺失值替换为指定值-3.8

# 确保孕酮值落在目标范围[-3.8, 1.4]内
progesterone <- pmin(pmax(progesterone, -3.8), 1.4)

# 生成factor1数据
factor1 <- rnorm(num_part * num_timepoints, mean = 0.03, sd = 1)
# 确保factor1值落在目标范围[-1.840773, 2.149765]内
factor1 <- pmin(pmax(factor1, -1.840773), 2.149765)

# 创建数据框
df <- data.frame(
  ID = rep(1:num_part, each = num_timepoints),
  timepoint = rep(1:num_timepoints, times = num_part),
  progesterone = progesterone,
  Factor1 = factor1
)

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