构建具有时间依存协变量的Elastic-net Cox模型

介绍

生存分析的目标是通过事件发生率来模拟过程。在应用统计设置中，这通常意味着识别与某些与基线相比可能导致更高或更低事件率的协变量。当事件是缓解或死亡时，较高的事件发生率对应于风险较高的患者。在生存建模中一个特别的挑战是存在删失（censoring），即患者在随访时间被记录但事件尚未（或可能永远不会）发生[1]。例如，可以对高风险癌症患者监察一整年，在结束时其中一些患者已经过世而一些患者仍然活着——后者在一年时间内有一个删失观察。

医学研究中最常见的模型是比例风险模型（proportional hazards model），该模型试图捕捉患者发生事件的相对顺序。在这样的框架中，删失患者可以提供信息来模拟事件发生率，因为一个删失的患者比一个经历某事件的患者的寿命更长，他可能是风险较低的患者[2]。在许多具有多个随访点的研究中，患者的协变量信息会随时间而变化。例如，在survival package中的heart数据集中，在心脏移植的等待列表中测量患者存活率。以下是数据片段。

library(survival)
head(heart)

## start stop event        age      year    surgery transplant id
## 1     0   50     1 -17.155373 0.1232033       0          0  1
## 2     0    6     1   3.835729 0.2546201       0          0  2
## 3     0    1     0   6.297057 0.2655715       0          0  3
## 4     1   16     1   6.297057 0.2655715       0          1  3
## 5     0   36     0  -7.737166 0.4900753       0          0  4
## 6    36   39     1  -7.737166 0.4900753       0          1  4

我们可以看到id号是3的患者（第3和4行数据）在等了一个月然后接受移植手术，移植后又过了15个月后死亡；id号是4的患者（第5和6行数据）在等了36个月然后接受移植手术，移植后又过了3个月后死亡。如上所示，时间相关数据集的格式必须用长格式（long-format），其中每个新的测量周期都被视为新观察。这个填充附加行可以被认为是添加左删失观察（除了上面讨论的标准右删失），因为这个新的“患者”在其开始时间的左边没有提供任何信息。

虽然在生存或时间 - 事件分析的许多数据集中，都具有这些与时间相关的属性，但目前使用glmnet或fastcox并不支持来估计这些模型的正则化（regularized version），这些基准库将用于拟合elastic-net模型到比例风险损失函数（proportional hazards loss function）。在本篇中，我们将展示如何构建一个自定义的近端梯度下降算法（proximal gradient descent algorithm），该算法可以包含时间依存协变量（time-dependent covariates）。

Elastic-net cox model

为了建立符号公式，将每个观察定义为具有以下元组（tuple）：(til,tiu,δi,xi)(t_i^l, t_i^u, \delta_i, x_i)(til,tiu,δi,xi)，在这里tilt_i^ltil和tiut_i^utiu是第i号病人的协变量的时间区间，δi\delta_iδi是病人在tiut_i^utiu时刻是否发生时间的指示器，xix_ixi是病人协变量测量值矢量。cox模型的部分似然（ partial likelihood）可以很容易地适用于处理这种时间依存协变量的情况。
Partial Likelihood:
L(β)=∏i=1N(exiTβ∑j∈R(ti)NexjTβ)δiL(\beta)=\prod\limits_{i=1}^N(\dfrac{e^{x_i^T \beta}}{\sum\limits_{j\in{R(t_i)}}^Ne^{x_j^T\beta}})^{\delta_i}L(β)=i=1∏N(j∈R(ti)∑NexjTβexiTβ)δi
Partial Log-Likelihood:
ℓ(β)=∑i=1Nδi(xiTβ−log[∑j∈R(ti)Nexp(xiTβ)])\ell(\beta)=\sum\limits_{i=1}^N\delta_i(x_i^T\beta-log[\sum\limits_{j\in R(t_i)}^N exp(x_i^T\beta)])ℓ(β)=i=1∑Nδi(xiTβ−log[j∈R(ti)∑Nexp(xiTβ)])

其中，xi=xi(t)x_i=x_i(t)xi=xi(t)是时间函数；R(ti)R(t_i)R(ti)是在事件时间tit_iti存活/未发生删失的病人的风险指标，特别的：
R(ti)={j:(tju≥ti)∧(tjl<ti)}R(t_i)=\left\{j:(t_j^u\ge t_i)\land (t_j^l<t_i) \right\}R(ti)={j:(tju≥ti)∧(tjl<ti)}
第一个条件确保患者要么经历事件，要么在稍后晚于tit_iti的时间点被审查（因此当i号患者经历事件时是活着的）；并且第二个条件确保开始时间发生在事件之前。请注意，在heart dataset中，[t3l,t3u]=(0,1][t_3^l,t_3^u]=(0,1][t3l,t3u]=(0,1]并且[t4l,t4u]=(1,16][t_4^l,t_4^u]=(1,16][t4l,t4u]=(1,16]保证了3号病人永远不会经历他自身的风险2遍。在yi(tj)=I[i∈R(tj)]时，y_i(t_j)=I[i\in R(t_j)]时，yi(tj)=I[i∈R(tj)]时，使用独热（注释1）^{（注释1）}（注释1）编码[Y]ij=yi(tj)[Y]_{ij}=y_i(t_j)[Y]ij=yi(tj)将极大降低计算量。

对于高维数据集和预测问题，其研究目标是找到一些β\betaβ能够平衡模型的拟合（例如使用数据集的全部信息）并且保证其高的泛化能力（例如忽略数据集特异性噪声）。正则化是一种返回系数向量的技术，该系数向量比原本返回的系数向量“更小”，从而减少了模型估计的方差并改善了泛化。此外，在具有特征数量多于观察数量的高维数据集中，正则化也是一种确保存在唯一解的方法。弹性网模型结合了系数向量的加权L1和L2惩罚项，前者可以导致稀疏性（即严格为零的系数），后者确保平滑的系数收缩。弹性网优化如下：
Elastic-net loss for the Cox model
β~=argminβ∑i=1Nδi{xiTβ−log[∑j∈R(ti)Nexp(xjTβ)]}+λ(α∣∣β∣∣1+0.5(1−α)∣∣β∣∣22)\tilde{\beta}=arg min_\beta \sum\limits_{i=1}^N \delta_i \left\{x_i^T\beta-log[\sum\limits_{j \in R(t_i)}^Nexp(x_j^T\beta)] \right\}+\lambda(\alpha||\beta||_1+0.5(1-\alpha)||\beta||_2^2)β~=argminβi=1∑Nδi{xiTβ−log[j∈R(ti)∑Nexp(xjTβ)]}+λ(α∣∣β∣∣1+0.5(1−α)∣∣β∣∣22)
其中，arg min 就是使后面这个式子达到最小值时的变量的取值，这里是β\betaβ的取值。超参数（hyperparameter）λ\lambdaλ定义为整个水平的正则化，而α\alphaα则定义为平衡稀疏解（α=1\alpha=1α=1，即Lasso模型）和零稀疏方法（α=0\alpha=0α=0，即Ridge模型）。超参数可以通过交叉验证的方法（cross-validation）最终确定。

基础代码

在survival package中，Surv()对象用于存储时间/事件信息，并可以使用以下函数转换为矩阵YYY。

risksets <- function(So) {n <- nrow(So)Y <- matrix(0,nrow=n, ncol=n)if (ncol(So) == 2) {endtime <- So[,1]event <- So[,2]for (i in seq(n)) {Y[i,] <- endtime[i] >= endtime}} else {starttime <- So[,1]endtime <- So[,2]event <- So[,3]for (i in seq(n)) {Y[i,] <- (endtime[i] >= endtime) & (starttime[i] < endtime)}}return(Y)
}

以下是两个不同时间过程的示例，一个是时间不变的（即开始时间始终为零），另一个是时间相关的。

So.ti <- Surv(time=c(1,2,3), event=c(0,1,1))
risksets(So.ti)

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    1    1    0
## [3,]    1    1    1

So.td <- Surv(time=c(0,1,0), time2=c(1,10,8), event=c(0,1,1))
risksets(So.td)

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    1
## [3,]    1    0    1

因为在~~文章~~ 中已经概述了使用近端梯度下降的细节，所以我将简要概括我们感兴趣的梯度更新目标。
Elastic-net Cox proximal update
β(k)=Sγαλ(β(k−1)+γ(k)[1NXT(δ−Pδ)−λ(1−α)β(k−1)])\beta^{(k)}=S_{\gamma \alpha \lambda}(\beta^{(k-1)}+\gamma^{(k)}[\cfrac{1}{N}X^T (\delta-P\delta)-\lambda(1-\alpha)\beta^{(k-1)}])β(k)=Sγαλ(β(k−1)+γ(k)[N1XT(δ−Pδ)−λ(1−α)β(k−1)])
其中γ(k)\gamma^{(k)}γ(k)是是每次迭代时的梯度步长。下面我们概述了此更新步骤的每个组件所需的代码：
S():S():S():

softhresh <- function(x,t) { sign(x) * pmax(abs(x) - t, 0) }

P:P:P:

Pfun <- function(Y,tY,eta) {rsk <- exp(eta)haz <- as.vector( tY %*% rsk )Pmat <- outer(rsk,haz,'/') * Yreturn(Pmat)
}

δ−Pδ:\delta-P\delta :δ−Pδ:

resfun <- function(X, b, Y, tY, l, a, P, d, ll) {eta <- as.vector(X %*% b)Phat <- P(Y,tY,eta)nll <- ll(Phat, d, b, l, a)res <- d - Phat %*% dreturn(list(res=res, nll=nll))
}

−1NXT(δ−Pδ)+λ(1−α)β(k−1)-\cfrac{1}{N}X^T (\delta-P\delta)+\lambda(1-\alpha)\beta^{(k-1)}−N1XT(δ−Pδ)+λ(1−α)β(k−1)

gradfun <- function(X, r, b, l, a) {      grad <- -t(X) %*% r / nrow(X) + l*(1-a)*breturn(grad)
}

等式（4）：

proxstep <- function(b, g, s, a, l) {btilde <- b - s * gb2 <- softhresh(btilde, a*l*s)  return(b2)
}

等式（3）：

llfun <- function(P,d,b,l,a) {    -mean(log(diag(P)[d==1])) + l*( a*sum(abs(b)) + (1-a)*sum(b^2)/2 )
}

（注释1）什么是独热编码（One-Hot）？
—————————————————————————————————————
One-Hot编码，又称为一位有效编码，主要是采用N位状态寄存器来对N个状态进行编码，每个状态都由他独立的寄存器位，并且在任意时候只有一位有效。
One-Hot编码是分类变量作为二进制向量的表示。这首先要求将分类值映射到整数值。然后，每个整数值被表示为二进制向量，除了整数的索引之外，它都是零值，它被标记为1。