使用R计算基于linear mixed effects model的Longitudinal repeated measures study的样本量

阿华AIGC实验室

2026-3-31

Hey 你好呀！我完全懂你的困扰——找了好多R的样本量计算例子都是有分组的，偏偏你的研究是单组重复测量，要验证患者在3次随访里能不能稳住基线得分，而且你对R还不太熟。刚好我可以给你一步步演示怎么用线性混合效应模型（LMM）来搞定这个计算，都是实操的代码，你跟着跑就行~

先理清楚你的核心参数：

单组，每个患者有4个时间点（基线+3次随访）
基线得分均值20，标准差4.5
临床意义上的得分下降是6分（也就是如果下降超过6分，就认为有临床影响）
研究目的：确认患者随访时的得分是否维持在基线水平（要么检测是否出现临床相关下降，要么证明下降没到这个程度）

第一步：安装并加载需要的包

我们要用lme4来拟合混合模型，simr来做样本量的模拟和power分析——模拟法是单组重复测量LMM最直观的样本量计算方式，因为没有现成的简单公式，而且灵活度高，适合你这种场景。

# 安装包（第一次用的时候跑）
install.packages(c("lme4", "simr"))

# 加载包
library(lme4)
library(simr)

第二步：构建模拟数据集（先假设一个初始样本量）

我们先从一个小的初始样本量（比如20人）开始，构建符合你研究场景的模拟数据，这样才能拟合模型做后续的power分析。这里我会把你给的参数都加进去，还会考虑个体间的变异（毕竟不同患者的基线得分肯定有差异）：

# 设定核心参数，你可以根据你的研究调整！
n_initial <- 20  # 先假设20个样本，后面再找最优值
time_points <- c(0,1,2,3)  # 0=基线，1/2/3=三次随访
baseline_mean <- 20  # 基线均值
clin_drop <- 6  # 临床相关的下降幅度
sd_individual <- 3  # 个体间的变异SD（不同患者基线的差异）
sd_residual <- sqrt(4.5^2 - sd_individual^2)  # 个体内的残差SD（同一个患者不同时间点的波动）

# 生成模拟数据
set.seed(123)  # 固定随机种子，让结果可重复
dat <- expand.grid(id = 1:n_initial, time = time_points)  # 每个患者对应4个时间点
dat$rand_int <- rnorm(n_initial, 0, sd_individual)[dat$id]  # 给每个患者加个体随机截距
# 计算得分：基线得分=基线均值+个体截距；随访时下降clin_drop，再加上残差波动
dat$score <- baseline_mean + dat$rand_int + ifelse(dat$time == 0, 0, -clin_drop) + rnorm(nrow(dat), 0, sd_residual)

第三步：拟合线性混合效应模型

我们拟合一个带随机截距的LMM，把时间作为固定效应（因为我们要检测时间对得分的影响）：

# 拟合模型
model <- lmer(score ~ time + (1|id), data = dat)
# 看看模型结果，确认参数符合我们的设定
summary(model)

第四步：用simr做Power分析和样本量估算

现在我们要计算，当真实存在6分的下降时，不同样本量下我们有多大的概率能检测到这个差异（也就是power），然后找到power达到80%（行业通用标准）的最小样本量。

# 先看初始样本量（20人）的power
powerSim(model, test = fixed("time"))

# 模拟不同样本量的power变化（从10到100，每次加10）
power_curve <- powerCurve(model, test = fixed("time"), along = "id", breaks = seq(10, 100, 10))
# 可视化power曲线，更直观
plot(power_curve)

# 直接计算达到80% power的最小样本量
required_n <- findPower(model, test = fixed("time"), target = 0.8, along = "id")
print(required_n)

一些你可能需要调整的细节

个体间变异的设定：我这里设的sd_individual=3，如果你不确定你的研究里个体差异有多大，也可以把它设为0（这时候残差SD直接用4.5），不过现实中不同患者的基线肯定有差异，所以尽量根据类似研究的结果调整这个值。

等效性检验vs优效性检验：上面的代码是检测“是否存在临床相关的下降”（优效性检验），如果你的研究目的是证明患者能维持基线（即下降不超过6分），这属于等效性检验，代码要稍微调整：

# 自定义等效性检验函数：检验时间效应的系数是否大于-6（即下降不超6分）
eq_test <- function(mod) {
  time_coef <- fixef(mod)["time"]
  time_se <- se.fixef(mod)["time"]
  t_stat <- (time_coef + 6) / time_se  # 检验H0: time <= -6 vs H1: time > -6
  p_val <- 1 - pt(t_stat, df = df.residual(mod))
  return(p_val)
}

# 假设真实情况是没有下降（time系数为0），计算能证明等效的power
fixef(model)["time"] <- 0
powerSim(model, test = eq_test)