一个 R 包搞定医学统计全流程（基线表/Cox/RCS/森林图/ROC），公卫/临床科研人必看！

开源R包medstats：搞定医学统计全流程（附 10 个代码示例）

🤯 临床研究者的"统计焦虑"

做过临床研究的小伙伴一定深有体会——数据分析这条路，坑太多了：

📦 medstats

medstats 是一个专为临床医学统计和流行病学研究设计的 R 包，由开发者（GitHub: @shanjiayu1）开发并开源。

它将临床研究中最常用的统计分析流程封装为 15 个核心函数，覆盖从数据清洗、统计建模到可视化、结果导出的全链路工作流。

一句话概括：装一个包，干完一整篇论文的统计分析。

✨ 核心功能与亮点

medstats 的 15 个导出函数可分为 四大模块：

1️⃣ 数据处理模块

• long_to_surv_data() —— 将纵向临床记录一键转换为"一人一行"的生存分析数据格式

• merge_duplicate_records() —— 合并重复记录，自动取每列第一个非缺失值

• make_table1() —— 生成基线特征表，支持分组比较

2️⃣ 统计建模模块

• run_glm_auto() —— 自动批量运行单因素 + 多因素广义线性模型（线性回归、Logistic 回归）

• run_cox_auto() —— 自动批量运行单因素 + 多因素 Cox 回归

• longdata_analysis() —— 重复测量数据综合分析（描述统计 + 组间比较 + GEE 模型）

3️⃣ 可视化模块

• plot_km() —— K-M 累积事件曲线 + 风险表

• plot_forest() —— 发表级森林图

• plot_rcs() —— 限制性立方样条图（支持线性模型和 Cox 模型）

• plot_roc() —— ROC 曲线 + AUC + 最佳截断值

• plot_sankey() —— 桑基图（状态流转可视化）

• plot_meanse() —— 均值 ± 标准误折线图

• plot_stacked() —— 堆叠百分比柱状图

4️⃣ 表格与导出模块

• format_flextable() —— 一行代码将数据框/gtsummary 对象格式化为发表级三线表

• export_word() —— 将多个表格和图片一键导出到同一个 Word 文档

🌟 亮点总结

📥 安装方式

medstats 目前在 GitHub 上开发维护，使用 remotes 包安装即可：

# 安装 remotes（如已安装可跳过）
install.packages("remotes")

# 安装 medstats
remotes::install_github(
  "shanjiayu1/medstats",
  dependencies = TRUE
)

# 加载包
library(medstats)

💡 使用 dependencies = TRUE 会同时安装示例、测试和 vignette 所需的依赖包，推荐使用。

📝 使用示例

示例 1：一键生成发表级三线表

library(medstats)

# 直接格式化数据框
ft_data <- format_flextable(head(mtcars[, 1:5]))

table1 <- make_table1(
  data = gtsummary::trial,
  vars = c("age", "marker", "stage", "grade"),
  specific_vars = "marker",  # Report as median (P25, P75)
  group_var = "trt"
)

table1

示例 2：批量跑 Logistic /COX 回归生成整理后的结果

# 一行代码搞定单因素 + 多因素 Logistic 回归
logistic_results <- run_glm_auto(
  data = gtsummary::trial,
  vars = c("age", "stage"),
  outcome_var = "response",
  family = "binomial"
)

# 格式化为发表级表格
format_flextable(logistic_results)

回归表自动使用上标模型标记，脚注标注 ¹ Univariable analysis. 和 ² Multivariable analysis.，无需手动添加。

示例 3：K-M 生存曲线 + 风险表

library(survival)

# 准备数据
lung_data <- survival::lung
lung_data$status_event <- as.integer(lung_data$status == 2)
lung_data$sex <- factor(
  lung_data$sex,
  levels = c(1, 2),
  labels = c("Male", "Female")
)

# 一行画 K-M 曲线
plot_km(
  data = lung_data,
  group_var = "sex",
  time_var = "time",
  status_var = "status_event",
  legend_labs = c("Male", "Female"),
  legend_title = "Sex",
  xlab = "Follow-up time (days)",
  ylab = "Cumulative mortality (%)",
  xlim = c(0, 1000),
  break_time = 200,
  show_risk_table = TRUE,    # 自动添加风险表
  save_filename = "Lung_KM.png"
)

示例 4：限制性立方样条（RCS）

# 线性模型 RCS
rcs_linear <- plot_rcs(
  data = mtcars,
  exposure = "wt",
  outcome = "mpg",
  covars = c("hp", "disp"),
  nk = 4,
  model_type = "linear",
  xlab = "Weight",
  ylab = "Predicted MPG"
)
rcs_linear$plot

# Cox 模型 RCS
lung_rcs <- survival::lung
lung_rcs$status_event <- as.integer(lung_rcs$status == 2)

rcs_cox <- plot_rcs(
  data = lung_rcs,
  exposure = "age",
  outcome = "Surv(time, status_event)",
  covars = c("sex", "ph.ecog"),
  model_type = "cox",
  xlab = "Age",
  ylab = "Hazard Ratio"
)
rcs_cox$plot

示例 5：森林图

forest_data <- data.frame(
  Variable = c("Age", "Stage II", "Stage III"),
  `OR (95% CI)` = c(
    "1.02 (0.99, 1.05)",
    "1.45 (0.82, 2.56)",
    "2.10 (1.12, 3.94)"
  ),
  `P value` = c("0.180", "0.200", "0.021"),
  check.names = FALSE
)

plot_forest(
  data = forest_data,
  ci_column = "OR (95% CI)",
  x_ticks = c(0, 0.5, 1, 2, 4),
  output_name = "Forest_plot.png"
)

示例 6：ROC 曲线 + AUC + 最佳截断值 📈

诊断试验评价利器！一行代码画出 ROC 曲线，自动标注 AUC 值和最佳截断点，再也不用手动算 Youden 指数了。

library(medstats)

# 使用 survival 包的 lung 数据集演示
# 假设我们想用 age 预测患者是否在随访期内死亡
library(survival)
lung_data <- survival::lung
lung_data$status_event <- as.integer(lung_data$status == 2)

# 一行代码绘制 ROC 曲线
roc_result <- plot_roc(
  data = lung_data,
  marker_var = "age",         # 预测变量（连续变量）
  outcome_var = "status_event", # 二分类结局（0/1）
  positive_label = 1,          # 阳性事件的编码值
  xlab = "1 - Specificity",
  ylab = "Sensitivity",
  save_filename = "ROC_age.png"
)

# 查看结果：AUC、最佳截断值、灵敏度、特异度
roc_result$auc          # AUC 值
roc_result$best_cutoff  # 最佳截断值
roc_result$sensitivity  # 对应灵敏度
roc_result$specificity  # 对应特异度

📌 函数自动使用 Youden 指数（灵敏度 + 特异度 - 1）确定最佳截断值，并标注在曲线上。适合诊断试验、标志物评价等场景。

示例 7：桑基图（状态流转可视化）🌊

随访数据中，患者状态在不同时间点之间如何流转？桑基图一图胜千言，直观展示患者从"轻症"到"重症"再到"康复"的动态变化轨迹。

library(medstats)

# 使用 ChickWeight 数据集演示患者状态流转
sankey_data <- datasets::ChickWeight |>
  dplyr::filter(Time %in% c(0, 10, 20)) |>
  dplyr::mutate(
    # 将连续变量离散化为状态类别
    visit = factor(
      paste0("Day ", Time),
      levels = c("Day 0", "Day 10", "Day 20")
    ),
    weight_status = dplyr::case_when(
      weight < 50 ~ "Light",     # 轻量
      weight < 150 ~ "Normal",   # 正常
      TRUE ~ "Heavy"             # 重量级
    ),
    weight_status = factor(
      weight_status,
      levels = c("Light", "Normal", "Heavy")
    )
  )

# 一行代码绘制桑基图
sankey_plot <- plot_sankey(
  data = sankey_data,
  id_var = "Chick",           # 个体 ID
  time_var = "visit",         # 时间点变量
  state_var = "weight_status", # 状态变量
  na_strategy = "show",       # 缺失值的处理策略：显示为单独的流
  missing_label = "Drop-out"  # 缺失值的标签名
)

sankey_plot

📌 非常适合展示纵向随访研究中患者状态的流转路径，如疾病分期变化、治疗方案的切换等。na_strategy = "show" 可将脱访患者单独显示，方便评估失访情况。

示例 8：均值 ± 标准误折线图 📉

纵向随访研究中，各组的均值随时间如何变化？折线图 + 误差棒 + 组间差异检验，一气呵成。

library(medstats)

# 使用 ChickWeight 数据集，提取关键时间点
growth_data <- datasets::ChickWeight |>
  dplyr::filter(Time %in% c(0, 4, 10, 14, 21)) |>
  dplyr::mutate(
    time_label = paste0("Day ", Time),
    diet_label = paste0("Diet ", Diet)
  )

# 绘制各饮食组的均值 ± 标准误折线图
meanse_result <- plot_meanse(
  data = growth_data,
  target_var = "weight",      # 数值结局变量
  time_var = "time_label",    # 时间点变量
  group_var = "diet_label",   # 分组变量
  xlab = "Growth time (days)",
  ylab = "Mean weight (g)",
  legend_title = "Diet"
)

meanse_result$plot

# 🎯 两两组比较：自动在每个时间点进行组间差异检验
two_diet_result <- plot_meanse(
  data = dplyr::filter(growth_data, diet_label %in% c("Diet 1", "Diet 2")),
  target_var = "weight",
  time_var = "time_label",
  group_var = "diet_label",
  test_method = "wilcox",     # 使用 Wilcoxon 秩和检验
  legend_title = "Diet"
)

# 查看每个时间点的 p 值
two_diet_result$test_data

📌 设置 test_method = "t" 或 "wilcox" 后，函数会在每个时间点自动进行两组比较，显著结果直接标注在图上对应位置，p 值可通过 $test_data 获取。非常适合纵向随访的组间趋势对比！

示例 9：堆叠百分比柱状图 📊

想看各分类在不同时间点的占比变化？堆叠百分比柱状图让趋势一目了然。

library(medstats)

# 使用 ChickWeight 数据集，将体重离散化
stacked_data <- datasets::ChickWeight
stacked_data$Time <- factor(
  stacked_data$Time,
  levels = sort(unique(stacked_data$Time))
)

# 绘制堆叠百分比柱状图
stacked_result <- plot_stacked(
  data = stacked_data,
  target_var = "weight",     # 数值变量（将被离散化）
  time_var = "Time",         # 时间点变量
  group_var = "Diet",        # 分组变量（可选）
  breaks = c(-Inf, 100, 200, 300, Inf),  # 离散化分界点
  labels = c("≤100 g", "101–200 g", "201–300 g", ">300 g"),  # 分段标签
  colors = c("#B5D1E8", "#A3D9A5", "#F2C68F", "#EB938F"),   # 自定义配色
  legend_title = "Weight range",
  label_size = 4.5           # 柱内百分比标签字号
)

stacked_result$plot

📌 group_var = NULL（默认）时，每个时间点画一个堆叠柱；设置 group_var 后，百分比按"时间×组别"计算，各组并排展示。用 label_size 调整柱内百分比字号，colors 自定义配色方案，出图即发表级！

示例 10：表格 + 图片一键导出 Word

# 准备表格和图片
table1 <- head(mtcars)
p1 <- ggplot2::ggplot(mtcars, ggplot2::aes(wt, mpg)) +
  ggplot2::geom_point()

# 详细方式：自定义标题
export_word(
  data_list = list(table1, p1),
  table_titles = c(
    "Table 1. mtcars dataset",
    "Figure 1. MPG and weight"
  ),
  output_file = "tables_and_plots.docx",
  figure_width = 6,
  figure_height = 5
)

# 简化方式：直接传入对象，标题自动生成
export_word(table1, p1, "tables_and_plots.docx")

📌 表格标题在表格上方，图片标题在图片下方，符合学术规范。还可以混入 PNG 文件路径，自动按比例缩放。

🎯 应用场景

medstats 适用于但不限于以下场景：

无论你是：

• 🏥 临床医生——想快速完成论文统计部分

• 🎓 研究生/博士生——正在为毕业论文的数据分析发愁

• 📊 流行病学研究者——需要处理大量随访数据

• 💊 药企统计师——需要标准化、可复现的统计流程

medstats 都能帮你大幅提升效率！

🏗️ 技术架构

medstats 的代码结构清晰，按功能模块组织：

medstats-clinical-r/
├── R/
│   ├── data-process.R       # 数据处理（生存数据转换、重复记录合并、基线表）
│   ├── regression.R          # 统计建模（GLM、Cox、GEE）
│   ├── plots-survival.R     # 生存可视化（K-M 曲线、桑基图）
│   ├── plots-diagnostic.R   # 诊断可视化（ROC、RCS、森林图）
│   ├── plots-line.R         # 趋势可视化（均值折线图、堆叠柱状图）
│   ├── format.R             # 表格格式化与 Word 导出
│   └── utils.R              # 工具函数
├── man/                      # 函数文档
├── vignettes/                # 教程文档
├── tests/                    # 单元测试
├── DESCRIPTION               # 包描述文件
├── NAMESPACE                 # 导出函数声明
├── README.md                 # 英文文档
└── README.zh-CN.md           # 中文文档

核心依赖包

medstats 站在巨人的肩膀上，整合了 R 生态中最优秀的包：

🙌 总结

medstats 是一个"小而美"的 R 包，它不试图重新发明轮子，而是将临床研究中最常用的统计分析流程串联起来，用最简洁的接口完成最复杂的工作：

1. ✅ 15 个核心函数覆盖临床统计全流程

2. ✅ 一行代码完成批量回归 + 格式化

3. ✅ 发表级三线表、K-M 曲线、森林图、RCS 图

4. ✅ 一键导出 Word，表格图片混合排版

5. ✅ 中文文档支持，上手零门槛

如果你正在做临床研究、写论文、处理随访数据，medstats 绝对值得一试！

🔗 项目地址：https://github.com/shanjiayu1/medstats-clinical-r

📦 安装命令：remotes::install_github("shanjiayu1/medstats", dependencies = TRUE)

觉得好用的话，别忘了给个 Star ⭐ 支持一下开发者！也欢迎在 Issues 中提出建议和反馈。

本文基于 medstats v0.1.0 版本的仓库内容撰写，所有代码示例均来自项目 README。