landmark/竞争风险/加权KM绘制的一体化包-jskm

教程地址：https://github.com/jinseob2kim/jskm

背景知识

landmark分析

生存分析中的 landmark 分析是一种特殊的生存分析方法，它的主要目的是在指定的一些固定时间点上评估患者的生存状况。这些固定时间点被称为 'landmark' 时间点。与传统的生存分析方法不同，landmark 分析可以帮助解决一些特定问题，并提供更加详细的信息。

Landmark 分析的意义和应用包括：

1. 处理时间依赖性问题：在生存分析中，某些因素可能在初始时期对生存率产生更大影响，随着时间的推移逐渐减弱。Landmark 分析可以帮助解决时间依赖性问题，通过在不同时间点进行分析，更准确地捕捉因素的影响。

2. 评估中间事件影响：在某些研究中，中间事件（如治疗改变、疾病进展）可能影响到生存分析的结果。通过在 landmark 时间点上进行分析，可以在考虑中间事件的情况下评估生存情况。

3. 动态监测治疗效果：对于临床试验或长期随访研究，landmark 分析可以用来监测不同时间点的治疗效果，从而帮助医生更好地决定治疗策略。

4. 预测患者生存情况：通过在多个 landmark 时间点上进行分析，可以为患者提供更加个体化的生存预测，从而更好地指导临床决策。

5. 验证模型鲁棒性：landmark 分析可以用来验证生存预测模型在不同时间点上的鲁棒性，从而评估模型的预测能力。

总之，landmark 分析在生存分析中具有重要的意义，可以帮助研究人员更全面地理解和评估患者的生存情况，并为医疗决策和临床研究提供更详细的信息。

Weighted Kaplan-Meier plot

Weighted Kaplan-Meier plot 是一种基于加权的 Kaplan-Meier 生存曲线，用于分析不同组之间的生存差异，但每组的样本具有不同的权重。这种方法对于处理样本不平衡或者某些样本比例更具有代表性的情况非常有用。

竞争风险模型

竞争风险模型（Competing Risks Model）是一种在生存分析中用于处理多个相互竞争的事件的统计模型。在这种模型中，每个事件都被视为一个“风险”，而多个事件可能会同时发生，但只有一个事件被视为“失败”。竞争风险模型的主要目的是在考虑多个竞争事件的情况下，分析每个事件的发生概率以及它们对生存时间的影响。

竞争风险模型的意义和应用包括：

1. 更真实的数据建模：在某些研究中，一个人或一个物体可能面临多种竞争性风险，例如在医学研究中，患者可能会死于不同原因，如疾病复发、其他疾病或自然因素。竞争风险模型可以更准确地建模这些情况，避免简单地将其他事件视为“丧失随访”。

2. 事件之间的相互影响：竞争风险模型允许分析不同事件之间的相互影响，例如一个事件的发生可能会影响其他事件的发生概率。这种模型可以提供更全面的信息，以更好地理解各个事件之间的关系。

3. 更精细的分析：竞争风险模型可以分析每个事件的发生概率，以及事件发生对生存时间的影响。这有助于更精细地理解每个事件对研究对象的影响，而不仅仅关注主要事件。

4. 预测患者风险：在医学和生物统计学中，竞争风险模型可以用于预测患者在面临多种竞争性风险时的生存概率。这对于制定治疗策略和做出临床决策非常有价值。

总之，竞争风险模型在生存分析中的意义在于提供一种更全面、更准确地分析多个竞争性事件的统计工具，帮助研究人员更好地理解和预测在这些竞争性风险下的生存情况。

代码示例教程

常规生存及landmark分析

install.packages('jskm')

## From github: latest version
install.packages('remotes')
remotes::install_github('jinseob2kim/jskm')

library(jskm)

#基本用法
library(survival)
data(colon)
fit <- survfit(Surv(time,status)~rx, data=colon)

#Plot the data
jskm(fit)


初步图形

# Customize the jskm plot
jskm(
  fit,
  table = TRUE, pval = TRUE, label.nrisk = 'No. at risk', size.label.nrisk = 8,
  xlabs = 'Time (Day)', ylabs = 'Survival', ystratalabs = c('Obs', 'Lev', 'Lev   5FU'), ystrataname = 'rx',
  marks = FALSE, timeby = 365, xlims = c(0, 3000), ylims = c(0.25, 1), showpercent = TRUE
)

以下是一些参数的解释：

• table = TRUE：在绘图中显示生存表格。
• pval = TRUE：在绘图中显示统计显著性 p 值。
• label.nrisk = 'No. at risk'：设置风险人数标签的文本。
• size.label.nrisk = 8：设置风险人数标签的字体大小。
• xlabs = 'Time (Day)'：设置 x 轴标签的文本。
• ylabs = 'Survival'：设置 y 轴标签的文本。
• ystratalabs 和 ystrataname：用于自定义不同治疗组的标签。
• marks = FALSE：禁用标记（mark）在曲线上的显示。
• timeby = 365：设置 x 轴的刻度间隔。
• xlims = c(0, 3000)：设置 x 轴的范围。
• ylims = c(0.25, 1)：设置 y 轴的范围。
• showpercent = TRUE：在 y 轴上显示百分比。

定制化图

# Fit the survival model
fit <- survfit(Surv(time, status) ~ rx, data = colon)

# Plot the Kaplan-Meier survival curves with customized features
jskm(
  fit,
  ci = TRUE,         # Display confidence intervals
  cumhaz = TRUE,     # Display cumulative hazard function
  mark = FALSE,      # Do not display marks on the plot
  ylab = 'Cumulative Incidence (%)',   # Customize y-axis label
  surv.scale = 'percent',   # Set survival scale to percentage
  pval = TRUE,       # Display p-values
  pval.size = 6,     # Set p-value font size
  pval.coord = c(300, 0.7)   # Set coordinates for displaying p-value
)


在这段代码中，我们调用了 jskm() 函数，并使用参数来定制绘图的外观和内容。以下是每个参数的详细解释：

• ci = TRUE：在绘图中显示置信区间。
• cumhaz = TRUE：在绘图中显示累积危险函数。
• mark = FALSE：不在曲线上显示标记（mark）。
• ylab = 'Cumulative Incidence (%)'：自定义 y 轴标签。
• surv.scale = 'percent'：将生存刻度设置为百分比。
• pval = TRUE：显示 p 值。
• pval.size = 6：设置 p 值的字体大小。
• pval.coord = c(300, 0.7)：设置在绘图中显示 p 值的坐标。

累计风险绘图

# Plot the Kaplan-Meier survival curves with customized features
jskm(
  fit,
  mark = FALSE,       # Do not display marks on the plot
  surv.scale = 'percent',   # Set survival scale to percentage
  pval = TRUE,         # Display p-values
  table = TRUE,        # Display the survival table
  cut.landmark = 500   # Set the landmark time point to 500
)

• mark = FALSE：不在曲线上显示标记（mark）。
• surv.scale = 'percent'：将生存刻度设置为百分比。
• pval = TRUE：显示 p 值。
• table = TRUE：显示生存表格。
• cut.landmark = 500：设置 landmrak 时间点为 500。

固定时间点landmark分析


# Plot the Kaplan-Meier survival curves with customized features
jskm(
  fit,
  mark = FALSE,        # Do not display marks on the plot
  surv.scale = 'percent',    # Set survival scale to percentage
  pval = TRUE,          # Display p-values
  table = TRUE,         # Display the survival table
  cut.landmark = 500,   # Set the landmark time point to 500
  showpercent = TRUE    # Display percentages on the plot
)


在这段代码中，我们调用了 jskm() 函数，并使用参数来定制绘图的外观和内容。以下是每个参数的详细解释：

• mark = FALSE：不在曲线上显示标记（mark）。
• surv.scale = 'percent'：将生存刻度设置为百分比。
• pval = TRUE：显示 p 值。
• table = TRUE：显示生存表格。
• cut.landmark = 500：设置 landmark 时间点为 500。
• showpercent = TRUE：在绘图上显示百分比。

这些参数可以帮助你更好地控制绘制出的图形的外观和内容，使其更符合你的需求。根据你的实际需要，你可以进一步调整这些参数来获得所期望的绘图效果。

增加table后的landmark分析

Competing risk analysis（竞争风险模型）

status2 variable: 0 - censoring, 1 - event, 2 - competing risk 以下是详细注释的代码：

# 假设你有一个数据框 colon 包含时间、状态和治疗组信息
# Make competing risk variable (not real)
# 在这个示例中，我们创建一个名为 status2 的新变量，用于模拟竞争风险变量
colon$status2 <- colon$status
colon$status2[1:400] <- 2  # 假设前 400 个观察属于竞争事件 2
colon$status2 <- factor(colon$status2)  # 将 status2 转换为因子变量

# 使用 survfit 计算 competing risk 情况下的 Kaplan-Meier 生存曲线
fit2 <- survfit(Surv(time, status2) ~ rx, data = colon)

# 绘制竞争风险的 Kaplan-Meier 生存曲线，同时显示生存表格
jskm(
  fit2,
  mark = FALSE,          # 不显示标记
  surv.scale = 'percent',  # 将生存刻度设置为百分比
  table = TRUE,           # 显示生存表格
  status.cmprsk = '1'     # 设置竞争风险事件为 1
)

在这段代码中，我们创建了一个模拟的竞争风险变量 status2，将前 400 个观察标记为竞争事件 2，然后使用 survfit() 计算针对这个竞争风险变量的 Kaplan-Meier 生存曲线。接着，我们使用 jskm() 函数绘制了竞争风险的 Kaplan-Meier 生存曲线，并设置了一些参数来控制绘图的外观和内容，包括不显示标记、将生存刻度设置为百分比、显示生存表格以及设置竞争风险事件为 1。请注意，这里使用的竞争风险变量 status2 是一个模拟的示例，实际应用中需要根据你的数据和研究问题来定义和计算竞争风险变量。

竞争风险模型的KM曲线

下面是详细注释的代码：


# 绘制竞争风险的 Kaplan-Meier 生存曲线，同时显示生存表格
jskm(
  fit2,
  mark = FALSE,              # 不显示标记
  surv.scale = 'percent',    # 将生存刻度设置为百分比
  table = TRUE,              # 显示生存表格
  status.cmprsk = '1',       # 设置竞争风险事件为 1
  showpercent = TRUE,        # 在绘图上显示百分比
  cut.landmark = 500         # 设置 landmark 时间点为 500
)


竞争风险显示risk table

Weighted Kaplan-Meier plot - svykm.object in survey package（加权KM）

# 加载所需的包
library(survey)
data(pbc, package = 'survival')

# 在数据框 pbc 中添加一个随机化变量 randomized，用于指示是否随机分组
pbc$randomized <- with(pbc, !is.na(trt) & trt > 0)

# 使用广义线性模型 glm() 拟合随机化模型，其中预测变量是 age 和 edema
biasmodel <- glm(randomized ~ age * edema, data = pbc)

# 在数据框 pbc 中添加一个随机化概率变量 randprob，使用拟合的模型的预测值
pbc$randprob <- fitted(biasmodel)

# 创建一个复杂抽样设计对象 dpbc，用于进行加权分析
dpbc <- svydesign(
  id = ~1,           # 指定一个虚拟的标识符变量，因为我们只需要一个观察一个群组
  prob = ~randprob,  # 指定权重变量为 randprob
  strata = ~edema,   # 指定分层变量为 edema
  data = subset(pbc, randomized)  # 仅使用随机分组的观察
)

# 使用 svykm() 函数计算加权的 Kaplan-Meier 生存曲线
s1 <- svykm(Surv(time, status > 0) ~ 1, design = dpbc)

# 使用 svykm() 函数计算加权的 Kaplan-Meier 生存曲线，按性别分组
s2 <- svykm(Surv(time, status > 0) ~ sex, design = dpbc)

# 使用 svyjskm() 函数绘制加权的 Kaplan-Meier 生存曲线
svyjskm(s1)

在这段代码中，我们首先加载了 survey 包，并使用 data() 函数加载了一个示例数据集 pbc。然后，我们通过拟合广义线性模型来创建一个随机化模型，计算随机分组的概率值，并将其添加到数据框中。接着，我们使用 svydesign() 函数创建了一个复杂抽样设计对象 dpbc，其中指定了权重变量、分层变量和数据。然后，我们使用 svykm() 函数计算了两组加权的 Kaplan-Meier 生存曲线。最后，我们使用 svyjskm() 函数绘制了加权的 Kaplan-Meier 生存曲线。

一个分组的KM曲线

svyjskm(s2, pval = T,  table = T, design = dpbc)


两个分组的加权KM

s3 <- svykm(Surv(time,status>0) ~ sex, design=dpbc, se = T)
svyjskm(s3)

带可信区间的加权KM

svyjskm(s3, ci = F)


SVY对象带SE的同时可以通过绘图参数指定不显示

svyjskm(s3, ci = F,  surv.scale = 'percent', pval =T, table = T, cut.landmark = 1000, showpercent = T)

加权KM的landark分析