医学研究中的生存数据建模
生存函数
\[
S(t)=P(T>t)=1-F(t)=\int_{t}^{+\infty}f(x)dx
\]
其中S(t)是累计生存概率或生存率，量化了生存时间大于t的概率。f(x)是密度函数，呈右偏态分布，反映了任意时间点 t 终点事件的瞬时发生率。F(t)=P(T<t)是f(t)在区间[0,t]的累计形式，也称为分布函数或累积函数。
风险函数
\[
\begin{aligned}
h(t)=&\lim_{\Delta t\to 0}\frac{P(t\le T<t+\Delta t|T\ge t)}{\Delta t}\\
=&\lim_{\Delta t\to 0}\frac{P(t\le T<t+\Delta t\ \&\ T\ge t)}{\Delta t·P(T\ge t) }\\
=&\lim_{\Delta t\to 0}\frac{S(t)- S(t+\Delta t)}{\Delta t·S(t)}\\
=&-\frac{d(\ln S(t))}{dt}
\end{aligned}
\]
\[S(t)=e^{-\int_0^t h(u)du} \]
\[
\begin{aligned}
h(t)\Delta t=&P(t\le T<t+\Delta t|T\ge t)
=\frac{P(t\le T<t+\Delta t\ \&\ T\ge t)}{P(T\ge t) }\\
=&\frac{P(t\le T<t+\Delta t)}{P(T\ge t) }\\
=&\frac{f(t)\Delta t}{S(t)}
\end{aligned}
\]
\(f(t)=h(t)S(t)\)
乘积极限法
product limit method 也称为Kaplan-Meier 法。
\(t_1<t_2<t_3<...<t_n\)，样本量大小n，ti 代表个体i发生终点事件或右删失的时间。由于一些个体有相同的生存时间，它们被称为 tied 观测时间，生存时间的个数小于样本量n。
点估计S(t)
\(n_1>n_2>n_3>...>n_n\) ,ni d代表在时间点ti 暴露于特定事件风险的幸存者数量。
\(d_i\) 代表在时间点ti 发生终点事件的数量。（如果没有 tie，di=1或0）
生存率的KM估计计算公式：
\[
\hat S(t)=\prod_{t_i\le t}\frac{n_i-d_i}{n_i}
\tag{21.1}\]
Equation 21.1 包括了删失情况，如果从ti-1 到ti 发生了删失，但没有终点事件，di =0，条件概率等于1。
区间估计S(t)
(1-α)×100% CI \([\hat S(t)-z_{1-\alpha/2}\sqrt{Var\left [\hat S(t) \right]},\hat S(t)+z_{1-\alpha/2}\sqrt{Var\left [\hat S(t) \right]}]\)
其中\(Var\left [\hat S(t) \right]=\hat S(t)^2\sum_{t_i\le t}\frac{d_i}{n_i(n_i-d_i)}\) (Greenwood method )
Code
df_raw <- tibble(
id=1:5,
sex=factor(c("Male","Male","Female","Male","Female")),
age_years=c(53,62,75,73,61),
outcome=c("Loss to follow-up","Death","Death","Relapse","Survival"),
time=c("37.5+",44.3,25.3,18.1,"56.7+")
)
df <- df_raw |> arrange(time) |> select(id,outcome,time) |>
mutate(
d_i=if_else(outcome=="Relapse"|outcome=="Death",1,0),
time=parse_number(time),
)
df
| 4 |
Relapse |
18.1 |
1 |
| 3 |
Death |
25.3 |
1 |
| 1 |
Loss to follow-up |
37.5 |
0 |
| 2 |
Death |
44.3 |
1 |
| 5 |
Survival |
56.7 |
0 |
Code
library(survminer)
library(survival)
km_fit<-survfit(Surv(time,d_i)~1,data=df)
# t_i
km_fit$time
#> [1] 18.1 25.3 37.5 44.3 56.7
# d_i
km_fit$n.event
#> [1] 1 1 0 1 0
# cnesored
km_fit$n.censor
#> [1] 0 0 1 0 1
# n_i
km_fit$n.risk
#> [1] 5 4 3 2 1
# 生存率
km_fit$surv
#> [1] 0.8 0.6 0.6 0.3 0.3
summary(km_fit)
#> Call: survfit(formula = Surv(time, d_i) ~ 1, data = df)
#>
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 18.1 5 1 0.8 0.179 0.5161 1
#> 25.3 4 1 0.6 0.219 0.2933 1
#> 44.3 2 1 0.3 0.239 0.0631 1
单因素分组生存曲线的比较
| log-rank test |
1 |
| Wilcoxon test |
nj
|
| Tarone-Ware test |
\(\sqrt{n_j}\) |
| Peto test |
\(\hat S(t_j)\) |
\[
\chi^2=\frac{\left(\sum_jw(t_j)(m_{ij}-e_{ij})\right)^2}{\hat {Var}\left(\sum_jw(t_j)(m_{ij}-e_{ij})\right)}
\]
log-rank test
H0 :两总体的生存曲线是相同的。
-
计算当第j次发生终点事件各组终点事件的期望值（e1j ,e2j ）
\(e_{1j}=\left ( \frac{n_{1j}}{n_{1j}+n_{2j}}\right)\times (m_{1j}+m_{2j})\)
\(e_{2j}=\left ( \frac{n_{2j}}{n_{1j}+n_{2j}}\right)\times (m_{1j}+m_{2j})\)
其中mij 表示在第 j 个时间点第 i 组终点事件的数量，nij 表示在第 j 个时间点第 i 组初始观测的数量
-
对所有时间点对终点事件的观测值和期望值的差异求和
\(O_i-Ei=\sum_j(m_{ij}-e_{ij})\ \ \ (i=1,2)\)
计算其方差估计值
\(\hat{Var}=\sum_j\frac{n_{1j}n_{2j}(m_{1j}+m_{2j})(n_{1j}+n_{2j}-m_{1j}-m_{2j})}{(n_{1j}+n_{2j})^2(n_{1j}+n_{2j}-1)}\ \ \ (i=1,2)\)
-
计算log-rank test 的检验统计量
\[
\chi^2=\frac{(O_1-E_1)^2}{\hat{Var}(O_1-E_1)} \ 或者 \ \chi^2=\frac{(O_2-E_2)^2}{\hat{Var}(O_2-E_2)}
\]
也可以近似估计为
\[
\chi^2=\sum_{i=1}^2\frac{(O_i-E_i)^2}{E_i} \sim \chi^2(\nu=1)
\]
Code
df <- read_csv("data/log-rank-survival.csv")
# 使用 Kaplan-Meier 方法创建一个Surv对象
surv_formula <- Surv(Days,status)~treatment
km<-survfit(surv_formula,data=df)
summary(km)
#> Call: survfit(formula = surv_formula, data = df)
#>
#> treatment=CON
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 11 14 1 0.9286 0.0688 0.8030 1.000
#> 13 13 1 0.8571 0.0935 0.6921 1.000
#> 14 12 2 0.7143 0.1207 0.5129 0.995
#> 15 10 1 0.6429 0.1281 0.4351 0.950
#> 17 9 3 0.4286 0.1323 0.2341 0.785
#> 20 6 2 0.2857 0.1207 0.1248 0.654
#> 21 4 2 0.1429 0.0935 0.0396 0.515
#> 25 2 1 0.0714 0.0688 0.0108 0.472
#> 27 1 1 0.0000 NaN NA NA
#>
#> treatment=DPVB
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 20 14 1 0.929 0.0688 0.803 1.000
#> 23 13 1 0.857 0.0935 0.692 1.000
#> 27 12 1 0.786 0.1097 0.598 1.000
#> 28 11 1 0.714 0.1207 0.513 0.995
#> 30 10 1 0.643 0.1281 0.435 0.950
#> 32 9 1 0.571 0.1323 0.363 0.899
#> 38 8 1 0.500 0.1336 0.296 0.844
#> 39 7 1 0.429 0.1323 0.234 0.785
#> 45 6 1 0.357 0.1281 0.177 0.721
#>
#> treatment=LDRT
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 13 14 2 0.8571 0.0935 0.6921 1.000
#> 15 12 1 0.7857 0.1097 0.5977 1.000
#> 16 11 1 0.7143 0.1207 0.5129 0.995
#> 18 10 1 0.6429 0.1281 0.4351 0.950
#> 19 9 2 0.5000 0.1336 0.2961 0.844
#> 20 7 3 0.2857 0.1207 0.1248 0.654
#> 24 4 1 0.2143 0.1097 0.0786 0.584
#> 25 3 1 0.1429 0.0935 0.0396 0.515
#> 27 2 1 0.0714 0.0688 0.0108 0.472
#> 30 1 1 0.0000 NaN NA NA
#>
#> treatment=LR_DPVB
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 30 14 1 0.929 0.0688 0.803 1
#> 40 13 1 0.857 0.0935 0.692 1
# 执行Log-rank检验
logrank_test <- survdiff(surv_formula,data = df,subset = T,na.action = "na.omit")
logrank_test$chisq
#> [1] 58.43627
logrank_test$pvalue
#> [1] 1.268397e-12