生存分析

1. 基本概念

1.1 生存数据

当我们想分析受试者发生某一事件的时长时，需要2个因变量：是否发生结局事件，以及所经历的时间。我们把发生事件称为死亡，则发生了死亡status = 1，未发生死亡status = 0；从开始随访到发生死亡经历的时长称为生存时间。

生存数据的一个特点是生存时间存在删失，即并不是每个受试者都能知道死亡的时间。如某个研究时长3年，患者A在0时刻入组，第2年死亡，则生存时间是2年；患者B在0时刻入组，第1.5年脱落或失访，则生存时间是1.5+年；患者C在第1年入组，在第3年研究结束出组，则生存时间是2+年；患者D在第2年入组，在第2.5年因其他原因死亡，则生存时间是0.5+年。以上，患者B、C、D的生存时间都是（右）删失的，不知道具体值。

因此，只有在研究期间发生了死亡事件的受试者才是未删失数据，这些人知道具体生存时间；而其余的人（始终存活、脱落、因其他原因死亡）都是删失数据，虽然也有生存时间，但真实的时间应大于这个值。所以，对于删失的生存时间，不能简单地用删失时间点的值来代替，这样会低估生存时间。

综上，生存数据的特点是：（1）不仅要研究生存时间，还要考虑是否发生了事件；（2）事件是否发生导致生存时间出现了删失，需要对删失数据进行处理。

1.2 生存函数

用随机变量\(T\)表示患者的生存时间，其概率密度函数是\(f(t), t\geq 0\)，分布函数是

\[F(t)=P(T \leq t)=\int_{0}^{t}f(x)dx \]

生存分析的核心是研究生存时间，主要目的是对以下函数进行估计和建模：

生存函数\(S(t)\)：患者活到时间\(t\)的概率

\[S(t) = P(T>t) = \int_{t}^{\infty}f(x)dx=1-F(t), t \geq 0 \]

风险函数\(h(t)\)：患者在\(t\)时刻的瞬时死亡率，即患者活过了\(t\)却没有活过\(t+\Delta t\)

\[P(T<t+\Delta t | T>t)=\frac{P(t<T<t+\Delta t)}{P(T>t)}=\frac{f(t)dt}{S(t)}=h(t)dt \]

则

\[h(t)=\frac{f(t)}{S(t)} \]

累积风险函数\(H(t)\)：

\[H(t)=\int_{0}^{t}h(x)dx \]

由上可知，\(f(t)\)、\(F(t)\)、\(S(t)\)、\(h(t)\)、\(H(t)\)可以互相转化，因此，只要知道生存时间\(T\)的分布类型，就可以计算出这些值。由此，衍生出了参数法生存分析。

除了估计生存时间的分布外，还可以直接根据数据而不假定分布类型来估计生存函数，称为非参数法生存分析。

还有一种半参数法生存分析。

2. 非参数法

2.1 用Kaplan-meier法估计生存函数

使用KM法估计\(S(t)\)时无需假定任何的分布类型，需要使用KM估计量，又称乘积极限估计量。有\(k\)个不同的生存时间，按从小到大排序：\(t_1<t_2<...<t_k\)，在\(t_i\)时刻有\(n_i\)个人存活（number at risk），则KM估计量为

\[\hat{S}(t)=\prod_{i\ : \ t_i<t}(1-\frac{d_i}{n_i}), t \geq 0 \]

示例：
有10个受试者，生存时间分别是：24+、16+、8、19、10、8+、5、17、20、10，则有8个不同的时间点，并按大小排序，得：5、8、10、16、17、19、20、24，得到这8个时间点的KM估计值为

R语言实现：

time <- c(24, 16, 8, 19, 10, 8, 5, 17, 20, 10)
status<- c(0, 0, 1, 1, 1, 0, 1, 1, 1, 1)
data <- data.frame(ID, time, status)

ggsurvplot(survfit(Surv(time, status) ~ 1), data = data,   
           risk.table = T, censor = T, conf.int = F)

2.2 用log-rank检验比较生存函数

检验2个组别的生存函数是否有显著差异，假设形式如下：

\[H_0:S_1(t)=S_2(t) \ \ \ for \ all \ t \]

\[H_1:S_1(t)\neq S_2(t) \ \ \ for \ some \ t \]

最常用的是log-rank检验。将未删失的生存时间排序：\(t_1<t_2<...<t_k\)，\(n_i\)表示\(t_i\)时刻仍然存活的人数，\(d_i\)表示恰好在\(t_i\)死亡的人数，则生存或死亡以及不同的组别构成了2*2的表格：

上表中，任意一个象限的值确定后，则其余三个象限的值也随之确定。考虑\(d_{1i}\)，根据超几何分布，可以得到随机变量\(d_{1i}\)的期望：

\[{\rm E} (d_{1i})=\frac{n_{1i} \cdot d_i}{n_i} \]

将所有的\(i\)求和，得正态分布统计量

\[U = \sum_{i=1}^{k}(d_{1i}-{\rm E}(d_{1i})) \]

将\(U\)标准化后得到标准正态分布统计量

\[z = \frac{U}{\sqrt{{\rm Var}(U)}} \]

\(z^2\)即为log-rank检验的检验统计量，服从自由度为1的卡方分布。

示例：
治疗组5例受试者的生存时间是：3.4、3.6+、4.1、4.9、5.8+；对照组的4例受试者生存时间是：2、3.7+、4.3、4.9。未删失的时间是：2、3.4、4.1、4.3，各时间点的\(d_{1i}\)分别是：

则\(U = (0-5/9)+(1-5/8)+(1-3/5)+(0-1/2)=-0.2806\)，\({\rm Var}(U)=20/81+15/64+6/25+1/4=0.9713\)；\(z^2 = (-0.2806/\sqrt{0.9713})^2=0.081\)，P值为0.7759，未能拒绝原假设，两组没有显著差异。

R语言实现：

time <- c(3.4, 3.6, 4.1, 4.9, 5.8, 2, 3.7, 4.3, 4.9)
status<- c(1, 0, 1, 0, 0, 1, 0, 1, 0)
group <- c("treatment", "treatment", "treatment", "treatment", "treatment",
           "control", "control", "control", "control")
data <- data.frame(ID, time, status, group)

survdiff(Surv(time, status) ~ group, data = data)

ggsurvplot(survfit(Surv(time, status) ~ group), data = data, 
           conf.int = F,  pval = T, risk.table = T, censor = T)

# Call:
#   survdiff(formula = Surv(time, status) ~ group, data = data)
# 
#                 N Observed Expected (O-E)^2/E (O-E)^2/V
# group=control   4        2     1.72    0.0458     0.081
# group=treatment 5        2     2.28    0.0345     0.081
# 
# Chisq= 0.1  on 1 degrees of freedom, p= 0.8 

3. 参数法

3.1 指数分布

我们可以为生存时间指定一种已知的分布类型，估计相应分布函数中的参数，从而得到生存函数，这种方法称为参数法。

以指数分布为例。指数分布是最简单的生存时间分布模型，是一种纯随机死亡模型，在任何时间的风险函数都是一个常数，具有“无记忆性”。指数分布的概率密度函数是

\[f(t)=\lambda e^{-\lambda t} \]

分布函数是

\[F(t)=1-e^{-\lambda t} \]

则生存函数是

\[S(t)=1-F(t)=e^{-\lambda t} \]

风险函数是

\[h(t)=\frac{f(t)}{S(t)}=\lambda \]

指数分布中只有一个参数\(\lambda\)，称为风险率，是一个尺度参数。它决定了生存时间的长短，\(\lambda\)越大，生存率下降越快，生存时间越短。可使用极大似然估计法来估计\(\lambda\)的值。

R语言实现

time <- c(1, 1, 2, 3, 3, 5, 8, 10, 16, 18)
status<- c(1, 1, 1, 0, 1, 1, 1, 1, 0, 1)
data1 <- data.frame(time, status)

KM <- ggsurvplot(survfit(Surv(time, status) ~ 1, data = data1), conf.int = F) 

fit <- survreg(Surv(time, status) ~ 1, data = data1, dist = "exponential")

prob <- seq(1, 0, by = -0.01)
pred <- data.frame(t = predict(fit, type = "quantile", p = 1 - prob)[1,], prob)

KM$plot + geom_line(data = pred, aes(t, prob))

3.2 回归模型

与非参数法相比，参数法的优势是可以进行协变量分析。设生存时间\(T\)有\(m\)个影响因素，可以写成如下的回归模型：

\[\ln(T) = \beta_0 + \beta_1x_1+ ... +\beta_mx_m + \varepsilon \]

即

\[T = e^{\beta_0 + \beta_1x_1+ ... +\beta_mx_m + \varepsilon} \]

将\(T\)取对数可以将其值域从\([0, +\infty)\)扩展到\((-\infty, +\infty)\)。

以指数分布为例，生存时间的概率密度函数为

\[f(t)=\lambda e^{-\lambda t}, t \geq 0 \]

生存函数为

\[S(t)=e^{-\lambda t}, t\geq 0 \]

其中\(\lambda = e^{-(\beta_0 + \beta_1x_1+ ... +\beta_mx_m)}\)，各影响因素的回归系数\(\beta_0, \beta_1, ..., \beta_m\)可以通过极大似然估计得到。

如何理解回归系数

\(T\)的期望是\({\rm E}(T)=\frac{1}{\lambda}=e^{\beta_0 + \beta_1x_1+ ... +\beta_mx_m}\)。
对于连续型协变量\(x_i\)，每增加一个单位，与原来相比

\[\frac{{\rm E}(T_{x_i+1)}}{{\rm E}(T_{x_i})} = \frac{e^{\beta_0 +\beta_i(x_i+1)}}{e^{\beta_0 +\beta_ix_i}}=e^{\beta_i} \]

对于分类型协变量，类别\(i\)的回归系数是\(\beta_i\)，类别\(j\)的回归系数是\(\beta_j\)，则

\[\frac{{\rm E}(T_{y_i)}}{{\rm E}(T_{y_j})} = \frac{e^{\beta_0 +\beta_i}}{e^{\beta_0 + \beta_j}}=e^{\beta_i-\beta_j} \]

如果与baseline相比，则

\[\frac{{\rm E}(T_{y_i)}}{{\rm E}(T_{y_0})} = \frac{e^{\beta_0 +\beta_i}}{e^{\beta_0}}=e^{\beta_i} \]

R语言实现

# K-M法
km <- ggsurvplot(survfit(Surv(time, status) ~ sex, data = lung)) 

# 指数分布
fit <- survreg(Surv(time, status) ~ sex, data = lung, dist = "exponential")
prob <- seq(1, 0, by = -0.01)
pred1 <- data.frame(t = predict(fit, newdata = list(sex = 1), type = "quantile", p = 1 - prob), prob)
pred2 <- data.frame(t = predict(fit, newdata = list(sex = 2), type = "quantile", p = 1 - prob), prob)

km$plot + geom_line(data = pred1, aes(t, prob)) + geom_line(data = pred2, aes(t, prob))

3.3 其他分布

Weibull分布是指数分布的推广形式，概率密度函数为

\[f(t)=\lambda \gamma t^{\alpha-1} e^{-\lambda t^{\gamma}} \]

分布函数为

\[F(t)=1- e^{-\lambda t^{\gamma}} \]

生存函数为

\[S(t)=e^{-\lambda t^{\gamma}} \]

风险函数为

\[h(t)=\lambda \gamma t^{\gamma-1} \]

其中，\(\lambda\)称为尺度参数，决定分布的分散度；\(\gamma\)称为形状参数，决定分布的形状。当\(\gamma <1\)时，风险函数随时间递减；当\(\gamma >1\)时，风险函数随时间递增；当\(\gamma =1\)时，Weibull分布简化为指数分布。

R语言survival包的survreg函数支持的分布有：extreme, logistic, gaussian, weibull, exponential, rayleigh, loggaussian, lognormal, loglogistic, t。

R语言flexsurv包的flexsurvreg函数支持的分布有：genf, genf.orig, gengamma, gengamma.orig, exp, weibull, weibullph, lnorm, gamma, gompertz, llogis, exponential, lognormal。

library(flexsurv)

# K-M法
km <- ggsurvplot(survfit(Surv(time, status) ~ sex, data = lung)) 

# gamma分布
fit <- flexsurvreg(Surv(time, status) ~ sex, data = lung, dist = "gamma")

prob <- seq(0.999, 0.001, by = -0.01)
pred <- predict(fit, newdata = data.frame(sex = c(1, 2)), type = "quantile", p = prob)
pred1 <- pred[[1]][[1]] 
pred2 <- pred[[1]][[2]]

km$plot + 
  geom_line(data = pred1, aes(`.pred_quantile`, 1 - `.quantile`)) + 
  geom_line(data = pred2, aes(`.pred_quantile`, 1 - `.quantile`))

4. 半参数法

4.1 Cox比例风险模型

参数法考虑的是生存函数\(S(t)\)与协变量\(x_i\)的关系，风险函数\(h(t)\)取决于时间\(t\)和依赖于时间的协变量\(x_i(t)\)。Cox提出了一种简化模型，建立\(h(t)\)与\(x_i\)的关系，称为Cox比例风险模型：

\[h(t, x_i, \beta_i)=h_0(t)e^{\sum x_i \beta_i},\ \ i= 1,...,m \]

其中，\(h_0(t)\)称为基线风险函数。由上式可知，风险函数中只有\(h_0(t)\)与时间有关，而后面协变量的部分与时间无关。任意两个受试者的风险函数之比都是一个不取决于时间的常数，由此得名“比例风险”：

\[\frac{h(t,x_{i1},...,x_{im},\beta_1,...,\beta_m)}{h(t,x_{j1},...,x_{jm},\beta_1,...,\beta_m)}=e^{\beta_1(x_{i1}-x_{j1})+...+\beta_m(x_{im}-x_{jm})} \]

在Cox模型中，\(h_0(t)\)和回归系数\(\beta_i\)都是未知的。参数法需要假定生存时间的分布类型，而Cox模型未对\(h_0(t)\)作任何假定。在许多场景（协变量考察）下，我们只需要知道回归系数\(\beta_i\)，即使\(h_0(t)\)未知也不会影响\(\beta_i\)的估计。因此，Cox模型尽管含有参数\(\beta_i\)，但也含有\(h_0(t)\)这一部分，不是完全的参数模型，属于半参数模型。

4.2 风险比

两个风险函数之比称为风险比（\(HR\)）。对于连续型协变量\(x_i\)，每增加一个单位，与原来相比

\[HR=\frac{h(t,x_i+1,\beta_i)}{h(t,x_i,\beta_i)} = \frac{h_0(t) e^{\beta_i(x_i+1)}}{h_0(t) e^{\beta_i x_i}} =e^{\beta_i} \]

对于分类型协变量，类别\(i\)的回归系数是\(\beta_i\)，类别\(j\)的回归系数是\(\beta_j\)，则

\[HR=\frac{h(t,\beta_i)}{h(t,\beta_j)} = \frac{h_0(t) e^{\beta_i}}{h_0(t) e^{\beta_j}} =e^{\beta_i-\beta_j} \]

如果与baseline相比，则

\[HR=\frac{h(t,\beta_i)}{h(t)} = \frac{h_0(t) e^{\beta_i}}{h_0(t)} =e^{\beta_i} \]

当\(\beta_i>0\)时，风险增加，\(\beta_i<\)时，风险降低；因此\(HR>1\)称为危险因素，\(HR<1\)称为保护因素。

R语言实现

# 指数分布
fit_exp <- survreg(Surv(time, status) ~ sex, data = lung, dist = "exponential")

1/exp(fit_exp$coefficients) # sex2: 0.6062

exp1 = predict(fit_exp, newdata = list(sex = 1), type = "quantile", p = seq(0, 1, 0.01))
exp2 = predict(fit_exp, newdata = list(sex = 2), type = "quantile", p = seq(0, 1, 0.01))
pred_exp <- data.frame(S = seq(1, 0, -0.01), sex1 = exp1, sex2 = exp2)
pred_exp_long <- gather(pred_exp, key = "sex", value = "time", -S) 

# Cox模型
fit_cox <- coxph(Surv(time, status) ~ sex, data = lung)

exp(fit_cox$coefficients) # sex2: 0.5880

cox1 <- survfit(fit_cox, newdata = data.frame(sex = 1))
pred_cox1 <- data.frame(S = cox1$surv, time = cox1$time, sex = "sex1")
cox2 <- survfit(fit_cox, newdata = data.frame(sex = 2))
pred_cox2 <- data.frame(S = cox2$surv, time = cox2$time, sex = "sex2")

# 画图
km <- survfit(Surv(time, status) ~ sex, data = lung)
kmplot <- ggsurvplot(km, risk.table = T)

kmplot$plot <- kmplot$plot + 
  geom_line(data = pred_exp_long, aes(x = time, y = S, group = sex), linetype = 2) +
  geom_line(data = pred_cox1, aes(x = time, y = S), col = "blue") +
  geom_line(data = pred_cox2, aes(x = time, y = S), col = "blue") 

kmplot