Processing math: 100%

Proportional Hazards

We will explore the relationship between survival and explanatory variables by modeling. In this class, we consider a broad class of regression models: Proportional Hazards (PH) models. logλ(t;Z)=logλ0(t)+βZ Suppose Z=1 for treated subjects and Z=0 for untreated subjects. Then this model says that the hazard is increased by a factor of eβ for treated subjects versus untreated subjects (eβ might be < 1).

This group of PH models divides further into:

  • Parameteric PH models: Assume parametric form for λ0(t), e.g., Weibull distribution.
  • Semi-parametric PH models: No assumptions on λ0(t): Cox PH model.
+/- r Code
library(survival)
library(KMsurv)

Cox Proportional Hazards model

λ(t;Z)=λ0(t)exp{βZ}

Why do we call it proportional hazards? Think of the first example, where Z=1 for treated and Z=0 for control. Then if we think of λ1(t) as the hazard rate for the treated group, and λ0(t) as the hazard for control, then we can write: λ1(t)=λ(t;Z=1)=λ0(t)exp{βZ}=λ0(t)exp{β} This implies that the ratio of the two hazards is a constant, ϕ, which does NOT depend on time, t. In other words, the hazards of the two groups remain proportional over time. ϕ=λ1(t)λ0(t)=eβ ϕ is referred to as the hazard ratio.

Baseline Hazard Function

In the example of comparing two treatment groups, λ0(t) is the hazard rate for the control group.

In general, λ0(t) is called the baseline hazard function, and reflects the underlying hazard for subjects with all covariates Z1,,Zp equal to 0 (i.e., the 'reference group').

The general form is: λ(t;Z)=λ0(t)exp{β1Z1++βpZp} So when we substitute all of the Zj equal to 0, we get λ(t;Z)=λ0(t).

In the general case, we think of the i-th individual having a set of covariates Zi=(Z1i,Z2i,,Zpi), and we model their hazard rate as some multiple of the baseline hazard rate: λ(t;Zi)=λ0(t)exp{β1Z1i++βpZpi} This means we can write the log of the hazard ratio for the i-th individual to the reference group as: logλ(t;Zi)λ0(t)=β1Z1i++βpZpi The Cox Proportional Hazards model is a linear model for the log of the hazard ratio.

One of the biggest advantages of the framework of the Cox PH model is that we can estimate the parameters β which reflect the effects of treatment and other covariates without having to make any assumptions about the form of λ0(t).

In other words, we don???t have to assume that λ0(t) follows an exponential model, or a Weibull model, or any other particular parametric model. That???s what makes the model semi-parametric.

Questions:

  1. Why don't we just model the hazard ratio, ϕ=λi(t)λ0(t), directly as a linear function of the covariates Z?
  2. Why doesn't the model have an intercept?

Parameters of interest, β, are estimated via partial likelihood. λ0() is treated as a nuisance parameter.

+/- r Code
leuk = read.table("http://www.karlin.mff.cuni.cz/~pesta/prednasky/NMFM404/Data/leuk.dat", sep=" ", head=T, skip=7) # read data
head(leuk)
+/- Output
##   weeks remiss trtmt
## 1     6      0     1
## 2     6      1     1
## 3     6      1     1
## 4     6      1     1
## 5     7      1     1
## 6     9      0     1
+/- r Code
# Cox PH Model
leuk.ph = coxph(Surv(weeks,remiss) ~ trtmt, data=leuk)
# Note: default = Efron method for handling ties
summary(leuk.ph)
+/- Output
## Call:
## coxph(formula = Surv(weeks, remiss) ~ trtmt, data = leuk)
## 
##   n= 42, number of events= 30 
## 
##          coef exp(coef) se(coef)      z Pr(>|z|)    
## trtmt -1.5721    0.2076   0.4124 -3.812 0.000138 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## trtmt    0.2076      4.817   0.09251    0.4659
## 
## Concordance= 0.69  (se = 0.053 )
## Rsquare= 0.322   (max possible= 0.988 )
## Likelihood ratio test= 16.35  on 1 df,   p=5.261e-05
## Wald test            = 14.53  on 1 df,   p=0.0001378
## Score (logrank) test = 17.25  on 1 df,   p=3.283e-05
+/- r Code
# Breslow handling of ties
leuk.phb = coxph(Surv(weeks, remiss) ~ trtmt, data=leuk, method="breslow")
summary(leuk.phb)
+/- Output
## Call:
## coxph(formula = Surv(weeks, remiss) ~ trtmt, data = leuk, method = "breslow")
## 
##   n= 42, number of events= 30 
## 
##          coef exp(coef) se(coef)      z Pr(>|z|)    
## trtmt -1.5092    0.2211   0.4096 -3.685 0.000229 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## trtmt    0.2211      4.523   0.09907    0.4934
## 
## Concordance= 0.69  (se = 0.053 )
## Rsquare= 0.304   (max possible= 0.989 )
## Likelihood ratio test= 15.21  on 1 df,   p=9.615e-05
## Wald test            = 13.58  on 1 df,   p=0.0002288
## Score (logrank) test = 15.93  on 1 df,   p=6.571e-05
+/- r Code
# Exact handling of ties
leuk.phe = coxph(Surv(weeks, remiss) ~ trtmt, data=leuk, method="exact")
summary(leuk.phe)
+/- Output
## Call:
## coxph(formula = Surv(weeks, remiss) ~ trtmt, data = leuk, method = "exact")
## 
##   n= 42, number of events= 30 
## 
##          coef exp(coef) se(coef)      z Pr(>|z|)    
## trtmt -1.6282    0.1963   0.4331 -3.759  0.00017 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## trtmt    0.1963      5.095   0.08398    0.4587
## 
## Rsquare= 0.321   (max possible= 0.98 )
## Likelihood ratio test= 16.25  on 1 df,   p=5.544e-05
## Wald test            = 14.13  on 1 df,   p=0.0001704
## Score (logrank) test = 16.79  on 1 df,   p=4.169e-05

Adjustments for ties

The proportional hazards model assumes a continuous hazard ??? ties are not possible. There are four proposed modifications to the likelihood to adjust for ties.

Bottom Line: Implications of Ties (See Allison (1995), p.127-137)

  • When there are no ties, all options give exactly the same results.
  • When there are only a few ties, it won???t make much difference which method is used. However, since the exact methods won't take much extra computing time, you might as well use one of them.
  • When there are many ties (relative to the number at risk), the Breslow option (default) performs poorly (Farewell & Prentice, 1980; Hsieh, 1995). Both of the approximate methods, Breslow and Efron, yield coefficients that are attenuated (biased toward 0).
  • The choice of which exact method to use should be based on substantive grounds - are the tied event times truly tied? ... or are they the result of imprecise measurement?
  • Computing time of exact methods is much longer than that of the approximate methods. However, in most cases it will still be less than 30 seconds even for the exact methods.
  • Best approximate method - the Efron approximation nearly always works better than the Breslow method, with no increase in computing time, so use this option if exact methods are too computer-intensive.

Confidence intervals and hypothesis tests

Many software packages provide estimates of β, but the hazard ratio HR=exp{β} is usually the parameter of interest.

For each covariate of interest, the null hypothesis is H0:HRj=1βj=0 Wald test is used. For nested models, a likelihood ratio test is constructed.

MAC prevention clinical trial data

ACTG 196 was a randomized clinical trial to study the effects of combination regimens on prevention of MAC (mycobacterium avium complex), one of the most common opportunistic infections in AIDS patients.

The treatment regimens were:

  • clarithromycin (new)
  • rifabutin (standard)
  • clarithromycin plus rifabutin
+/- r Code
mac = read.csv("http://www.karlin.mff.cuni.cz/~pesta/prednasky/NMFM404/Data/mac.csv", head=F, col.names =
    c("patid", "age", "agecat", "sex", "karnof", "ivdrug",
    "antiret", "cd4", "cd4cat", "ctg", "dthstat", "dthtime",
    "macstat", "mactime", "swdrstat", "swdrtime", "rif",
    "clari", "toxstat", "toxtime"))
# mactime = time to mac
# macstat = event indicator
head(mac)     # Check out a few data rows:
+/- Output
##   patid age agecat sex karnof ivdrug antiret cd4 cd4cat ctg dthstat
## 1     1  42      1   0     90      0       1   8      0   1       0
## 2     2  33      0   0     90      0       0  30      1   1       0
## 3     3  39      1   0    100      0       1  80      1   1       1
## 4     4  35      0   0     80      0       0  58      1   1       0
## 5     5  42      1   0     90      0       0  59      1   1       0
## 6     6  45      1   0     90      0       0  18      0   1       1
##   dthtime macstat mactime swdrstat swdrtime rif clari toxstat toxtime
## 1     623       1     560        1      106   1     0       0     623
## 2     651       0     651        1       86   1     0       0     651
## 3     464       0      26        0       26   1     0       0      26
## 4     622       0     622        1      111   0     1       0     501
## 5     643       0     643        1      119   0     1       0     643
## 6     216       0     171        0      171   0     1       1     171
+/- r Code
# Treatment arms:
#  rif + clari (N=389)
#  clari (N=397)
#  rif (N=391)
table(mac$clari, mac$rif)
+/- Output
##    
##       0   1
##   0 389 391
##   1 397   0
+/- r Code
# Fit and plot KM survival curves
fit.km = survfit(Surv(mactime, macstat) ~ clari + rif, data=mac)
plot of chunk ph-mac-km
+/- r Code
plot(fit.km, mark.time=F, lty=1:3, xlab="Days", ylab="MAC-Free Survival")
legend(lty=1:3, x=100, y=.4, legend=c("Rif + Clari", "Rifabutin", "Clarithromycin"), cex=0.8)

Model 1

+/- r Code
# Fit Cox PH model
fit.mac1 = coxph(Surv(mactime,macstat) ~ karnof + clari + rif, data=mac)
summary(fit.mac1)
+/- Output
## Call:
## coxph(formula = Surv(mactime, macstat) ~ karnof + clari + rif, 
##     data = mac)
## 
##   n= 1177, number of events= 121 
## 
##            coef exp(coef) se(coef)      z Pr(>|z|)    
## karnof -0.04485   0.95614  0.01064 -4.217 2.48e-05 ***
## clari   0.27557   1.31728  0.25801  1.068 0.285509    
## rif     0.87197   2.39161  0.23694  3.680 0.000233 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## karnof    0.9561     1.0459    0.9364    0.9763
## clari     1.3173     0.7591    0.7944    2.1842
## rif       2.3916     0.4181    1.5032    3.8051
## 
## Concordance= 0.649  (se = 0.028 )
## Rsquare= 0.027   (max possible= 0.73 )
## Likelihood ratio test= 32.02  on 3 df,   p=5.193e-07
## Wald test            = 32.29  on 3 df,   p=4.548e-07
## Score (logrank) test = 33.16  on 3 df,   p=2.977e-07

Model 2

+/- r Code
# Fit Cox PH model adding baseline CD4 variable
fit.mac2 = coxph(Surv(mactime,macstat) ~ karnof + clari + rif + cd4, data=mac)
summary(fit.mac2)
+/- Output
## Call:
## coxph(formula = Surv(mactime, macstat) ~ karnof + clari + rif + 
##     cd4, data = mac)
## 
##   n= 1177, number of events= 121 
## 
##             coef exp(coef)  se(coef)      z Pr(>|z|)    
## karnof -0.036874  0.963798  0.010665 -3.457 0.000546 ***
## clari   0.252345  1.287041  0.258337  0.977 0.328664    
## rif     0.879749  2.410294  0.237092  3.711 0.000207 ***
## cd4    -0.018360  0.981807  0.003684 -4.984 6.23e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## karnof    0.9638     1.0376    0.9439    0.9842
## clari     1.2870     0.7770    0.7757    2.1354
## rif       2.4103     0.4149    1.5145    3.8360
## cd4       0.9818     1.0185    0.9747    0.9889
## 
## Concordance= 0.716  (se = 0.028 )
## Rsquare= 0.053   (max possible= 0.73 )
## Likelihood ratio test= 63.77  on 4 df,   p=4.682e-13
## Wald test            = 55.59  on 4 df,   p=2.449e-11
## Score (logrank) test = 56.22  on 4 df,   p=1.806e-11

We can also compute the hazard ratio ourselves, by exponentiating the β coefficients: HRcd4=exp{???0.01835}=0.98 Why is this HR so close to 1, and yet still highly significant? What is the interpretation of this HR?

The likelihood ratio test for the effect of CD4 is twice the difference in minus log-likelihoods between the two models (on 1 df): ??2LR=2???(63.8???32.0)=31.8 How does this test statistic compare to the Wald χ2 test?

In the mac study, there were three treatment arms (rif, clari, and the rif+clari combination (for rif=clari=0). The combination therapy is the reference group, so the coefficients of rif and clari are in comparison to the combination group.

How do we compare the Rifa arm against the Clari arm? That is equivalent to testing H0:βclari=βrifa For that, we need the matrix varˆβ.

+/- r Code
names(fit.mac2)   # check components of fit object
+/- Output
##  [1] "coefficients"      "var"               "loglik"           
##  [4] "score"             "iter"              "linear.predictors"
##  [7] "residuals"         "means"             "concordance"      
## [10] "method"            "n"                 "nevent"           
## [13] "terms"             "assign"            "wald.test"        
## [16] "y"                 "formula"           "call"
+/- r Code
beta_hat = fit.mac2$coef
Var_beta = fit.mac2$var   # Var(beta_hat)
A = c(0,1,-1,0)   # contrast vector
beta_dif = t(A) %*% beta_hat     # beta_clari - beta_rif
var_dif = t(A) %*% Var_beta %*% A   # Var(beta_clari-beta_rif)
se_dif = sqrt(var_dif)        # SE(beta_clari-beta_rif)
t_dif = beta_dif/se_dif
t_dif     # Wald test t-statistic, df=n-p
+/- Output
##           [,1]
## [1,] -2.959871
+/- r Code
2* pt( -abs(t_dif), df=1177-4) # p-value for difference in arms
+/- Output
##             [,1]
## [1,] 0.003139506

Testing Model 2 for overall treatment effect Two ways:

  1. Wald test. βrif=0&βclari=0. Null distribution is Chi-square with df=2 (why?).
    +/- r Code
    A = matrix(0, ncol=4, nrow=2)
    A[1,2] = A[2,3] = 1
    beta_trt = A %% beta_hat
    var_trt = A %% Var_beta %% t(A)
    (Wald_trt = t(beta_trt) %% solve(var_trt) %*% beta_trt)
    
    +/- Output
    ##          [,1]
    ## [1,] 17.00252
    
    +/- r Code
    1-pchisq(Wald_trt, df=2)  # p-value for Wald test for treatment
    
    +/- Output
    ##              [,1]
    ## [1,] 0.0002032125
    
  2. Likelihood ratio test. Fit model with and without rif, clari and compare log-likelihood, against Chi-square with df=2.
    +/- r Code
    fit.mac3 = update(fit.mac2, .~. -rif -clari)
    anova(fit.mac3, fit.mac2, test="Chisq")
    
    +/- Output
    ## Analysis of Deviance Table
    ##  Cox model: response is  Surv(mactime, macstat)
    ##  Model 1: ~ karnof + cd4
    ##  Model 2: ~ karnof + clari + rif + cd4
    ##    loglik  Chisq Df P(>|Chi|)
    ## 1 -747.07
    ## 2 -738.62 16.914 2 0.0002125 *** ## --- ## Signif. codes: 0 '**' 0.001 '' 0.01 '' 0.05 '.' 0.1 ' ' 1

Predicted survival

For the i-th patient with covariates Zi, we have: Si(t)=[S0(t)]exp{βZi} Say we are interested in the survival pattern for single males in the nursing home study. Based on the previous formula, if we had an estimate for the survival function in the reference group, i.e., ˆS0(t), we could get estimates of the survival function for any set of covariates Zi.

How can we estimate the survival function, S0(t)? We could use the KM estimator, but there are a few disadvantages of that approach:

  • It would only use the survival times for observations contained in the reference group, and not all the rest of the survival times.
  • It would tend to be somewhat choppy, since it would reflect the smaller sample size of the reference group.
  • It???s possible that there are no subjects in the dataset who are in the 'reference' group (ex. say covariates are age and sex; there is no one of age=0 in our dataset).

Instead, we will use a baseline hazard estimator which takes advantage of the proportional hazards assumption to get a smoother estimate. ˆSi(t)=[ˆS0(t)]exp{ˆβZi} Using the above formula, we substitute ˆβ based on fitting the Cox PH model, and calculate ˆS0(t) by one of the following approaches:

  • Breslow estimator (R) ... extending the concept of the Nelson-Aalen estimator to the proportional hazards model
  • Kalbfleisch/Prentice estimator (R, SAS) ... analogous to the Kaplan-Meier estimator

Consider the Nursing Home study. We fit the model λ(t)=λ0(t)exp{β1[married]+β2[health status]} Let's predict the survival for single individuals with health status = 5 (worst), and for those with health status = 4 (second worst).

+/- r Code
nurshome = read.table("http://www.karlin.mff.cuni.cz/~pesta/prednasky/NMFM404/Data/NursingHome.dat", head=F, skip=14, col.names = c("los", "age", "rx", "gender", "married", "health", "censor"))
nurshome$discharged = 1 - nurshome$censor # event indicator
head(nurshome)
+/- Output
##    los age rx gender married health censor discharged
## 1   37  86  1      0       0      2      0          1
## 2   61  77  1      0       0      4      0          1
## 3 1084  75  1      0       0      3      1          0
## 4 1092  77  1      0       1      2      1          0
## 5   23  86  1      0       0      4      0          1
## 6 1091  71  1      1       0      3      1          0
+/- r Code
table(nurshome$health)
+/- Output
## 
##   2   3   4   5 
## 343 576 513 169
+/- r Code
# Fit Cox PH model
fit.nurs = coxph(Surv(los,discharged) ~ married + health, data=nurshome)
# Predict survival for single patients with health = 5
newdat = data.frame(married=0, health=5)
newdat    # data frame with same predictors as fit.nurs
+/- Output
##   married health
## 1       0      5
+/- r Code
ps5 = survfit(fit.nurs, newdata=newdat)    # predicted survival
# Compare with Kaplan-Meier
nursh5 = nurshome[nurshome$health==5 & nurshome$married==0,]
fit.km5 = survfit(Surv(los, discharged) ~ 1, data = nursh5)
# Predict survival for single patients, health =4
newdat[1,] = c(0,4)
ps4 =  survfit(fit.nurs, newdata=newdat)
# Compare with Kaplan-Meier
nursh4 = nurshome[nurshome$health==4 & nurshome$married==0,]
fit.km4 = survfit(Surv(los, discharged) ~ 1, data = nursh4)
plot of chunk ph-nurs-predict
+/- r Code
plot(ps5, xlab= "Length of stay in nursing home (days)", ylab = "Proportion not discharged", col=2, conf.int=F)
lines(fit.km5, mark.time=F, col=2)     # add line to existing plot
lines(ps4, xlab= "Length of stay in nursing home (days)", ylab = "Proportion not discharged", col=4, conf.int=F)
lines(fit.km4, mark.time=F, col=4)     # add line to existing plot

Continuing the Nursing Home example, let's get and plot the baseline survival S0(t) and cumulative hazard Λ0(t). They correspond to married=0, health=0.

+/- r Code
newdat[1,] = c(0,0)
ps0 = survfit(fit.nurs, newdata=newdat)   # S_0 = ps0$surv
bh = basehaz(fit.nurs, centered=F)
Lambda0 = bh$hazard   # it???s *cummulative* hazard!
# same as Lambda0 = -log(ps0$surv)
plot of chunk ph-nurs-baseline
+/- r Code
plot(bh$time, Lambda0, type="l", xlab="Length of stay in nursing home (days)", ylab="Cummulative hazard", col=4)

Model selection

Collett (Section 3.6) has an excellent discussion of various approaches for model selection. In practice, model selection proceeds through a combination of

  • knowledge of the science
  • trial and error, common sense
  • automatic variable selection procedures
    • forward selection
    • backward selection
    • stepwise selection

Many advocate the approach of first doing a univariate analysis to 'screen' out potentially significant variables for consideration in the multivariate model (see Collett).

Moreover, typically univariate analysis is a part of a larger analysis, in order to identify the importance of each predictor taken in itself.

Atlantic halibut survival times data

One conservation measure suggested for trawl fishing is a minimum size limit for halibut (32 inches). However, this size limit would only be effective if captured fish below the limit survived until the time of their release. An experiment was conducted to evaluate the survival rates of halibut caught by trawls or longlines, and to assess other factors which might contribute to survival (duration of trawling, maximum depth fished, size of fish, and handling time).

+/- r Code
halibut = read.csv("http://www.karlin.mff.cuni.cz/~pesta/prednasky/NMFM404/Data/halibut.csv", head=T)
dim(halibut)
+/- Output
## [1] 294   8
+/- r Code
head(halibut)
+/- Output
##   No Time Death Towd Deldepth Length Handtime Logcat
## 1  1  209     1   30       13     41        8  6.992
## 2  2  209     1   30       13     44        8  6.992
## 3  3  209     1   30       13     47       10  6.992
## 4  4  209     1   30       13     34       10  6.992
## 5  5   38     1   30       13     40       11  6.992
## 6  6  209     1   30       13     42       11  6.992
+/- r Code
# Time = Survival time, in minutes
# Death = Event indicator
# Length =  length of fish (cm)
# Handtime = handling time (min)
# Logcat = log total catch
+/- r Code
### Univariate analysis
# Towing duration:
table(halibut$Towd)
+/- Output
## 
##  30 100 
## 107 187
+/- r Code
kmtow = survfit(Surv(Time, Death) ~ Towd, data=halibut)
# Length of fish:
summary(halibut$Length)
+/- Output
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   29.00   38.00   43.00   42.51   46.75   57.00
+/- r Code
halibut$Lencat = (halibut$Length > 42) + 0
table(halibut$Lencat)
+/- Output
## 
##   0   1 
## 134 160
+/- r Code
kmlen = survfit(Surv(Time, Death) ~ Lencat, data=halibut)
# Depth range of towing:
table(halibut$Deldepth)
+/- Output
## 
##  1  2  3  4  5  6  8 10 13 15 23 49 58 
##  6 17  8  6 45 46 40 68  8 12 32  3  3
+/- r Code
summary(halibut$Deldepth)
+/- Output
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   5.000   8.000   9.966  10.000  58.000
+/- r Code
halibut$Depthcat = cut(halibut$Deldepth, breaks=c(0,5,9,11,60))
kmdep = survfit(Surv(Time, Death) ~ Depthcat, data=halibut)
# Handling time:
summary(halibut$Handtime)
+/- Output
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    6.00   10.00   12.51   17.00   38.00
+/- r Code
halibut$Handcat = cut(halibut$Handtime, breaks = c(0,6,10,17,40))
kmhan = survfit(Surv(Time, Death) ~ Handcat, data=halibut)
# Log total catch:
summary(halibut$Logcat)
+/- Output
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.904   4.059   4.203   4.886   5.323   8.690
+/- r Code
halibut$Catchcat = cut(halibut$Logcat, breaks=c(0,4.06,4.89,5.33,9))
kmlgc = survfit(Surv(Time, Death) ~ Catchcat, data = halibut)
plot of chunk ph-halibut
+/- r Code
par(mfrow=c(3,2))
plot(kmtow, lty=1:2, col=2:3, xlab="Minutes", ylab="Survival", main="Towing Duration")
plot(kmlen, lty=1:2, col=2:3, xlab="Minutes", ylab="Survival", main="Fish Length")
plot(kmdep, lty=1:4, col=1:4, xlab="Minutes", ylab="Survival", main="Depth of Towing")
plot(kmhan, lty=1:4, col=1:4, xlab="Minutes", ylab="Survival", main="Handling Time")
plot(kmlgc, lty=1:4, col=1:4, xlab="Minutes", ylab="Survival", main="Total Catch")
par(mfrow=c(1,1))

Model selection approach

  1. Fit a univariate model for each covariate, and identify the predictors significant at some level p1, say 0.20 (Hosmer and Lemeshow recommend p1=0.25).
  2. Fit a multivariate model with all significant univariate predictors, and use backward selection to eliminate nonsignificant variables at some level p2, say 0.10.
  3. Starting with final step 2. model, consider each of the non-significant variables from step 1. using forward selection, with significance level p3, say 0.10.
  4. Do final pruning of main-effects model (omit variables that are non-significant, add any that are significant), using stepwise regression with significance level p4. At this stage, you may also consider adding interactions between any of the main effects currently in the model, under the hierarchical principle.

Collett recommends using a likelihood ratio test for all variable inclusion/exclusion decisions.

Each step uses a 'greedy' approach:

  • Backward step: among several candidates for elimi- nation from the model, the one with the smallest effect on the criterion (AIC, LR), is eliminated None is eliminated if AIC increases (AIC) or LR significant (LR)
  • Forward step: among several candidates for addition to the model, the one with the largest effect on the criterion (AIC, LR), is added None is added if AIC increases (AIC) or LR not significant (LR)
+/- r Code
library(MASS)
args(stepAIC)
+/- Output
## function (object, scope, scale = 0, direction = c("both", "backward", 
##     "forward"), trace = 1, keep = NULL, steps = 1000, use.start = FALSE, 
##     k = 2, ...) 
## NULL
+/- r Code
# Model selection for Halibut data
fit = coxph(Surv(Time, Death) ~ Towd + Length + Deldepth + Handtime + Logcat, data=halibut)
# Backward selection using AIC
fitb = stepAIC(fit, direction="backward", k=2)
+/- Output
## Start:  AIC=2520.15
## Surv(Time, Death) ~ Towd + Length + Deldepth + Handtime + Logcat
## 
##            Df    AIC
## - Deldepth  1 2519.7
## <none>        2520.2
## - Towd      1 2531.7
## - Logcat    1 2532.8
## - Length    1 2533.0
## - Handtime  1 2544.9
## 
## Step:  AIC=2519.7
## Surv(Time, Death) ~ Towd + Length + Handtime + Logcat
## 
##            Df    AIC
## <none>        2519.7
## - Logcat    1 2530.8
## - Length    1 2531.2
## - Towd      1 2532.7
## - Handtime  1 2547.7
+/- r Code
summary(fitb)
+/- Output
## Call:
## coxph(formula = Surv(Time, Death) ~ Towd + Length + Handtime + 
##     Logcat, data = halibut)
## 
##   n= 294, number of events= 273 
## 
##               coef exp(coef)  se(coef)      z Pr(>|z|)    
## Towd      0.007744  1.007774  0.002017  3.839 0.000123 ***
## Length   -0.036960  0.963715  0.010028 -3.686 0.000228 ***
## Handtime  0.054947  1.056485  0.009870  5.567 2.59e-08 ***
## Logcat   -0.183373  0.832458  0.050982 -3.597 0.000322 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## Towd        1.0078     0.9923    1.0038    1.0118
## Length      0.9637     1.0377    0.9450    0.9828
## Handtime    1.0565     0.9465    1.0362    1.0771
## Logcat      0.8325     1.2013    0.7533    0.9199
## 
## Concordance= 0.683  (se = 0.02 )
## Rsquare= 0.25   (max possible= 1 )
## Likelihood ratio test= 84.5  on 4 df,   p=0
## Wald test            = 90.71  on 4 df,   p=0
## Score (logrank) test = 94.51  on 4 df,   p=0
+/- r Code
# Forward selection using AIC
fit0 = coxph(Surv(Time, Death) ~ 1, data=halibut)
  # (starting model)
fitf = stepAIC(fit0, scope=formula(fit), direction="forward", k=2)
+/- Warning
## Warning in is.na(fit$coefficients): is.na() applied to non-(list or
## vector) of type 'NULL'
+/- Output
## Start:  AIC=2596.2
## Surv(Time, Death) ~ 1
+/- Warning
## Warning in is.na(fit$coefficients): is.na() applied to non-(list or
## vector) of type 'NULL'
+/- Warning
## Warning in is.na(fit$coefficients): is.na() applied to non-(list or
## vector) of type 'NULL'
+/- Output
##            Df    AIC
## + Handtime  1 2556.8
## + Towd      1 2568.0
## + Length    1 2586.1
## + Deldepth  1 2594.5
## <none>        2596.2
## + Logcat    1 2596.8
## 
## Step:  AIC=2556.85
## Surv(Time, Death) ~ Handtime
## 
##            Df    AIC
## + Logcat    1 2540.3
## + Towd      1 2543.6
## + Length    1 2549.9
## <none>        2556.8
## + Deldepth  1 2558.8
## 
## Step:  AIC=2540.3
## Surv(Time, Death) ~ Handtime + Logcat
## 
##            Df    AIC
## + Towd      1 2531.2
## + Length    1 2532.7
## <none>        2540.3
## + Deldepth  1 2541.4
## 
## Step:  AIC=2531.23
## Surv(Time, Death) ~ Handtime + Logcat + Towd
## 
##            Df    AIC
## + Length    1 2519.7
## <none>        2531.2
## + Deldepth  1 2533.0
## 
## Step:  AIC=2519.7
## Surv(Time, Death) ~ Handtime + Logcat + Towd + Length
## 
##            Df    AIC
## <none>        2519.7
## + Deldepth  1 2520.2
+/- r Code
summary(fitf)
+/- Output
## Call:
## coxph(formula = Surv(Time, Death) ~ Handtime + Logcat + Towd + 
##     Length, data = halibut)
## 
##   n= 294, number of events= 273 
## 
##               coef exp(coef)  se(coef)      z Pr(>|z|)    
## Handtime  0.054947  1.056485  0.009870  5.567 2.59e-08 ***
## Logcat   -0.183373  0.832458  0.050982 -3.597 0.000322 ***
## Towd      0.007744  1.007774  0.002017  3.839 0.000123 ***
## Length   -0.036960  0.963715  0.010028 -3.686 0.000228 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## Handtime    1.0565     0.9465    1.0362    1.0771
## Logcat      0.8325     1.2013    0.7533    0.9199
## Towd        1.0078     0.9923    1.0038    1.0118
## Length      0.9637     1.0377    0.9450    0.9828
## 
## Concordance= 0.683  (se = 0.02 )
## Rsquare= 0.25   (max possible= 1 )
## Likelihood ratio test= 84.5  on 4 df,   p=0
## Wald test            = 90.71  on 4 df,   p=0
## Score (logrank) test = 94.51  on 4 df,   p=0
+/- r Code
# same final model as backward selection!
# Stepwise model selection (backward and forward)
fits = stepAIC(fit, direction="both", k=2)
+/- Output
## Start:  AIC=2520.15
## Surv(Time, Death) ~ Towd + Length + Deldepth + Handtime + Logcat
## 
##            Df    AIC
## - Deldepth  1 2519.7
## <none>        2520.2
## - Towd      1 2531.7
## - Logcat    1 2532.8
## - Length    1 2533.0
## - Handtime  1 2544.9
## 
## Step:  AIC=2519.7
## Surv(Time, Death) ~ Towd + Length + Handtime + Logcat
## 
##            Df    AIC
## <none>        2519.7
## + Deldepth  1 2520.2
## - Logcat    1 2530.8
## - Length    1 2531.2
## - Towd      1 2532.7
## - Handtime  1 2547.7
+/- r Code
summary(fits)
+/- Output
## Call:
## coxph(formula = Surv(Time, Death) ~ Towd + Length + Handtime + 
##     Logcat, data = halibut)
## 
##   n= 294, number of events= 273 
## 
##               coef exp(coef)  se(coef)      z Pr(>|z|)    
## Towd      0.007744  1.007774  0.002017  3.839 0.000123 ***
## Length   -0.036960  0.963715  0.010028 -3.686 0.000228 ***
## Handtime  0.054947  1.056485  0.009870  5.567 2.59e-08 ***
## Logcat   -0.183373  0.832458  0.050982 -3.597 0.000322 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## Towd        1.0078     0.9923    1.0038    1.0118
## Length      0.9637     1.0377    0.9450    0.9828
## Handtime    1.0565     0.9465    1.0362    1.0771
## Logcat      0.8325     1.2013    0.7533    0.9199
## 
## Concordance= 0.683  (se = 0.02 )
## Rsquare= 0.25   (max possible= 1 )
## Likelihood ratio test= 84.5  on 4 df,   p=0
## Wald test            = 90.71  on 4 df,   p=0
## Score (logrank) test = 94.51  on 4 df,   p=0
+/- r Code
# Stepwise selection with alpha=k=3
fitk3 = stepAIC(fit, direction="both", k=3)
+/- Output
## Start:  AIC=2525.15
## Surv(Time, Death) ~ Towd + Length + Deldepth + Handtime + Logcat
## 
##            Df    AIC
## - Deldepth  1 2523.7
## <none>        2525.2
## - Towd      1 2535.7
## - Logcat    1 2536.8
## - Length    1 2537.0
## - Handtime  1 2548.9
## 
## Step:  AIC=2523.7
## Surv(Time, Death) ~ Towd + Length + Handtime + Logcat
## 
##            Df    AIC
## <none>        2523.7
## + Deldepth  1 2525.2
## - Logcat    1 2533.8
## - Length    1 2534.2
## - Towd      1 2535.7
## - Handtime  1 2550.7
+/- r Code
summary(fitk3)
+/- Output
## Call:
## coxph(formula = Surv(Time, Death) ~ Towd + Length + Handtime + 
##     Logcat, data = halibut)
## 
##   n= 294, number of events= 273 
## 
##               coef exp(coef)  se(coef)      z Pr(>|z|)    
## Towd      0.007744  1.007774  0.002017  3.839 0.000123 ***
## Length   -0.036960  0.963715  0.010028 -3.686 0.000228 ***
## Handtime  0.054947  1.056485  0.009870  5.567 2.59e-08 ***
## Logcat   -0.183373  0.832458  0.050982 -3.597 0.000322 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## Towd        1.0078     0.9923    1.0038    1.0118
## Length      0.9637     1.0377    0.9450    0.9828
## Handtime    1.0565     0.9465    1.0362    1.0771
## Logcat      0.8325     1.2013    0.7533    0.9199
## 
## Concordance= 0.683  (se = 0.02 )
## Rsquare= 0.25   (max possible= 1 )
## Likelihood ratio test= 84.5  on 4 df,   p=0
## Wald test            = 90.71  on 4 df,   p=0
## Score (logrank) test = 94.51  on 4 df,   p=0

Notes:

  • When the halibut data was analyzed with the forward, backward and stepwise options, the same final model was reached. However, this will not always (in fact, rarely) be the case.
  • Variables can be forced into the model using the include option in SAS. Any variables that you want to force in- clusion of must be listed first in your model statement. In stepAIC, these are included in the scope part.
  • SAS uses the score test to decide what variables to add and the Wald test for what variables to remove.
  • When might we want to force certain variables into the model?
    1. to examine interactions
    2. to keep main effects in the model
    3. to include scientifically important variables, e.g. treatment, gender, race in a medical study, or other covariates which are known from previous studies or from theory to be important.
  • Model selection is an imperfect art!
    • No method is universally agreed upon
    • Choosing between two models using LR or AIC has an information theory grounding
    • Forward/backward/stepwise selection have no theo- retical bases
    • Leads to potential biases, nominal p-values smaller than true p-values
    • The model selection step is usually not taken into account, and inference is done as if the chosen model were 'true'

Model diagnostics - Assessing overall model fit

How do we know if the model fits well?

  • Always look at univariate plots (Kaplan-Meiers). Construct a Kaplan-Meier survival plot for each of the important predictors.
  • Check proportional hazards assumption.
  • Check residuals.
    1. generalized (Cox-Snell)
    2. martingale
    3. deviance
    4. Schoenfeld
    5. weighted Schoenfeld

Generalized residuals

+/- r Code
# Analysis of halibut data, continued

# 1. Compute cummulative hazard Lambda_0:
Lambda0 = basehaz(fitb, centered=F)
Lambda0[1:3,]
+/- Output
##        hazard time
## 1 0.009323417  1.1
## 2 0.018739005  1.2
## 3 0.057162782  1.3
+/- r Code
# 2. Compute Lambda_0(X_i) for all i:
predLambda0 = function(Lambda0, t)
{
    u = length(Lambda0$time[Lambda0$time0, Lambda0$hazard[u], 0) )
}
n = length(halibut$Time)
L0 = rep(NA, n)   # for storing Lambda_0(X_i)
for(i in 1:n) L0[i] = predLambda0(Lambda0, halibut$Time[i])
# 3. Compute exp(beta Z_i)
ebetaZ =  exp(as.matrix(halibut[,c(4,6,7,8)]) %*% fitb$coef)
# 4. Compute e_i for all i:
e = L0 * ebetaZ
# 5. Run KM for (e, delta)
resfit = survfit(Surv(e,halibut$Death) ~ 1)
plot of chunk ph-generalized-res
+/- r Code
# 6.1 Plot log survival function against time
plot(resfit$time, log(resfit$surv), type="l")
plot of chunk ph-generalized-res
+/- r Code
# 6.2 Plot log-log surv against log time
plot(log(resfit$time), log(-log(resfit$surv)), type="l", xlab = "log time", ylab = "log-log Survival", main = "Halibut: Cox-Snell residuals")

Martingale residuals

plot of chunk unnamed-chunk-16
+/- r Code
# Halibut data, continued

# Calculates martingale, deviance, score or Schoenfeld residuals for a Cox proportional hazards model.
martres = residuals(fitb, type="martingale")  # Martingale residuals
martres[1:3]
+/- Output
##          1          2          3 
## -0.2384461 -0.1084664 -0.1073743
+/- r Code
boxplot(martres ~ halibut$Towd, xlab="Towing duration", ylab="Martingale Residuals", boxwex=0.4)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-16
+/- r Code
plot(halibut$Length, martres, xlab="Halibut Length", ylab="Martingale Residuals", col=2)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-16
+/- r Code
plot(halibut$Handtime, martres, xlab="Handling Time", ylab="Martingale Residuals", col=3)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-16
+/- r Code
plot(halibut$Logcat, martres, xlab="Log Catch Size", ylab="Martingale Residuals", col=4)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-16
+/- r Code
plot(fitb$linear.predictors, martres, xlab="Linear Predictor", ylab="Martingale Residuals", col=6)
abline(h=0, lty=2)

Deviance residuals

plot of chunk unnamed-chunk-17
+/- r Code
# Halibut data, continued
devres = residuals(fitb, type="deviance")  # Deviance residuals
devres[1:3]
+/- Output
##          1          2          3 
## -0.2217595 -0.1047756 -0.1037553
+/- r Code
boxplot(devres ~ halibut$Towd, xlab="Towing duration", ylab="Deviance Residuals", boxwex=0.4)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-17
+/- r Code
plot(halibut$Length, devres, xlab="Halibut Length", ylab="Deviance Residuals", col=2)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-17
+/- r Code
plot(halibut$Handtime, devres, xlab="Handling Time", ylab="Deviance Residuals", col=3)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-17
+/- r Code
plot(halibut$Logcat, devres, xlab="Log Catch Size", ylab="Deviance Residuals", col=4)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-17
+/- r Code
plot(fitb$linear.predictors, devres, xlab="Linear Predictor", ylab="Deviance Residuals", col=6)
abline(h=0, lty=2)

Schoenfeld residuals

plot of chunk unnamed-chunk-18
+/- r Code
# Halibut data, continued
schoeres = residuals(fitb, type="schoenfeld")  # Schoenfeld residuals
schoeres[1:3,]
+/- Output
##         Towd     Length Handtime     Logcat
## 1.1 13.59661 -10.584914 4.165367 -0.7233896
## 1.2 13.73103  -6.689555 3.206545 -0.7305410
## 1.3 14.00863   1.193347 2.279487 -0.7453105
+/- r Code
evtime = sort(halibut$Time[halibut$Death==1]) # event times: 1.1,1.2,1.3 etc
plot(log(evtime), schoeres[,"Towd"], xlab="Log Survival time", ylab="Schoenfeld Residuals, Towing Duration")
abline(h=0, lty=2)
plot of chunk unnamed-chunk-18
+/- r Code
plot(log(evtime), schoeres[,"Length"], xlab="Log Survival time", ylab="Schoenfeld Residuals, Fish Length", col=2)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-18
+/- r Code
plot(log(evtime), schoeres[,"Handtime"], xlab="Log Survival time", ylab="Schoenfeld Res, Handling Time", col=3)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-18
+/- r Code
plot(log(evtime), schoeres[,"Logcat"], xlab="Log Survival time", ylab="Schoenfeld Res, Log Catch Size", col=4)
abline(h=0, lty=2)

Scaled (weighted) Schoenfeld residuals

plot of chunk unnamed-chunk-19
+/- r Code
# Halibut data, continued
scschres = residuals(fitb, type="scaledsch")  # Schoenfeld residuals
scschres[1:3,]
+/- Output
##           [,1]          [,2]      [,3]       [,4]
## 1.1 0.02160255 -0.3136602583 0.1220311 -0.8067589
## 1.2 0.02019807 -0.2117188889 0.1166763 -0.7854618
## 1.3 0.01532368 -0.0002137067 0.1327947 -0.7959172
+/- r Code
evtime = sort(halibut$Time[halibut$Death==1]) # event times: 1.1,1.2,1.3 etc
plot(log(evtime), scschres[,1], xlab="Log Survival time", ylab="Scaled Sch Residuals, Towing Duration")
abline(h=0, lty=2)
plot of chunk unnamed-chunk-19
+/- r Code
plot(log(evtime), scschres[,2], xlab="Log Survival time", ylab="Scaled Sch Residuals, Fish Length", col=2)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-19
+/- r Code
plot(log(evtime), scschres[,3], xlab="Log Survival time", ylab="Scaled Sch Res, Handling Time", col=3)
abline(h=0, lty=2)
plot of chunk unnamed-chunk-19
+/- r Code
plot(log(evtime), scschres[,4], xlab="Log Survival time", ylab="Scaled Sch Res, Log Catch Size", col=4)
abline(h=0, lty=2)

Using Residual plots to explore relationships

If you calculate martingale or deviance residuals without any covariates in the model and then plot them against covariates, you obtain a graphical impression of the relationship between the covariate and the hazard.

plot of chunk unnamed-chunk-20
+/- r Code
fit = coxph(Surv(Time, Death) ~ 1, data=halibut)
martres = residuals(fit, type="martingale")
plot(halibut$Deldepth, martres, xlab="Depth Range", ylab="Martingale Residuals")
lines(lowess(halibut$Deldepth, martres), col=2)
plot of chunk unnamed-chunk-20
+/- r Code
plot(halibut$Length, martres, xlab="Fish Length", ylab= "Martingale Residuals")
lines(lowess(halibut$Length, martres), col=2)
plot of chunk unnamed-chunk-20
+/- r Code
plot(halibut$Handtime, martres, xlab="Handling Time", ylab= "Martingale Residuals")
lines(lowess(halibut$Handtime, martres), col=2)
plot of chunk unnamed-chunk-20
+/- r Code
plot(halibut$Logcat, martres, xlab="Log Catch Size", ylab= "Martingale Residuals")
lines(lowess(halibut$Logcat, martres), col=2)

Assessing the PH assumption

There are several options for checking the assumption of pro- portional hazards:

  1. Graphical
    • Plots of logˆΛ(t)=log[logˆS(t)] vs logt for two subgroups
    • Plots of weighted Schoenfeld residuals vs time
  2. Tests of time-varying coefficient β=β(t)
    • Tests based on Schoenfeld residuals (Grambsch and Thereneau, 1994)
    • Including interaction terms between covariates Zj and time, or logt.
  3. Overall goodness of fit tests I(b), II(b), III available in R via cox.zph.

How do we interpret the above? Kleinbaum (and other texts) suggest a strategy of assuming that PH holds unless there is very strong evidence to counter this assumption:

  • estimated survival curves are fairly separated, then cross
  • estimated log cumulative hazard curves cross, or look very unparallel over time
  • weighted Schoenfeld residuals clearly increase or decrease over time (you could fit a OLS regression line and see if the slope is significant)
  • test for log(time) x covariate interaction term is significant (this relates to time-dependent covariates)

If PH doesn't exactly hold for a particular covariate but we fit the PH model anyway, then what we are getting is sort of an average HR, averaged over the event times.

In most cases, this is not such a bad estimate. Allison claims that too much emphasis is put on testing the PH assumption, and not enough to other important aspects of the model.

Example: Nursing Home - gender and marital status

+/- r Code
# Fit KM curves for men and women separately
fitgen = survfit(Surv(los, discharged) ~ gender, data=nurshome)
names(fitgen)
+/- Output
##  [1] "n"         "time"      "n.risk"    "n.event"   "n.censor" 
##  [6] "surv"      "type"      "strata"    "std.err"   "upper"    
## [11] "lower"     "conf.type" "conf.int"  "call"
+/- r Code
fitgen$time
+/- Output
##   [1]    0    1    2    3    4    5    6    7    8    9   10   11   12   13
##  [15]   14   15   16   17   18   19   20   21   22   23   24   25   26   27
##  [29]   28   29   30   31   32   33   34   35   36   37   38   39   40   41
##  [43]   42   43   44   45   46   47   48   49   50   51   52   53   54   55
##  [57]   56   57   59   60   61   62   63   64   66   67   68   69   70   71
##  [71]   72   74   75   76   77   78   80   81   82   83   84   85   86   88
##  [85]   89   90   91   92   93   94   95   96   97   98   99  100  101  102
##  [99]  104  105  106  107  108  109  110  111  113  114  115  116  118  119
## [113]  120  121  122  123  125  126  127  128  129  130  131  132  133  134
## [127]  135  137  139  140  141  142  143  144  146  147  148  149  150  152
## [141]  153  155  156  159  160  161  162  163  164  165  166  168  169  170
## [155]  172  174  175  176  178  180  182  184  185  187  189  190  191  192
## [169]  193  194  195  196  197  198  199  201  202  203  204  205  207  208
## [183]  209  211  214  215  217  218  221  222  223  224  225  226  227  230
## [197]  231  232  233  234  237  238  240  241  243  244  250  252  253  256
## [211]  258  259  261  262  263  265  268  269  270  271  273  274  276  277
## [225]  279  281  282  283  285  288  290  291  293  295  296  297  298  301
## [239]  302  305  306  307  309  310  311  312  314  315  317  318  322  326
## [253]  328  330  332  337  340  342  344  350  351  352  355  358  360  362
## [267]  363  364  365  366  367  369  370  372  373  374  375  376  377  378
## [281]  379  380  381  382  383  384  386  387  388  389  390  393  394  395
## [295]  396  397  401  402  403  405  407  408  409  410  411  413  415  416
## [309]  417  418  420  421  422  424  425  428  429  430  433  434  435  437
## [323]  439  442  444  446  449  450  451  452  453  456  457  461  462  463
## [337]  465  466  467  470  471  472  474  475  476  477  478  480  481  483
## [351]  487  488  489  492  494  495  497  498  500  505  506  507  508  512
## [365]  514  519  521  522  525  526  527  530  533  537  540  545  547  549
## [379]  553  554  556  561  562  563  569  570  572  576  579  584  585  589
## [393]  592  597  598  599  605  611  613  617  618  619  620  624  625  626
## [407]  627  628  634  635  639  642  645  648  653  654  660  661  664  665
## [421]  669  670  673  676  679  680  681  683  684  686  687  690  691  695
## [435]  696  697  698  702  707  708  709  710  711  717  719  721  722  724
## [449]  726  732  733  736  741  747  749  754  756  757  758  759  770  774
## [463]  775  780  783  788  790  792  803  810  811  812  817  818  823  832
## [477]  833  838  851  852  853  858  864  865  866  870  875  879  881  895
## [491]  896  900  903  906  909  911  913  921  924  929  936  941  942  945
## [505]  962  963  964  965  970  982  986  992  993  994  999 1005 1006 1008
## [519] 1010 1012 1022 1026 1028 1033 1034 1041 1043 1047 1053 1054 1055 1057
## [533] 1062 1074 1084 1088 1092    0    1    2    3    4    5    6    7    8
## [547]    9   10   11   12   13   14   15   16   17   18   19   20   21   22
## [561]   23   24   25   26   27   29   30   31   32   33   34   35   37   38
## [575]   39   40   41   42   43   44   45   46   47   48   49   50   51   54
## [589]   55   57   58   60   61   62   63   65   66   68   69   71   72   73
## [603]   76   77   80   81   82   83   84   85   86   88   89   90   91   92
## [617]   95   96   97  100  103  105  106  108  109  112  116  118  119  120
## [631]  121  122  123  126  131  133  138  143  144  146  147  148  151  152
## [645]  153  156  159  160  162  164  165  167  169  170  171  174  176  178
## [659]  182  185  189  190  193  195  202  205  206  207  212  216  227  234
## [673]  237  241  250  260  262  267  269  270  273  274  275  276  277  278
## [687]  279  290  294  297  302  303  306  310  326  335  365  370  374  375
## [701]  381  389  393  396  397  399  400  404  408  412  416  448  449  456
## [715]  463  465  467  476  487  492  495  501  541  548  561  565  568  574
## [729]  578  582  584  586  597  599  602  605  607  620  624  638  646  649
## [743]  651  663  669  674  709  725  738  740  743  752  756  773  797  816
## [757]  864  881  899  934  938  948  971  973  985 1060 1074 1091
+/- r Code
fitgen$surv
+/- Output
##   [1] 0.99575552 0.97962649 0.97113752 0.96859083 0.95925297 0.94906621
##   [7] 0.93972835 0.92954160 0.92359932 0.91341256 0.90747029 0.89983022
##  [13] 0.89303905 0.88285229 0.87351443 0.86332767 0.85483871 0.85059423
##  [19] 0.84295416 0.83786078 0.83276740 0.82258065 0.81578947 0.80730051
##  [25] 0.80050934 0.79286927 0.79032258 0.78438031 0.77843803 0.77079796
##  [31] 0.76740238 0.76570458 0.75976231 0.75551783 0.75127334 0.74702886
##  [37] 0.74617997 0.74023769 0.73769100 0.73344652 0.73005093 0.72495756
##  [43] 0.72241087 0.71816638 0.71222411 0.70797963 0.70373514 0.69779287
##  [49] 0.69185059 0.69100170 0.69015280 0.68675722 0.68505942 0.68251273
##  [55] 0.67741935 0.67487267 0.66977929 0.66638370 0.66298812 0.66044143
##  [61] 0.65619694 0.65025467 0.64601019 0.64346350 0.64091681 0.64006791
##  [67] 0.63837012 0.63497453 0.63327674 0.62818336 0.62563667 0.62139219
##  [73] 0.61884550 0.61460102 0.61290323 0.61035654 0.60950764 0.60696095
##  [79] 0.60441426 0.60271647 0.59847199 0.59592530 0.59422750 0.59337861
##  [85] 0.58998302 0.58913413 0.58488964 0.58064516 0.57979626 0.57640068
##  [91] 0.57470289 0.56876061 0.56791171 0.56536503 0.56451613 0.56112054
##  [97] 0.56027165 0.55772496 0.55687606 0.55432937 0.55348048 0.55263158
## [103] 0.54838710 0.54753820 0.54584041 0.54414261 0.54244482 0.54074703
## [109] 0.53735144 0.53650255 0.53480475 0.53395586 0.53225806 0.53140917
## [115] 0.52886248 0.52546689 0.52376910 0.52292020 0.52122241 0.51867572
## [121] 0.51697793 0.51443124 0.51358234 0.51188455 0.51018676 0.50764007
## [127] 0.50679117 0.50509338 0.50254669 0.50169779 0.50084890 0.50000000
## [133] 0.49830221 0.49405772 0.49235993 0.49151104 0.48811545 0.48556876
## [139] 0.48471986 0.48387097 0.48217317 0.47962649 0.47707980 0.47538200
## [145] 0.47453311 0.47198642 0.46943973 0.46859083 0.46519525 0.46434635
## [151] 0.46349745 0.45925297 0.45840407 0.45755518 0.45415959 0.45246180
## [157] 0.45161290 0.45076401 0.44991511 0.44821732 0.44736842 0.44651952
## [163] 0.44567063 0.44397284 0.44312394 0.44142615 0.44057725 0.43972835
## [169] 0.43887946 0.43803056 0.43633277 0.43548387 0.43463497 0.43208829
## [175] 0.43123939 0.42954160 0.42784380 0.42614601 0.42359932 0.42275042
## [181] 0.42020374 0.41850594 0.41595925 0.41341256 0.41171477 0.41001698
## [187] 0.40916808 0.40831919 0.40747029 0.40577250 0.40322581 0.40067912
## [193] 0.39983022 0.39898132 0.39813243 0.39728353 0.39643463 0.39473684
## [199] 0.39388795 0.39219015 0.39134126 0.39049236 0.38879457 0.38794567
## [205] 0.38709677 0.38624788 0.38539898 0.38370119 0.38285229 0.38200340
## [211] 0.38115450 0.38030560 0.37775891 0.37691002 0.37606112 0.37521222
## [217] 0.37436333 0.37351443 0.37266553 0.37181664 0.37096774 0.36842105
## [223] 0.36757216 0.36502547 0.36417657 0.36332767 0.36247878 0.36078098
## [229] 0.35993209 0.35908319 0.35823430 0.35568761 0.35314092 0.35229202
## [235] 0.35059423 0.34974533 0.34804754 0.34719864 0.34634975 0.34550085
## [241] 0.34465195 0.34380306 0.34295416 0.34210526 0.34125637 0.34040747
## [247] 0.33955857 0.33786078 0.33531409 0.33361630 0.33276740 0.33191851
## [253] 0.33022071 0.32937182 0.32852292 0.32767402 0.32682513 0.32597623
## [259] 0.32427844 0.32342954 0.32258065 0.32173175 0.32088285 0.31918506
## [265] 0.31833616 0.31748727 0.31663837 0.31578947 0.31578947 0.31493368
## [271] 0.31493368 0.31407084 0.31320563 0.31234043 0.31234043 0.31060520
## [277] 0.30973759 0.30886754 0.30712253 0.30712253 0.30537254 0.30449755
## [283] 0.30449755 0.30449755 0.30449755 0.30360980 0.30272205 0.30183430
## [289] 0.30183430 0.30183430 0.30093865 0.30093865 0.29913662 0.29913662
## [295] 0.29823289 0.29732640 0.29732640 0.29641436 0.29549382 0.29457040
## [301] 0.29457040 0.29457040 0.29364116 0.29270896 0.29177379 0.29083862
## [307] 0.29083862 0.29083862 0.29083862 0.28989434 0.28895006 0.28895006
## [313] 0.28800268 0.28800268 0.28800268 0.28703946 0.28703946 0.28703946
## [319] 0.28606973 0.28510000 0.28413027 0.28413027 0.28023808 0.28023808
## [325] 0.27827837 0.27827837 0.27728806 0.27629774 0.27530387 0.27430639
## [331] 0.27430639 0.27430639 0.27430639 0.27329791 0.27329791 0.27329791
## [337] 0.27227814 0.27125454 0.27125454 0.26919177 0.26919177 0.26814839
## [343] 0.26814839 0.26709683 0.26604526 0.26604526 0.26393379 0.26393379
## [349] 0.26179668 0.26072812 0.26072812 0.26072812 0.26072812 0.25963721
## [355] 0.25854629 0.25854629 0.25631745 0.25631745 0.25631745 0.25519325
## [361] 0.25519325 0.25519325 0.25519325 0.25519325 0.25519325 0.25402798
## [367] 0.25402798 0.25402798 0.25285193 0.25167587 0.25049982 0.24932376
## [373] 0.24932376 0.24932376 0.24932376 0.24812509 0.24692642 0.24452907
## [379] 0.24333040 0.24333040 0.24333040 0.24333040 0.24333040 0.24333040
## [385] 0.24209522 0.24086004 0.24086004 0.24086004 0.24086004 0.23959899
## [391] 0.23833127 0.23833127 0.23833127 0.23833127 0.23833127 0.23702176
## [397] 0.23702176 0.23702176 0.23702176 0.23702176 0.23702176 0.23702176
## [403] 0.23702176 0.23702176 0.23563567 0.23563567 0.23563567 0.23423308
## [409] 0.23423308 0.23423308 0.23423308 0.23423308 0.23423308 0.23423308
## [415] 0.23275059 0.23125860 0.23125860 0.23125860 0.22974710 0.22823561
## [421] 0.22823561 0.22823561 0.22823561 0.22823561 0.22668299 0.22668299
## [427] 0.22668299 0.22668299 0.22668299 0.22507531 0.22507531 0.22507531
## [433] 0.22507531 0.22507531 0.22507531 0.22338301 0.22338301 0.22338301
## [439] 0.22338301 0.22158153 0.21978006 0.21978006 0.21978006 0.21978006
## [445] 0.21978006 0.21978006 0.21978006 0.21978006 0.21778206 0.21778206
## [451] 0.21778206 0.21570794 0.21363383 0.21363383 0.21363383 0.21363383
## [457] 0.21363383 0.21136113 0.20908843 0.20908843 0.20908843 0.20908843
## [463] 0.20908843 0.20908843 0.20908843 0.20908843 0.20908843 0.20650709
## [469] 0.20650709 0.20650709 0.20382518 0.20382518 0.20382518 0.20382518
## [475] 0.20382518 0.20382518 0.20382518 0.20382518 0.20382518 0.20382518
## [481] 0.20382518 0.20382518 0.20382518 0.20382518 0.20382518 0.20042809
## [487] 0.20042809 0.20042809 0.20042809 0.20042809 0.20042809 0.19657371
## [493] 0.19657371 0.19657371 0.19657371 0.19239129 0.19239129 0.19239129
## [499] 0.19239129 0.19239129 0.19239129 0.18769882 0.18288603 0.18288603
## [505] 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603
## [511] 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603
## [517] 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603
## [523] 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603
## [529] 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603 0.18288603
## [535] 0.18288603 0.18288603 0.18288603 0.98817967 0.98108747 0.96453901
## [541] 0.95035461 0.93617021 0.91252955 0.90307329 0.88888889 0.87234043
## [547] 0.86288416 0.85106383 0.83451537 0.82269504 0.81087470 0.79432624
## [553] 0.78486998 0.78014184 0.76122931 0.75177305 0.73995272 0.72813239
## [559] 0.72576832 0.71631206 0.70921986 0.70212766 0.68321513 0.67375887
## [565] 0.65957447 0.65011820 0.64302600 0.63593381 0.62884161 0.62174941
## [571] 0.61938534 0.61465721 0.60992908 0.60283688 0.59338061 0.58628842
## [577] 0.58392435 0.58156028 0.57919622 0.57210402 0.56973995 0.56501182
## [583] 0.56028369 0.55791962 0.55082742 0.54846336 0.54137116 0.53900709
## [589] 0.53664303 0.53427896 0.53191489 0.52955083 0.52009456 0.51773050
## [595] 0.51536643 0.50591017 0.50118203 0.49881797 0.49408983 0.48699764
## [601] 0.48226950 0.47517730 0.47281324 0.47044917 0.46335697 0.46099291
## [607] 0.45862884 0.45626478 0.45153664 0.44917258 0.44444444 0.44208038
## [613] 0.43498818 0.43262411 0.42553191 0.42316785 0.41843972 0.41607565
## [619] 0.41371158 0.40661939 0.39952719 0.39716312 0.39479905 0.39243499
## [625] 0.39007092 0.38770686 0.38534279 0.38297872 0.37825059 0.37352246
## [631] 0.37115839 0.36643026 0.36406619 0.35697400 0.35224586 0.34988180
## [637] 0.34751773 0.34278960 0.34042553 0.33806147 0.33333333 0.33096927
## [643] 0.32624113 0.32387707 0.32151300 0.31914894 0.31678487 0.31442080
## [649] 0.31205674 0.30969267 0.30732861 0.30260047 0.29787234 0.29550827
## [655] 0.29314421 0.29078014 0.28841608 0.28605201 0.28368794 0.28132388
## [661] 0.27659574 0.27423168 0.27186761 0.26950355 0.26713948 0.26477541
## [667] 0.26241135 0.26004728 0.25768322 0.25531915 0.25295508 0.25059102
## [673] 0.24822695 0.24586288 0.24349882 0.23877069 0.23640662 0.23404255
## [679] 0.23167849 0.22931442 0.22695035 0.22458629 0.22222222 0.21985816
## [685] 0.21749409 0.21513002 0.21276596 0.21040189 0.20803783 0.20567376
## [691] 0.20330969 0.20094563 0.19621749 0.19385343 0.19148936 0.18912530
## [697] 0.18676123 0.18439716 0.18439716 0.18200239 0.17960763 0.17721286
## [703] 0.17481809 0.17481809 0.16996203 0.16753400 0.16510597 0.16510597
## [709] 0.16510597 0.16260437 0.16010276 0.15760116 0.15760116 0.15760116
## [715] 0.15760116 0.15760116 0.15760116 0.15483622 0.15483622 0.15202102
## [721] 0.15202102 0.15202102 0.15202102 0.15202102 0.15202102 0.15202102
## [727] 0.15202102 0.14871622 0.14541141 0.14210661 0.14210661 0.13872312
## [733] 0.13525504 0.13178696 0.12831888 0.12831888 0.12831888 0.12831888
## [739] 0.12831888 0.12831888 0.12831888 0.12831888 0.12404159 0.11976429
## [745] 0.11976429 0.11532857 0.11089286 0.11089286 0.11089286 0.11089286
## [751] 0.10534822 0.10534822 0.10534822 0.10534822 0.10534822 0.10534822
## [757] 0.10534822 0.08779018 0.08779018 0.08779018 0.08779018 0.08779018
## [763] 0.08779018 0.08779018 0.08779018 0.08779018 0.08779018 0.08779018
+/- r Code
(n = fitgen$strata)  # nr fail times for each level of "gender"
+/- Output
## gender=0 gender=1 
##      537      231
plot of chunk ph-ph-assump
+/- r Code
# Compute and plot log cumulative hazard for each stratum
Lambda = -log(fitgen$surv)
plot( log(fitgen$time), log(Lambda), type="n", xlab = "Log Length of Stay", ylab = "Log Cum Hazard")
lines( log(fitgen$time)[1:n[1]], log(Lambda)[1:n[1]], col=2 )
lines( log(fitgen$time)[(n[1]+1):sum(n)], log(Lambda)[(n[1]+1):sum(n)], col=3, lty=2)
legend(4,-3, legend=c("Females", "Males"), lty=1:2, col=2:3)
+/- r Code
# Fit KM curves for single and married separately
fitmar = survfit(Surv(los, discharged) ~ married, data=nurshome)
names(fitmar)
+/- Output
##  [1] "n"         "time"      "n.risk"    "n.event"   "n.censor" 
##  [6] "surv"      "type"      "strata"    "std.err"   "upper"    
## [11] "lower"     "conf.type" "conf.int"  "call"
+/- r Code
(n = fitmar$strata)  # nr fail times for each level of marital status
+/- Output
## married=0 married=1 
##       563       174
plot of chunk ph-ph-assump2
+/- r Code
# Compute and plot log cumulative hazard for each stratum
Lambda = -log(fitmar$surv)
plot( log(fitmar$time), log(Lambda), type="n", xlab = "Log Length of Stay", ylab = "Log Cum Hazard")
lines( log(fitmar$time)[1:n[1]], log(Lambda)[1:n[1]], col=2 )
lines( log(fitmar$time)[(n[1]+1):sum(n)], log(Lambda)[(n[1]+1):sum(n)], col=3, lty=2)
legend(4,-3, legend=c("Not Married", "Married"), lty=1:2, col=2:3)
+/- r Code
# Check PH for model including both gender and marital status
fitboth = survfit(Surv(los, discharged) ~ gender + married, data=nurshome)
fitboth
+/- Output
## Call: survfit(formula = Surv(los, discharged) ~ gender + married, data = nurshome)
## 
##                     records n.max n.start events median 0.95LCL 0.95UCL
## gender=0, married=0    1065  1065    1065    811  148.0     128     168
## gender=0, married=1     113   113     113     96   96.0      63     149
## gender=1, married=0     261   261     261    228   66.0      49     100
## gender=1, married=1     162   162     162    144   68.5      39      89

Testing PH using cox.zph

+/- r Code
# Nursing Home example, continued
fitnh = coxph(Surv(los, discharged) ~ married + gender + age, data=nurshome)
(testnh = cox.zph(fitnh))
+/- Output
##              rho  chisq     p
## married  0.00476 0.0299 0.863
## gender  -0.03721 1.8414 0.175
## age      0.00281 0.0104 0.919
## GLOBAL        NA 2.0614 0.560
plot of chunk ph-nurs-assess
+/- r Code
par(mfrow=c(2,2))
plot(testnh)  # Schoenfeld residuals - see plot
+/- r Code
# Halibut data, revisited
fithal = coxph(Surv(Time, Death) ~ Towd + Length + Deldepth + Handtime + Logcat, data=halibut)
(testhal = cox.zph(fithal))
+/- Output
##               rho   chisq        p
## Towd     -0.12120  3.6198 0.057096
## Length    0.00719  0.0123 0.911639
## Deldepth -0.06451  0.7290 0.393202
## Handtime -0.03716  0.3819 0.536596
## Logcat   -0.16300  6.4651 0.011002
## GLOBAL         NA 20.8640 0.000859
plot of chunk ph-halibut-assess
+/- r Code
par(mfrow=c(3,2))
plot(testhal) # see plot

What if proportional hazards fails?

  • do a stratified analysis
  • include a time-varying covariate to allow changing hazard ratios over time
  • include interactions with time

The second two options relate to time-dependent covariates.

We will focus on the first alternative, and then the second two options will be briefly described.

Example: PH assumption for gender (nursing home data):

  • include married and health as covariates in a Cox PH model, but stratify by gender.
  • calculate the baseline survival function for each level of the variable gender (i.e., males and females)
  • plot the log-cumulative hazards for males and females and evaluate whether the lines (curves) are parallel

In the above example, we make the PH assumption for married and health, but not for gender.

This is like getting a KM survival estimate for each gen- der without assuming PH, but is more flexible since we can control for other covariates.

We would repeat the stratification for each variable for which we wanted to check the PH assumption.

+/- r Code
# Nursing Home analysis, continued
fitnh = coxph(Surv(los, discharged) ~ married + health + gender, data=nurshome)
fitnh
+/- Output
## Call:
## coxph(formula = Surv(los, discharged) ~ married + health + gender, 
##     data = nurshome)
## 
## 
##          coef exp(coef) se(coef)    z       p
## married 0.161      1.17   0.0767 2.10 3.6e-02
## health  0.169      1.18   0.0311 5.42 6.0e-08
## gender  0.355      1.43   0.0660 5.37 7.7e-08
## 
## Likelihood ratio test=74.1  on 3 df, p=5.55e-16  n= 1601, number of events= 1279
+/- r Code
fitnhs = coxph(Surv(los, discharged) ~ married + health + strata(gender), data = nurshome)
fitnhs
+/- Output
## Call:
## coxph(formula = Surv(los, discharged) ~ married + health + strata(gender), 
##     data = nurshome)
## 
## 
##          coef exp(coef) se(coef)    z       p
## married 0.163      1.18   0.0768 2.12 3.4e-02
## health  0.169      1.18   0.0312 5.43 5.6e-08
## 
## Likelihood ratio test=34.5  on 2 df, p=3.3e-08  n= 1601, number of events= 1279
plot of chunk ph-stratify
+/- r Code
# Plot log cum. base. haz. for strata
bh = basehaz(fitnhs)  # cum. base. hazards
plot(log(bh$time), log(bh$hazard), type="n", xlab = "log Time", ylab="log Cumulative Hazard")
lines(log(bh$time)[bh$strata=="gender=0"], log(bh$hazard)[bh$strata=="gender=0"], lty=2, col=2)
lines(log(bh$time)[bh$strata=="gender=1"], log(bh$hazard)[bh$strata=="gender=1"], lty=3, col=3)
legend(4,-3, legend=c("Females", "Males"), lty=2:3, col=2:3)

References

  • Collett: Modelling Survival Data in Medical Research
  • Cox & Oakes: Analysis of Survival Data
  • Kalbfleisch & Prentice: The Statistical Analysis of Failure Time Data
  • Lee: Statistical Methods for Survival Data Analysis
  • Fleming & Harrington: Counting Processes and Survival Analysis
  • Hosmer & Lemeshow: Applied Survival Analysis
  • Kleinbaum: Survival Analysis: A self-learning text
  • Klein & Moeschberger: Survival Analysis: Techniques for censored and truncated data
  • Cantor: Extending SAS Survival Analysis Techniques for Medical Research
  • Allison: Survival Analysis Using the SAS System
  • Jennison & Turnbull: Group Sequential Methods with Applications to Clinical Trials