050. Cox Regression

Cox Regression

Proportional Hazards Model

Proposed by Cox (1972, JRSS-B), primarily to model the relationship between hazard function and covariates. most cited paper in statistics ( 41; 000 as of April 2016), one of the most cited in science.

Several extensions to more complex data structures, e.g., clustered failure time data, or recurrent event data, etc.

※ Data Structure

Observed data: $\{X_i=T_i\wedge C_i,{\Delta }_i=I(T_i<C_i),{\mathbf{Z}}_i(\cdot )\}\stackrel{iid}{\sim }$

추가로 $N_i=I(X_i\le t,{\Delta }_i=1)$, $Z_i(t)$ = covariate vector (possibly time-dependent).

Cox PH Model

$${\lambda }_i(t)=\lambda (t|Z_i)={\lambda }_0(t)\exp ({\beta }^{\prime }{Z}_i)\text{\hspace{1em}(Cox Model)}$$

semiparametric model:

• $\exp ({\beta }^{\prime }{Z}_i)$, parametric assumption on covariate effects
• multiplicative model
• $\beta$ : $p\times 1$ vector, $p<\infty$
• ${\lambda }_0(t)$, nonparametric; is $\infty$ dimensional
• shape of hazard function is unspecified

Due to nonparametric component, standard maximum likelihood theory does not apply

Let $Z_{ij}$ be the $j$-th element of $Z_i$

• ${\beta }_j$ = difference in log hazards

• $\exp ({\beta }_j)$ = ratio of hazards; assumed constant for all $t$

• ${\lambda }_0(t)$: baseline hazard; common to all subjects, ${\lambda }_0(t)={\lambda }_i(t|Z_i=0)$

The hazard ratio, $\exp ({\beta }_j)$, is sometimes referred to as a relative risk

• risk = probability, not a rate
• hazard is a rate, not a probability
• in ratio of hazards, time dimension cancels out

Direction of effect:

$$\begin{array}{rlrlr}{\beta }_j>0:& & \uparrow {\lambda }_i& & \downarrow S_i(t)\\ {\beta }_j<0:& & \downarrow {\lambda }_i& & \uparrow S_i(t)\end{array}$$

Magnitude of effect is easy to interpret w.r.t. ${\lambda }_i(t)$

Cumulative hazard function:

$$\begin{array}{rl}{\lambda }_i(t)& ={\lambda }_0(t)\exp (\beta Z_i)\\ {\Lambda }_i(t)& ={\int }_0^t{\lambda }_0(s)\exp (\beta Z_i)ds\\ & ={\Lambda }_0(t)\exp (\beta Z_i)\end{array}$$

Survival function:

$$\begin{array}{rl}{S}_i(t)& =\exp \{-{\Lambda }_i(t)\}\\ & =\exp \{-{\Lambda }_0(t)\exp ({\beta }^{\prime }{Z}_i)\}\\ & =S_0(t{)}^{\exp \left\{{\beta }^{\prime }{Z}_i\right\}}\end{array}$$

By fitting a Cox model, one can readily interpret the multiplicative effect on the hazard:

• ex) randomized trial: treatment ($Z_i=1$) versus placebo ($Z_i=0$); $\hat{\beta }=0.405$ ($\exp (\hat{\beta })=1.5$)
• ${\lambda }_i(t)$ for treated patients is 50% more of that of the controls.
• irrespective of ${\lambda }_0(t)$

Nevertheless, ${\Lambda }_0(t)$ is required in order to determine $Z_i$‘s effect on $S_i(t)$, e.g.,

$$\begin{array}{rlrlr}S(t|Z_i=0)=0.95& & vs.& & S(t|Z_i=1)=0.93\\ S(t|Z_i=0)=0.70& & vs.& & S(t|Z_i=1)=0.59\end{array}$$

Cox Model: Independent Censoring

Independent censoring assumption is less stringent than in nonparametric estimation.

Assumption is often written as $T_i\perp C_i|Z_i$:

$$\begin{array}{rlrl}& \underset{\delta \to 0}{\lim }\frac{1}{\delta }P(t\le T_i<t+\delta |T_i\ge t,C_i\ge t,& & Z_i)\\ =& \underset{\delta \to 0}{\lim }\frac{1}{\delta }P(t\le T_i<t+\delta |T_i\ge t,& & Z_i)\end{array}$$

※ Note: $C_i$ is allowed to depend on $Z_i$

Semiparametric PH Model: General

• General expression for multiplicative proportional hazards model:

$${\lambda }_i(t)={\lambda }_0(t)g({\beta }^{\prime }{Z}_i)$$

$g(x)$ is link function, specified. $\forall x:g(x)\ge 0$, $\exists g^{\prime \prime }(x)$, and in special case, $g(x)=\exp (x)$.

• Other choices for link function (e.g., Self & Prentice, 1983): $g(x)=1+x=(1+x{)}^{-1}=\log (1+x)$

※ Notes:

• not all choices of $g(x)$ lead to clear interpretation of ${\beta }_j$
• certain choices of $g(x)$ lead to numerical issues; e.g., likelihood is flat; local maxima, etc.
• $g(x)\ne exp(x)$ has received little attention in the literature

Multiplicative Model

Cox model is a multiplicative model, i.e., covariates assumed to affect survival probability by multiplying the baseline hazard.

• Additive models also been proposed

Proportional Hazards Regression and Multiplicative Intensity Model

• Recall Counting process: martingale representation

$$\begin{array}{rl}N(t)& =I(X\le t,\Delta =1)\\ Y(t)& =I(X\ge t)\\ M(t)& =N(t)-{\int }_0^{t}Y(u){\lambda }_0(u)e^{{\beta }^{\prime }Z}du\\ {\mathcal{F}}_t& =\sigma \{N(u),Y(u+),Z:0\le u\le t\}\end{array}\text{\hspace{1em}(1)}$$

1. intensity $l(u)=Y(u){\lambda }_0(u)e^{{\beta }^{\prime }Z}$, therefore integrated form is cumulative intensity $A(t)$.

• Multiplicative Intensity Model:

$$l(t)=Y(t){\lambda }_0(t)e^{{\beta }^{\prime }Z(t)}$$

• Counting process: $N(t)$ = Number of events of a specified type that have occurred by time $t$

  • $N(t)$ may take more than one jump
  • multiple infections, repeated breakdowns, hospital admissions
  • $EN(t)<\infty$

• At-risk process: $Y(t)$, left-continuous process, $1$ if failure can be observed at time $t$, otherwise $0$.

  • $Y(t)$ can be used to represent situation in which a subject enter and exit risk sets several times
  • $Y(t)$ may be $1$ even after an observed failure

• Covariate process: $Z(t)$ = (bounded) predictable process

  • time-dependent treatment, risk factors
  • model checking and relaxing PH assumption

• Baseline hazard function: ${\lambda }_0(\cdot )$ = an arbitrary deterministic function

• Filtration: ${\mathcal{F}}_t=\sigma \{N(u),Y(u+),Z(u):0\le u\le t\}$

• Martingale: $M(t)=N(t)-{\int }_0^{t}l(u)du$

• Intensity function: E \Big \{ dN(t) \Big \| \mathcal F_{t-} \Big} = l(t) dt

• Data: $n$ independent observations on $\{N(\cdot ),Y(\cdot ),Z(\cdot )\}$

Likelihood; conditional, marginal and partial likelihoods

• $X=$ vector of observations; $f_X(x,\theta )=$ density of $X$

• $\theta =$ vector parameter; $\theta =({\beta }^{\prime },{\varphi }^{\prime }{)}^{\prime }$

• $\beta =$ parameter of interest; $\varphi =$ nuisance parameter

• likelihood: $f_X(x,\theta )=f_{W|V}(w|v,\theta )f_V(v,\theta )$

  • $X=(V^{\prime },W^{\prime }{)}^{\prime }$
  • infinite-dimensional $\varphi$
  • $f_{W|V}(w|v,\theta )$ does not involve $\varphi$ $\Rightarrow$ use $f_{W|V}(w|v,\beta )$ (conditional likelihood)
  • $f_V(v,\theta )$ does not involve $\varphi$ $\Rightarrow$ use $f_V(v,\beta )$ (marginal likelihood)

$$X=(V_1,W_1,\cdots ,V_K,W_K)$$ $$\begin{array}{rl}{f}_X(x,\theta )& =f_{V_1,W_1,\cdots ,V_K,W_K}(v_1,w_1,\cdots ,v_K,w_K;\theta )\\ & =f_{V_1}(v_1;\theta )f_{W_1|V_1}(w_1|v_1;\theta )f_{V_2|V_1,W_1}(v_2|v_1,w_1;\theta )\times \cdots \\ & =\left\{\prod _{i=1}^K{f}_{W_i|Q_i}(w_i|q_i;\theta )\right\}\left\{\prod _{i=1}^K{f}_{V_i|P_i}(v_i|p_i;\theta )\right\}\end{array}$$ $$\begin{array}{rlrl}{P}_1=\varphi ,& & & P_i=(V_1,W_1,\cdots ,V_{i-1},W_{i-1})\\ Q_1=V1,& & & Q_i=(V_1,W_1,\cdots ,W_{i-1},V_i)\end{array}$$

${\prod }_{i=1}^K{f}_{W_i|Q_i}(w_i|q_i;\theta )$ is free of $\varphi$ $\Rightarrow$ use ${\prod }_{i=1}^K{f}_{W_i|Q_i}(w_i|q_i;\beta )$ (partial likelihood)

Partial & Marginal Likelihoods

Focus on Proportional Hazards Model: i.e., $(X_i,{\delta }_i,Z_i),i=1,\cdots ,n$ ($n$ independent triplets)

$$\begin{array}{rlrl}& \lambda (t|Z)={\lambda }_0(t)e^{{\beta }^{\prime }Z}& & S(t|Z)=\{S_0(t){\}}^{e^{{\beta }^{\prime }Z}}\end{array}\text{\hspace{1em}(1)}$$

위에서 ${\lambda }_0(t)$는 unspecified.

• Partial Likelihood: assume no ties, absolutely continuous failure distribution

Suppose there are L observed failures at ${\tau }_1<\cdots <{\tau }_L$ (set ${\tau }_0\equiv 0$ & ${\tau }_{L+1}\equiv \infty$)

16.png

Let (i) be the label for individual failing at ${\tau }_i$ (set $(L+1)\equiv n+1$). Note $t_{(i)}={\tau }_i$

Covariates for $L$ failures: $(Z_{(1)},\cdots ,Z_{(L)})$. (Hereafter, condition on $\{Z_i:i=1,\cdots ,n\}$)

Censorship times in $[{\tau }_i;{\tau }_{i+1})$: $({\tau }_{i1},\cdots ,{\tau }_{i,m_i})$ with covariates $(Z_{(i,1)},\cdots ,Z_{(i,m_i)})$, i.e., $(i,j)$ is label for item censored at ${\tau }_{ij}$

17.png

The data can be divided into sets

$$(V_1,W_1,\cdots ,V_{L+1},W_{L+1})$$

where, for $i=1,\cdots ,L,L+1$,

$$\begin{array}{rl}{V}_i& =\{{\tau }_i,{\tau }_{i-1,j};(i-1,j):j=1,\cdots ,m_{i-1}\}\\ and{W}_i& =\{(i)\}\end{array}$$

18.png

19.png

GOAL: Build a likelihood on a subset of the full data set

• carrying most of the information about $\beta$
• carrying no information on nuisance parameters $\{{\lambda }_0(t):t\ge 0\}$

PROPOSAL: Generate likelihood of $\{W_1,\cdots ,W_L\}$

JUSTIFICATION, WHY?:

• Timing of events $\{{\tau }_1,\cdots ,{\tau }_L\}$ can be explained by ${\lambda }_0(\cdot )$.
• Censoring times and labels can be ignored if we assume non-informative censorship (independent censoring).

So this is a partial likelihood in the sense that it is only part of the likelihood of the observed data.

If $Q_i\equiv (V_1,W_1,\cdots ,V_{i-1},W_{i-1},V_i)$ and ${\mathcal{F}}_{{\tau }_i}\equiv (Q_i,Z)$, the partial likelihood is ${\prod }_{i=1}^{L}P(W_i=(i)|{\mathcal{F}}_{{\tau }_i})$, i.e., given the risk set at ${\tau }_i$, and given event occurs at ${\tau }_i$.

Denote $R_i\equiv \{j:X_j\ge {\tau }_i\}$ as risk set at ${\tau }_i$. Then, by the assumption of independent censoring,

$$\begin{array}{rl}P(W_i=(i)|{\mathcal{F}}_{{\tau }_i})& =\frac{P\{t_{(i)}\in [{\tau }_i,{\tau }_i+d\tau \left)\right|{\mathcal{F}}_{{\tau }_i}\}\cdot \prod _{j\in R_i-(i)}P\{t_j\notin [{\tau }_i,{\tau }_i+d\tau )\left|{\mathcal{F}}_{{\tau }_i}\right\}}{\sum _{l\in R_i}\left[P\{t_l\in [{\tau }_i,{\tau }_i+d\tau \left)\right|{\mathcal{F}}_{{\tau }_i}\}\cdot \prod _{j\in R_i-l}P\{t_j\notin [{\tau }_i,{\tau }_i+d\tau )\left|{\mathcal{F}}_{{\tau }_i}\right\}\right]}\\ & \\ & \\ & =\frac{d\Lambda \left({\tau }_i\right|Z_{(i)}\left)\prod _{j\in R_i-(i)}\right\{1-d\Lambda \left({\tau }_i\right|Z_j\left)\right\}}{\sum _{l\in R_i}\left[d\Lambda \left({\tau }_i\right|Z_l\left)\prod _{j\in R_i-l}\right\{1-d\Lambda \left({\tau }_i\right|Z_j\left)\right\}\right]}÷\frac{d{\tau }_i}{d{\tau }_i}\\ & \\ & \\ & =\frac{\lambda \left({\tau }_i\right|Z_{(i)})}{\frac{P\{T\in [t,t+dt)|T\ge t,Z\}}{dt}=\sum _{l\in R_i}\left[\lambda \left({\tau }_i\right|Z_l)\right]}\stackrel{(1)}{=}\frac{\exp ({\beta }^{\prime }{Z}_{(i)})}{\sum _{l\in R_i}\exp ({\beta }^{\prime }{Z}_l)}\end{array}\text{\hspace{1em}(a)(2)}$$

• at (a), $P\{t_j\notin [{\tau }_i,{\tau }_i+d\tau \left)\right|{\mathcal{F}}_{{\tau }_i}\}=1-P\{t_j\in [{\tau }_i,{\tau }_i+d\tau )\left|{\mathcal{F}}_{{\tau }_i}\right\}$
• at (2), $d\Lambda \left({\tau }_i\right|Z_j)=0$

Thus, the Partial Likelihood is

$$\prod _{i=1}^L\frac{\exp ({\beta }^{\prime }{Z}_{(i)})}{\sum _{l\in R_i}\exp ({\beta }^{\prime }{Z}_l)}=L(\beta )\text{\hspace{1em}(3)}$$

Note: unspecified ${\lambda }_0(\cdot )$ + noninformative censoring $\Rightarrow$ $\prod _{i=1}^L{f}_{V_i|P_i}(v_i|p_i;\theta )$ contains little or no information about $\beta$.

• Counting process notation:

$$\begin{array}{rlr}L(\beta )=\prod _{i=1}^n\prod _{t\ge 0}{\left\{\frac{\exp ({\beta }^{\prime }{Z}_i)}{\sum _{j=1}^n{Y}_j(t)\exp ({\beta }^{\prime }{Z}_j)}\right\}}^{d{N}_i(t)},& & d{N}_i(t)=\{\begin{array}{ll}1& N_i(t)-N_i(t-)=1\\ 0& o.w.\end{array}\end{array}$$

• Maximum partial likelihood estimator (MPLE): $L(\hat{\beta })={\max }_{\beta }L(\beta )$ (using Newton-Raphson (NR) algorithm)

  • Specifically, the log partial likelihood is then
  $$l(\beta )=\sum _{i=1}^n{\int }_0^{\infty }\left[Y_i(t)Z_i\beta -\log \left(\sum _{j=1}^n{Y}_j(t)\exp ({\beta }^{\prime }{Z}_j)\right)\right]d{N}_i(t)$$
  • The score vector, $U(\beta )$, can be obtained by differentiating $l(\beta )$ w.r.t. $\beta$:
  $$\begin{array}{rlrl}U(\beta )& =\sum _{i=1}^n{\int }_0^{\infty }\{Z_i-\bar{Z}(\beta ,t)\}& & d{N}_i(t)\\ & =\sum _{i=1}^n{\int }_0^{\infty }\left\{Z_i-\frac{{\sum }_{i=1}^n{Y}_i(t)Z_i\exp ({\beta }^{\prime }{Z}_i)}{{\sum }_{i=1}^n{Y}_i(t)\exp ({\beta }^{\prime }{Z}_i)}\right\}& & d{N}_i(t)\end{array}$$

  • where $\bar{Z}(\beta ,t)$ is a weighted mean of $Z$ over those observations still at risk at time $t$.

  • The information matrix, $\mathcal{I}(\beta )$, is the negative second derivative where

  $$\begin{array}{rl}\mathcal{I}(\beta )& =\sum _{i=1}^n{\int }_0^{\infty }V(\beta ,t)d{N}_i(s)\\ & \\ V(\beta ,t)& =\frac{\sum _{i=1}^n{Y}_i(t)\exp ({\beta }^{\prime }{Z}_i)\{Z_i-\hat{Z}(\beta ,t){\}}^{\prime }\{Z_i-\hat{Z}(\beta ,t)\}}{\sum _{i=1}^n{Y}_i(t)\exp ({\beta }^{\prime }{Z}_i)}\end{array}$$
  • and $V(\beta ,t)$ is the weighted variance of $Z$ at time $t$.

Then, the MPLE, $\hat{\beta }$, is found by solving the partial likelihood equation: $U(\hat{\beta })=0$.

Under some regularity conditions, $\hat{\beta }$ is consistent and asymptotically normally distributed with mean $\beta$ and variance $E\{\mathcal{I}(\beta ){\}}^{-1}$ (will be shown later.)

The NR algorithm to solve the partial likelihood equation: Compute iteratively until convergence (requires an initial value ${\hat{\beta }}^{(0)}$).

$${\hat{\beta }}^{(n+1)}={\hat{\beta }}^{(n)}+{\mathcal{I}}^{-1}\left({\hat{\beta }}^{(n)}\right)\cdot U\left({\hat{\beta }}^{(n)}\right)$$

※ Note:

1. (incredibly) Robust algorithm!
2. ${\hat{\beta }}^{(0)}=0$ usually works.

Cox Proportional Hazards Model

Cox model:

$$\begin{array}{rl}{\lambda }_i(t)=\lambda (t|Z_i)& ={\lambda }_0(t)\exp ({\beta }^{\prime }{Z}_i)\\ & ={\lambda }_0(t)\exp ({\beta }_1{Z}_{i1}+\cdots +{\beta }_k{Z}_{ik})\\ & \Updownarrow \\ \log \lambda (t|Z_i)& =\log [{\lambda }_0(t)]+{\beta }_1{Z}_{i1}+\cdots +{\beta }_k{Z}_{ik}\\ S(t|Z_i)& =[S_0(t){]}^{\exp ({\beta }_1{Z}_{i1}+\cdots +{\beta }_k{Z}_{ik})}\end{array}$$

※ Note:

$$\begin{array}{rl}{\lambda }_0(t)& =\lambda (t|Z_1=\cdots =Z_k=0)\\ & \\ \exp ({\beta }_1{Z}_1+\cdots +{\beta }_k{Z}_k)& =RR\\ & =\frac{\lambda (t|Z_1,\cdots ,Z_k)}{\lambda (t|Z_1=\cdots =Z_k=0)}\end{array}\text{\hspace{1em}(1)}$$

• (1) is relative risk of hazard of death comparing covariates values $Z_1,\cdots ,Z_k$ to $Z_1=\cdots =Z_k=0$

Interpreting Cox Model Coeffcients: ${\beta }_k$ is the log RR (hazard ratio) for a unit change in $Z_k$, given all other covariates remain constant, i.e.,

$$\begin{array}{rl}\frac{\lambda \left[t\right|Z_1,\cdots ,(Z_{k^{\prime }}+1),\cdots ,Z_k]}{\lambda \left[t\right|Z_1,\cdots ,Z_{k^{\prime }},\cdots ,Z_k]}& =\exp ({\beta }_1\cdot 0+\cdots +{\beta }_{k^{\prime }}\cdot (Z_{k^{\prime }}+1-Z_{k^{\prime }})+\cdots +{\beta }_k\cdot 0)\\ & =\exp ({\beta }_{k^{\prime }})\end{array}$$

The RR comparing 2 sets of values for the covariates $(Z_1,\cdots ,Z_k)$ vs. $(Z_1^{\prime },\cdots ,Z_k^{\prime })$:

$$RR=\frac{\lambda (t|Z_1,\cdots ,Z_k)}{\lambda (t|Z_1^{\prime },\cdots ,Z_k^{\prime })}=\exp \{{\beta }_1(Z_1-Z_1^{\prime })+\cdots +{\beta }_k(Z_k-Z_k^{\prime })\}$$

20.png

Comparison of Nested Models

• Nested Models:

$$\begin{array}{rl}\lambda (t)& ={\lambda }_0(t)\exp ({\beta }_1{Z}_1+\cdots {\beta }_p{Z}_p+{\beta }_{p+1}{Z}_{p+1}+\cdots +{\beta }_k{Z}_k)\\ & ={\lambda }_0(t)\exp ({\beta }_1{Z}_1+\cdots {\beta }_p{Z}_p)\end{array}\text{\hspace{1em}(Full Model)(Reduced Model)}$$

To test:

• Nested Models:

$$\begin{array}{rlrlrlrl}& H_0:& & RM& & \iff & & H_0:{\beta }_{p+1}=\cdots ={\beta }_k=0\\ & H_A:& & RM& & \iff & & H_A:\ne \text{ somewhere}\end{array}$$

Use the partial likelihood ratio statistic, $X_{Cox}^2=-2[\log PL(RM)-\log PL(FM)]$.

Under $H_0$: Reduced model, and when $n$ is large:

\begin{align} X^2_{Cox} \sim \chi^2_{k-p} && k-p \text{ is the ## of parameters set to 0 by }H_0 \end{align}

20.png, 21.png

Stratification

Two Ways to Stratify. Suppose a confounder $C$ has 3 levels on which we would like to stratify when comparin g $\lambda (t|E)$ and $\lambda (t|\bar{E})$. How? $X_E=\{\begin{array}{llc}1& E& \text{(exposed)}\\ 0& \bar{E}& \text{(not exposed)}\end{array}$

22.png

• Which Way to Stratify?

1. Under dummy variable stratification model, the proportional stratum-to-stratum hazards assumption may not be correct. If not, the con-founder $C$ may be inadequately controlled.
2. Proportionality assumption can be checked using time-dependent covariates.
3. True stratification is a more thorough adjustment, as long as observations within each level are homogeneous. If $C$ can be measured continuously and the strata were formed by grouping values of it, better control for $C$ might be achieved with continuous (could be time-dependent) covariate adjustment.
4. If $C$ is controlled using the true stratification there is no way to estimate one summary relative risk comparing two levels of $C$. However, we can estimate ${\lambda }_{0i}(t)$ for each stratum then we can estimate a RR function.
5. True stratification generally requires more data to obtain the same precision in coefficient estimates.

23.png

24.png

Test statistics

The standard asymptotic likelihood inference tests, Wald, score, and likelihood ratio (LR), still can be applied for the Cox partial likelihood.

25.png

Their finite sample properties may differ; in general, the LRT is the most reliable, the Wald test is the least.

26.png

When $p=1$ and the single covariate is categorical, the score test is identical to the log-rank test.

27.png

Handling ties

Real data sets often contain tied event times.

• When do we have ties?

1. Continuous event times are grouped into intervals.
2. Event time scale is discrete.

Four commonly used ways of handling ties: 1) Breslow approximation, 2) Efron approximation, 3) Exact partial likelihood, and 4) Averaged likelihood.

When the underlying time is continuous but ties are generated due to a grouping, the contribution to the partial likelihood for the $i$-th event at time $t_i$ is $\frac{\exp ({\beta }^{\prime }{Z}_i)}{\sum _{j\in R_i}{Y}_j(t_i)\exp ({\beta }^{\prime }{Z}_j)}$

Two commonly used methods are

1. Breslow approximation
2. Efron approximation

Example: Assume 5 subjects are at risk of dying at time $t$ and two die at the same time $t$ (because of grouping of time) If the time data had been more precise, then the first two terms in the likelihood would be either

28.png

29.png

30.png