生存分析期末复习（理论部分）

Survival Analysis

一、基本概念：

$T^*:$ survival time 生存时间
复发事件：住院 $\rightarrow$ 好了 $\rightarrow$ 住院 $\rightarrow$ 好了 $\rightarrow$ 住院 $\rightarrow$ 挂了
竞争风险：有多种原因中的一个导致了时间的结束（死亡）如并发症
Multistate Model：多状态模型，

如 illness-death model：自己画一下！

censoring 删失数据：知道生存时间（survival time）的部分信息，但不知道确切的生存时间
Right-censored ： Survival time > observed time
left-censored ： Survival time < observed time
interval-censored 生存时间在某个区间内

左截断：生存时间需要足够长才能造成左截断！此时的样本并非随机且无偏的，因为：例如对艾滋病的研究，可能会遗漏哪些感染了艾滋病并很快死亡的样本。

这被称为： biased sample 幸存者偏差

两个重要的函数

对连续的 $T^*$ : 生存函数（单调不增）

$S(t) = P(T^*>t) = 1-P(T^*\leq t)$

风险函数：

$\lambda(t) = \lim_{h\rightarrow0+}\frac{P(t\leq T^* \leq t+h|T^*\geq t)}{h} \qquad [t,t+h)$

即： $\lambda(t) =\lim_{h\rightarrow0+}\frac{P(t\leq T^* \leq t+h)}{h}/S(t) = \frac{S(t)-S(t+h)}{h} /S(t) = -\frac{dS(t)}{dt}/S(t) = -\frac{dlnS(t)}{dt}$

反过来， $S(t) = exp\{-\int_{0}^{t}\lambda(u)du\}$

对删失的假设：用 $C^*$ 表示删失

1、 Random censoring : 生存时间和删失独立

2、独立删失假设：
$\lim_{h\rightarrow0+}\frac{P(t\leq T^* \leq t+h|T^*\geq t,C^*\geq t,X)}{h} = \lim_{h\rightarrow0+}\frac{P(t\leq T^* \leq t+h|T^*\geq t,x)}{h}$
对得到的数据， $T_i = min(T_i^*,C_i^*)$ 是个体i的观测值

令： $\delta_i = I(T_i^*\leq C_i^*) = \left\{\begin{matrix}1,见后面~ \\0,\qquad censored \end{matrix}\right.$ event occurred, $T_i = T_i^*$

右删失的数据长这样： $\{(T_i,\delta_i,X_i(t));0\leq t \leq T_i,i=1,2,\cdots,n\}$

二、对生存函数和风险函数的非参数估计：

预备数学知识

注：累积风险函数： $\Lambda(t) = \int_{0}^{t}\lambda(u)du$

$S(t) = exp\{-\int_0^t\lambda(u)du\} = exp\{-\Lambda(t)\}$ 生存函数可以由累积风险函数表示
$E(T^*) = \int_{0}^{\infty} S(t)dt$ why?

$T^*$ 非连续：

$f(a_j) = P(T^* = a_j)$
$S(t) = P(T^*>t) = \sum_{j|a_j>t}f(a_j)$ 生存函数
风险函数：一个个体活到 $a_j$ 那么 $(0,a_j)$ 未发生事件，而在 $a_j$ 这个时间点发生事件的概率：

$\lambda_j = P(T^* = a_j|T^*\geq a_j) = \frac{P(T^* = a_j)}{P(T^*\geq a_j)} = \frac{f(a_j)}{S(a_{j}-)}=\frac{f(a_j)}{S(a_{j-1})}$

注： $P(T^*\geq a_j) \neq S(a_j) ;\quad P(T^*\geq a_j) = S(a_{j-1})$

这里的 $S(a_j-)$
$S(a_j-) = \lim_{t\rightarrow a_j-}S(t) = P(T^*\geq a_j) = P(T^*> a_{j-1}) = S(a_{j-1})$
此外，
$S(t) = \Pi_{j|a_j\leq t}(1-\lambda_j) = \Pi_{j|a_j\leq t}P(T^* \geq a_{j+1}|T^*\geq a_j)$
注： $1-\lambda_j = 1-P(T^*=a_j|T\geq a_j) = P(T^*>a_j|T^*\geq a_j) = P(T^*\geq a_{j+1}|T^* \geq a_j)$

注： $f(a_j) = P(T^*=a_j) = \lambda_j \Pi_{k=1}^{j-1}(1-\lambda_k)$

离散生存函数的证明.png

Kaplan-Meier Estimator of S(t)

$d_j$ : 在时间点 $t_j$ 发生的时间

$m_j:$ 在时间段 $[t_j,t_{j+1})$ 删失的数据

$n_j: = (d_j+m_j)+\cdots+(d_k+m_k)$ 为 $t_j$ 时间段之前处于风险中的事件个数

风险集： $\{i:T_i^*\geq t_j,C_i^*\geq t_j\}$

当处于时间原点的时候， $n_0=n,d_0=0$ 用似然函数的方法推导K-M estimator：

Nelson-Aalen Estimator of $\Lambda(t)$

首先，风险的估计： $\hat{\lambda_j} = \frac{d_j}{n_j}$ 在时间点 $t_j$ 而： $d\Lambda(t) = \left\{\begin{matrix}\lambda_j,\qquad t=t_j,j=1,\cdots,k \\0,\qquad \qquad\quad\quad otherwise \end{matrix}\right.$

$\therefore d\hat{\Lambda}(t) = \left\{\begin{matrix}\hat{\lambda_j},\qquad t=t_j,j=1,\cdots,k \\0,\qquad \qquad\quad\quad otherwise \end{matrix}\right.$

$\therefore \hat{\Lambda}(t) = \int_0^td\hat{\Lambda}(u) = \sum_{j|t_j\leq t}\frac{d_j}{n_j}$

这就是 Nelson-Aalen 估计

variance Estimator of $\hat{S}(t)$ and $\hat{\Lambda}(t)$

给定风险集， $\{i:t_i\geq t_j,i=1,\cdots,n\}$

$\lambda_j = P(T^*=t_j|T^*\geq t_j)$ 它的点估计是： $\hat{\lambda_j} = \frac{d_j}{n_j}$

and $d_j\sim Binomial(n_j,\lambda_j)$ 所以： $E(d_j) = n_j\lambda_j$ $var(d_j) = n_j\lambda_j(1-\lambda_j)$

从而： $E(\hat{\lambda_j}) = \lambda_j$ $var(\hat{\lambda_j}) = \frac{1}{n_j}\lambda_j(1-\lambda_j)$

$var(\hat{\lambda_j}) = \frac{1}{n_j}\hat{\lambda_j}(1-\hat{\lambda_j}) = \frac{d_j(n_j-d_j)}{n_j^3}$

当 $n_j$ 很大， $n_jp_j$ 很大时，渐进正态，即： $\hat{\lambda_j}-\lambda_j\sim N(0,\frac{d_j(n_j-d_j)}{n_j^3})$
$Log\hat{S}(t) = \sum_{j|t_j\leq t} log(1-\hat{\lambda_t}) \qquad 连乘变连加$
从而 $var\{log\hat{S}(t)\}\approx \sum_{j|t_j\leq t} var\{log(1-\hat{\lambda_j})\} \approx \sum_{j|t_j\leq t}\frac{1}{(1-\hat{\lambda_j})^2}var(1-\hat{\lambda_j})$

$= \sum_{j|t_j\leq t}\frac{1}{(1-\hat{\lambda_j})^2}var\{\hat{\lambda_j}\} \approx \sum_{j|t_j\leq t}\frac{d_j}{n_j(n_j-d_j)}$

上式中有用到 $delta's$ $method:$ $var[f(Y)] \approx [f^{'}(Y)]^2var(Y)$

$\therefore \hat{var}\{\hat{S}(t)\} = \hat{S}^2(t)\hat{var}\{log\hat{S}(t)\} = \{\Pi_{j|t_j\leq t}(1-\frac{d_j}{n_j})^2\}\times \sum_{j|t_j\leq t}\frac{d_j}{n_j(n_j-d_j)}$

$and\quad 95\% \quad CI : \hat{S}(t)\pm 1.96\sqrt{\hat{var}\{\hat{S}(t)\}}$

example k-M估计

例：数据： 24+，16+，8,19,10,8+，5,17,20,10

重排： 5,8,10,16,17,19,20,24

Time	At risk	Died	censored	Survival Rate	Estimator
$t_i$	$n_i$	$d_i$	是否删失	1- $\frac{d_i}{n_i}$	$\hat{S}(t),t_i<t<t_{i+1}$
0	10	0	0	1-0=1	1
5	10	1	0	$1-\frac{1}{10}=0.9$	$1\times 0.9=0.9$
8	9	1	1	$1-\frac{1}{9} = 0.89$	$0.9\times 0.89=0.80$
10	7	2	0	$1-\frac{2}{7} = 0.71$	$0.80\times 0.71=0.57$
16	5	0	1	1-0=1	$0.57\times1=0.57$
17	4	1	0	···	···
···				···	···

客观比较两条或多条生存曲线的差异：是否统计显著？ log-rank test 对数秩检验

Note： Median Survival time ： kk2012 P30 这个也挺重要！

R code: kk2012 P625-629

三、 Nonparametric Test of Survival Differences:

Log-Rank Test (非参数方法)

用以检验分布函数是否相同，而我们知道： $S(t) = 1-F(t)$ 所以等价于生存函数是否相同

$H_0: S_1(t) = S_2(t)=\cdots = S_p(t),\forall t\geq 0$ $H_{\alpha}:\exist j,k,1\leq j\neq k\leq p,t_0\geq 0,s.t.\quad S_j(t_0)\neq S_k(t_0)$

上课虽然确实讲了手算公式，但是应该不会真的让手算吧hhh：

推广： stratified Log-Rank Test 不同组别的S(t)不同

新的检验统计量：

如何分层:

weighted Log-Rank Test 一般为两个生存曲线有交叉时使用

例如：某新药在服用初期较显著，但过了一段时间药效一般，这时候用weighted Log-Rank Test，p=1 能检测出该药比安慰剂有效，如果采取标准 Log-Rank Test则可能得出该药与安慰剂差不多的结论~

四、混淆(Confounding)和交互项(Interaction) 效应

混杂: 例如辛普森悖论中的性别

In summary， we should stratify on "gender" in stratified analysis; we should adjust for "gender" in regression analysis.

特别的，当C是一个类别变量时，例如：（two groups : A=1,A=0) ,分层/回归的方法用来消除混杂效应是必要的。

例子3：白血病 logWBC kk2012P30

一个不需要考虑混杂的 Randomization 随机化 A和X独立

其实，混杂问题应该在试验设计时被处理，而非在数据分析阶段（被处理）例如： a double-blinded randomized clinical trial

NOte：并非任何时候都能做到 Randomization 例如 e.g. kk2012P30 不能随机分配到吸烟/非吸烟组吸不吸烟不由我们决定。

Interaction Effect

交叉项(e.g. $x_i\times x_j$ )的存在是否有必要？

1、 KM-plot in all strata

2、stratified log-rank test ; log-rank test in all strata

3、 cox regression

4、stratified cox regression

注： stratification uses less assumption than regression, but requires more data to perform well.

五、Survival Model 的简单介绍：

指数模型 Exponential Model

$T^*\sim exponential \quad distribution ,$ i.e. $f(t) = \lambda e^{-\lambda t}$

survival function & hazard function :
$S(t) = e^{-\lambda t} \qquad, t\geq 0$

$\lambda(t|x) = exp(\alpha+\beta^TX)$

Weibull Model

$T^*\sim weibull \quad distribution$
$S(t) = e^{-\lambda t^p},\qquad t>0$

$\lambda(t) = \lambda pt^{p-1} ,\quad \lambda>0,\quad p>0$

其中 $\lambda = exp(\alpha+\beta^TX)$

weibull model:
$\lambda(t|X) = exp(\alpha+\beta^T X)pt^{p-1}$
可以看出，p=1时是指数模型

Log-logistic Model

$T^*\sim Log-logistic \quad distribution$ then:
$S(t) = \frac{1}{1+\lambda t^p}$

$\lambda(t) = \frac{\lambda pt^{p-1}}{1+\lambda t^p}$

以后讨论 $\lambda$

Log-normal Model (印象里没讲？)

$T^* \sim log-normal$ i.e. $logT^*\sim N(\mu,\sigma^2)$ $f(t) = \frac{\Phi(\frac{logt-\mu}{\sigma})}{\sigma-t}$ ? ; $\Phi(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$

$\lambda(t) = \frac{f(t)}{s(t)}$ we let $\mu = \alpha+\beta^TX$ :
$f(t|X) = \frac{\Phi(\frac{logt-\alpha-\beta^TX}{\sigma})}{t}$
?

令 $Y = logT^*$ 那么： $Y = \alpha+\beta^TX+\sigma W$ where $W\sim N(0,1)$

这也是一种加速失效模型，详见下：

Accelerated Failure Time Model 加速失效模型

model： $Y = \alpha+\beta^TX+\sigma W$ 其中 W 是随机误差项； $Y = logT^*$

对数生存时间 $logT^*$ ，(if $x\rightarrow x+1,Y\rightarrow Y+\beta$ )

weibull model is a special case of AFT models: √ 指数模型，log-normal 也是特殊的AFT

而log-logitic 属于 proportional odds 等比例比值模型。

为什么叫加速失效？加速因子是什么？

简单分析，令 $\alpha=0,\sigma=1$ 且 $X是协变量$ X=1代表接受治疗，treatment，X=0 则为安慰剂placebo

$S(t|trt) = P(T^*>t|trt) = P(logT^*-\beta>logt-\beta|x=1) = P(w>logt-\beta)$

$S(t|placebo) = P(w>logt)$

$\therefore S(te^\beta|trt) = P(w>log(te^\beta)-\beta) = S(t|placebo)$

例如： $\beta=log0.5,那么：S_{trt}(0.5t) = S_{placebo}(t)$ 从而trt相比placebo会加速失效，且中位数生存时间trt组是placebo组的一半！其实是一个道理： $\beta<0,$ ↑trt组的病人会更快的发生事件。

我们把 $e^\beta$ 叫做加速因子√

Cox Model

模型如下：
$\lambda(t|X) = \lambda_0(t)e^{\beta^TX}$
其中 $\lambda_0(t)$ 称为基准风险函数，即： $\lambda_0(t) = \lambda(t|X=0)$

Hazard ratio(HR) 风险比： $\frac{\lambda(t|x_1,\cdots,x_{j-1},x_{j},x_{j+1},\cdots,x_p)}{\lambda(t|x_1,\cdots,x_{j-1},x_{j+1},\cdots,x_p)} = exp(\beta_j)$

so： one unit increase of $X_j$ will change the Hazard function to $exp(\beta_j)\times hazard\quad function$

i.e. $if\quad x_j\rightarrow x_j+1,\lambda(t|x)\rightarrow \lambda(t|x)\times exp(\beta_j)$

同样的，对两个数据 $X,X^*$ $HR = \frac{\lambda(t|X^*)}{\lambda(t|X)} = exp(\beta^T(X^*-X))$ 它与时间t无关，这也是为什么cox模型被称为等比例风险模型（cox proportional hazards model)

不同于参数化的模型，cox模型不限制风险函数的形状。cox模型没有截距项（相比于其他模型 $\alpha+\beta^TX$ )

Connection of Cox Model and AFT Model

connection of AFT and cox.png

六、生存分析模型的似然估计：*

生存分析： $T^*$ 密度函数 $t(t|x;\theta)$ 风险函数： $\lambda(t|x;\theta)$ 生存函数： $S(t|X;\theta)$

in weibull model: $\theta = (\alpha,\beta,p)$

in AFT model: $\theta = (\alpha,\beta,\sigma)$

in Cox model: $\theta=(\lambda_0(t),t\geq 0;\beta)$

我们要做的就是估计 $\hat{\theta}$ go！！

延用之前的符号： $T_i = min(T_i^*,C_i^*);\delta_i = I(T_i^*\leq C_i^*)$

$P(T_i \in [t,t+dt),\delta_i=1|x_i;\theta) = P(T_i^* \in[t,t+dt),C_i^*>t|x_i;\theta )$ random censoring 假设

$=P(T_i^*\in[t,t+dt)|x_i;\theta) \times P(C_i^*>t|x_i;\theta)$

$=f(t|x_i;\theta)dtP(C_i^*>t|x_i)$

从而： $f(T_i,\delta_i=1|X_i;\theta) = f(T_i|X_i;\theta)P(C_i^*>T_i|X_i)$ (1)

Note: 后面那个P是删失时间 $C_i^*$ 的生存函数，其中 $C_i^*$ 是随机变量，而此时 $T_i$ 是观测值

同理： $P(T_i \in [t,t+dt),\delta_i=0|x_i;\theta) = P(C_i^* \in[t,t+dt),T_i^*>t|x_i;\theta )$

$=P(T_i^*>t|x_i;\theta)\times P(C_i^*\in[t,t+dt)|x_i;\theta)$

$=S(T_i|X_i;\theta) f(C_i^*=T_i|X_i)$ (2)

接下来是似然函数：*

对右删失数据（我们目前主要就关心右删失） $\{(T_i,\delta_i,X_i),i=1,\cdots,n\}$ 其似然函数：

$L(\theta) = \Pi_{i=1}^{n}f(T_i,\delta_i,X_i) = \Pi_{i=1}^{n}f(T_i,\delta_i|X_i)f(X_i)$

$这里仔细想想！精髓=\Pi_{i=1}^n\{f(T_i,\delta_i=1|X_i)\}^{\delta_i}\times \{f(T_i,\delta_i=0|X_i)\}^{1-\delta_i}\times f(X_i)$

由(1)(2) 两式子：

$=\Pi_{i=1}^n\{f(T_i|X_i,\theta)\}^{\delta_i}\times\{S(T_i|X_i;\theta)\}^{1-\delta_i}\times \{P(C_i^*>T_i|X_i)\}^{\delta_i}\times\{f(C_i^*=T_i|X_i)\}^{1-\delta_i}\times f(X_i)$

后三者是不含有 $\theta$ 的，因此，当我们取对数似然函数并求导研究其极值时，会消失，可以认为：

$\propto \Pi_{i=1}^n\{f(T_i|X_i,\theta)\}^{\delta_i}\times\{S(T_i|X_i;\theta)\}^{1-\delta_i}$ (3)

更神奇的是，还记得： $\lambda(t) = lim_{t\rightarrow 0} \frac{P(t\leq T^*<t+h)}{h}/S(t) = \frac{f(t)}{S(t)}$

而（3）式 $= \Pi_{i=1}^n[\frac{f(T_i|X_i,\theta)}{S(T_i|X_i;\theta)}]^{\delta_i}\times S(T_i|X_i;\theta) = \{\lambda(T_i|X_i;\theta)\}^{\delta_i}\times S(T_i|X_i;\theta)$ （4）

$= [\lambda(T_i|X_i;\theta)]^{\delta_i}\times exp(-\Lambda(T_i|X_i;\theta)) =[\lambda(T_i|X_i;\theta)]^{\delta_i}\times exp(-\int_0^{\infty}\sum_{j\in R(u)}\lambda(u|X_j;\theta)du)$

其中 $R(u) = \{i;T_i\geq u,i=1,2,\cdots n\}$ 是时间u的风险集；似然函数可以由风险函数写出↑

极大似然估计及其性质：

对（4）取对数并求 $\frac{\part logL(\theta)}{\part \theta}=0$ ; 因为 $L(\theta)$ 对 $\theta$ 是凸函数，肯定存在唯一确定的全局最大值 $\hat{\theta}$ ,由牛顿迭代法可以求得数值解： $\theta^{(t+1)} = \theta^{(t)}-[\frac{\part^2 logL(\theta)}{\part \theta \part \theta^T}]^{-1}|_{\theta=\theta^{(t)}}\times \frac{\part logL(\theta)}{\part \theta}|_{\theta=\theta^{(t)}}$

迭代停止条件： $||\theta^{(t+1)}-\theta^{(t)}||_2< \varepsilon$ so,少年，选择合适的初值，开始迭代吧（不是）

理论上，e.g. Weibull model： $\lambda(t|X;\theta) = e^{\alpha+\beta^TX}pt^{p-1}$ $S(t|X,\theta) = e^{(-e^\alpha+\beta^TX)t^p}$

由（4） $L(\theta) \propto \{\lambda(T_i|X_i;\theta)\}^{\delta_i}\times S(T_i|X_i;\theta) = \Pi_{i=1}^n(e^{\alpha+\beta^TX}pT_i^{p-1})^{\delta_i}\times e^{(-e^\alpha+\beta^TX)T_i^p}$

$\therefore logL(\theta) \propto \sum_{i=1}^n[\delta_i(\alpha+\beta^Tx_i+logP+(p-1)logT_i)-e^{\alpha+\beta^Tx_i}T_i^p]$

然后你就分别对 $\alpha,\beta,p$ 求导让他们等于0呗~ 数值解还是需要牛顿迭代法

MLE的性质： $\hat{\theta}$

$\sqrt{n}(\hat{\theta}-\theta) \overset{d}{\rightarrow} N(0,I^{-1}(\theta))$

$I(\theta) = E[\frac{\part logL_i(\theta)}{\part \theta}\times\frac{\part logL_i(\theta)}{\part \theta^T}] = -E[\frac{\part^2 logL_i(\theta)}{\part \theta \part \theta^T}] = -\frac{1}{n}E[\frac{\part^2 logL(\theta)}{\part \theta \part \theta^T}]$ 是 Fisher信息矩阵

！ under $H_0: \theta=\theta_0$ $n \rightarrow \infty$ $\sqrt{n}(\hat{\theta}-\theta) \overset{d}{\rightarrow} N(0,I^{-1}(\theta))$ wald test
$n(\hat{\theta}-\theta)^TI(\theta_0)(\hat{\theta}-\theta) \overset{d}{\rightarrow} \chi^2_r$ score test r是 $\theta$ 的维数

$-2log\frac{L(\theta_0)}{L(\hat{\theta})}\overset{d}{\rightarrow} \chi^2_r$ likelihood ratio test

Cox 的偏似然函数方法*

只关心 $\beta$ 而不管 $\lambda_0(t)$

$L_p(\beta) = \Pi_{j=1}^k\frac{\lambda(t_j|X_{(j);\theta})}{\sum_{l\in R(t_j)}\lambda(t_j|X_l;\theta)} = \Pi_{j=1}^k \frac{e^{\beta^TX_{(j)}}}{\sum_{l\in R(t_j)}e^{\beta^TX_l}}$ 后面那个等式是因为cox模型

其中： $R(t) = \{i; T_i\geq t,i=1,\cdots ,n\}$ 是时间点t的风险集。