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       0

First, 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.66

plot(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 original 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-16

Later, we will compare these results to those obtained from a joint regression model.

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.455 3814.91 3882.284
## 
## 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.4429 <0.0001
## year               0.1824  0.0182 10.0038 <0.0001
## year:drugTreatment 0.0046  0.0244  0.1897  0.8495
## 
## Event Process
##                 Value Std.Err  z-value p-value
## drugD-penicil -0.0119  0.1863  -0.0637  0.9492
## sexfemale      0.1416  0.2509   0.5645  0.5724
## I(age - 50)    0.0653  0.0093   7.0042 <0.0001
## Assoct         1.3584  0.1025  13.2531 <0.0001
## log(xi.1)     -4.9206  0.3783 -13.0087        
## log(xi.2)     -4.6684  0.3822 -12.2152        
## log(xi.3)     -4.9028  0.4151 -11.8115        
## log(xi.4)     -4.8548  0.4545 -10.6812        
## log(xi.5)     -4.4194  0.4221 -10.4692        
## log(xi.6)     -4.1235  0.4468  -9.2296        
## log(xi.7)     -4.6982  0.5696  -8.2485        
## 
## 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.086 3816.172 3891.032
## 
## 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.4918  0.0583  8.4416 <0.0001
## year               0.1829  0.0183 10.0060 <0.0001
## year:drugTreatment 0.0048  0.0245  0.1957  0.8448
## 
## Event Process
##                 Value Std.Err z-value p-value
## drugD-penicil  0.0036  0.1861  0.0195  0.9845
## sexfemale      0.1395  0.2520  0.5534  0.5800
## I(age - 50)    0.0642  0.0093  6.9076 <0.0001
## Assoct         1.3538  0.1028 13.1655 <0.0001
## bs1           -4.6126  0.5778 -7.9835 <0.0001
## bs2           -5.2674  0.6700 -7.8618 <0.0001
## bs3           -4.4891  0.6822 -6.5807 <0.0001
## bs4           -4.9337  0.5847 -8.4388 <0.0001
## bs5           -4.7577  0.6195 -7.6796 <0.0001
## bs6           -4.7085  0.6658 -7.0717 <0.0001
## bs7           -2.5598  1.1680 -2.1916  0.0284
## bs8           -8.1073  2.2146 -3.6609  0.0003
## bs9           -3.3354  1.8687 -1.7849  0.0743
## 
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5 
## 
## Optimization:
## Convergence: 0

The 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 survival 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.9558

There 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 longitudinal 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: 0

There 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: 0

Note 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.1420 <0.0001
## drugTreatment           -0.1510  0.0977 -1.5459  0.1221
## year                     0.1764  0.0315  5.6074 <0.0001
## I(year^2)                0.0022  0.0033  0.6492  0.5162
## drugTreatment:year       0.0018  0.0421  0.0429  0.9657
## drugTreatment:I(year^2) -0.0016  0.0042 -0.3766  0.7065
## 
## Event Process
##                 Value Std.Err  z-value p-value
## (Intercept)   -6.3716  0.5640 -11.2974 <0.0001
## drugD-penicil -0.0198  0.2008  -0.0987  0.9214
## sexfemale      0.2308  0.2656   0.8691  0.3848
## I(age - 50)    0.0682  0.0101   6.7461 <0.0001
## Assoct         1.4068  0.1267  11.1019 <0.0001
## Assoct.s       2.6830  0.6489   4.1347 <0.0001
## log(shape)     0.2758  0.0924   2.9859  0.0028
## 
## Scale: 1.3176 
## 
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5 
## 
## Optimization:
## Convergence: 0

anova(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.0001

The 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 computationally 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.