survminer customizable plots of survival functions (requires ggplot2)
broom: storing output tables as data frames.
library(survival)library(survminer) # for customizable graphs of survival functionlibrary(broom) # for tidy output library(ggplot2) # for graphing (actually loaded by survminer)
Outline
Review of survival and hazard functions
Data setup
Kaplan-Meier estimator
Cox proportional hazards model
Assessing the proportional hazards assumption
Time-varying covariates
A very quick review of survival analysis (POLL)
What is survival analysis?
Survival analysis models survival time, the length of time to an event.
Other names for survival time are failure time or time to event.
Example events:
death from disease
divorce
malfunctioning of a machine
first job
The survival time clock starts when a subject becomes at risk for the event, at which point they enter the risk set.
Survival function
Let \(T\) represent survival time.
The survival function, \(S(t)\) expresses the probability that a subject survives beyond time \(T\).
\[S(t)=Pr(T>t)\]
For the survival curve above,
\(S(100)=.577\), the probability that a subject survives beyond 100 days is 0.577.
\(S(200)=.122\), the probability that a subject survives beyond 200 days is 0.122.
Typically, we also assume:
\(S(0)=1\), all subjects survive the very first moment
\(S(\infty)=0\), all subjects fail after infinite time
Median survival time
Median survival time: time \(t\) at which 50% of population expected to be still surviving:
Hazard function
Some methods, e.g. Cox model, focus on the hazard function, \(h(t)\), which is inversely related to the \(S(t)\).
\(h(t)\), is the rate of events, or hazard at time \(t\), given that the subject has survived until time \(t\).
As the hazard increases, more events are expected per unit time, survival will decrease.
Below we see three examples of hazard functions.
In the green hazard function above, \(h(200)=\frac{\text{.0204 events}}{day}\).
Cumulative hazard
The cumulative hazard function, \(H(t)\), expresses how much hazard has accumulated over time up to time \(t\).
\[H(t)=\int_0^t h(u)du\]
The probability that a subject fails increases as the hazard accumulates.
Because the hazard function \(h(t)\) is never negative, the cumulative hazard \(H(t)\) never decreases with time.
Relationship between the hazard and survival functions
The survival function is inversely related to the cumulative hazard function:
\[S(t)=exp(-H(t))\]
Therefore, by modeling either the survival function or the hazard function, we can infer the other.
Censored survival times
Right Censoring
Often the exact time the event occurs is unknown or censored.
Right censoring means that a subject’s actual survival time is greater than their observed time
study ends before event occurs
subject is lost to follow-up
Left and Interval Censoring
Left-censoring means that a subject’s actual survival time is less than their observed time.
positive tests for the infection (the event) may be considerably delayed
Interval censoring means that a subject’s survival time is only known to lie between 2 time points.
We will only discuss methods that handle right-censoring in this workshop.
Many standard methods, such as linear regression, are not equipped to deal with censored outcomes.
*Image adapted from Kleinbaum and Klein, Survival Analysis: A Self-Learning Text, Third Edition, Springer, 2012.
Assumption of non-informative censoring
Most survival analysis methods assume non-informative censoring: a subject’s censoring time should not be related to the unobserved survival time
Below, we see 2 graphs diagrams where subjects have the same true survival times (red squares), but censoring times (black circles) are informative on the left and non-informative on the right.
Bias from informative censoring
Failing to account for informative censoring may result in biased estimates of survival.
Below are estimated survival functions from the data diagrammed on the previous slide with either noninformative (purple) or informative censoring (green).
Examples of possible informative censoring and resulting bias if not addressed:
oldest, most vulnerable, subjects drop out of study of time to death after surgery, biasing survival estimates upward
Travel-loving subjects, who delay marriage to travel, drop out of study of time to first marriage, biasing survival estimates downward
Data set up for survival analysis
Data for survival analysis
A data record for subject in a survival analysis typically includes:
survival time variable measuring time to event (or censoring time, time of last observation)
a status variable coding whether the subject experienced the event or not (censored)
covariates
The aml dataset
We’ll start with the aml dataset in the survival package.
These data come from a study looking at time to death for patients with acute myelogenous leukemia, comparing “maintained” chemotherapy treatment to “nonmaintained”.
Variables:
time survival or censoring time
status 0=censored, 1=death
x chemotherapy “maintained” or “nonmaintained”
time
status
x
9
1
Maintained
13
1
Maintained
13
0
Maintained
18
1
Maintained
23
1
Maintained
28
0
Maintained
31
1
Maintained
34
1
Maintained
45
0
Maintained
48
1
Maintained
161
0
Maintained
5
1
Nonmaintained
5
1
Nonmaintained
8
1
Nonmaintained
8
1
Nonmaintained
12
1
Nonmaintained
16
0
Nonmaintained
23
1
Nonmaintained
27
1
Nonmaintained
30
1
Nonmaintained
33
1
Nonmaintained
43
1
Nonmaintained
45
1
Nonmaintained
The Surv() function for survival outcomes
Use Surv() to specify the survival outcome variables.
Surv(time, status)
time survival/censoring time variable
status status variable. To code for censored/event use:
0/1
1/2
FALSE/TRUE
Censoring is assumed to be right-censored unless otherwise specified with the type argument.
Surv() specification for start-stop format
Some survival analyses require time to be recorded in 2 variables that mark the beginning and end of time intervals. We need this format to model:
time-varying covariates
interval censoring
recurrent events data
In this format, some or all subjects may have multiple rows of data.
This format is sometimes called start-stop format.
The jasa1 data set has this setup, where start and stop are the time variables, and event is the status variable:
plot() produces separate KM survival curves when supplied a stratified survfit object.
However, confidence intervals are no longer produced by default so we request them with conf.int=TRUE.\(\dagger\)
We also use col to specify two colors.
# stratified KM curves with 95% CI, 2 colorsplot(KM.x, ylab="survival probability", xlab="months",conf.int=TRUE, col=c("red", "blue"))
\(\dagger\)plot() is a generic base R function and calls plot.survfit() when supplied a survfit object. The option conf.int= is an option of plot.survfit().
Customized survival plots with survminer
plot.survfit() has limited options to customize a plot of survival functions.
ggsurvplot() from the survminer package leverages the graphical power of the ggplot2 package to allow extensive customization of survival function plots.
ggsurvplot() by default adds + symbols to denote censored observations.
ggsurvplot(KM.x, conf.int=T)
ggsurvplot() makes adding a risk table (number at risk at selected timepoints) very easy.
ggsurvplot(KM.x, conf.int=T,risk.table=T)
You can pass most arguments to various functions in ggplot2 through ggsurvplot().
ggsurvplot(KM.x, conf.int=T, risk.table=T,palette="Accent", # argument to scale_color_brewer()size=2, # argument to geom_line()ggtheme =theme_minimal()) # changing ggtheme
You can also use traditional ggplot2 syntax by extracting the ggplot object as ggsurvplot()$plot.
g <-ggsurvplot(KM.x, conf.int=T,risk.table=T)$plot # this is the ggplot objectg +scale_fill_grey() +scale_color_grey()
See our ggplot2 seminar to learn how to use the ggplot2 package.
Testing equivalence of survival functions with survdiff()
A common hypothesis of interest:
\(H_0\): survival curves across 2 or more groups are equivalent \(H_A\): survival curves across 2 or more groups are not equivalent
A popular method to evaluate this hypothesis is the log-rank statistic, which is \(\chi^2\) distributed with \(g-1\) degrees of freedom under the null (\(g\) is the number of groups).
The function survdiff() performs the log-rank test.
Call:
survdiff(formula = Surv(time, status) ~ x, data = aml)
N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained 11 7 10.69 1.27 3.4
x=Nonmaintained 12 11 7.31 1.86 3.4
Chisq= 3.4 on 1 degrees of freedom, p= 0.07
We see weak evidence that the survival curves may be different.
Weighted log-rank tests
For some problems, early differences in survival curves might be more important than later differences, such as for death after surgery.
Early differences can be emphasized by applying weights that decrease as time advances.
One popular weighting method is to weight the survival times by the survival probabilities \(\hat{S}(t)\) themselves, which decrease with time.
survdiff() allows weights in calculation of the \(\chi^2\) statistic with the rho= argument. The weights are defined as \(\hat{S}(t)^\rho\), where \(0 \le \rho \le 1\).
rho=1 weights are the survival estimates, \(\hat{S}(t)^1=\hat{S}(t)\), equivalent to the Gehan-Wilcoxon test
rho=0 equal weights, \(\hat{S}(t)^0=1\), the traditional log-rank test and the default
0<rho<1 values closer to 1 weight earlier time points more
# rho=1 specifies Peto & Peto modification of Gehan-Wilcoxon,# more weight put on earlier time pointssurvdiff(Surv(time, status) ~ x, data=aml, rho=1)
Call:
survdiff(formula = Surv(time, status) ~ x, data = aml, rho = 1)
N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained 11 3.85 6.14 0.859 2.78
x=Nonmaintained 12 7.18 4.88 1.081 2.78
Chisq= 2.8 on 1 degrees of freedom, p= 0.1
Exercise 1
The veteran data set describes survival times for veterans with lung cancer.
Variables:
time: survival time
status: status, 0=censored, 1=dead
trt: treatment, 1=standard, 2=test
Create a graph and table of the Kaplan-Meier estimated survival function for the entire data set. What is the median survival time?
Create a graph of KM estimated survival functions stratified by treatment, adding 95% confidence intervals and coloring the 2 functions. Does survival appear different for the 2 treatment groups?
Use the log-rank test to provide more evidence for your assessment of survival of the 2 groups.
The Cox proportional hazards model
Why the Cox model is so widely used
The Cox model:
is by far the most popular survival analysis method
estimates changes to the hazard function, \(h(t)\), rather than the survival function \(S(t)\)
accommodates multiple predictors (covariates)
does not assume any distribution for survival times
naturally accommodates right-censoring and time-varying covariates
can be extended in many ways:
time-varying coefficients
random effect frailties for recurrent events or clustering
competing risks modeling
The Cox proportional hazards model
For simplicity, we begin with a Cox model with a single time-constant predictor, \(X_1\):
\[h(t|X_1=x_1) = h_0(t)exp(b_1x_1)\]
\(h(t|X_1=x_1)\): the hazard at time \(t\) for a subject with predictor \(X_1\) equal to the value \(x_1\)
\(h_0(t)\): the baseline hazard at time \(t\), the hazard for a subject with all predictors equal to zero
\(exp(b_1x_1)\): the hazard ratio comparing the hazard for a subject with \(X_1=x_1\) to a subject with \(X_1=0\)
Hazard ratio
For example, imagine \(X_1=1\) for treatment and \(X_1=0\) for control.
\(HR=.25\) means treatment has \(\frac{1}{4}\) (25%) the hazard of control, or a 75% decrease. Treatment would lead to better survival with a lower hazard.
\(HR=2\) here means that treatment has twice the hazard of control, or a 100% increase, and thus worse survival.
In general, \(exp(b_1)\) expresses the hazard ratio for a 1-unit increase in the predictor.
\(b_1\) itself is the log-hazard ratio.
Proportional hazards
The standard Cox model assumes proportional hazards, which means that covariate effects (expressed as hazard ratios) are constant over time.
Notice that the expression for the hazard ratio for the treatment effect does not contain time, so is constant:
\[HR=exp(b_1)\]
Note: Under proportional hazards, a subject’s hazard can change over time, but the hazard ratio comparing 2 subjects with different covariate values cannot.
Visually, this is represented by “parallel” survival curves that should not cross:
Proportional hazards is easier to evaluate if we plot \(-log(-log(S(t)))\), sometimes called a log minus log survival plot. Parallel curves suggest proportional hazards holds.
Non-proportional hazards means that a predictor’s effect changes over time.
Hazard and survival functions will often cross:
As will graphs of \(-log(-log(S(t)))\):
Failing to account for non-constant hazard ratios threatens validity of Cox model estimates.
The Cox model can be extended in various ways to accommodate non-proportional hazards.
Baseline hazard function \(h_0(t)\)
The Cox model is popular because it does not require specification of the baseline hazard function, \(h_0(t)\), the hazard function for a subject with zero on all covariates.
Essentially, this means we can use the Cox model without knowing the shape of the underlying hazard functions.
Because no parameters related to a distribution (e.g. mean) are estimated, the Cox model is semiparametric.
Because \(h_0(t)\) is left unspecified, the Cox model only estimates covariate effects on the hazard function, not the hazard function itself.
No constant/intercept in Cox models, which prevents estimation of the hazard function (without additional methods).
Cox model with multiple predictors
The Cox is easily extended to accommodate multiple predictors, each of whose effects is assumed to be proportional over time.
\[h(t|X_1, X_2, ...X_p) = h_0(t)exp(b_1X_1 + b_2X_2 + ... + b_pX_p)\] Each coefficient \(b_i\) can be exponentiated to calculate a hazard ratio.
Fitting the Cox model with the survival package
New data set for Cox modeling
We will use the lung data set from the survival package for Cox modeling.
The data describe survival of patients with advanced lung cancer.
Variables:
time: survival time in days
status: 1=censored, 2=dead
age: age in years (assessed at beginning)
sex: 1=male, 2=female\(^\dagger\)
wt.loss: weight loss (pounds) in last 6 months
\(^\dagger\)Normally binary variables should be coded 0/1 so that the intercept is interpretable. However, with no intercept in the Cox model 1/2 coding is equivalent.
Fitting a Cox model
Use coxph(formula,data=)
The formula is a typical R regression formula:
Surv(time, event) ~ x1 + x2…
where x1 + x2… is a list of one or more predictor variables separated by +.
# fit cox model and save resultslung.cox <-coxph(Surv(time, status) ~ age + sex + wt.loss, data=lung)# summary of resultssummary(lung.cox)
Call:
coxph(formula = Surv(time, status) ~ age + sex + wt.loss, data = lung)
n= 214, number of events= 152
(14 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.0200882 1.0202913 0.0096644 2.079 0.0377 *
sex -0.5210319 0.5939074 0.1743541 -2.988 0.0028 **
wt.loss 0.0007596 1.0007599 0.0061934 0.123 0.9024
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0203 0.9801 1.0011 1.0398
sex 0.5939 1.6838 0.4220 0.8359
wt.loss 1.0008 0.9992 0.9887 1.0130
Concordance= 0.612 (se = 0.027 )
Likelihood ratio test= 14.67 on 3 df, p=0.002
Wald test = 13.98 on 3 df, p=0.003
Score (logrank) test = 14.24 on 3 df, p=0.003
Output interpretation:
sample size=214, 152 events, 14 observations dropped due to missingness
coef: log hazard ratio coefficients
a positive coefficient means the hazard increases (and survival decreases) as the predictor increases
We then pass the data frame of model results to ggplot2 to plot the hazard ratios and their confidence intervals.
# plot of hazard ratios and 95% CIsggplot(lung.cox.tab, aes(y=term, x=estimate, xmin=conf.low, xmax=conf.high)) +geom_pointrange() +# plots center point (x) and range (xmin, xmax)geom_vline(xintercept=1, color="red") +# vertical line at HR=1labs(x="hazard ratio", title="Hazard ratios and 95% CIs") +theme_classic()
Predicting survival with Cox estimates
We often want to estimate the survival functions for a subject with a particular set of covariate values.
The semiparametric Cox model does not estimate the survival function (or hazard function) directly.
However, we can use non-parameteric methods (similar to KM) to estimate the baseline survival function, \(S_0(t)\), the survival function for a subject with zero on all covariates.
\[{S}_0(t)={S}(t|X_1=0,X_2=0,...,X_p=0)\]
The Cox model coefficients then allow estimation of the survival function for a subject with any set of covariate values:
Because \(\hat{S}_0(t)\) is estimated non-parametrically, survival functions will again be step functions that change values only at observed event times.
Predicting survival after coxph() with survfit()
survfit(), when supplied a coxph model object, predicts survival functions.
If no covariate values are supplied to survfit(), the survival function is estimated for a subject with covariates at their means.
tidy() (or summary()) produces a table of the survival function estimated by survfit():
# predict survival function for subject with means on all covariatessurv.at.means <-survfit(lung.cox) #table of survival functiontidy(surv.at.means)
plot() plots the predicted survival curve and CI estimated by survfit().
# plot of predicted survival for subject at means of covariatesplot(surv.at.means, xlab="days", ylab="survival probability")
Predicting survival at specific covariate values
Specify a data frame of covariate values in the newdata= argument of survfit() to predict a survival function for each row of the data frame:
First we create the new data frame of 2 rows of covariate values.
# create new data for plotting: 1 row for each sex# and mean age and wt.loss for both rowsplotdata <-data.frame(age=mean(lung$age),sex=1:2,wt.loss=mean(lung$wt.loss, na.rm=TRUE))# look at new dataplotdata
age sex wt.loss
1 62.44737 1 9.831776
2 62.44737 2 9.831776
Then we supply the new data to survfit to predict separate survival functions for males and females, both at mean age and mean weight loss.
In the table of predicted survival functions produced by tidy(), the column suffixes .1 and .2 differentiate estimates between the 2 sexes.
e.g., estimate.1 is the survival function estimates for sex=1, males, while estimate.2 is for females.
# get survival function estimates for each sexsurv.by.sex <-survfit(lung.cox, newdata=plotdata) # one function for each sex# columns with .1 are for males, .1 for femalestidy(surv.by.sex)
Use ggsurvplot() from survminer for more customized plots.
Note: if a dataset was specified as newdata= in survfit() to generate survival function estimates, that dataset must be specified as data= in ggsurvplot() as well.
# data= is the same data used in survfit()# censor=F removes censoring symbolsggsurvplot(surv.by.sex, data=plotdata, censor=F,legend.labs=c("male", "female"))
We see that females have overall better survival, reflected by \(\widehat{HR}_{sex}=.594\).
Assessing the proportional hazards (PH) assumption
The plausibility of a model’s assumptions affects our confidence in the estimates.
Several methods have been developed to assess the PH assumption of the Cox model (including the log minus log plots discussed earlier).
Here we discuss 2 tools developed by Grambsch and Therneau (1994) based on Schoenfeld residuals1 that are available in the survival package.
cox.zph() performs a \(\chi^2\) test based on Schoenfeld residuals to test the hypothesis:
\(H_0\): covariate effect is constant (proportional) over time \(H_A\): covariate effect changes over time
The null hypothesis of PH is tested for each covariate individually and jointly as well.
cox.zph(lung.cox)
chisq df p
age 0.5077 1 0.48
sex 2.5489 1 0.11
wt.loss 0.0144 1 0.90
GLOBAL 3.0051 3 0.39
No strong evidence of violation of proportional hazards for any covariate, but some weak evidence for sex perhaps.
Another tool used to assess PH is a plot of a smoothed curve over the Schoenfeld residuals.
To create this plot, plot() the object returned by cox.zph().
These plots depict an estimate of the coefficient (labelled “Beta(t)” in the plot) (y-axis) across time (x-axis). PH is indicated by flat smoothed curve
plot(cox.zph(lung.cox))
Again we see some evidence of non-proportional hazards for sex, as the effect of sex seems to increase with time.
The effect of sex seems to be strongest at the beginning of follow-up, but then trends toward zero as time passes.
1Schoenfeld residuals are defined for each covariate and for each subject at the subject’s event time. The residual measures the difference between that subject’s covariate values and the average covariate values of those still at risk. For example, a positive residual for age means that a subject who failed is older than those still risk. If the residuals for age increase on average as time passes, that means that older patients are failing at higher rates than expected by the Cox model as time passes, suggesting that the effect of age is increasing with time.
Dealing with proportional hazards violations
Several strategies have been proposed to account for violation of PH. We discuss two here:
stratify by the non-PH variable
add an interaction of the non-PH variable with time to the model
If the change in the coefficient is not large enough to be clinically meaningfully, it can perhaps be ignored as well.
Imagine we are concerned that possible PH violation by sex in our Cox models may bias all of the model estimates.
We can address this concern with a sensitivity analysis to show how model inferences change depending on whether we ignore or address the PH violation by sex.
In the Cox model assuming PH for sex, we find evidence that the hazard increases with age but the effect of wt.loss is practically null or undeterminable.
# reprinting original model resultssummary(lung.cox)
Call:
coxph(formula = Surv(time, status) ~ age + sex + wt.loss, data = lung)
n= 214, number of events= 152
(14 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.0200882 1.0202913 0.0096644 2.079 0.0377 *
sex -0.5210319 0.5939074 0.1743541 -2.988 0.0028 **
wt.loss 0.0007596 1.0007599 0.0061934 0.123 0.9024
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0203 0.9801 1.0011 1.0398
sex 0.5939 1.6838 0.4220 0.8359
wt.loss 1.0008 0.9992 0.9887 1.0130
Concordance= 0.612 (se = 0.027 )
Likelihood ratio test= 14.67 on 3 df, p=0.002
Wald test = 13.98 on 3 df, p=0.003
Score (logrank) test = 14.24 on 3 df, p=0.003
Stratified Cox model
One approach to addressing a PH-violating variable is to estimate the model stratified by that variable.
Drawback: no coefficient will be estimated for the stratification variable.
In the stratified Cox model:
the model is estimated separately by stratum
the baseline hazard function, \(h_0(t)\), is allowed to be different across strata
this can accommodate the non-proportional effects of the stratification variable
parameter estimates are averaged across strata to generate one final set of estimates
Put the stratification variable inside strata() within the coxph() model formula:
lung.strat.sex <-coxph(Surv(time, status) ~ age + wt.loss +strata(sex), data=lung)summary(lung.strat.sex)
Call:
coxph(formula = Surv(time, status) ~ age + wt.loss + strata(sex),
data = lung)
n= 214, number of events= 152
(14 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.0192190 1.0194049 0.0096226 1.997 0.0458 *
wt.loss 0.0001412 1.0001412 0.0062509 0.023 0.9820
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.019 0.9810 1.000 1.039
wt.loss 1.000 0.9999 0.988 1.012
Concordance= 0.561 (se = 0.027 )
Likelihood ratio test= 4.09 on 2 df, p=0.1
Wald test = 3.99 on 2 df, p=0.1
Score (logrank) test = 4 on 2 df, p=0.1
The coefficients for age and wt.loss, averaged across males and females, are similar to those from the model where we assumed PH for sex, so our overall inferences are unchanged.
However, if the effect of sex is of primary interest, we should not stratify by sex.
Modeling time-varying coefficients
The Cox model can be extended to allow the coefficient representing a covariate effect to change over time, by interacting that covariate with some function of time.
However, unlike in other regression models, we cannot create this interaction term by simply multiplying the covariate by the time variable (unless the data are in a special structure, see survSplit()).
Instead, we strongly recommend the use of the time-transform function, tt(), to avoid easily-made mistakes.
Specify the nonPH covariate both by itself and inside of tt() in the model formula
Define the function tt() within coxph()
always start with tt = function(x,t,...)
then define the relationship between the covariate x and time t. For example:
x*t will allow coefficient to change with time linearly
x*log(t) allows coefficient to change with log (magnitudes) of time
x*(t>100) allows coefficient to take on 2 different values, one value when \(t<=100\) and another value \(t>100\).
Below we specify an interaction of sex with time itself, so that the effect of sex is allowed to change linearly with time.
# notice sex and tt(sex) in model formulalung.sex.by.time <-coxph(Surv(time, status) ~ age + wt.loss + sex +tt(sex), # sex and tt(sex) in formuladata=lung,tt=function(x,t,...) x*t) # linear change in effect of sexsummary(lung.sex.by.time)
Call:
coxph(formula = Surv(time, status) ~ age + wt.loss + sex + tt(sex),
data = lung, tt = function(x, t, ...) x * t)
n= 214, number of events= 152
(14 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.0194343 1.0196244 0.0096522 2.013 0.0441 *
wt.loss 0.0001260 1.0001261 0.0062502 0.020 0.9839
sex -0.9417444 0.3899470 0.3224791 -2.920 0.0035 **
tt(sex) 0.0013778 1.0013787 0.0008581 1.606 0.1084
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0196 0.9808 1.0005 1.0391
wt.loss 1.0001 0.9999 0.9879 1.0125
sex 0.3899 2.5645 0.2073 0.7337
tt(sex) 1.0014 0.9986 0.9997 1.0031
Concordance= 0.613 (se = 0.027 )
Likelihood ratio test= 17.23 on 4 df, p=0.002
Wald test = 15.86 on 4 df, p=0.003
Score (logrank) test = 16.44 on 4 df, p=0.002
The coef estimated for sex is the log-hazard-ratio at day \(t=0\), \(\hat{b}_{sex}=-.942\), corresponding to a hazard ratio of \(\hat{HR}_{sex}=exp(\hat{b}_{sex})=.39\).
The coef for tt(sex) is the change in the log-hazard-ratio for each additional day that passes.
Again, our inferences regarding age and wt.loss are mostly unchanged.
These sex coefficient estimates match the graph of the smoothed Schoenfeld residuals for sex. At the beginning of follow-up, the coefficient is close to -1, and it increases gradually over time.
plot(cox.zph(lung.cox), var="sex")
Time-varying covariates: data set up
The Cox PH model accommodates time-varying covariates, but the data should be structured in start-stop time format:
each subject can have multiple rows of data
2 time variables record the beginning and end of time intervals
status variable records the event status at the end of each interval
for intervals where the subject remains in the risk set, the status should be censored
for single event data, if the event occurs, it will be in the subject’s last interval
with no time gaps, the start time of intervals is the stop time of the previous interval (except for the first interval)
the time when a covariate changes value should be recorded as the beginning of a new interval
The survival package provides a function tmerge() to help get your data in this format. See vignette(timedep) for guidance on its usage.
The jasa1 dataset, which looks at survival for patients on a waiting list for heart transplant, is already set up in this format:
In this data set, transplant is a time-varying covariate:
Call:
coxph(formula = Surv(start, stop, event) ~ transplant + age +
surgery, data = jasa1)
n= 170, number of events= 75
coef exp(coef) se(coef) z Pr(>|z|)
transplant 0.01405 1.01415 0.30822 0.046 0.9636
age 0.03055 1.03103 0.01389 2.199 0.0279 *
surgery -0.77326 0.46150 0.35966 -2.150 0.0316 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
transplant 1.0142 0.9860 0.5543 1.8555
age 1.0310 0.9699 1.0033 1.0595
surgery 0.4615 2.1668 0.2280 0.9339
Concordance= 0.599 (se = 0.036 )
Likelihood ratio test= 10.72 on 3 df, p=0.01
Wald test = 9.68 on 3 df, p=0.02
Score (logrank) test = 10 on 3 df, p=0.02
We can follow with the same procedures as before, checking PH assumptions and plotting predicted survival curves.
# check PH assumptionscox.zph(jasa1.cox)
chisq df p
transplant 0.118 1 0.73
age 0.897 1 0.34
surgery 0.097 1 0.76
GLOBAL 1.363 3 0.71
No strong evidence that PH is violated.
# plot predicted survival by transplant group at mean age and surgery=0plotdata <-data.frame(transplant=0:1, age=-2.48, surgery=0)surv.by.transplant <-survfit(jasa1.cox, newdata=plotdata)ggsurvplot(surv.by.transplant, data=plotdata, # remember to supply data to ggsurvplot() for predicted survival after coxph()legend.labs=c("0", "1"), legend.title="transplant") # give proper labels and title to legend
Exercise 2
Now we will fit a Cox proprotional hazards model to the veteran data set.
time: survival time in weeks
status: status, 0=censored, 1=dead
trt: treatment, 1=standard, 2=test
age: age
Use a Cox proportional hazards model to estimate hazard ratios for the effects of treatment and age. Interpret the estimated hazard ratios.
Assess the proportional hazards assumption for both covariates using a chi-square test and graphs.
Fit a Cox model that allows the effect of trt to change over time (even though the original model was not very useful).
More about survival analysis
More R stuff
The coxph() function has some additional flexibility, including:
recurrent events modeling (see vignette("survival"))
competing risks modeling (see vignette("survival"))
frailty random effects, see ?frailty for more, but in general, use the coxme package
robust and cluster robust variance estimation, see the robust and cluster options in ?coxph
weights that can be treated as replication/frequency or sampling weights, see the weights argument in ?coxph
The survival package has the survreg() function for parametric regression models.
Additional useful survival analysis packages
coxme package, frailty random effects models, also written by Therneau
icenreg package, for interval censored data (note: data set up is different from survival)
References
References for the survival package
The survival package has some of the best vignettes (tutorials) of any R package.
Key vignettes for getting started (use vignette(package="survival") to see a full list):
vignette("survival") for a general introduction of the usage and capabilities of the survival package, as well as some theoretical survival analysis background
vignette("timedep") for guidance modeling time-varying covariates and time-varying coefficients
References for survival analysis
Grambsch, P. & Therneau, T. (1994), Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81, 515-26.
Therneau, T. M. & Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model, Third Edition. New York: Springer.
Kleinbaum, D. G., & Klein, M. (2012). Survival analysis: a self-learning text (Vol. 3). New York: Springer.
