教程地址：https://genentech.github.io/survivalCCW/
survivalCCW包是一个轻量级的包，有助于构建R中时间到事件端点的克隆删失加权（CCW）分析数据。
背景
在观察性研究中估计非随机化治疗的效果可能具有挑战性，因为在患者的治疗策略启动之前，患者的治疗分配通常是未知的。天真地使用观察性治疗策略通常会导致不朽的时间偏差和其他偏差。例如，考虑治疗持续时间对总体生存期的影响。比较观察到接受治疗时间更长的患者的分析会遭受不朽的时间偏差，因为治疗持续时间更长与生存期更长有关。解决这种偏差的传统方法包括地标分析（landmark analysis），它受到选择偏差和统计学效能降低的影响，以及时间相关的cox模型，它阻止了对绝对风险的估计，导致可解释性更差。克隆删失加权（CCW）是一种稳健的方法，可以解决不朽的时间偏差，同时保持研究效能和可解释性，并允许绝对风险估计。
实施CCW可能很耗时，不仅因为这是一种大多数应用生物统计学家不熟悉的新方法，还因为它涉及细致入微的数据操作。所有患者都被克隆并分配到基线治疗或未治疗；然后他们被转换成“长”随访格式，每个patient-clone-period组合一行；最后，计算每个patient-clone-period的删失权重的反比概率。
{survivalCCW}
为了便于CCW的轻松实现并提供标准方法，我们引入了survivalCCW包，这是一个轻量级、开源的R包，旨在简化使用CCW进行生存分析的过程。该软件包提供了格式化数据、计算权重和轻松执行生存分析的功能。{survivalCCW包的工作流程直观且用户友好，便于快速数据操作和推理。
下面显示了一个说明性示例。
# 使用管道操作符处理数据框df
df |>
# 创建克隆数据
# create_clones函数用于生成克隆数据集,主要用于时间依赖性分析
# id: 指定个体ID列名
# event: 指定事件状态列名
# time_to_event: 指定到事件发生的时间列名
# exposure: 指定暴露状态列名
# time_to_exposure: 指定到暴露发生的时间列名
# ced_window: 设置累积暴露时间窗口为100个时间单位
create_clones(
id = 'id',
event = 'event',
time_to_event = 'timetoevent',
exposure = 'exposure',
time_to_exposure = 'timetoexposure',
ced_window = 100
) |>
# 将克隆数据转换为长格式
# cast_clones_to_long()函数将宽格式的克隆数据转换为长格式,便于后续分析
cast_clones_to_long() |>
# 生成累积暴露权重
# generate_ccw()函数计算累积暴露权重
# 参数为需要调整的协变量向量c('cov1', 'cov2')
generate_ccw(c('cov1', 'cov2'))
如何使用 survivalCCW 进行克隆删失加权生存分析
克隆删失加权
这个轻量级的包描述了如何进行克隆删失加权（CCW），以解决生存分析中不朽的时间偏差问题。这个小插曲将浏览马林格等人2020年发布的应用教程。
前言
CCW在观察性研究中存在不朽的人时偏差是有用的。例如，当比较非小细胞肺癌（NSCLC）中手术接受者和非接受者时，手术组将比非手术组有更长的生存时间，因为非手术组包括在接受手术之前死亡的患者。这是一种不朽时间偏差。
马林格发布的ccw数据集使用这个精确的设置作为示例。让我们探索一下数据集。
library(tidyverse)
library(DT) # For better printing of data.frames
data(dummy_data)
head(dummy_data) |>
DT::datatable(
rownames = FALSE,
options = list(
scrollX = TRUE,
paging = FALSE
)
)
glimpse(dummy_data)
数据说明
- id：患者标识符
- timetoexposure：暴露时间，连续，非暴露患者的NA
- exposure：二进制，1=暴露，0=未暴露
- event：事件发生，二进制，1=事件，0=删失
- timetoevent：事件之间的时间，连续，上限为365
- cov1：二元协变量——与暴露时间和事件时间正相关
- cov2：连续协变量——与暴露时间和事件时间正相关
请注意，此包解决了协变量都在基线定义的情况。
创建克隆
第一步是创建克隆。这可以使用create_clones的SurvivalCCW函数对任何时间到事件的结果完成。为了create_clones工作，我们需要传递一个包含以下列的每患者一行的data. frame：
- id变量（在本例中为id）
- 表示删失（0）或事件（1）（在本例中为event）的传统结果变量。请注意，尚不允许附加值。
- 第一个事件的时间（在本例中为timetoEvent）
- 暴露变量，在删失/事件之前的任何时间定义暴露（在这种情况下是exposure）。必须是（0）或（1）
- 暴露时间变量（在本例中为timetoexposure）
- 临床资格窗口（在这种情况下，我们将做100个时间单位或天）
让我们看看这在实践中是什么样子
# 创建克隆数据集
# create_clones()函数用于生成时间依赖性分析所需的克隆数据
clones <- create_clones(
# df: 输入的原始数据框
df = dummy_data,
# id: 指定个体ID列名
id = 'id',
# event: 指定事件状态列名(0/1)
event = 'event',
# time_to_event: 指定从基线到事件发生的时间列名
time_to_event = 'timetoevent',
# exposure: 指定暴露状态列名(0/1)
exposure = 'exposure',
# time_to_exposure: 指定从基线到暴露发生的时间列名
time_to_exposure = 'timetoexposure',
# ced_window: 设置累积暴露时间窗口长度(时间单位)
ced_window = 100
)
# 展示克隆数据集的前几行
# 使用DT包的datatable函数创建交互式数据表格
head(clones) |>
DT::datatable(
# rownames=FALSE: 不显示行名
rownames = FALSE,
options = list(
# scrollX=TRUE: 允许水平滚动
scrollX = TRUE,
# paging=FALSE: 不分页显示
paging = FALSE
)
)
请注意，此对象只是一个data. frame，带有一个额外的自定义类，未来的函数将评估该类：
class(clones)
根据文档内容，我来解释 create_clones 函数的数据转换过程:
# create_clones 函数的主要步骤:
# 1. 为每个病人创建两个克隆(clone=0和clone=1)
df_0 <- df_1 <- df # 复制原始数据框两份
# 2. 初始化新的列
df_0$outcome <- df_1$outcome <- rep(NA_integer_, NROW(df)) # 结局列
df_0$fup_outcome <- df_1$fup_outcome <- rep(NA_real_, NROW(df)) # 随访时间列
df_0$censor <- df_1$censor <- rep(NA_integer_, NROW(df)) # 删失状态列
df_0$fup_censor <- df_1$fup_censor <- rep(NA_real_, NROW(df)) # 删失时间列
# 3. 添加克隆标识
df_0$clone <- 0L # 未暴露组
df_1$clone <- 1L # 暴露组
# 4. 合并数据并添加类和属性
df_clones <- rbind(df_0, df_1) # 合并两个克隆组
class(df_clones) <- c("ccw_clones", class(df_clones)) # 添加ccw_clones类
# 5. 传递重要属性
attributes(df_clones)$id <- id # 病人ID
attributes(df_clones)$event <- event # 事件
attributes(df_clones)$time_to_event <- time_to_event # 事件时间
attributes(df_clones)$exposure <- exposure # 暴露状态
主要转换过程是:
对每个病人创建两个克隆(clone=0和clone=1),分别代表未暴露和暴露状态
-
为每个克隆添加以下新列:
- outcome: 结局状态
- fup_outcome: 随访时间
- censor: 删失状态
- fup_censor: 删失时间
合并两个克隆组的数据,并添加ccw_clones类和相关属性
-
这种转换允许我们:
- 比较同一病人在暴露和未暴露状态下的结局差异
- 计算因果权重
- 进行生存分析
这种克隆权重方法可以帮助我们评估暴露对结局的因果效应。
您可以在创建克隆后随着时间的推移可视化删失情况分布：
# 绘制随时间变化的审查情况图
# plot_censoring_over_time()函数用于可视化展示数据中随时间推移的审查模式
# 参数clones: 包含时间依赖性数据的克隆数据集对象
plot_censoring_over_time(clones)
转换为长格式
现在我们只需要将数据转换为长格式。cast_to_long函数将为我们执行此操作。不需要额外的参数（clones对象是一个工件，可以让您更好地查看和理解该方法）：
clones_long <- cast_clones_to_long(clones)
head(clones_long, row.names = FALSE) |>
DT::datatable(
rownames = FALSE,
options = list(
scrollX = TRUE,
paging = FALSE
)
)
让我详细解释 cast_clones_to_long 函数的数据转换过程：
# cast_clones_to_long 函数主要步骤:
# 1. 获取所有事件时间点
event_times <- sort(unique(c(
df[,time_to_event], # 结局事件时间
df[,time_to_exposure], # 暴露事件时间
ced_window # 临床资格窗口时间
)))
# 2. 分别处理未暴露组(clone=0)和暴露组(clone=1)的数据
# 以未暴露组为例:
df_0 <- df[df[, "clone"] == 0L, ] # 获取未暴露组数据
df_0$t_start <- 0 # 设置起始时间为0
# 3. 对每个结局类型进行迭代处理
for (i in 1:NROW(unique(df_0$outcome))) {
df_outcome <- df_0[df_0$outcome == unique(df_0$outcome)[i], ]
df_outcome$outcome <- 1L
# 使用 survival::survSplit 函数将数据展开成多个时间区间
df_0_long_outcome <- survival::survSplit(
data = df_outcome,
cut = event_times, # 在每个事件时间点切分
end = "fup_outcome", # 结局随访时间
start = "t_start", # 开始时间
event = "outcome" # 结局事件
)
}
# 4. 处理删失数据
df_0_long_censor <- survival::survSplit(
df_0,
cut = event_times,
end = "fup_outcome",
start = "t_start",
event = "censor"
)
# 5. 合并结局和删失数据
df_0_long_outcome$censor <- df_0_long_censor$censor
df_0_long_outcome$t_stop <- df_0_long_outcome$fup_outcome
# 6. 最后合并暴露组和未暴露组的数据
df_long <- rbind(
df_1_long, # 暴露组长格式数据
df_0_long # 未暴露组长格式数据
)
举个具体例子：
# 原始数据(宽格式)
df_original <- data.frame(
id = 1,
exposure = 1,
time_to_event = 10,
time_to_exposure = 5
)
# 转换后的数据(长格式)
df_long <- data.frame(
id = rep(1, 4),
clone = rep(c(0,1), each=2), # 0=未暴露,1=暴露
t_start = c(0,5, 0,5), # 时间区间起点
t_stop = c(5,10, 5,10), # 时间区间终点
outcome = c(0,1, 0,1), # 在每个时间区间的结局状态
censor = c(0,0, 0,0) # 在每个时间区间的删失状态
)
主要转换特点：
将每个病人的数据按时间点切分成多个区间
-
每个区间包含:
- 起始时间(t_start)
- 结束时间(t_stop)
- 该区间内的结局状态
- 该区间内的删失状态
-
这种转换允许我们:
- 精确计算每个时间点的风险
- 处理时依赖性协变量
- 进行复杂的生存分析
这种长格式数据结构更适合进行生存分析，因为它能够捕捉到随时间变化的信息。
让我详细解释时间点区间的设定过程：
# 假设我们有这样的原始数据
df <- data.frame(
id = c(1, 2),
time_to_event = c(10, 15), # 结局事件时间
time_to_exposure = c(5, 8), # 暴露事件时间
ced_window = 3 # 临床资格窗口时间
)
# 1. 首先收集所有关键时间点
event_times <- sort(unique(c(
df$time_to_event, # 结局事件时间: 10, 15
df$time_to_exposure, # 暴露事件时间: 5, 8
df$ced_window # 临床资格窗口: 3
)))
# event_times 结果为: c(3, 5, 8, 10, 15)
# 2. 对每个病人，根据这些时间点切分数据
# 以病人1为例，转换后的长格式数据：
df_long_patient1 <- data.frame(
id = 1,
# 时间区间
t_start = c(0, 3, 5, 8, 10), # 区间起点
t_stop = c(3, 5, 8, 10, 15), # 区间终点
# 在每个区间内的状态
exposure = c(0, 0, 1, 1, 1), # 0=未暴露, 1=暴露
outcome = c(0, 0, 0, 0, 1), # 0=无事件, 1=有事件
censor = c(0, 0, 0, 0, 0) # 0=未删失, 1=删失
)
时间区间设定的关键点：
-
关键时间点的来源：
- time_to_event: 结局事件发生时间
- time_to_exposure: 开始暴露的时间
- ced_window: 临床资格窗口时间（通常是观察开始的一段时间）
-
区间划分原则：
# 每个区间都是左闭右开的，例如： [0, 3) # 观察开始到临床资格窗口 [3, 5) # 临床资格窗口到暴露时间 [5, 8) # 暴露后的第一个区间 [8, 10) # 暴露后到结局事件前 [10, 15] # 结局事件发生的区间 -
状态记录规则：
# exposure（暴露状态） - 在暴露时间点之前 = 0 - 在暴露时间点之后 = 1 # outcome（结局状态） - 在结局事件时间之前 = 0 - 在结局事件发生时 = 1 # censor（删失状态） - 未删失期间 = 0 - 发生删失时 = 1 -
实际应用中的考虑：
# 时间区间的特点 - 每个区间的起点(t_start)是上一个区间的终点 - 区间长度不需要相等 - 区间数量取决于关键事件的发生次数 # 状态变化的记录 - 状态的改变只发生在区间的边界点 - 每个区间内的状态保持不变 - 可以精确捕捉状态改变的时间点
这种时间区间设定的优势：
- 可以精确追踪状态随时间的变化
- 便于计算时间依赖性的风险和概率
- 适合进行复杂的生存分析
- 能够处理时变协变量
这种设计使得我们可以：
- 准确估计在任意时间点的风险
- 计算累积风险和生存概率
- 评估时间依赖性暴露的影响
- 进行各种生存分析方法的建模
生成权重
现在，我们只需要生成权重。
clones_long_weights <- generate_ccw(clones_long, predvars = c("cov1", "cov2"))
head(clones_long_weights) |>
DT::datatable(
rownames = FALSE,
options = list(
scrollX = TRUE,
paging = FALSE
)
)
您还可以使用plot_ccw_over_time（）可视化权重，使用plot_var_mean_over_time（）可视化平均值。
clones_long_weights |>
plot_ccw_over_time()
clones_long_weights |>
plot_var_mean_over_time("cov1")
评估结果
我们现在拥有进行CCW分析所需的一切。例如，我们可以将东西组合在一起来评估手术与不手术的风险比：
# 加载survival包,用于生存分析
library(survival)
# 创建分析数据集
# 使用管道操作符依次处理数据
df <- dummy_data |>
# 使用create_clones()创建克隆数据集
# id: 个体ID变量名
# event: 事件变量名
# time_to_event: 到事件发生的时间变量名
# exposure: 暴露变量名
# time_to_exposure: 到暴露发生的时间变量名
# ced_window: 设置时间窗口为100
create_clones(
id = 'id',
event = 'event',
time_to_event = 'timetoevent',
exposure = 'exposure',
time_to_exposure = 'timetoexposure',
ced_window = 100
) |>
# 将克隆数据转换为长格式
cast_clones_to_long() |>
# 生成克隆权重,使用cov1和cov2作为预测变量
generate_ccw(c('cov1', 'cov2'))
# 使用Cox比例风险模型进行分析
# Surv()创建生存对象,使用t_start和t_stop定义时间区间,outcome为事件状态
# clone为主要预测变量
# weights = weight_cox使用生成的克隆权重
coxph(Surv(t_start, t_stop, outcome) ~ clone, data = df, weights = weight_cox)
请注意，我们在coxph（）模型中使用了结果而不是事件。尽管如此，这种分析当然有问题，因为克隆过程会使方差无效。解决这个问题的最简单方法是引导方差。我还没有创建一个函数来做到这一点，但是把下面的例子作为如何做到这一点的例子。
# 加载boot包,用于进行Bootstrap分析
library(boot)
# 定义Bootstrap分析的主函数
# data: 输入的原始数据集
# indices: Bootstrap重抽样的索引
boot_cox <- function(data, indices) {
# 使用Bootstrap样本创建带权重的长格式数据框
# data[indices,]: 根据indices抽取Bootstrap样本
ccw_df <- data[indices, ] |>
# 创建克隆数据集
create_clones(
id = 'id', # 个体ID变量
event = 'event', # 事件状态变量
time_to_event = 'timetoevent',# 到事件时间变量
exposure = 'exposure', # 暴露状态变量
time_to_exposure = 'timetoexposure', # 到暴露时间变量
ced_window = 100 # 时间窗口设为100
) |>
# 转换为长格式数据
cast_clones_to_long() |>
# 生成克隆权重
generate_ccw(c('cov1', 'cov2'))
# 使用Cox比例风险模型计算风险比(HR)
cox_ccw <- coxph(Surv(t_start, t_stop, outcome) ~ clone,
data = ccw_df,
weights = weight_cox)
# 提取并计算HR值
hr <- cox_ccw |>
coef() |> # 提取系数
exp() # 取指数得到HR
# 初始化输出向量
out <- c("hr" = hr)
# 分别为治疗组(clone=1)和对照组(clone=0)创建生存对象
# 使用survfit函数拟合KM曲线
surv_1 <- survfit(Surv(t_start, t_stop, outcome) ~ 1L,
data = ccw_df[ccw_df$clone == 1, ],
weights = weight_cox)
surv_0 <- survfit(Surv(t_start, t_stop, outcome) ~ 1L,
data = ccw_df[ccw_df$clone == 0, ],
weights = weight_cox)
# 计算限制平均生存时间(RMST)差异
# 设定365天为截止时间点
rmst_1 <- surv_1 |>
summary(rmean = 365) |> # 计算365天内的平均生存时间
(\(summary) summary$table)() |> # 提取结果表
(\(table) table["rmean"])() # 提取rmean值
rmst_0 <- surv_0 |>
summary(rmean = 365) |>
(\(summary) summary$table)() |>
(\(table) table["rmean"])()
# 计算RMST差异
rmst_diff <- rmst_1 - rmst_0
# 将RMST差异添加到输出向量
out <- c(out, "rmst_diff" = rmst_diff)
# 计算1年生存率差异
# 找到最接近365天的时间点索引
index_1yr_1 <- which.min(abs(surv_1$time - 365))
index_1yr_0 <- which.min(abs(surv_0$time - 365))
# 获取1年生存概率
surv_1_1yr <- surv_1$surv[index_1yr_1] # 治疗组1年生存率
surv_0_1yr <- surv_0$surv[index_1yr_0] # 对照组1年生存率
# 计算生存率差异
surv_diff_1yr <- surv_1_1yr - surv_0_1yr
# 将1年生存率差异添加到输出向量
out <- c(out, "surv_diff_1yr" = surv_diff_1yr)
}
# 执行Bootstrap分析
# data: 使用dummy_data数据集
# statistic: 使用boot_cox函数
# R: Bootstrap重复次数为10次
boot_out <- boot(data = dummy_data, statistic = boot_cox, R = 10)
# 计算Cox回归系数的Bootstrap置信区间
boot.ci(boot_out, # Bootstrap分析结果对象
type = "norm", # 使用正态分布方法计算置信区间
index = 1) # 提取第1个统计量(Cox回归系数)的置信区间
# 计算RMST差异的Bootstrap置信区间
boot.ci(boot_out, # Bootstrap分析结果对象
type = "norm", # 使用正态分布方法计算置信区间
index = 2) # 提取第2个统计量(RMST差异)的置信区间
# 计算1年生存率差异的Bootstrap置信区间
boot.ci(boot_out, # Bootstrap分析结果对象
type = "norm", # 使用正态分布方法计算置信区间
index = 3) # 提取第3个统计量(1年生存率差异)的置信区间
竞争风险
将包功能扩展到竞争风险问题是微不足道的。我们将结果变量保留为整数，但添加了附加值。
dummy_data$competing <- sample(0:2, nrow(dummy_data), replace = TRUE)
head(dummy_data) |>
DT::datatable(
rownames = FALSE,
options = list(
scrollX = TRUE,
paging = FALSE
)
)
compete_ccw <- dummy_data |>
create_clones(
id = 'id',
event = 'competing',
time_to_event = 'timetoevent',
exposure = 'exposure',
time_to_exposure = 'timetoexposure',
ced_window = 100
) |>
cast_clones_to_long() |>
generate_ccw(c('cov1', 'cov2'))
table(compete_ccw$outcome)
head(compete_ccw) |>
DT::datatable(
rownames = FALSE,
options = list(
scrollX = TRUE,
paging = FALSE
)
)
裁剪权重极端值
您可以使用winsorize_ccw_weights（）函数修剪/取消排序极端权重：
compete_ccw_trim <- compete_ccw |>
winsorize_ccw_weights(quantiles = c(0.10, 0.90))
max(compete_ccw$weight_cox)
max(compete_ccw_trim$weight_cox)
References
Hernán, Miguel A., and James M. Robins. “Using big data to emulate a target trial when a randomized trial is not available.” American journal of epidemiology 183.8 (2016): 758-764.
Gaber, Charles E., et al. “The Clone-Censor-Weight Method in Pharmacoepidemiologic Research: Foundations and Methodological Implementation.” Current Epidemiology Reports (2024): 1-11.
Maringe, Camille, et al. “Reflection on modern methods: trial emulation in the presence of immortal-time bias. Assessing the benefit of major surgery for elderly lung cancer patients using observational data.” International journal of epidemiology 49.5 (2020): 1719-1729.
本文由mdnice多平台发布
