SISCER 2021: Joint Modeling of Longitudinal and Survival Data
Benjamin French; b.french@vumc.org
This file provides an analysis in R of the data collected from a randomized controlled trial of treatment with D-penicillamine for primary biliary cirrhosis (PBC), a rare autoimmune liver disease.
Disclaimer: because this is a didactic exercise, we will consider a range of potential models. In practice, model fitting should be guided by a priori scientific knowledge and exploratory data analysis, and should be directly relevant to the scientific question of interest.
1 Preliminaries
First, install the JM, ggplot2, survminer, and dplyr extension packages.
install.packages(c('JM', 'ggplot2', 'survminer', 'dplyr'), repos='https://cloud.r-project.org')Next, load these packages. Note that the JM package automatically loads the nlme (for fitting mixed-effects models) and survival (for fitting Cox regression models) extension packages, among others. I also included some options for graphics as a user-specified function. (Some of the figures in this document are meant to be of publication quality, while others are only used for exploring the data and evaluating fitted models.)
library(JM)## Loading required package: MASS## Loading required package: nlme## Loading required package: splines## Loading required package: survivallibrary(ggplot2)
library(survminer)
library(dplyr)
theme_set(theme_bw())
mytheme <- function(...){
theme(panel.grid.minor=element_blank(),
plot.margin=unit(c(0.5, 0.5, 0.5, 0.5), 'lines'),
axis.title.x=element_text(size=12, vjust=0),
axis.title.y=element_text(size=12, vjust=1, angle=90),
legend.title=element_text(size=12, face='bold'),
legend.text=element_text(size=12),
legend.position='top', legend.direction='horizontal',
legend.key=element_blank(),
legend.key.width=unit(1.5,"cm"),
legend.key.height=unit(0.35,'cm'),
strip.text.y=element_text(size=12))}2 PBC data
The PBC data are provided as part of the JM extension package. Note that there is a help file for the dataset that can be accessed using ?pbc2. The dataset is in so-called long format, in which each participant has multiple rows, with each row corresponding to a measurement time.
data(pbc2)
head(pbc2, 11)## id years status drug age sex year ascites hepatomegaly spiders
## 1 1 1.09517 dead D-penicil 58.76684 female 0.0000000 Yes Yes Yes
## 2 1 1.09517 dead D-penicil 58.76684 female 0.5256817 Yes Yes Yes
## 3 2 14.15234 alive D-penicil 56.44782 female 0.0000000 No Yes Yes
## 4 2 14.15234 alive D-penicil 56.44782 female 0.4983025 No Yes Yes
## 5 2 14.15234 alive D-penicil 56.44782 female 0.9993429 No Yes Yes
## 6 2 14.15234 alive D-penicil 56.44782 female 2.1027270 No Yes Yes
## 7 2 14.15234 alive D-penicil 56.44782 female 4.9008871 Yes Yes Yes
## 8 2 14.15234 alive D-penicil 56.44782 female 5.8892783 Yes Yes Yes
## 9 2 14.15234 alive D-penicil 56.44782 female 6.8858833 Yes Yes Yes
## 10 2 14.15234 alive D-penicil 56.44782 female 7.8907020 Yes Yes Yes
## 11 2 14.15234 alive D-penicil 56.44782 female 8.8325485 Yes Yes Yes
## edema serBilir serChol albumin alkaline SGOT platelets prothrombin histologic status2
## 1 edema despite diuretics 14.5 261 2.60 1718 138.0 190 12.2 4 1
## 2 edema despite diuretics 21.3 NA 2.94 1612 6.2 183 11.2 4 1
## 3 No edema 1.1 302 4.14 7395 113.5 221 10.6 3 0
## 4 No edema 0.8 NA 3.60 2107 139.5 188 11.0 3 0
## 5 No edema 1.0 NA 3.55 1711 144.2 161 11.6 3 0
## 6 No edema 1.9 NA 3.92 1365 144.2 122 10.6 3 0
## 7 edema no diuretics 2.6 230 3.32 1110 131.8 135 11.3 3 0
## 8 edema despite diuretics 3.6 NA 2.92 996 131.8 100 11.5 3 0
## 9 edema despite diuretics 4.2 NA 2.73 860 145.7 103 11.5 3 0
## 10 edema despite diuretics 3.6 244 2.80 779 119.0 113 11.5 3 0
## 11 edema despite diuretics 4.6 237 2.67 669 88.0 100 11.5 3 0str(pbc2)## 'data.frame': 1945 obs. of 20 variables:
## $ id : Factor w/ 312 levels "1","2","3","4",..: 1 1 2 2 2 2 2 2 2 2 ...
## $ years : num 1.1 1.1 14.2 14.2 14.2 ...
## $ status : Factor w/ 3 levels "alive","transplanted",..: 3 3 1 1 1 1 1 1 1 1 ...
## $ drug : Factor w/ 2 levels "placebo","D-penicil": 2 2 2 2 2 2 2 2 2 2 ...
## $ age : num 58.8 58.8 56.4 56.4 56.4 ...
## $ sex : Factor w/ 2 levels "male","female": 2 2 2 2 2 2 2 2 2 2 ...
## $ year : num 0 0.526 0 0.498 0.999 ...
## $ ascites : Factor w/ 2 levels "No","Yes": 2 2 1 1 1 1 2 2 2 2 ...
## $ hepatomegaly: Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ spiders : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ edema : Factor w/ 3 levels "No edema","edema no diuretics",..: 3 3 1 1 1 1 2 3 3 3 ...
## $ serBilir : num 14.5 21.3 1.1 0.8 1 1.9 2.6 3.6 4.2 3.6 ...
## $ serChol : int 261 NA 302 NA NA NA 230 NA NA 244 ...
## $ albumin : num 2.6 2.94 4.14 3.6 3.55 3.92 3.32 2.92 2.73 2.8 ...
## $ alkaline : int 1718 1612 7395 2107 1711 1365 1110 996 860 779 ...
## $ SGOT : num 138 6.2 113.5 139.5 144.2 ...
## $ platelets : int 190 183 221 188 161 122 135 100 103 113 ...
## $ prothrombin : num 12.2 11.2 10.6 11 11.6 10.6 11.3 11.5 11.5 11.5 ...
## $ histologic : int 4 4 3 3 3 3 3 3 3 3 ...
## $ status2 : num 1 1 0 0 0 0 0 0 0 0 ...Data are available for 312 study participants, with an unequal number of observations per participant.
subset(pbc2, !is.na(serBilir)) %>%
count(id, name='m') -> num
with(num, table(m))## m
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 27 26 32 44 30 23 19 23 17 23 16 12 6 5 6 3Our analysis will focus on serum bilirubin levels as a measure of liver function. A histogram of bilirubin levels at baseline (year 0) does not exhibit a symmetric distribution.
pbc2 <- within(pbc2, {
levels(drug) <- c('Placebo', 'Treatment')})
ggplot(subset(pbc2, year==0), aes(x=serBilir)) +
geom_histogram(color='black', fill='white') +
facet_wrap(~ drug) +
scale_x_continuous(name='Serum bilirubin, mg/dl') +
scale_y_continuous(name='Count') +
mytheme()
Using ggplot, plot the observed serum bilirubin levels over time for the first 16 study participants, along with a linear regression line and an indicator of when the participant died or was censored. (Liver transplantation is included as a censoring event, which we could debate.) Note that because bilirubin levels are positively skewed, the y axis is displayed on the log (base 10) scale.
ggplot(subset(pbc2, is.element(id, 1:16)), aes(year, serBilir)) +
facet_wrap( ~ id) +
geom_vline(aes(xintercept=years, color=factor(status!='dead')), lwd=1) +
geom_point() +
geom_smooth(se=FALSE, method='lm') +
scale_x_continuous(name='Time, years', limits=c(0,15)) +
scale_y_continuous(trans='log10', name='Serum bilirubin, mg/dl') +
scale_colour_manual(name='Status', breaks=c(FALSE, TRUE), values=c('red', 'darkgray'),
labels=c('Death', 'Censored')) +
mytheme()## `geom_smooth()` using formula 'y ~ x'
Note that there are other options for scatterplot smoothers; see ?geom_smooth for details.
3 Longitudinal data analysis
First, using ggplot, plot the observed individual-level profiles (or, trajectories) in bilirubin level over time, along with a smoothing line for each treatment group to display the average trend.
ggplot(pbc2, aes(year, serBilir, group=id)) +
geom_line(aes(color=drug), alpha=0.4) +
geom_smooth(aes(year, serBilir, group=drug, color=drug), se=FALSE, size=1.5) +
scale_x_continuous(name='Time, years', limits=c(0,15)) +
scale_y_continuous(trans='log10', name='Serum bilirubin, mg/dl') +
scale_color_manual(name='Intervention', values=c('#00BFC4', '#F8766D')) +
mytheme()
Next, using lme, fit a series a linear mixed-effects models, considering alternative specifications for the mean model (i.e., main effects of treatment and time, treatment-by-time interaction) and alternative specifications for the random effects (i.e., random intercepts only, or random intercepts and random slopes). Note that because bilirubin levels are log-transformed, coefficient estimates correspond to differences on the log scale. Exponentiate those coefficients to obtain ratios in bilirubin levels. Math facts: \(\log(a)-\log(b)=\log(a/b)\) and \(\exp\{\log(a/b)\}=a/b\).
# Main effects of treatment and time, random intercepts
ml1 <- lme(log(serBilir) ~ drug + year, random = ~ 1 | id, data=pbc2)
summary(ml1)## Linear mixed-effects model fit by REML
## Data: pbc2
## AIC BIC logLik
## 3797.89 3825.747 -1893.945
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.09385 0.4920293
##
## Fixed effects: log(serBilir) ~ drug + year
## Value Std.Error DF t-value p-value
## (Intercept) 0.6265245 0.09080995 1632 6.899293 0.0000
## drugTreatment -0.1107256 0.12709756 310 -0.871186 0.3843
## year 0.0950894 0.00433025 1632 21.959307 0.0000
## Correlation:
## (Intr) drgTrt
## drugTreatment -0.707
## year -0.103 -0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.46753951 -0.54604155 -0.02275856 0.49525690 4.42743992
##
## Number of Observations: 1945
## Number of Groups: 312# Interaction between treatment and time, random intercepts
ml2 <- lme(log(serBilir) ~ drug*year, random = ~ 1 | id, data=pbc2)
summary(ml2)## Linear mixed-effects model fit by REML
## Data: pbc2
## AIC BIC logLik
## 3805.78 3839.206 -1896.89
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.09333 0.4919563
##
## Fixed effects: log(serBilir) ~ drug * year
## Value Std.Error DF t-value p-value
## (Intercept) 0.6392973 0.09127541 1631 7.004048 0.0000
## drugTreatment -0.1357919 0.12842957 310 -1.057326 0.2912
## year 0.0891536 0.00621618 1631 14.342180 0.0000
## drugTreatment:year 0.0115238 0.00866303 1631 1.330227 0.1836
## Correlation:
## (Intr) drgTrt year
## drugTreatment -0.711
## year -0.147 0.104
## drugTreatment:year 0.105 -0.147 -0.718
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.40293143 -0.53726172 -0.03042006 0.49145048 4.46892643
##
## Number of Observations: 1945
## Number of Groups: 312# Treatment-specific time effects, random intercepts
ml3 <- lme(log(serBilir) ~ year + drug:year, random = ~ 1 | id, data=pbc2)
summary(ml3)## Linear mixed-effects model fit by REML
## Data: pbc2
## AIC BIC logLik
## 3802.631 3830.489 -1896.316
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.093389 0.4919691
##
## Fixed effects: log(serBilir) ~ year + drug:year
## Value Std.Error DF t-value p-value
## (Intercept) 0.5707091 0.06421480 1631 8.887502 0.000
## year 0.0898390 0.00618248 1631 14.531219 0.000
## year:drugTreatment 0.0101796 0.00856947 1631 1.187895 0.235
## Correlation:
## (Intr) year
## year -0.104
## year:drugTreatment 0.001 -0.714
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.41229935 -0.53864027 -0.02809354 0.49249404 4.46134316
##
## Number of Observations: 1945
## Number of Groups: 312# Treatment-specific time effects, random intercepts + slopes
ml4 <- lme(log(serBilir) ~ year + drug:year, random = ~ year | id, data=pbc2)
summary(ml4)## Linear mixed-effects model fit by REML
## Data: pbc2
## AIC BIC logLik
## 3082.323 3121.323 -1534.161
##
## Random effects:
## Formula: ~year | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.9990381 (Intr)
## year 0.1722826 0.417
## Residual 0.3489259
##
## Fixed effects: log(serBilir) ~ year + drug:year
## Value Std.Error DF t-value p-value
## (Intercept) 0.4956686 0.05807539 1631 8.534915 0.0000
## year 0.1761726 0.01754759 1631 10.039704 0.0000
## year:drugTreatment 0.0027708 0.02411083 1631 0.114920 0.9085
## Correlation:
## (Intr) year
## year 0.177
## year:drugTreatment 0.002 -0.705
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.31683669 -0.49673826 -0.01978638 0.45207679 5.28538333
##
## Number of Observations: 1945
## Number of Groups: 312# 95% confidence intervals
intervals(ml4)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 0.38175835 0.495668558 0.60957876
## year 0.14175444 0.176172633 0.21059083
## year:drugTreatment -0.04452064 0.002770814 0.05006227
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd((Intercept)) 0.9204623 0.9990381 1.0843215
## sd(year) 0.1504235 0.1722826 0.1973183
## cor((Intercept),year) 0.2565228 0.4171400 0.5553377
##
## Within-group standard error:
## lower est. upper
## 0.3359261 0.3489259 0.3624289exp(intervals(ml4)$fixed)## lower est. upper
## (Intercept) 1.4648581 1.641595 1.839656
## year 1.1522937 1.192644 1.234407
## year:drugTreatment 0.9564559 1.002775 1.051337
## attr(,"label")
## [1] "Fixed effects:"# Likelihood ratio test for inclusion of random slopes
anova(ml4, ml3)## Model df AIC BIC logLik Test L.Ratio p-value
## ml4 1 7 3082.323 3121.323 -1534.161
## ml3 2 5 3802.631 3830.489 -1896.316 1 vs 2 724.3088 <.0001I prefer the model with treatment-specific time effects and random intercepts and slopes (ml4). At baseline, the mean bilirubin level is \(\exp(0.496)=1.642\) mg/dl, which is constrained to be equal between the treatment groups. There is substantial heterogeneity in bilirubin levels at baseline; compare \(\hat{\beta}_0=0.496\) to \(\hat{D}_{11}=0.999^2\). Note that \(\exp(0.496 \pm 1.96\times 0.999) = (0.232, 11.632)\) mg/dl quantifies variability across participants in bilirubin levels at baseline. NB: this is not a confidence interval for \(\beta_0\).
For the placebo group, a one-year increase in time is associated with a \(\exp(0.176)=1.193\)-factor increase in bilirubin levels. For the treatment group, a one-year increase in time is associated with a \(\exp(0.179)=1.196\)-factor increase in bilirubin levels. There is no evidence that the average trend in bilirubin levels differs between the placebo and treatment groups \((P = 0.91)\).
There is substantial heterogeneity in bilirubin trends over time; compare \(\hat{\beta}_1=0.176\) and \((\hat{\beta}_1 + \hat{\beta}_2)=0.179\) to \(\hat{D}_{22}=0.172^2\). Note that \(\exp(0.176 \pm 1.96\times 0.172) = ( 0.851, 1.672)\) and \(\exp(0.179 \pm 1.96\times 0.172) = (0.853, 1.676)\) quantify the variability across participants in the per-year increase in bilirubin levels among the placebo and treatment groups, respectively.
4 Survival data analysis
4.1 Standard Cox model
For fitting standard Cox regression models to the PBC data, a dataset is available that only includes a single row for each participant that contains the covariate values at baseline, survival time, and event indicator.
data(pbc2.id)
head(pbc2.id, 2)## id years status drug age sex year ascites hepatomegaly spiders edema
## 1 1 1.09517 dead D-penicil 58.76684 female 0 Yes Yes Yes edema despite diuretics
## 2 2 14.15234 alive D-penicil 56.44782 female 0 No Yes Yes No edema
## serBilir serChol albumin alkaline SGOT platelets prothrombin histologic status2
## 1 14.5 261 2.60 1718 138.0 190 12.2 4 1
## 2 1.1 302 4.14 7395 113.5 221 10.6 3 0First, using survminer, produce a Kaplan-Meier plot to display the survival probability over time for males and females. The figure includes confidence intervals for the survival probability, vertical dashes that indicate censoring events, and a table of the number of participants at risk below the figure. (A similar figure stratified by treatment group is less interesting.)
km <- survfit(Surv(years, status=='dead') ~ sex, data=pbc2.id)
ggsurvplot(km, data=pbc2.id, xlab='Time, years',
conf.int=TRUE, risk.table=TRUE, risk.table.y.text=FALSE,
ggtheme=theme(panel.grid.minor=element_blank()), tables.height=0.15,
tables.theme=theme(panel.border = element_blank(),
panel.grid.major=element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank()),
legend.title='Sex', legend.labs=c('Male', 'Female'),
palette=c('#F8766D', '#00BFC4'))
Next, use coxph to fit a series of Cox regression models, with and without adjustment for bilirubin level at baseline. Later, we will use ms1 for fitting joint regression models, which requires the x=TRUE argument. Use cox.zph to evaluate the proportional hazards assumption.
# Standard Cox model with treatment; note x=TRUE argument
ms1 <- coxph(Surv(years, status=='dead') ~ drug + sex + I(age-50), data=pbc2.id, x=TRUE)
summary(ms1)## Call:
## coxph(formula = Surv(years, status == "dead") ~ drug + sex +
## I(age - 50), data = pbc2.id, x = TRUE)
##
## n= 312, number of events= 140
##
## coef exp(coef) se(coef) z Pr(>|z|)
## drugD-penicil -0.146013 0.864146 0.172143 -0.848 0.3963
## sexfemale -0.470905 0.624437 0.221785 -2.123 0.0337 *
## I(age - 50) 0.042842 1.043773 0.008505 5.037 4.72e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## drugD-penicil 0.8641 1.1572 0.6167 1.2109
## sexfemale 0.6244 1.6014 0.4043 0.9644
## I(age - 50) 1.0438 0.9581 1.0265 1.0613
##
## Concordance= 0.629 (se = 0.024 )
## Likelihood ratio test= 33.25 on 3 df, p=3e-07
## Wald test = 34.87 on 3 df, p=1e-07
## Score (logrank) test = 35.31 on 3 df, p=1e-07# Standard Cox model with treatment and baseline bilirubin
ms2 <- coxph(Surv(years, status=='dead') ~ drug + sex + I(age-50) + log(serBilir), data=pbc2.id)
summary(ms2)## Call:
## coxph(formula = Surv(years, status == "dead") ~ drug + sex +
## I(age - 50) + log(serBilir), data = pbc2.id)
##
## n= 312, number of events= 140
##
## coef exp(coef) se(coef) z Pr(>|z|)
## drugD-penicil -0.127653 0.880159 0.176473 -0.723 0.469
## sexfemale -0.004180 0.995829 0.233030 -0.018 0.986
## I(age - 50) 0.046376 1.047468 0.008222 5.640 1.7e-08 ***
## log(serBilir) 1.085856 2.961975 0.092587 11.728 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## drugD-penicil 0.8802 1.1362 0.6228 1.244
## sexfemale 0.9958 1.0042 0.6307 1.572
## I(age - 50) 1.0475 0.9547 1.0307 1.064
## log(serBilir) 2.9620 0.3376 2.4704 3.551
##
## Concordance= 0.811 (se = 0.018 )
## Likelihood ratio test= 167.8 on 4 df, p=<2e-16
## Wald test = 174.7 on 4 df, p=<2e-16
## Score (logrank) test = 199.9 on 4 df, p=<2e-16# Check proportional hazards assumption
cox.zph(ms2)## chisq df p
## drug 0.0122 1 0.91
## sex 1.2882 1 0.26
## I(age - 50) 0.1800 1 0.67
## log(serBilir) 0.6665 1 0.41
## GLOBAL 2.4187 4 0.66plot(cox.zph(ms2)[4]) # Selects the fourth variable (bilirubin)
According to the model that includes baseline serum bilirubin (ms2):
- \(\exp(-0.128) = 0.88\) is the hazard ratio comparing the hazard of death between the treatment and placebo groups.
- \(\exp(-0.004) = 0.996\) is the hazard ratio comparing the hazard of death between females and males.
- \(\exp(0.046) = 1.047\) is the hazard ratio comparing the hazard of death for each one-year increase in age.
- \(\exp(1.086) = 2.962\) is the hazard ratio comparing the hazard of death for each one-unit increase in log bilirubin. Note that an additive difference on the log scale is a multiplicative difference on the orignal scale. Math fact: \(\exp(1)\approx2.718\).
There is no evidence for violation of the proportional hazards assumption. We fail to reject the null hypothesis of proportional hazards. The figure indicates approximately constant (log) hazard ratios over time for serum bilirubin.
4.2 Extended Cox model
Above, we plotted the observed individual-level profiles (or, trajectories) in bilirubin level over time since randomization, to explore how bilirubin level changed over time. To explore how bilirubin levels are associated with risk of death, plot the observed individual-level profiles in bilirubin level over time from the end of follow-up (subtract the event time years from the visit time year), stratified by follow-up status (died or censored), along with a smoothing line for each treatment group to display the average trend.
pbc2 <- within(pbc2, {
status.lab <- ifelse(status=='dead', 'Died', 'Censored')})
ggplot(pbc2, aes(year-years, serBilir, group=id)) +
facet_wrap( ~ status.lab) +
geom_line(aes(color=drug), alpha=0.4) +
geom_smooth(aes(year-years, serBilir, group=drug, color=drug), se=FALSE, size=1.5) +
scale_x_continuous(name='Time from end of follow-up, years', limits=c(-15,0)) +
scale_y_continuous(trans='log10', name='Serum bilirubin, mg/dl') +
scale_color_manual(name='Intervention', values=c('#00BFC4', '#F8766D')) +
mytheme()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
For fitting extended Cox models with serum bilirubin as a time-dependent covariate, a dataset is required that includes the start times and stop times for each time interval. This can be achieved using tmerge; see ?tmerge for details.
temp <- subset(pbc2.id, select=c(id, years, status, drug, sex, age, serBilir))
temp <- tmerge(temp, temp, id=id, event=event(years, status))
temp <- tmerge(temp, pbc2, id=id, serBilir=tdc(year, serBilir))## Warning in tmerge(temp, pbc2, id = id, serBilir = tdc(year, serBilir)): replacement of variable 'serBilir'head(temp, 11)## id years status drug sex age serBilir tstart tstop event
## 1 1 1.09517 dead D-penicil female 58.76684 14.5 0.0000000 0.5256817 alive
## 2 1 1.09517 dead D-penicil female 58.76684 21.3 0.5256817 1.0951703 dead
## 3 2 14.15234 alive D-penicil female 56.44782 1.1 0.0000000 0.4983025 alive
## 4 2 14.15234 alive D-penicil female 56.44782 0.8 0.4983025 0.9993429 alive
## 5 2 14.15234 alive D-penicil female 56.44782 1.0 0.9993429 2.1027270 alive
## 6 2 14.15234 alive D-penicil female 56.44782 1.9 2.1027270 4.9008871 alive
## 7 2 14.15234 alive D-penicil female 56.44782 2.6 4.9008871 5.8892783 alive
## 8 2 14.15234 alive D-penicil female 56.44782 3.6 5.8892783 6.8858833 alive
## 9 2 14.15234 alive D-penicil female 56.44782 4.2 6.8858833 7.8907020 alive
## 10 2 14.15234 alive D-penicil female 56.44782 3.6 7.8907020 8.8325485 alive
## 11 2 14.15234 alive D-penicil female 56.44782 4.6 8.8325485 14.1523382 alive# Extended Cox model with treatment and time-dependent serum bilirubin
ms3 <- coxph(Surv(tstart, tstop, event=='dead') ~ drug + sex + I(age-50) + log(serBilir), data=temp)
summary(ms3)## Call:
## coxph(formula = Surv(tstart, tstop, event == "dead") ~ drug +
## sex + I(age - 50) + log(serBilir), data = temp)
##
## n= 1945, number of events= 140
##
## coef exp(coef) se(coef) z Pr(>|z|)
## drugD-penicil -0.025800 0.974530 0.173150 -0.149 0.882
## sexfemale 0.179692 1.196849 0.235209 0.764 0.445
## I(age - 50) 0.068315 1.070702 0.008632 7.914 2.49e-15 ***
## log(serBilir) 1.464772 4.326556 0.094280 15.536 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## drugD-penicil 0.9745 1.0261 0.6941 1.368
## sexfemale 1.1968 0.8355 0.7548 1.898
## I(age - 50) 1.0707 0.9340 1.0527 1.089
## log(serBilir) 4.3266 0.2311 3.5966 5.205
##
## Concordance= 0.885 (se = 0.015 )
## Likelihood ratio test= 352.8 on 4 df, p=<2e-16
## Wald test = 258.8 on 4 df, p=<2e-16
## Score (logrank) test = 398.6 on 4 df, p=<2e-16Later, we will compare these results to those obtained from a joint regression model.
5 Joint regression models
5.1 Basic models
Joint regression models are fit using jointModel; see ?jointModel for details. The key inputs are: lmeObject, a linear mixed-effects model for the longitudinal marker; survObject, a standard Cox regression model for the survival outcome; timeVar, a time variable that corresponds to the time variable used in the mixed-effects model and that is on the same scale as the survival outcome; and a method for specifying the baseline hazard function, the parameterization of the relative risk model, and the procedure for numerical integration (weibull-PH-aGH is the default). Note that the underlying datasets to which the longitudinal and survival submodels are fit (e.g., pbc2 and pbc2.id, respectively) need to be in the same order (e.g., by id). See ?jointModelObject for a description of the objects returned by jointModel.
Using jointModel, fit a series of joint models with alternative specifications for the baseline hazard function: a parametric Weibull distribution, a piecewise constant step function, and regression splines. Recall that ml4 is the fitted mixed-effects model and ms1 is the fitted Cox regression model. NB: survObject is not the extended Cox model that includes serum bilirubin as a time-dependent covariate (e.g., ms3).
# Weibull hazard function
mj1 <- jointModel(lmeObject=ml4, survObject=ms1, timeVar='year', method='weibull-PH-aGH')
summary(mj1)##
## Call:
## jointModel(lmeObject = ml4, survObject = ms1, timeVar = "year",
## method = "weibull-PH-aGH")
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 1945 Number of Events: 140 (44.9%)
## Number of Groups: 312
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Weibull relative risk model
## Parameterization: Time-dependent
##
## log.Lik AIC BIC
## -1892.113 3810.227 3858.886
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.0018 (Intr)
## year 0.1793 0.4250
## Residual 0.3473
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 0.4920 0.0583 8.4422 <0.0001
## year 0.1815 0.0182 9.9620 <0.0001
## year:drugTreatment 0.0052 0.0243 0.2119 0.8322
##
## Event Process
## Value Std.Err z-value p-value
## (Intercept) -4.9792 0.3942 -12.6302 <0.0001
## drugD-penicil -0.0344 0.1848 -0.1860 0.8524
## sexfemale 0.1626 0.2500 0.6502 0.5156
## I(age - 50) 0.0650 0.0092 7.0399 <0.0001
## Assoct 1.3583 0.1014 13.3899 <0.0001
## log(shape) 0.1112 0.0800 1.3897 0.1646
##
## Scale: 1.1176
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0# Piecewise constant step function
mj2 <- jointModel(lmeObject=ml4, survObject=ms1, timeVar='year', method='piecewise-PH-aGH')
summary(mj2)##
## Call:
## jointModel(lmeObject = ml4, survObject = ms1, timeVar = "year",
## method = "piecewise-PH-aGH")
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 1945 Number of Events: 140 (44.9%)
## Number of Groups: 312
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Relative risk model with piecewise-constant
## baseline risk function
## Parameterization: Time-dependent
##
## log.Lik AIC BIC
## -1889.456 3814.911 3882.285
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.0009 (Intr)
## year 0.1797 0.4315
## Residual 0.3474
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 0.4919 0.0583 8.4433 <0.0001
## year 0.1824 0.0182 10.0036 <0.0001
## year:drugTreatment 0.0046 0.0244 0.1899 0.8494
##
## Event Process
## Value Std.Err z-value p-value
## drugD-penicil -0.0124 0.1863 -0.0664 0.9470
## sexfemale 0.1393 0.2508 0.5556 0.5785
## I(age - 50) 0.0653 0.0093 7.0020 <0.0001
## Assoct 1.3578 0.1025 13.2529 <0.0001
## log(xi.1) -4.9167 0.3780 -13.0078
## log(xi.2) -4.6649 0.3819 -12.2136
## log(xi.3) -4.8996 0.4149 -11.8095
## log(xi.4) -4.8511 0.4543 -10.6782
## log(xi.5) -4.4159 0.4219 -10.4662
## log(xi.6) -4.1187 0.4464 -9.2255
## log(xi.7) -4.6926 0.5691 -8.2454
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0# Regression splines
mj3 <- jointModel(lmeObject=ml4, survObject=ms1, timeVar='year', method='spline-PH-aGH')
summary(mj3)##
## Call:
## jointModel(lmeObject = ml4, survObject = ms1, timeVar = "year",
## method = "spline-PH-aGH")
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 1945 Number of Events: 140 (44.9%)
## Number of Groups: 312
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Relative risk model with spline-approximated
## baseline risk function
## Parameterization: Time-dependent
##
## log.Lik AIC BIC
## -1888.156 3816.312 3891.172
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.0003 (Intr)
## year 0.1800 0.4337
## Residual 0.3474
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 0.4920 0.0583 8.4440 <0.0001
## year 0.1828 0.0183 10.0033 <0.0001
## year:drugTreatment 0.0048 0.0245 0.1974 0.8435
##
## Event Process
## Value Std.Err z-value p-value
## drugD-penicil 0.0036 0.1861 0.0196 0.9844
## sexfemale 0.1403 0.2519 0.5569 0.5776
## I(age - 50) 0.0644 0.0093 6.9242 <0.0001
## Assoct 1.3543 0.1028 13.1786 <0.0001
## bs1 -4.6216 0.5783 -7.9915 <0.0001
## bs2 -5.2576 0.6700 -7.8475 <0.0001
## bs3 -4.5098 0.6822 -6.6105 <0.0001
## bs4 -4.9134 0.5843 -8.4092 <0.0001
## bs5 -4.7995 0.6199 -7.7426 <0.0001
## bs6 -4.6132 0.6620 -6.9688 <0.0001
## bs7 -2.8239 1.1676 -2.4185 0.0156
## bs8 -7.4645 2.1586 -3.4580 0.0005
## bs9 -3.8692 2.0582 -1.8799 0.0601
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0The estimates under Variance Components correspond to the variance components of the linear mixed-effects model for (log) serum bilirubin (i.e., \(D\) and \(\sigma\)). The coefficients under Longitudinal Process correspond to those in the linear-mixed effects model for (log) serum bilirubin (i.e., \(\beta\)). The coefficients under Event Process correspond to those in the standard Cox regression model for survival (i.e., \(\gamma\)), plus the coefficient that measures the association between the value of the subject-specific evolution (i.e., Assoct) and the event hazard (i.e., \(\alpha\)), and any parameters for the baseline hazard function.
According to the AIC, mj1 (Weibull relative risk model) provides a comparatively better fit to the data. According to this model, there is no difference in the average trend in bilirubin levels between the placebo and treatment groups \((P=0.83)\). There is a strong association between bilirubin levels and mortality: hazard ratio, \(\exp(1.358) = 3.889\). The lecture notes provide tables that compare the coefficients and standard errors obtained from separate and joint analyses.
Although likelihood-based methods require several assumptions (e.g., the likelihood function is correctly specified) and can be computationally demanding, the availability of a likelihood permits likelihood-based hypothesis tests of nested models. With joint regression models, we can conduct a hypothesis test for a combination of parameters in the longitudinal and/or survival submodels. For example, we can test whether there was any effect of treatment on either serum bilirubin or mortality by evaluating the null hypothesis \(H: \beta_2=\gamma_1=0\). This involves fitting the ‘null’ model (without drug in either the longitudinal or suvival submodels) and conducting a likelihood ratio test.
ml0 <- lme(log(serBilir) ~ year, random = ~ year | id, data=pbc2)
ms0 <- coxph(Surv(years, status=='dead') ~ sex + age, data=pbc2.id, x=TRUE)
mj0 <- jointModel(lmeObject=ml0, survObject=ms0, timeVar='year', method='weibull-PH-aGH')
anova(mj0, mj1)##
## AIC BIC log.Lik LRT df p.value
## mj0 3806.32 3847.49 -1892.16
## mj1 3810.23 3858.89 -1892.11 0.09 2 0.9558There is no evidence of a treatment effect on either serum bilirubin or mortality.
5.2 Model diagnostics
Residual analysis can be used to evaluate the adequacy of a fitted joint regression model. See ?residuals.jointModel and ?fitted.jointModel for details on extracting the residuals and fitted values, respectively, from a jointModel object. In particular, there are various options for extracting residuals and fitted values for the longituidinal and survival submodels.
plotResid <- function(x, y, col.lowess='red', ...){
plot(x, y, ...)
abline(h=0, col='gray', lwd=2)
lines(lowess(x, y), col=col.lowess, lwd=2)}
plotQQ <- function(x, ...){
qqnorm(x, main='')
qqline(x, col='gray', lwd=2)}
r_l_s <- residuals(mj1, process='Longitudinal', type='stand-Subject')
f_l_s <- fitted(mj1, process='Longitudinal', type='Subject')
r_l_m <- residuals(mj1, process='Longitudinal', type='stand-Marginal')
f_l_m <- fitted(mj1, process='Longitudinal', type='Marginal')
r_s <- residuals(mj1, process='Event', type='Martingale')
f_s <- fitted(mj1, process='Longitudinal', type='EventTime')
par(mfrow=c(2,2), mar=c(5,5,1,1))
plotResid(f_l_s, r_l_s, xlab='Fitted values (subject)', ylab='Residuals (subject)')
plotQQ(r_l_s)
plotResid(f_l_m, r_l_m, xlab='Fitted values (marginal)', ylab='Residuals (marginal)')
plotResid(f_s, r_s, xlab='Fitted values (longitudinal)', ylab='Residuals (Martingale)')
(Longitudinal) The plot of subject-specific residuals versus fitted values (upper left panel) indicates no systematic departure from constant variance. The quantile-quantile plot of subject-specific residuals (upper right panel) could indicate a departure from normality, but this is often difficult to judge; it’s certainly better than if we had not log-transformed the outcome. The plot of marginal residuals versus fitted values (lower left panel) shows a clear trend. This trend could indicate that the mean model is not correctly specified. However, serum bilirubin is positively associated with mortality, such that lack of fit could be due to the impact of missing data. This motivates the use of multiple imputation residuals (§6.3), in which the observed data are augmented with randomly imputed longitudinal outcomes under the complete-data model, corresponding to the longitudinal outcomes that would have been observed had the participant not died. The code to accomplish this is a bit dense, so is not displayed here.
The smoother based on the complete data (blue line) does not exhibit the systematic trend that was clearly exhibited by the smoother based on the observed data (red line). Therefore, the specification for the mean model appears to be adequate.
(Survival) The plot of Martingale residuals versus fitted values (lower right panel of the previous figure) does not indicate a systematic departure from linearity.
5.3 Advanced models
Use the interFact argument to specify an interaction between serum bilirubin and treatment in the survival submodel.
mj4 <- jointModel(lmeObject=ml4, survObject=ms1, timeVar='year', method='weibull-PH-aGH',
interFact=list(value = ~ drug, data=pbc2.id))
summary(mj4)##
## Call:
## jointModel(lmeObject = ml4, survObject = ms1, timeVar = "year",
## method = "weibull-PH-aGH", interFact = list(value = ~drug,
## data = pbc2.id))
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 1945 Number of Events: 140 (44.9%)
## Number of Groups: 312
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Weibull relative risk model
## Parameterization: Time-dependent
##
## log.Lik AIC BIC
## -1891.945 3811.89 3864.292
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.0016 (Intr)
## year 0.1794 0.4269
## Residual 0.3473
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 0.4917 0.0583 8.4387 <0.0001
## year 0.1814 0.0182 9.9515 <0.0001
## year:drugTreatment 0.0055 0.0243 0.2263 0.8210
##
## Event Process
## Value Std.Err z-value p-value
## (Intercept) -4.8725 0.4270 -11.4099 <0.0001
## drugD-penicil -0.2733 0.4412 -0.6193 0.5357
## sexfemale 0.1836 0.2525 0.7271 0.4672
## I(age - 50) 0.0660 0.0094 7.0237 <0.0001
## Assoct 1.3005 0.1381 9.4160 <0.0001
## Assoct:drugD-penicil 0.1186 0.1996 0.5939 0.5526
## log(shape) 0.1076 0.0803 1.3399 0.1803
##
## Scale: 1.1136
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0There is no evidence that the association between serum bilirubin and mortality differs between treatment groups.
So far, we have assumed that the association between serum bilirubin and mortality is fully captured by the value of the subject-specific evolution, which changes over time. (This is what is denoted by Parameterization: Time-dependent in the joint model output.) However, we can also allow the association between serum bilirubin and mortality to depend on the slope (or rate of change) of the subject-specific evolution, which could also change over time depending on how we parameterize the model.
Consider a model with time variable \(t\), group variable \(g\), and \(m_i(t) = \beta_0 + \beta_1 t + \beta_2 g_i\times t + b_{i0} + b_{i1}t\); this model corresponds to ml4, which is linear in time. Then \(m'_i(t) = \beta_1 + \beta_2 g_i + b_{i1}\). Math fact: if \(f(x)=ax+b\), then \(f'(x)=a\). To fit the model, specify parameterization=“both” and provide a derivForm argument, which includes formulas for the derivatives of the fixed-effects and random-effects components with respect to time, as well as numeric vectors to indicate which variables in lmeObject correspond to the derivatives.
# Level and slope models for serum bilirubin in survival submodel
# Linear time effects
dFormL <- list(fixed = ~ 1 + drug, indFixed=c(2,3), random = ~ 1, indRandom=2)
mj5 <- jointModel(lmeObject=ml4, survObject=ms1, timeVar='year', method='weibull-PH-aGH',
parameterization='both', derivForm=dFormL)
summary(mj5)##
## Call:
## jointModel(lmeObject = ml4, survObject = ms1, timeVar = "year",
## parameterization = "both", method = "weibull-PH-aGH", derivForm = dFormL)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 1945 Number of Events: 140 (44.9%)
## Number of Groups: 312
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Weibull relative risk model
## Parameterization: Time-dependent + time-dependent slope
##
## log.Lik AIC BIC
## -1889.902 3807.804 3860.206
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.0001 (Intr)
## year 0.1840 0.4612
## Residual 0.3471
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 0.4901 0.0580 8.4557 <0.0001
## year 0.1878 0.0186 10.0974 <0.0001
## year:drugTreatment 0.0053 0.0244 0.2163 0.8288
##
## Event Process
## Value Std.Err z-value p-value
## (Intercept) -5.3856 0.4743 -11.3552 <0.0001
## drugD-penicil -0.0458 0.1919 -0.2388 0.8112
## sexfemale 0.1547 0.2567 0.6025 0.5469
## I(age - 50) 0.0647 0.0096 6.7511 <0.0001
## Assoct 1.2308 0.1246 9.8805 <0.0001
## Assoct.s 1.8741 0.9484 1.9760 0.0482
## log(shape) 0.1814 0.0906 2.0019 0.0453
##
## Scale: 1.1988
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0Note that the model output now indicates Parameterization: Time-dependent + time-dependent slope. According to this model, there is strong evidence that both the value and slope of serum bilirubin are associated with mortality, as indicated by the coefficients for Assoct and Assoct.s, respectively. However, because the longitudinal submodel is linear in time, its slope is constant (i.e., does not change over time). We can obtain a non-constant slope by fitting a longitudinal submodel that is non-linear in time. For example, consider a linear-quadratic function for time. The model specification is a bit messy, but doable thanks to calculus. Math fact: if \(f(x)=ax^2 + bx + c\), then \(f'(x)=2ax + b\).
# Quadratic time effects
ml6 <- lme(log(serBilir) ~ drug*(year + I(year^2)), random = ~ year + I(year^2) | id, data=pbc2)
dFormQ <- list(fixed = ~ 1 + I(2*year) + drug + I(2*year):drug, indFixed=3:6,
random = ~ 1 + I(2*year), indRandom=2:3)
mj6 <- jointModel(lmeObject=ml6, survObject=ms1, timeVar='year', method='weibull-PH-aGH',
parameterization='both', derivForm=dFormQ)
summary(mj6)##
## Call:
## jointModel(lmeObject = ml6, survObject = ms1, timeVar = "year",
## parameterization = "both", method = "weibull-PH-aGH", derivForm = dFormQ)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 1945 Number of Events: 140 (44.9%)
## Number of Groups: 312
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Weibull relative risk model
## Parameterization: Time-dependent + time-dependent slope
##
## log.Lik AIC BIC
## -1777.869 3595.737 3670.597
##
## Variance Components:
## StdDev Corr
## (Intercept) 0.9952 (Intr) year
## year 0.3184 0.2607
## I(year^2) 0.0255 -0.0430 -0.8613
## Residual 0.3023
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 0.5833 0.0716 8.1419 <0.0001
## drugTreatment -0.1510 0.0977 -1.5459 0.1221
## year 0.1764 0.0315 5.6075 <0.0001
## I(year^2) 0.0022 0.0033 0.6492 0.5162
## drugTreatment:year 0.0018 0.0421 0.0429 0.9658
## drugTreatment:I(year^2) -0.0016 0.0042 -0.3765 0.7065
##
## Event Process
## Value Std.Err z-value p-value
## (Intercept) -6.3717 0.5640 -11.2973 <0.0001
## drugD-penicil -0.0198 0.2008 -0.0985 0.9215
## sexfemale 0.2309 0.2656 0.8692 0.3848
## I(age - 50) 0.0682 0.0101 6.7461 <0.0001
## Assoct 1.4068 0.1267 11.1018 <0.0001
## Assoct.s 2.6832 0.6489 4.1349 <0.0001
## log(shape) 0.2758 0.0924 2.9858 0.0028
##
## Scale: 1.3176
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0anova(mj5, mj6)##
## AIC BIC log.Lik LRT df p.value
## mj5 3807.80 3860.21 -1889.90
## mj6 3595.74 3670.60 -1777.87 224.07 6 <0.0001The quadratic terms substantially improve model fit, as indicated by AIC and the likelihood ratio test. There remains strong evidence that both the value and slope of serum bilirubin are associated with mortality.
NB: fitting more complex models requires more computation time. In addition, fitting models to larger datasets can be computatinally demanding.
5.4 Prediction
The development and evaluation of risk-prediction models is a complex topic. (In fact, there are other modules devoted to the topic at SISCER.) Below are some commands to obtain and display dynamic predictions of conditional survival from a joint regression model. See ?survfitJM and ?plot.survfitJM for details.
ND <- subset(pbc2, id==15)
survPreds <- vector('list', nrow(ND))
for(i in 1:nrow(ND)){
set.seed(4)
survPreds[[i]] <- survfitJM(mj6, newdata=ND[1:i,])}
par(mfrow=c(1,3))
for(i in c(3,6,9)){
plot(survPreds[[i]], estimator='median', conf.int=TRUE, include.y=TRUE, pch=19, lwd=2, col='blue',
xlab='', ylab='', ylab2='', main=paste(round(survPreds[[i]]$last.time, 1), 'years'))}
mtext('Time, years', side=1, line=-2, outer=TRUE)
mtext('log serum bilirubin', side=2, line=-1, outer=TRUE)
mtext('Survival probability', side=4, line=-1, outer=TRUE)
At year 1, serum bilirubin (value) remains low, with almost no change (slope); predicted survival is relatively favorable for this participant. However, at year 4.1, serum bilirubin begins to increase, with a corresponding decrease in predicted survival. At year 6.8, serum bilirubin increases further; predicted survival is relatively poor. The participant died at year 9.8.
This file was created from:
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## attached base packages:
## [1] splines stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.4 survminer_0.4.7 ggpubr_0.4.0 ggplot2_3.3.2 JM_1.4-8 survival_3.2-7
## [7] nlme_3.1-152 MASS_7.3-53
##
## loaded via a namespace (and not attached):
## [1] zoo_1.8-8 tidyselect_1.1.0 xfun_0.15 purrr_0.3.4 haven_2.3.1
## [6] lattice_0.20-41 carData_3.0-4 colorspace_1.4-1 vctrs_0.3.6 generics_0.0.2
## [11] htmltools_0.5.0 mgcv_1.8-33 yaml_2.2.1 survMisc_0.5.5 rlang_0.4.10
## [16] pillar_1.4.4 foreign_0.8-81 glue_1.4.1 withr_2.4.1 DBI_1.1.1
## [21] readxl_1.3.1 lifecycle_0.2.0 stringr_1.4.0 cellranger_1.1.0 munsell_0.5.0
## [26] ggsignif_0.6.0 gtable_0.3.0 zip_2.1.1 evaluate_0.14 labeling_0.3
## [31] knitr_1.29 rio_0.5.16 forcats_0.5.0 curl_4.3 broom_0.5.6
## [36] Rcpp_1.0.4.6 xtable_1.8-4 scales_1.1.1 backports_1.1.8 abind_1.4-5
## [41] farver_2.0.3 km.ci_0.5-2 gridExtra_2.3 hms_0.5.3 digest_0.6.25
## [46] stringi_1.4.6 openxlsx_4.1.5 rstatix_0.6.0 KMsurv_0.1-5 grid_4.0.4
## [51] tools_4.0.4 magrittr_1.5 tibble_3.0.1 crayon_1.3.4 tidyr_1.1.0
## [56] car_3.0-8 pkgconfig_2.0.3 ellipsis_0.3.1 Matrix_1.3-2 data.table_1.13.6
## [61] assertthat_0.2.1 rmarkdown_2.6 R6_2.4.1 compiler_4.0.4