纵向数据建模前的数据探索

医学和生信笔记

发布于 2026-03-17 18:07:02

570

文章被收录于专栏：医学和生信笔记医学和生信笔记

给大家介绍一下对于重复测量数据或者纵向数据可以使用哪些图形进行探索，因为重复测量数据可使用的方法很多，比如重复测量方差分析、广义估计方程、混合效应模型等，在进行正式的建模前非常有必要对数据进行一番探索。

unsetunset生成数据unsetunset

下面模拟一个数据，数据生成的过程如果看不懂可以不看。

60名患者接受两种治疗，分别在第0周、第1周、第2周、第3周记录患者的疼痛评分。

# ===================================================================
# 生成模拟重复测量数据（完全可复现！）
# ===================================================================
set.seed(123)  # 确保每次运行结果一致

n_subjects <- 60          # 60 名受试者
time_points <- c(0, 1, 2, 3)  # 4 个时间点：基线、1周、2周、3周

# 构建数据框：每个受试者在每个时间点都有一个观测
df <- expand.grid(
  subject_id = 1:n_subjects,
  time = time_points,
  stringsAsFactors = FALSE
)

# 分组：随机分配 Treatment / Control
df$group <- ifelse(df$subject_id <= n_subjects/2, "Treatment", "Control")

# 模拟结局变量 outcome（连续型，如疼痛评分）
# 假设：
# - 基线均值：20
# - Treatment 组随时间线性下降（每周 -1.5）
# - Control 组基本稳定（每周 -0.2）
# - 个体随机效应：每个受试者有自己的基线偏移（正态，SD=2）
# - 时间内相关：残差为 AR(1) 结构（相邻时间点相关）
# - 随机误差：N(0, 1.5)

library(dplyr)

df <- df %>%
  group_by(subject_id) %>%
  mutate(
    # 个体随机截距（基线偏移）
    random_intercept = rnorm(1, mean = 0, sd = 2),
    # 时间效应：Treatment 组下降更快
    time_effect = ifelse(group == "Treatment", -1.5 * time, -0.2 * time),
    # 残差（自相关结构：相邻时间点相关）
    # 使用简单方法：用 AR(1) 模拟，ε_t = ρ * ε_{t-1} + noise
    # 用递归方式模拟相关残差
    residual = {
      eps <- numeric(length(time))
      eps[1] <- rnorm(1, 0, 1.5)
      for(i in2:length(time)) {
        eps[i] <- 0.6 * eps[i-1] + rnorm(1, 0, sqrt(1.5^2 * (1 - 0.6^2)))
      }
      eps
    },
    # 最终 outcome
    outcome = 20 + random_intercept + time_effect + residual
  ) %>%
  ungroup()

# 查看前10行
head(df, 10)
## # A tibble: 10 × 7
##    subject_id  time group     random_intercept time_effect residual outcome
##         <int> <dbl> <chr>                <dbl>       <dbl>    <dbl>   <dbl>
##  1          1     0 Treatment           -1.12            0  0.569      19.4
##  2          2     0 Treatment           -0.460           0 -1.61       17.9
##  3          3     0 Treatment            3.12            0  1.38       24.5
##  4          4     0 Treatment            0.141           0  1.51       21.6
##  5          5     0 Treatment            0.259           0 -0.427      19.8
##  6          6     0 Treatment            3.43            0  0.00865    23.4
##  7          7     0 Treatment            0.922           0 -0.331      20.6
##  8          8     0 Treatment           -2.53            0 -0.489      17.0
##  9          9     0 Treatment           -1.37            0  0.358      19.0
## 10         10     0 Treatment           -0.891           0  3.28       22.4

unsetunset个体轨迹图unsetunset

首先看一下每个患者在4个时间点的疼痛评分。这张图展示了每一个受试者随时间变化的完整路径。每一条线代表一个患者的观测轨迹。

library(ggplot2)
ggplot(df, aes(x = time, y = outcome, group = subject_id, color = group)) +
  geom_line(alpha = 0.5, linewidth = 0.7) +
  labs(
    title = "图1：个体轨迹图（Spaghetti Plot）",
    x = "时间（周）", y = "疼痛评分",
    color = "分组"
  ) +
  theme_minimal() +
  scale_color_manual(values = c("Treatment" = "#E41A1C", "Control" = "#377EB8"))

这张图能展示很多问题，比如：

Treatment组个体是否都呈下降趋势？Control是否稳定？如果是，则说明“组别×时间”存在交互效应，提示应在模型中包含group*time项
个体间变异大吗？变异大提示存在显著的个体异质性（heterogeneity），如果用GEE建模则不用过度关注，因为GEE关注的是群体平均效应（又叫边际效应），但是如果混合效应模型则需要加入随机效应模型
是否有异常轨迹？可能是离群值，应检查数据录入错误，或考虑稳健标准误
Treatment和Control轨迹频繁交叉？表明表明个体响应差异大，但仍可能有总体趋势差异（GEE适合处理这种情况）

unsetunset平均趋势图unsetunset

然后是各组的平均趋势图，且带95%的置信区间。这张图展示的是每组在每个时间点的平均结局及其置信区间。它把个体复杂性“压缩”成群体趋势，突出整体模式。

# 先分组计算一下每组在每个时间点的平均疼痛评分，并计算置信区间
summary_df <- df %>%
  group_by(group, time) %>%
  summarise(
    mean_outcome = mean(outcome),
    se = sd(outcome) / sqrt(n()),
    .groups = "drop"
  ) %>%
  mutate(
    lower = mean_outcome - 1.96 * se,
    upper = mean_outcome + 1.96 * se
  )

# 然后画图即可
ggplot(summary_df, aes(x = time, y = mean_outcome, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = group), alpha = 0.2) +
  labs(
    title = "图2：组平均趋势（含95%置信区间）",
    x = "时间（周）", y = "平均结局",
    color = "分组", fill = "分组"
  ) +
  theme_minimal() +
  scale_color_manual(values = c("Treatment" = "#E41A1C", "Control" = "#377EB8")) +
  scale_fill_manual(values = c("Treatment" = "#E41A1C", "Control" = "#377EB8"))

这张图可以展示的问题：

两组基线是否可比？（0周时Treatment和Control的CI是否重叠，若不重叠，说明随机化可能失败，可能需要调整基线）
Treatment组是否随时间下降更快？如果是，则提示需要在模型中加入group:time交互项
下降是线性的吗？若非线性，可考虑将time设为分类变量，或加入二次项
误差条宽度，CI宽表示变异性大，标准误较大，检验效能降低

unsetunset箱线图unsetunset

这张图可以展示每个时间点、每组的分布特征：中位数、四分位距、偏态、离群值、分布形状。可以用于判断：

数据是否符合正态假设？是否需要转换？
是否存在极端个体影响结果？
干预是否不仅改变均值，还改变了变异性？

ggplot(df, aes(x = factor(time), y = outcome, fill = group)) +
  geom_violin(alpha = 0.7, position = position_dodge(0.8)) +
  geom_boxplot(width = 0.2, position = position_dodge(0.8), alpha = 0.8) +
  labs(
    title = "图3：各时间点结局分布（小提琴 + 箱线）",
    x = "时间点（周）", y = "结局变量",
    fill = "分组"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c("Treatment" = "#E41A1C", "Control" = "#377EB8"))

这张图可以看出以下问题：

分布是否对称？若严重偏态可考虑进行数据变换，比如log
是否有离群值？箱线图外的点可能是离群值，分析时可能需要剔除
组内变化vs组间差异，如Control组内部波动大，但Treatment组稳定下降，提示干预可能减少变异。

unsetunset时间点相关性热图unsetunset

这张图可以展示不同时间点之间的皮尔逊相关系数。颜色越深红，相关性越高。

这张图对于广义估计方程比较有用，因为可以直观展示不同时间点之间的相关性结构。比如：

所有时间点间相关性接近相同，选用”exchangeable”
相邻强相关，越远越弱，选用自回归结构
各对时间点相关差异大，无规律，选择”unstructured”
几乎无相关，选择”independence”

# 转为宽格式，计算各时间点之间的相关
df_wide <- df %>%
  select(subject_id, group, time, outcome) %>%
  tidyr::pivot_wider(names_from = time, values_from = outcome, names_prefix = "time_") %>%
  select(-group)  # 只保留时间点变量

cor_matrix <- cor(df_wide[, c("time_0", "time_1", "time_2", "time_3")], use = "complete.obs")
print("时间点间相关系数矩阵：")
## [1] "时间点间相关系数矩阵："
print(round(cor_matrix, 3))
##        time_0 time_1 time_2 time_3
## time_0  1.000  0.829  0.625  0.535
## time_1  0.829  1.000  0.835  0.677
## time_2  0.625  0.835  1.000  0.860
## time_3  0.535  0.677  0.860  1.000

# 可视化相关矩阵（热力图）
library(reshape2)
library(pheatmap)

cor_melt <- melt(cor_matrix)
pheatmap::pheatmap(
  cor_matrix,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  main = "图4：时间点间结局相关性热图"
)