Cox模型方法简要总结

参考书目 Applied Survival Analysis: Regression Modeling of Time-to-Event Data, 2nd. (2015)
以常见问题和理论为主,有空再整代码

1. Hazard function and hazard ratio

Hazard function $h(t,x,\beta)=h_0(t)e^{\beta x}$
Hazard Ratio = $HR =\cfrac{h_0(t)e^{X_1\beta}}{h_0(t)e^{X_2\beta}} = e^{(X_1 -X_2)\beta}$

Interpertation of HR: $\beta$ = log2 =0.69 > 0 -> people in group A are dying at twice the rate of people in group B.

Cox模型的前提:

Hazard Ratio 恒定, 不随时间变化
Cox model 只能获得一个HR, 当HR随时间变化时当然不能用. 网上随便找个图举例. 肉眼可见,前期HR接近1不显著, 后期survival experience 有较大差异. 建议使用Cox extension model加入时间变量. 因为前期的曲线交叉,非得使用Cox model 也不能用log-rank test 计算P值 (用Renyi test)

image.png

2. $X\beta$ 线性,并且X和时间无关.
反例1, 非线性. 30岁 vs 20岁 HR = 1.5, 40岁 vs 30 岁 HR = 3,风险非线性增长,但是把年龄放在同一个变量X里,导致只获得一个HR.
反例2: X 和时间相关. X = age,实际上年龄随着时间变化, 当随访时间较长时不适合该模型.

3. 每个人的生存时间都不相同
理论上是这样,但是由于记录的精度不同很难做到. 解决方法见 tied survival time.

Cox 模型无法得到baseline hazard function $h_0(t)$ , 但是能获得 $e^{\beta x}$ 的参数估计:
当不考虑censoring时: $l_p(\beta)=\prod \cfrac{e^{x_i \beta }}{\sum e^{x_j\beta}}$
当两组差异悬殊时likelihood无法收敛,结果HR很大或趋近0. (此同logistic regression. 在生存分析里例子: A组几乎全死的比B组早). 此时可使用Firth方法(基于Bayesian approximation)进行参数估计.
如果没有censored数据,直接用logistic model更合适. 生存数据一般是右倾分布,记得使用log transform. 最后得到的是生存时间的Geometric mean ratio.

1.1 Cumulative hazard function

Recall that
$H(t)=\int_0^th(u)du=\int_0^t\cfrac{p(u)}{S(u)}du=-log(S(t))$
$S(t)=e^{-H(t)}$ and $p(t)=\frac{d}{dt}P(t)=h(t)S(t)=h(t)e^{-H(t)}$

Thus,
$H(t,x,β)=\int^t_0h(u,X,\beta)du=\int^t_0h_0(u)exp(X'\beta)du\\=exp(X'\beta)\int^t_0h_0(u)du=exp(X'\beta)H_0(t)$

and $S(t,X,β) = e^{-exp(X'\beta)H_0(t)}$

$S(t,X,β) = [S_0 (t)]^{exp(\beta_1x_1+\beta_2x_x+...+\beta_px_p)}$
( $S_0(t)$ is the baseline survaival function)

The cumulative hazard is measured by the cumulative baseline risk $H_0(t)$ and modified by effect X

1.2 Partial Likelihood

Distribution of survival time --> cumulative distribution function F(t,β,x)
$S(t,β,x) = 1 − F(t,β,x)$

Example:
$F(1,β,80)$ : proportion of 80-year-old patients expected to die in less than or equal to 1 year
$f(t,β,x)$ : the probability that an x-year-old patient survives exactly t years

The actual likelihood function is constructed by considering the
contribution of the triplets (t,1,x) and (t,0,x) separately.
0 --> survive at least t (including censored)
1 --> die exactly in t

In general, the contribution of each triplet to the likelihood is
$[f(t,β,X)]^c [S(t,β,X)]^{1−c}$ , where c = 0 or 1 (所有censored的数据将不会对likelihood 产生影响)

$l(\beta)=\prod^n_{i=1}[f(t_i,\beta,X_i)]^{c_i}[S(t_i,\beta,X_i)]^{1-c_i}$

maximize the likelihood function
$L(\beta)=\sum^n_{i=1}c_ilog\ f(t_i,\beta,X_i)+(1-c_i)S(t_i,\beta,X_i)$
-> $L(\beta)=\sum^n_{i=1}c_ilog\ h_0(t_i)+\beta c_ix_i+e^{\beta x_i}logS_0(t_i)$

1.3 Estimates and 95%CI of Cox model

Estimates of $\beta$
$l_p(\beta)=\prod^n_{i=1}\bigg[\cfrac{e^{\beta x_i}}{\sum_{j\in R(t_i)}e^{\beta x_j}}\bigg]^{c_i}$

$L_p(\beta) = \sum^m_{i=1}\bigg[\beta x_i-log\sum_{j\in R(t_i)}e^{\beta x_j}\bigg]$
Varivance:

$I(\beta)=-\cfrac{d^2L_p(\beta)}{d\beta^2}$ , $\hat{Var}[\hat\beta]=I(\hat \beta)^{-1}$

显著性检验的三种方法

Partial likelihood ratio test $G = 2[L_p(\hat β) − L_p (0)]$
$L_p(0)=-\sum log\ n_i$ , $G \sim \chi^2_{df=1}$
P-value = $P[\chi^2_{df=1}\ge G]$
Wald test: $z=\cfrac{\hat \beta}{\hat{SE}(\hat \beta)}$ , $z\sim N(0,1)$
P-value = $2P[Z>|z|]$
Score test: $z^*=\cfrac{\frac{dL_p}{d\beta}}{\sqrt{I(\beta)}}\bigg|_{\beta=0}$ , $z^*\sim N(0,1)$
P-value = $2P[Z>|z^*|]$
In SAS, z follows chi-sq distribution with df = 1.

95% Wald confidence interval 和Wald test假设相同
95% CI of $\beta$ = $\hat \beta \pm 1.96\hat{SE}$

多变量检验
1. Maximum likelihood $G\sim{\chi^2_{df=p}}$ , p = number of predictors
2. Score test
First partial derivatives = $u(\beta)$
Statistic = $u\prime(0)[I(0)]^{-1}u(0)$ ;
Information matrix = covariance matrix (0)
Under the null hypothesis (变量全都不显著), the vector of scores u(0) will be distributed as multivariate normal with a mean vector equal to 0
3.Multivariable Wald test
Wald statistic = $\hat \beta^\prime I(\hat \beta)\hat\beta$
Asymptotically as chi-square with df = p

1.4 Tied survival times

生存时间重合解决方法

Exact method (Kalbfleisch and Prentice, 2002)
Breslow’s (1974) approximation
Efron’s (1977) approximation

Exact method

此方法认为同时发生的事件实际上存在先后顺序,于是计算所有可能性. 假设同一时间记录了5个事件, 排序有5!=120种, 每种可能记为Oi, 取并集...
缺点: 同一时间有很多事件时计算量太大

Breslow’s and Efron’s approximations

Breslow方法在ties比较少时好用, Efron在ties比较多时好用(算得快)

1.5 log-hazard function

$g(t,x,\beta)=log[h_0(t)]+\beta X$
$g(t,x=a,\beta)=g(t,x=b,\beta)=(a-b)\beta$
计算g(x)直接用加减法很方便,然后计算HR = exp(g(x1)-g(x2))
$HR(t,x=a\ vs\ x=b \beta) = \frac{h(t,x=a,\beta)}{h(t,x=b,\beta)} = e^{(a-b)\beta}$
同样,censored的数据不会改变HR

2. Model selection

变量选择首先要注意不收敛的情况,有时候共线性也会导致不收敛. 和logistic regression一样使用VIF检查共线性.

模型选择的三种方法

Stepwise and best subset selection
Purposeful selection
Fractional polynomials

2.1 Stepwise

Stepwise以p值为标准比较两个模型
$G = -2[L(加了新变量的模型) - L(没加)]$ ,然后计算p值. 选择加入后p值最小的变量, 如果 $p < P_E$ (可选15%, 也可选更高纳入更多变量,然后进一步分析) 就加入这个变量.
在加入了新变量后还要检查模型里的其它变量的显著性是否发生改变. 若有变量不再显著( $p > P_R, 这里P_R >P_E$ )则要删去这个变量.
不再显著的原因: 相关, 共线性等 (待补充说明)
Stepwise需要仔细检查所选出的变量. 如在使用dummy variable时,所有分类必须共同纳入(待补充说明)

2.2 Best-subsets selection

假设有5个备择变量,所有组合有2^5种,从中选出最好的. 比较标准包括R-Square, adjusted R-Square, Mallow’s C. 一般分别从单变量,二变量,三变量...每种中选出几个最优模型,然后进行比较. 变量过多可能过拟合(引入过多噪声), 太少可能欠拟合. 可以以每10个events加入一个变量为标准避免过拟合.
Mallow’s C
$C = W_q + (p − 2q)$
$W_q$ : Wald statistic, 可用 score test的值代替
p: 变量数目

Purposeful selection

根据研究目的选择变量.
研究目的可分为两类: 使用模型进行预测; 使用模型评估某因素的风险. 后者典型为临床试验,评估治疗的效果.
预测模型追求模型预测效果,共线性等问题可以忽略.所以难以评估每个因子具体影响大,或者说effect不准确.
风险评估模型以某因素或者某些因素为研究目的
(待续)

Cox模型方法简要总结