--- title: "Competing risk analysis" author: "Maxime Turgeon" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: fig_height: 6 fig_width: 8 toc: yes toc_depth: 4 vignette: > %\VignetteIndexEntry{Competing risk analysis} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} bibliography: reference.bib editor_options: chunk_output_type: console --- ## Data We will use the same data that was used in Scrucca *et al* [-@scrucca2010regression]. The data was downloaded from the main author's website; it is also available as part of this package. ```{r echo=TRUE, eval=TRUE} set.seed(12345) library(casebase) data(bmtcrr) head(bmtcrr) ``` We will perform a competing risk analysis on data from 177 patients who received a stem cell transplant for acute leukemia. The event of interest in relapse, but other competing causes (e.g. transplant-related death) need to be taken into account. We also want to take into account the effect of several covariates such as Sex, Disease (lymphoblastic or myeloblastic leukemia, abbreviated as ALL and AML, respectively), Phase at transplant (Relapse, CR1, CR2, CR3), Source of stem cells (bone marrow and peripheral blood, coded as BM+PB, or peripheral blood, coded as PB), and Age. Below, we reproduce their Table 1: | Variable | Description | Statistical summary | | -------- | ---------------------- | ---------------------- | | Sex | Sex | M=Male (100)
F=Female (77) | | D | Disease | ALL (73)
AML (104) | | Phase | Phase | CR1 (47)
CR2 (45)
CR3 (12)
Relapse (73) | | Source | Type of transplant | BM+PB (21)
PB (156) | | Age | Age of patient (years) | 4–62
30.47 (13.04) | | Ftime | Failure time (months) | 0.13–131.77
20.28 (30.78) | | Status | Status indicator | 0=censored (46)
1=relapse (56)
2=competing event (75) | The statistical summary is generated differently for continuous and categorical variables: - For continuous variables, we are given the range, followed by the mean and standard deviation. - For categorical variables, we are given the counts for each category. Note that failure time can also correspond to censoring. ## Population-time plots In order to try and visualize the incidence density of relapse, we can look at a population-time plot: on the X-axis we have time, and on the Y-axis we have the size of the risk set at a particular time point. Failure times associated to the event of interest can then be highlighted on the plot using red dots. ```{r poptime1, eval=TRUE} nobs <- nrow(bmtcrr) ftime <- bmtcrr$ftime ord <- order(ftime, decreasing = FALSE) # We split the person-moments in four categories: # 1) at-risk # 2) main event # 3) competing event # 4) censored yCoords <- cbind(cumsum(bmtcrr[ord, "Status"] == 2), cumsum(bmtcrr[ord, "Status"] == 1), cumsum(bmtcrr[ord, "Status"] == 0)) yCoords <- cbind(yCoords, nobs - rowSums(yCoords)) # Plot only at-risk plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs), xlab = 'Follow-up time', ylab = 'Population') polygon(c(0, 0, ftime[ord], max(ftime), 0), c(0, nobs, yCoords[,4], 0, 0), col = "grey90") cases <- bmtcrr[, "Status"] == 1 # randomly move the cases vertically moved_cases <- yCoords[cases[ord], 4] * runif(sum(cases)) points((ftime[ord])[cases[ord]], moved_cases, pch = 20, col = "red", cex = 1) ``` We can right away draw a few conclusions from this plot: first of all, we get a sense of how quickly the size of the risk set changes over time. We also see that the incidence density is non-constant: most relapses occur before 15 months. Finally, we also see that the risk set keeps shrinking after the last event has occurred; this could be due to either censoring or the competing event. To get an idea of whether only relapse is responsible for the shrinking of the risk set in the first few months of follow-up, we can also keep track of how many events have occurred at each time point: ```{r poptime2, eval=TRUE} # Plot at-risk and events plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs), xlab = 'Follow-up time', ylab = 'Population') polygon(x = c(0,ftime[ord], max(ftime), 0), y = c(0, yCoords[,2], 0, 0), col = "firebrick3") polygon(x = c(0, ftime[ord], ftime[rev(ord)], 0, 0), y = c(0, yCoords[,2], rev(yCoords[,2] + yCoords[,4]), nobs, 0), col = "grey90") # randomly move the cases vertically moved_cases <- yCoords[cases[ord], 2] + yCoords[cases[ord], 4] * runif(sum(cases)) points((ftime[ord])[cases[ord]], moved_cases, pch = 20, col = "red", cex = 1) legend("topright", legend = c("Relapse", "At-risk"), col = c("firebrick3", "grey90"), pch = 15) ``` Therefore, there is also censoring and loss due to competing events happening in the first few months. However, with this plot, we can't differentiate between the two contributions. For this reason we can also keep track of the number of competing events at each time point: ```{r poptime3, eval=TRUE} plot(0, type = 'n', xlim = c(0, max(ftime)), ylim = c(0, nobs), xlab = 'Follow-up time', ylab = 'Population') polygon(x = c(0, max(ftime), max(ftime), 0), y = c(0, 0, nobs, nobs), col = "white") # Event of interest polygon(x = c(0,ftime[ord], max(ftime), 0), y = c(0, yCoords[,2], 0, 0), col = "firebrick3") # Risk set polygon(x = c(0, ftime[ord], ftime[rev(ord)], 0, 0), y = c(0, yCoords[,2], rev(yCoords[,2] + yCoords[,4]), nobs, 0), col = "grey90") # Competing event polygon(x = c(0, ftime[ord], max(ftime), 0), y = c(nobs, nobs - yCoords[,1], nobs, nobs), col = "dodgerblue2") # randomly move the cases vertically moved_cases <- yCoords[cases[ord], 2] + yCoords[cases[ord], 4] * runif(sum(cases)) points((ftime[ord])[cases[ord]], moved_cases, pch = 20, col = "red", cex = 1) legend("topright", legend = c("Relapse", "Competing event", "At-risk"), col = c("firebrick3", "dodgerblue2", "grey90"), pch = 15) ``` From this last plot, we can see that there is no censoring during the first 10 months. Moreover, we see that the last competing event occurs around 20 months. Putting all this information together, we have evidence of two types of patients: very sick patients who either relapse or have a competing event early on, and healthier patients who are eventually lost to follow-up. ## Analysis We now turn to the analysis of this dataset. The population-time plots above give evidence of non-constant hazard; therefore, we will explicitly include time in the model. Note that we also include all other variables as possible confounders. First, we include time as a linear term: ```{r eval=TRUE, warning=FALSE} model1 <- fitSmoothHazard(Status ~ ftime + Sex + D + Phase + Source + Age, data = bmtcrr, ratio = 100, time = "ftime") summary(model1) ``` Because of the results in Turgeon *et al* [-@turgeonCompRisk], the standard errors we obtain from the multinomial logit fit are asymptotically correct, and therefore can be used to construct asymptotic confidence intervals. From this summary, we see that time is indeed significant, as is Phase (only relapse vs. CR1). Interestingly, we see that the type of disease is only significant for the event of interest, whereas the type of transplant and the age of the patient are only significant for the competing event. Next, we include the logarithm of time in the model (which leads to a Weibull hazard): ```{r eval=TRUE, warning=FALSE} model2 <- fitSmoothHazard(Status ~ log(ftime) + Sex + D + Phase + Source + Age, data = bmtcrr, ratio = 100, time = "ftime") summary(model2) ``` As we can see, the results are similar to the ones with a Gompertz hazard, although Sex is now significant for the competing event. Finally, using splines, we can be quite flexible about the way the hazard depends on time: ```{r eval=TRUE, warning=FALSE} model3 <- fitSmoothHazard( Status ~ splines::bs(ftime) + Sex + D + Phase + Source + Age, data = bmtcrr, ratio = 100, time = "ftime") summary(model3) ``` Again, we see that the results are quite similar for this third model. ### Absolute risk We now look at the 2-year risk of relapse: ```{r absRisk, eval=TRUE, warning = FALSE} linearRisk <- absoluteRisk(object = model1, time = 24, newdata = bmtcrr[1:10,]) logRisk <- absoluteRisk(object = model2, time = 24, newdata = bmtcrr[1:10,]) splineRisk <- absoluteRisk(object = model3, time = 24, newdata = bmtcrr[1:10,]) ``` ```{r absRiskPlot, eval=TRUE} plot(linearRisk, logRisk, xlab = "Linear", ylab = "Log/Spline", pch = 19, xlim = c(0,1), ylim = c(0,1), col = 'red') points(linearRisk, splineRisk, col = 'blue', pch = 19) abline(a = 0, b = 1, lty = 2, lwd = 2) legend("topleft", legend = c("Log", "Spline"), pch = 19, col = c("red", "blue")) ``` ## Session information ```{r echo=FALSE, eval=TRUE} print(sessionInfo(), locale = FALSE) ``` ## References

  1. Efron, Bradley. 1977. "The Efficiency of Cox's Likelihood Function for Censored Data." Journal of the American Statistical Association 72 (359). Taylor & Francis Group: 557–65.

  2. Hanley, James A, and Olli S Miettinen. 2009. "Fitting Smooth-in-Time Prognostic Risk Functions via Logistic Regression." The International Journal of Biostatistics 5 (1).

  3. Mantel, Nathan. 1973. "Synthetic Retrospective Studies and Related Topics." Biometrics. JSTOR, 479–86.

  4. Saarela, Olli. 2015. "A Case-Base Sampling Method for Estimating Recurrent Event Intensities." Lifetime Data Analysis. Springer, 1–17.

  5. Saarela, Olli, and Elja Arjas. 2015. "Non-Parametric Bayesian Hazard Regression for Chronic Disease Risk Assessment." Scandinavian Journal of Statistics 42 (2). Wiley Online Library: 609–26.

  6. Scrucca, L, A Santucci, and F Aversa. 2010. "Regression Modeling of Competing Risk Using R: An in Depth Guide for Clinicians." Bone Marrow Transplantation 45 (9). Nature Publishing Group: 1388–95.

  7. Kalbfleisch, John D., and Ross L. Prentice. The statistical analysis of failure time data. Vol. 360. John Wiley & Sons, 2011.

  8. Cox, D. R. "Regression models and life tables." Journal of the Royal Statistical Society 34 (1972): 187-220.