Exercise 1 (Nelson-Aalen, Kaplan-Meier, and Cox regression)

The exercise may contain more material than you will be able to cover during the tutorial session. If so, you should skip questions e, h, and j.

From 1962 to 1969 a number of patients with liver cirrhosis at several hospitals in Copenhagen were included in a randomized clinical trial. The purpose of the study was to investigate whether patients treated with the hormone prednisone had a better survival than patients who got an inactive placebo treatment. In this exercise, we will restrict attention to the 386 patients in the study who had no excess fluid in the abdomen at entry. (Source: Andersen, Borgan, Gill & Keiding: Statistical Models Based on Counting Processes. Springer, 1993.)

We will study the effect of treatment with prednisone. We will also investigate the effect on survival of the covariates sex, age, and prothrombin index (a measurement based on a blood test of some coagulation factors produced by the liver, in percent of the normal value).

You may read the data into R by the command:

liver=read.table("https://www.mn.uio.no/math/english/people/aca/borgan/iwsm2026/liver.txt", header=T)

(You may copy this command and the other commands below and paste them into the source editor of RStudio, possibly for further editing.)

The data are organized with one line for each of the 386 patients, and with the following variables in the six columns:

status: indicator for death/censoring (1=dead; 0=censored)

time: time in days from start of treatment to death/censoring

treat: treatment (0=prednisone; 1=placebo)

sex: gender (0=female; 1=male)

age: age in years at start of treatment

prot: prothrombin index

We will start out by investigating the effect of one covariate at a time using Nelson-Aalen and Kaplan-Meier plots and univariate Cox regressions. To this end we will use the survival library, so this has to be loaded.

a) We first consider the effect of treatment. We may obtain Nelson-Aalen plots for the two treatment groups by the following commands (when time is converted to years).

fit.naa=survfit(Surv(time/365.25,status)~factor(treat), data=liver, stype=2)

plot(fit.naa, fun="cumhaz", xlab="Years after randomization", lty=1:2, main="Treatment")

legend("topleft",c("Prednisone","Placebo"),lty=1:2)

Comment: The option stype=2 in the survfit command instructs R to estimate the survival function as the exponential of the Nelson-Aalen estimator, and then the plot command with option fun="cumhaz" gives the Nelson-Aalen plots.

Perform the commands and interpret the Nelson-Aalen plots.

b) To obtain Kaplan-Meier plots for the two treatment groups we may give the commands

fit.km=survfit(Surv(time/365.25,status)~factor(treat), data=liver)

plot(fit.km, xlab="Years after randomization", lty=1:2, main="Treatment")

legend("bottomleft",c("Prednisone","Placebo"),lty=1:2)

Perform the commands and interpret the Kaplan-Meier plots.

c) We will then take a closer look at the prednisone treated patients. The commands below give a Kaplan-Meier plot with 95 % pointwise confidence limits for these patients.

fit.km=survfit(Surv(time/365.25,status)~1, data=liver, subset=(treat==0), conf.type="plain")

plot(fit.km, xlab="Years after randomization", lty=1:2, main="Prednisone treatment")

Perform the commands and use the plot to make a rough estimate of the median survival time for the prednisone treated patients. Also make rough estimates of the 95 % confidence limits for the median survival time. Then give the command quantile(fit.km) to get a precise estimate of the median survival time with confidence limits. Also find an estimate of the median survival time for the placebo treated patients with confidence limits.

d) We will then fit a Cox regression model with treatment as the only covariate. This is obtained by the commands:

fit.cox=coxph(Surv(time/365.25,status)~factor(treat), data=liver)

summary(fit.cox)

Perform the commands and interpret the results of the Cox regression.

e) Perform similar analyses as in a-d for the covariate sex.

We will also consider the numeric covariates age and prothrombin index. In order to make Nelson-Aalen and Kaplan-Meier plots for a numeric covariate, we first have to group the individuals based on their values of the covariates. For age we will consider the groups (i) 49 years or younger, (ii) 50-59 years, (iii) 60-69 years, and (iv) 70 years or older. To this end we create a new categorial covariate agegroup by the command:

liver$agegroup=cut(liver$age, breaks=c(0,49,59,69,100), labels=1:4)

f) Make a grouping of age by the command given above. Then compute Nelson-Aalen and Kaplan-Meier estimates for the four age groups and interpret the plots. What are the estimated median survival times for the four age groups?

g) Fit a Cox regression model using age as a numeric covariate and interpret the results.

Comment: For a numeric covariate like age, one should check if it is reasonable to assume a log-linear effect. There are several methods do so. However, we do not have time to discuss model check in the short course. But we will mention that one method for checking log-linearity, is to see if a smoothing spline for age gives a better fit than a model which assumes a log-linear effect of age. This may be achieved by the commands

cox.sage=coxph(Surv(time/365.25,status)~pspline(age),data=liver)

print(cox.sage)

termplot(cox.sage,se=T)

If you perform these computations, you will see that it is reasonable to assume a log-linear effect of age.

h) Perform similar analyses as in f and g for the prothrombin index. For the Nelson-Aalen and Kaplan-Meier plots, you may use the prothrombin groups (i) 49 or less, (ii) 50-69, (iii) 70-89, and (iv) 90 or more.

i) Perform a Cox regression using all four covariates and interpret the results.

Comment: A Cox regression model assumes proportional hazards (when there are no time-varying covariates). One method for checking proportional hazards, is to consider a model that for each time-fixed covariate x, also includes a time-varying covariate of the form g(t)x for a given function g(t), and then use the score test to see if the regression coefficients for these time dependent covariate are significantly different from zero. This is obtained by the commands

cox.all.zph=cox.zph(cox.all)

print(cox.all.zph)

where cox.all is the coxph-object we obtain for the model with all four covariates in question i. If you perform these tests, you will see that the assumption of proportional hazards seems ok.

j) Use the fitted model from question i to provide estimated survival functions for the following combination of the covariates:

A prednisone treated woman aged 60 years who has prothrombin index 100

A prednisone treated woman aged 60 years who has prothrombin index 40.

A placebo treated man aged 60 years who has prothrombin index 100

A placebo treated man aged 60 years who has prothrombin index 40

See the lectures for R commands.

If time allows, you may also consider other combination of the covariates.