1 AID data set, competing risk
The first analysis uses the AID data set and a competing risks transition shown in the figure below.
oldpar <- par(mar=c(0,0,0,0))
states <- c("Event free", "AIDS", "SI")
smat <- matrix(0, 3, 3, dimnames=list(states, states))
smat[1,2] <- smat[1,3] <- 1
statefig(1:2, smat)
The statefig routine is designed primarily for ease of use and creates figures
that are “good enough” for most uses.
1.1 Aalen-Johansen curves
We first create a multistate status variable and use it to plot the competing risk curves for the outcome. A key tool for dealing with multistate outcomes is replacement of the usual “status” variable of 0= censored, 1=event (or FALSE=censored TRUE = event) with a factor variable in the Surv function. This allows us to specify not just that an event occurred, but what type of event.
aidssi$event <- factor(aidssi$status, 0:2, c("censored", "AIDS", "SI"))
# The correct Aalen-Johansen curves
ajfit <- survfit(Surv(time, event) ~1, data = aidssi)
ajfit$transitions
#> to
#> from AIDS SI (censored)
#> (s0) 114 108 107
#> AIDS 0 0 0
#> SI 0 0 0
plot(ajfit, xmax = 13, col = 1:2, lwd = 2,
xlab = "Years from HIV infection", ylab = "Probability")
legend(8, .2, c("AIDS", "SI"), lty = 1, lwd = 2, col = 1:2, bty = 'n')
Figure 1.1: Correct Aalen-Johansen curves for the cumulative incidence of AIDS and SI
Since an initial state was not specified in the data or the survfit call, the function assumes that all subjects started in a common state (s0) = “state 0”. Like the (Intercept) term in a linear model fit, the name created by the survfit function is placed in parenthesis to avoid overlap with any variable names in the data. The transitions matrix shows that 114 subjects transitioned from this initial state to AIDS, 108 transitioned to SI, and 107 were censored.
A small footnote: The survfit routine in R produces Aalen-Johansen (AJ) estimates, which are applicable to any state space diagram (an arrangement of boxes and arrows). For a simple two state model such as alive & dead, the AJ estimate reduces to a Kaplan-Meier (KM). For a competing risk model such as this, the AJ estimate produces the same values as the cumulative incidence estimator. Put another way, the KM and CI are special cases of the AJ. The tutorial uses all three labels of KM, CI, and AJ.
We will use “Txx” to stand for figures or page numbers in the tutorial. Figure 1.2 (T2) shows the two Kaplan-Meier curves, with one going uphill and the other downhill. The estimated fraction with AIDS is the area above the red curve, the fraction with SI the area below the blue one, and the middle part is the fraction with neither. The fact that they cross is used to emphasize the inconsistency of the two estimates, i.e., that they add to more than 1.0.
# re-create figure T2
# KM curves that censor the other endpoint (a bad idea)
bad1 <- survfit(Surv(time, event=="AIDS") ~ 1, data = aidssi)
bad2 <- survfit(Surv(time, event=="SI") ~1, data = aidssi)
plot(bad1, conf.int = FALSE, xmax = 13,
xlab = "Years from HIV infection", ylab = "Probability")
lines(bad2, conf.int = FALSE, fun = 'event', xmax = 13)
text(c(8,8), c(.8, .22), c("AIDS", "SI"))
Figure 1.2: This figure (T2) shows the estimated survival curve for AIDS and probability of SI appearance, based on the naive Kaplan-Meier estimate.
Figure 1.3 (T3) shows the Aalen-Johansen curves in the same form. The
default in the survival package is to plot each curve on the natural axis
\(p_k(t)\) = probability of being in state \(k\) at time \(t\), which is the pstate
component of the survfit object. The authors of the tutorial like to use a
stacked display: the distance between the horizontal axis and the first curve is
the probability of being in state 1, the distance between the first and second
lines is the probability of being in state 2, etc. Since \(\sum_k p_k(t)=1\)
(everyone has to be somewhere), the final curve is a horizontal line at 1. The
following helper function pstack for stacked curves draws the plots in this
form. At time 0 the two lines are at y= 0 and 1: everyone is in the “neither
AIDS or SI” group.
pstack <- function(fit, top=FALSE, ...) {
temp <- survfit0(fit) # add the point at time 0
if (is.matrix(temp$pstate)) # usual case
temp$pstate <- t(apply(temp$pstate, 1, cumsum))
else if (is.array(temp$pstate))
temp$pstate <- aperm(apply(temp$pstate, 1:2, cumsum), c(2,3,1))
# this works because we don't change any other aspect of the survfit
# object, but only modify the probabilities.
if (top) plot(temp, noplot="", ...)
else plot(temp, noplot=temp$states[length(temp$states)], ...)
}# re-create figure T3
pstack(ajfit[c(2,1,3)], col=1, xmax=13, lwd=2,
xlab="Years from HIV infection", ylab="Probability")
lines(bad1, conf.int=FALSE, col="lightgray")
lines(bad2, conf.int=FALSE, fun='event', col='lightgray')
text(c(4, 8,8), c(.5, .85, .15), c("Event free", "AIDS", "SI"), col=1)
Figure 1.3: This figure (T3) shows estimates of probabilities of AIDS and SI appearance, based on the naive Kaplan-Meier (grey) and on cumulative incidence/Aalen-Johansen functions (black)
Figure 1.4 (T4) reorders the states so the event free is the top group. This author prefers the unstacked version (1.1), which shows more clearly that the probabilities of the two outcomes are very nearly the same.
pstack(ajfit[c(2,3,1)], xmax=13, lwd=2, col=1, ylim=c(0,1),
xlab="Years from HIV infection", ylab="Probability")
text(c(11, 11, 11), c(.2, .55, .9), c("AIDS", "SI", "Event free"))
Figure 1.4: This figure (T4) shows cumulative incidence curves of AIDS and SI appearance. The cumulative incidence functions are stacked; the distance between two curves represent the probability of the different events.
A last point is to note that for cumulative hazard functions, you can do the estimates separately for each endpoint, censoring the other. In the figure below, the estimates from the joint fit and those from the “bad” fits completely overlay each other.
plot(ajfit, cumhaz=TRUE, xmax=13, col=1:2, lty=2,
xlab="Years from HIV infection", ylab="Cumulative incidence")
lines(bad1, cumhaz=TRUE, conf.int=FALSE)
lines(bad2, cumhaz=TRUE, col=2, conf.int=FALSE)
Figure 1.5: Cumulative incidence fit using the cumlative hazard function
1.2 Proportional hazards models
The code below first fits a joint model for the two endpoints, followed by individual models for the two rates, each of which treats the other endpoint as censored.
cfit0 <- coxph(Surv(time, event) ~ ccr5, data = aidssi, id = patnr)
print(cfit0, digits=2)
#> Call:
#> coxph(formula = Surv(time, event) ~ ccr5, data = aidssi, id = patnr)
#>
#>
#> 1:2 coef exp(coef) se(coef) robust se z p
#> ccr5WM -1.24 0.29 0.31 0.30 -4.2 3e-05
#>
#>
#> 1:3 coef exp(coef) se(coef) robust se z p
#> ccr5WM -0.25 0.78 0.24 0.23 -1.1 0.3
#>
#> States: 1= (s0), 2= AIDS, 3= SI
#>
#> Likelihood ratio test=23 on 2 df, p=9.3e-06
#> n= 324, number of events= 220
#> (5 observations deleted due to missingness)
cfit1 <- coxph(Surv(time, event=="AIDS") ~ ccr5, data = aidssi)
print(cfit1, digits=2)
#> Call:
#> coxph(formula = Surv(time, event == "AIDS") ~ ccr5, data = aidssi)
#>
#> coef exp(coef) se(coef) z p
#> ccr5WM -1.24 0.29 0.31 -4 6e-05
#>
#> Likelihood ratio test=22 on 1 df, p=2.8e-06
#> n= 324, number of events= 113
#> (5 observations deleted due to missingness)
cfit2 <- coxph(Surv(time, event=="SI") ~ ccr5, data = aidssi)
print(cfit2, digits=2)
#> Call:
#> coxph(formula = Surv(time, event == "SI") ~ ccr5, data = aidssi)
#>
#> coef exp(coef) se(coef) z p
#> ccr5WM -0.25 0.78 0.24 -1.1 0.3
#>
#> Likelihood ratio test=1.2 on 1 df, p=0.27
#> n= 324, number of events= 107
#> (5 observations deleted due to missingness)Notice that the coefficients for the joint fit are identical to those where each endpoint is fit separately. This highlights a basic fact of multistate models
- Item hazards and cumulative hazards can be estimated one by one.
- Probability in state (absolute risk) must be estimated jointly.
The Cox model is a model for the hazards, and the separability allows for a lot
of freedom in how code and data sets are constructed. (It also gives more
opportunity for error, and for this reason the authors prefer the joint approach
of cfit0). The tutorial fits separate Cox models, page T2404, of the form
found in cfit1 and cfit2. We can also fit the joint model ‘by hand’ using a
stacked data set, which will have 329 rows = number of subjects for the AIDS
endpoint, followed by 329 rows for the SI endpoint. We had to be a bit cautious
since the tutorial uses “cause” for the event type and the data set aidsii
already has a variable by that name; hence the initial subset call.
temp <- subset(aidssi, select= c(patnr, time, ccr5))
temp1 <- data.frame(temp, status= 1*(aidssi$event=="AIDS"), cause="AIDS")
temp2 <- data.frame(temp, status= 1*(aidssi$event=="SI"), cause="SI")
stack <- rbind(temp1, temp2)
cfit3 <- coxph(Surv(time, status) ~ ccr5 * strata(cause), data=stack)
print(cfit3, digits=2)
#> Call:
#> coxph(formula = Surv(time, status) ~ ccr5 * strata(cause), data = stack)
#>
#> coef exp(coef) se(coef) z p
#> ccr5WM -1.24 0.29 0.31 -4.0 6e-05
#> ccr5WM:strata(cause)SI 0.98 2.67 0.39 2.5 0.01
#>
#> Likelihood ratio test=23 on 2 df, p=9.3e-06
#> n= 648, number of events= 220
#> (10 observations deleted due to missingness)The use of an interaction term gives a different form for the coefficients; the
second is now the difference in CCR-5 effect between the two endpoints. Which
form one prefers is a matter of taste. In the tutorial they used the equation
Surv(time, status) ~ ccr5*cause + strata(cause), which leads to a redundant
variable in the \(X\) matrix of the regression and a consequent NA coefficient in
the data set, but does not otherwise affect the results. We can also add
individual indicator variables to the stacked data set for ccr within type,
which gives yet another way of writing the same model. Last, we verify that the
partial likelihoods for our three versions are all identical.
stack$ccr5.1 <- (stack$ccr5=="WM") * (stack$cause == "AIDS")
stack$ccr5.2 <- (stack$ccr5=="WM") * (stack$cause == "SI")
cfit3b <- coxph(Surv(time, status) ~ ccr5.1 + ccr5.2 + strata(cause), data = stack)
cfit3b$coef
#> ccr5.1 ccr5.2
#> -1.2358147 -0.2542042
temp <- cbind(cfit0=cfit0$loglik, cfit3=cfit3$loglik, cfit3b=cfit3b$loglik)
rownames(temp) <- c("beta=0", "beta=final")
temp
#> cfit0 cfit3 cfit3b
#> beta=0 -1116.689 -1116.689 -1116.689
#> beta=final -1105.102 -1105.102 -1105.102We can also fit a models where the effect of ccr5 on the two types of outcome is
assumed to be equal. (We agree with the tutorial that there is not good medical
reason for such an assumption, the model is simply for illustration.) Not
surprisingly, the realized coefficient is midway between the estimates of the
ccr effect on the two separate endpoints. The second fit uses the joint model
approach by adding a constraint. In this case the formula argument for coxph
is a list. The first element of the list is a standard formula containing the
response and a set of covariates, and later elements, and the second, third,
etc. elements of the list are of the form state1:state2 ~ covariates. These
later element modify the formula for selected pairs of states. In this case the
second element specifies transitions from state 1 to 2 and 1:3 should share a
common ccr5 coefficient.
common1 <- coxph(Surv(time, status) ~ ccr5 + strata(cause), data=stack)
print(common1, digits=2)
#> Call:
#> coxph(formula = Surv(time, status) ~ ccr5 + strata(cause), data = stack)
#>
#> coef exp(coef) se(coef) z p
#> ccr5WM -0.70 0.50 0.19 -3.8 2e-04
#>
#> Likelihood ratio test=16 on 1 df, p=5e-05
#> n= 648, number of events= 220
#> (10 observations deleted due to missingness)
common1b <- coxph(list( Surv(time, event) ~ 1,
1:2 + 1:3 ~ ccr5/common ),
data=aidssi, id=patnr)At this point the tutorial explores an approach that we find problematic, which
is to fit models to the stacked data set without including the stratum. The
partial likelihood for the Cox model has a term for each event time, each term
is a ratio that compares the risk score of the event (numerator) to the sum of
risk scores for all subjects who were at risk for the event (denominator). When
the stack data set is fit without a strata statement, like below, then at each
event time the “risk set” will have 2 clones of each subject, one labeled with
covariate cause = AIDS and the other as SI. If we look closely, the estimated
coefficient from this second fit is almost identical to the stratified fit
common1, however.
common2 <- coxph(Surv(time, status) ~ ccr5, data = stack)
all.equal(common2$coef, common1$coef)
#> [1] "Mean relative difference: 8.40399e-05"In fact, if the Breslow approximation is used for ties, one can show that the
partial likelihood (PL) values for the two fits will satisfy the identity
PL(common2) = PL(common1) - d log(2), where \(d\) is the total number of events.
Since the two partial likelihoods differ by a constant, they will maximize at
the same location, i.e., give exactly the same coefficient estimates. One can
further show that if cause is added to the second model as a covariate, that
this will not change the ccr5 coefficient, while adding an estimate of the
relative proportion of events of each type.
# reprise common1 and common2, using the breslow option
test1 <- coxph(Surv(time, status) ~ ccr5 + strata(cause), stack,
ties='breslow')
test2 <- coxph(Surv(time, status) ~ ccr5, stack, ties='breslow')
all.equal(test2$loglik + test2$nevent * log(2), test1$loglik)
#> [1] TRUE
all.equal(test2$coef, test1$coef)
#> [1] TRUE
test3 <- coxph(Surv(time, status) ~ ccr5 + cause, stack, ties='breslow')
test3
#> Call:
#> coxph(formula = Surv(time, status) ~ ccr5 + cause, data = stack,
#> ties = "breslow")
#>
#> coef exp(coef) se(coef) z p
#> ccr5WM -0.70115 0.49601 0.18597 -3.770 0.000163
#> causeSI -0.05456 0.94690 0.13489 -0.404 0.685867
#>
#> Likelihood ratio test=16.62 on 2 df, p=0.0002459
#> n= 648, number of events= 220
#> (10 observations deleted due to missingness)
all.equal(test3$coef[1], test1$coef)
#> [1] TRUEThese identities do not assure the author that this pseudo risk set approach,
where subjects are duplicated, is a valid way to estimate the ccr5 effect under
the assumption of a common baseline hazard. The first model common1 can be
directly fit in the multistate framework by adding the constraint of a common
ccr5 effect for the two transitions; this is found above as common1b. One can
not directly fit a version of test2 using the multistate model, however, as
the underlying code for multistate fits rigorously enforces a “one copy”
principle: during the entire period of time that a subject is at risk, there
should be exactly one copy of that subject present in the data set. See the
survcheck routine for a more detailed discussion.
1.3 Predicted curves
We can now generate predicted Aalen-Johansen curves from the Cox model fits. As with any Cox model, this starts by deciding who to predict, i.e. the set of covariate values at which to obtain a prediction. For a model with a single binary variable this is an easy task.
# re-create figure T5 in a single panel
dummy <- data.frame(ccr5=c("WW", "WM"))
pred.aj <- survfit(cfit0, newdata=dummy)
dim(pred.aj)
#> data states
#> 2 3
pred.aj$states
#> [1] "(s0)" "AIDS" "SI"The resulting curves have an apparent dimension of (number of strata, number of covariate patterns, number of states). We plot subsets of the curves by using subscripts. (When there are no strata in the coxph fit (1 stratum) the code allows one to omit the first subscript.)
oldpar <- par(mfrow=c(1,2))
plot(pred.aj[,,"AIDS"], lwd=2, col=c("black", "gray"),
xmax=13, ylim=c(0,.5),
xlab="Years from HIV infection", ylab="Probability of AIDS")
text(c(9.5, 10), c(.3, .1), c("WW", "WM"))
plot(pred.aj[,,"SI"], lwd=2, col=c("black", "gray"),
xmax=13, ylim=c(0,.5),
xlab="Years from HIV infection", ylab="Probability of SI")
text(c(8.5, 9), c(.33, .25), c("WW", "WM"))
Figure 1.6: Figure T5 showing predicted curves for AIDS and SI stratified by subjects with CCR5 wild-type (WW) and mutant (WM).
Predicted survival curves from the two fits to individual endpoints
suffer from the same issue as the individual Kaplan-Meier curves
bad1 and bad2: the predicted risk risk of having either
AIDS or SI will be greater than 1 for some time points, which is clearly
impossible. Absolute risk estimates must be done jointly.
The tutorial at this point uses simulation data to further elucidate
the underlying issues between per-endpoint and joint estimates,
which we will not replicate.
1.4 Fine-Gray fits
We can also fit Fine-Gray models for AIDS and SI appearance. In the survival package this is done by creating a special data set - one for each endpoint. Ordinary Cox model code can then be applied to those data sets.
fdata1 <- finegray(Surv(time, event) ~ ., data = aidssi, etype = 'AIDS')
fgfit1 <- coxph(Surv(fgstart, fgstop, fgstatus) ~ ccr5, data = fdata1,
weight = fgwt)
fgfit1
#> Call:
#> coxph(formula = Surv(fgstart, fgstop, fgstatus) ~ ccr5, data = fdata1,
#> weights = fgwt)
#>
#> coef exp(coef) se(coef) robust se z p
#> ccr5WM -1.0043 0.3663 0.3059 0.2830 -3.549 0.000387
#>
#> Likelihood ratio test=13.94 on 1 df, p=0.000189
#> n= 2668, number of events= 113
#> (21 observations deleted due to missingness)
fdata2 <- finegray(Surv(time, event) ~., aidssi, etype="SI")
fgfit2 <- coxph(Surv(fgstart, fgstop, fgstatus) ~ ccr5, fdata2,
weight = fgwt)
fgfit2
#> Call:
#> coxph(formula = Surv(fgstart, fgstop, fgstatus) ~ ccr5, data = fdata2,
#> weights = fgwt)
#>
#> coef exp(coef) se(coef) robust se z p
#> ccr5WM 0.02366 1.02394 0.23554 0.21832 0.108 0.914
#>
#> Likelihood ratio test=0.01 on 1 df, p=0.9202
#> n= 2195, number of events= 107
#> (31 observations deleted due to missingness)The predicted curves based on the Fine-Gray model 1.7 (T8) use the ordinary survival tools (not Aalen-Johansen), since they are ordinary Cox models on a special data set.
# re-create figure T8: Fine-Gray curves
fgsurv1<-survfit(fgfit1, newdata=dummy)
fgsurv2<-survfit(fgfit2, newdata=dummy)
oldpar <- par(mfrow=c(1,2), mar=c(4.1, 3.1, 3.1, 1)) #leave room for title
plot(fgsurv1, col=1:2, lty=c(1,1,2,2), lwd=2, xmax=13,
ylim=c(0, .5),fun='event',
xlab="Years from HIV infection", ylab="Probability")
title("AIDS")
plot(fgsurv2, col=1:2, lty=c(1,1,2,2), lwd=2, xmax=13,
ylim=c(0, .5), fun='event',
xlab="Years from HIV infection", ylab="Probability")
title("SI appearance") 
Figure 1.7: Reproduction of Figure T8 showing cumulative incidence functions for AIDS (left) and SI appearance (right), for CCR5 wild-type (WW) and mutant (WM), based on the Fine and Gray model.
The last plot 1.8 (T9) in this section of the tutorial contains the Aalen-Johansen non-parametric fits stratified by CCR5 status.
# re-create figure T9: curves by CCR type
aj2 <- survfit(Surv(time, event) ~ ccr5, data = aidssi)
oldpar <- par(mfrow=c(1,2))
plot(aj2[,"AIDS"], xmax=13, col=1:2, lwd=2, ylim=c(0, .5),
xlab="Years from HIV infection", ylab="Probability of AIDS")
text(c(10, 10), c(.35, .07), c("WW", "WM"))
plot(aj2[,"SI"], xmax=13, col=1:2, lwd=2, ylim=c(0, .5),
xlab="Years from HIV infection", ylab="Probability of SI")
text(c(8, 8), c(.34, .18), c("WW", "WM"))
Figure 1.8: Reproduction of Figure T9 showing non-parametric cumulative incidence functions for AIDS (left) and SI appearance (right), for CCR5 wild-type (WW) and mutant (WM)
