如何在R中生成重复测量模型用的临床试验数据集？

核心参数设定 #

先明确所有预设条件，统一定义参数方便后续调整：

• 每组样本量：n_per_group <- 309
• 试验分组：groups <- c("a", "b")
• 评估次数（基线+26次随访）：time_points <- 0:26 # 0代表基线，1-26对应每两周一次，共27次
• 结局变量标准差：sd_outcome <- 11
• 失访比例：dropout_rate <- 0.3
• 单评估随机缺失概率：miss_single_rate <- 0.05 # 可按需调整

步骤1：构建受试者-时间基础框架 #

用tidyr生成每个受试者在所有时间点的记录，并分配分组：

library(tidyverse)

# 生成受试者ID与分组信息
subjects <- tibble(
  subject_id = 1:(2*n_per_group),
  group = rep(groups, each = n_per_group)
)

# 扩展为长格式的受试者-时间点数据集
data_long <- subjects %>%
  crossing(time = time_points) %>%
  arrange(subject_id, time)

步骤2：定义时间依赖均值趋势+基线关联 #

假设两组有不同的时间趋势，且受试者得分与自身基线值强相关（通过随机截距实现关联）：

# 设定组水平的时间趋势参数
# 组a：基线均值50，每两周上升0.1；组b：基线均值52，每两周上升0.2
mean_params <- tibble(
  group = groups,
  baseline_mean = c(50, 52),
  time_slope = c(0.1, 0.2)
)

# 合并趋势参数到数据集
data_long <- data_long %>%
  left_join(mean_params, by = "group")

# 生成受试者随机截距（围绕组基线均值波动，控制基线值与后续测量的关联强度）
set.seed(123) # 设种子保证结果可重复
subjects <- subjects %>%
  left_join(mean_params, by = "group") %>%
  mutate(subject_intercept = rnorm(n(), mean = baseline_mean, sd = sqrt(8)))

# 合并随机截距到长数据
data_long <- data_long %>%
  left_join(subjects %>% select(subject_id, subject_intercept), by = "subject_id")

# 计算每个时间点的预期得分
data_long <- data_long %>%
  mutate(expected_value = subject_intercept + time_slope * time)

步骤3：加入随机误差并截断到0-100范围 #

set.seed(456)
data_long <- data_long %>%
  mutate(
    outcome = rnorm(n(), mean = expected_value, sd = sd_outcome),
    # 将结局变量限制在0-100范围内
    outcome = pmax(pmin(outcome, 100), 0)
  )

步骤4：处理失访与随机缺失 #

4.1 处理30%受试者提前失访

随机选择30%的受试者，设定他们的最后观察时间：

set.seed(789)
dropout_subjects <- subjects %>%
  sample_frac(dropout_rate) %>%
  mutate(
    # 失访时间设为基线后1-25之间的随机值（确保至少完成1次随访）
    last_time = sample(1:25, size = n(), replace = TRUE)
  )

# 标记失访后的记录为缺失
data_long <- data_long %>%
  left_join(dropout_subjects %>% select(subject_id, last_time), by = "subject_id") %>%
  mutate(
    outcome = ifelse(!is.na(last_time) & time > last_time, NA, outcome)
  ) %>%
  select(-last_time)

4.2 处理随机缺失个别评估

在未失访的记录中，随机标记部分评估为缺失：

set.seed(101)
data_long <- data_long %>%
  mutate(
    outcome = ifelse(
      !is.na(outcome) & runif(n()) < miss_single_rate,
      NA,
      outcome
    )
  )

步骤5：整理最终数据集 #

# 转换为宽格式（可选，适配部分重复测量分析方法）
data_wide <- data_long %>%
  select(subject_id, group, time, outcome) %>%
  pivot_wider(names_from = time, values_from = outcome, names_prefix = "time_")

# 查看长格式数据结构
glimpse(data_long %>% select(subject_id, group, time, outcome))

可调参数说明 #

• 时间趋势可替换为二次项（time_slope2 * time^2）或分段线性趋势，适配不同研究假设
• 随机截距的标准差可调整，控制基线值与后续测量的关联强度
• 单评估随机缺失概率miss_single_rate可根据实际场景修改

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