In drizopoulos/JMbayes2: Extended Joint Models for Longitudinal and Time-to-Event Data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library("JMbayes2")

Functional Forms

Simplified syntax

We have previously seen that function jm() via its functional_forms argument allows the specification of different functional forms to link the longitudinal and event time outcomes. This argument accepts either a single formula or a list of formulas per longitudinal outcome with the terms we wish to include.

We will illustrate some of these possibilities using the PBC dataset. We start by fitting a Cox model for the composite event transplantation or death, including sex as a baseline covariate:

pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive')
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)

Our aim is to assess the strength of the association between the risk of the composite event and whether the patients experienced hepatomegaly during follow-up. We will describe the patient-specific profiles over time for this biomarker using a mixed-effects logistic model, where both in the fixed- and random-effects we include an intercept and the time effect. The syntax to fit this model with mixed_model() is:

fm <- mixed_model(hepatomegaly ~ year, data = pbc2, random = ~ year | id, 
                  family = binomial())

The default call to jm() adds the subject-specific linear predictor of the mixed-effects logistic regression as a time-varying covariate in the survival relative risk model:

jointFit1 <- jm(CoxFit, fm, time_var = "year")
summary(jointFit1)

In the output this is named as value(hepatomegaly) to denote that the current value functional form is used. That is, we assume that the risk at a specific time $t$ is associated with the value of the linear predictor of the longitudinal outcome at the same time point $t$. In this particular case, the subject-specific linear predictor denotes the log odds of experiencing hepatomegaly at time $t$.

Transformation functions

The fact that the default version of the current value functional form is on the linear predictor scale of the mixed model may be problematic to interpret when this linear predictor is connected with a nonlinear link function to mean of the longitudinal outcome. In these situations we may want to transform the subject-specific linear predictor back to the scale of the outcome. To achieve this we can use a transformation function. Continuing on the previous example, we update jointFit1 by now linking the expit transformation of the linear predictor (i.e., $\mbox{expit}(x) = \exp(x) / {1 + \exp(x)}$) with the risk of an event. This is done using the vexpit() function:

jointFit2 <- update(jointFit1, functional_forms = ~ vexpit(value(hepatomegaly)))
summary(jointFit2)

Other available functions to use in the definition of the functional_forms argument are vexp() for calculating the exponent, vlog() for calculating the natural logarithm, and vsqrt() for calculating the square root.

If we want to include the time-varying slope of the transformed linear predictor, we also have the Dexpit() and Dexp() functions available. As an example, we extend jointFit2 by including the derivative of the $\mbox{expit}()$ transformation:

forms <- ~ vexpit(value(hepatomegaly)) + Dexpit(slope(hepatomegaly))
jointFit3 <- update(jointFit1, functional_forms = forms)
summary(jointFit3)

The call to Dexpit(slope(hepatomegaly)) is internally transformed to Dexpit(value(hepatomegaly)):slope(hepatomegaly), which calculates the derivative of the $\mbox{expit}()$ evaluated at the linear predictor times the derivative of the linear predictor. This is because $$\frac{d}{dt} \mbox{expit}{\eta(t)} = \mbox{expit}{\eta(t)} \, \Bigl [ 1 - \mbox{expit}{\eta(t)} \Bigr ] \times \frac{d}{dt}\eta(t)$$

The Slope functional form

As we have seen in previous examples, the slope() function is used to specify the slope functional form $d \eta(t)/dt$. According to the definition of the derivative, this corresponds to the change in the longitudinal profile ${\eta(t + \varepsilon) - \eta(t)}/ \varepsilon$ as $\varepsilon$ approaches zero. However, the interpretation of this term may be challenging in some settings. A possible alternative would be to increase the value of $\varepsilon$, e.g., $\varepsilon = 1$. For example, if the time scale is years, this would quantify the change of the longitudinal profile in the last year before $t$.

To fit a joint model with such a term, we can use the eps and direction arguments of the slope() function. We illustrate this in the following example, in which we use the serum bilirubin

gm <- lme(log(serBilir) ~ ns(year, 2), data = pbc2, random = ~ ns(year, 2) | id,
          control = lmeControl(opt = "optim"))

We first fit the joint model with time-varying slope term:

jFit1 <- jm(CoxFit, gm, time_var = "year",
            functional_forms = ~ value(log(serBilir)) + slope(log(serBilir)))

To specify that we want to include the change in the log serum bilirubin levels during the last year before $t$, we specify eps = 1 and direction = "back" in the call to slope(). This calculates the term ${\eta(t) - \eta(t - \varepsilon)} / \varepsilon$ for $\varepsilon$ set equal to eps = 1:

jFit2 <- jm(CoxFit, gm, time_var = "year",
            functional_forms = ~ value(log(serBilir)) + 
              slope(log(serBilir), eps = 1, direction = "back"))

We compare the two fits

summary(jFit1)$Survival

summary(jFit2)$Survival

drizopoulos/JMbayes2 documentation built on April 15, 2024, 3:44 a.m.