Population-level marginal summaries from smooth two-slope mixed models with random changepoints

Divan A. Burger, Sean van der Merwe, Emmanuel Lesaffre

Introduction

In general,

\[E[f(X)]\ne f(E[X])\]

and

\[\hat{f}(\boldsymbol\theta)\ne f(\hat{\boldsymbol\theta})\]

yet statisticians keep assuming equality.

Source and authors

This presentation is drawn from our paper that is under review at the Journal of Biopharmaceutical Statistics.

Divan A. Burger\(^{1,2}\), Sean van der Merwe\(^{2}\), Emmanuel Lesaffre\(^{3,4}\)

  1. Cytel Inc., Waltham, MA, USA
  2. Department of Mathematical Statistics and Actuarial Science, University of the Free State, Bloemfontein, South Africa
  3. I-BioStat, KU Leuven, Leuven, Belgium
  4. Department of Statistics and Actuarial Science, University of Stellenbosch, Stellenbosch, South Africa

Outline

  • Motivation
  • Dataset to be used as example
  • Modelling aspects
  • More modelling aspects
  • Conclusions

Motivation

  • We want people to get better from illness
  • Understand the process of getting better
  • Understand the factors that affect the process
  • Understand the statistical uncertainty and variation
    • Better modelling of the variation \(\longrightarrow\) better modelling of the process

Spoiler

Consider the rate at which people get better from HIV when taking ARVs

The expected rate of healing for a random future person is not the same as for the typical person (doesn’t exist)

Data set

The ACTG 315 dataset, available in the ushr R package (Morris et al. 2020), includes longitudinal measurements of HIV viral load (log\(_{10}\) RNA copies/mL) over time. It features data on 46 patients, with the longest measurement recorded on Day 196 after baseline (Day 0).

Data set

RNA by time and subject

Censoring

We observe two kinds of censoring:

  1. Subjects stop taking part in the study at different lengths of time
  2. A proportion of observations censored below the lower detection limit

We address these in different ways:

  1. By storing the data in long form and fitting a mixed effects model (each patient has their own curve)
  2. By putting censored observations into the likelihood via their CDF (not PDF), i.e. \(P[X <= LDL]\) instead of \(P[X = x]\); e.g. Lachos et al. (2011)

The modelling

We build up the model step by step.

Mixed effects models

  • Mixed effects models include fixed effects and random effects (e.g. random intercepts or random slopes per subject)
  • Mixed effects models are used when we have sampling of nominal observation groups from a population of possible groups
  • Typically a random sample of people from a population of people, but with multiple observations per person
    • Can also be a random sample of animals, random sample of fields, random sample of classes or schools

Random effects are important to help address the fact that our residuals are not independent.

Non-linear regression

  • Looking at the log RNA plot we see a non-linear pattern.
    • First a fast curve down
    • then a slow curve down/flat/up
    • with a smooth transition

Each person has their own curve

Each person has an intercept, an initial slope, a turning point, an ultimate slope, and their own turning smoothness. Some turn gradually while others are more kinky (literally, not figuratively).

Examples of papers implementing this idea in various fields include: Perelson et al. (1996), Heerspink et al. (2021), Bacon and Watts (1971), Donald and Diacon (2008), Burger and Schall (2015), Hall et al. (2000), Ghosh and Vaida (2007), van den Hout et al. (2011), etc.

Non-linear curve illustration

While we did formally use the DHARMa approach in the paper to check the goodness-of-fit for the two models (Hartig 2021), I want to encourage all researchers to try to plot the model and data together if at all possible (sometimes some creativity and some conditioning is needed).

Model fit procedure

  • The tool we used for fitting this particular model is JAGS
    • JAGS stands for Just Another Gibbs Sampler
  • The fitting was done using the R package runjags (Denwood 2016)
  • JAGS uses clever Gibbs sampling to arrive at posterior simulations for all the parameters
  • Once you have posterior simulations you can calculate anything you want
    • Fits and predictions for a particular subject
    • Fits for the average subject
    • Predictions for a random future subject, with full uncertainty

Model code

model {
  for (r in 1:3) {
    L_chol[r, r] ~ dt(0, 0.5, 3)T(0, )
    for (c in 1:(r - 1)) {
      L_chol[r, c] ~ dnorm(0, 0.001)
    }
    for (c in (r + 1):3) {
      L_chol[r, c] <- 0
    }
  }

  for (r in 1:3) {
    for (c in 1:3) {
      Prec_RE[r, c] <- inprod(L_chol[r, 1:3], L_chol[c, 1:3])
    }
  }

  for (j in 1:3) {
    Beta_RE_bar[j] ~ dnorm(0, 0.001)
  }

  Sigma_resid ~ dt(0, 0.5, 3)T(0, )
  Tau_resid <- pow(Sigma_resid, -2)

  Lambda ~ dunif(-1, 1)
  Eps ~ dexp(1)
  NuMinus2 ~ dgamma(2, Eps)
  Nu <- NuMinus2 + 2

  LogC <- loggam((Nu + 1)/2) - 0.5*log(pi*(Nu - 2)) - loggam(Nu/2)
  Cc <- exp(LogC)
  Aconst <- 4*Lambda*Cc*(Nu - 2)/(Nu - 1)
  B2 <- 1 + 3*pow(Lambda, 2) - pow(Aconst, 2)
  Bconst <- sqrt(B2)
  Kfac <- sqrt(Nu/(Nu - 2))
  Knot <- -Aconst/Bconst

  Mu_tau ~ dnorm(0, 0.001)
  Sigma_tau ~ dt(0, 0.5, 3)T(0, )
  Tau_tau <- pow(Sigma_tau, -2)

  Mu_kappa ~ dnorm(0, 0.001)
  Sigma_kappa ~ dt(0, 0.5, 3)T(0, )
  Tau_kappa <- pow(Sigma_kappa, -2)

  Phi0 <- 1/sqrt(2*pi)
  KappaSwitch <- 1.0E-4

  for (i in 1:N) {
    v[i, 1:3] ~ dmnorm(Beta_RE_bar[1:3], Prec_RE[,])
    Alpha0[i] <- v[i, 1]
    Beta1[i] <- v[i, 2]
    Beta2[i] <- v[i, 3]

    Tau[i] ~ dnorm(Mu_tau, Tau_tau)T(0, TauMax)
    Kappa[i] ~ dnorm(Mu_kappa, Tau_kappa)T(0, )
  }

  for (n in 1:Ntotal) {
    DeltaT[n] <- Month[n] - Tau[Id[n]]
    Delta0[n] <- -Tau[Id[n]]
    
    Zk[n] <- Kappa[Id[n]]*DeltaT[n]
    Zk0[n] <- Kappa[Id[n]]*Delta0[n]

    PhiZ[n] <- pnorm(Zk[n], 0, 1)
    PhiZ0[n] <- pnorm(Zk0[n], 0, 1)
    phiZ[n] <- dnorm(Zk[n], 0, 1)
    phiZ0[n] <- dnorm(Zk0[n], 0, 1)

    HingeExact[n] <- DeltaT[n]*PhiZ[n] + phiZ[n]/Kappa[Id[n]]
    Hinge0Exact[n] <- Delta0[n]*PhiZ0[n] + phiZ0[n]/Kappa[Id[n]]
    CenteredHingeExact[n] <- HingeExact[n] - Hinge0Exact[n]

    CenteredHingeApprox[n] <- 0.5*Month[n] +
      0.5*Phi0*Kappa[Id[n]]*(pow(DeltaT[n], 2) - pow(Delta0[n], 2))

    UseApprox[n] <- step(KappaSwitch - Kappa[Id[n]])
    CenteredHinge[n] <- UseApprox[n]*CenteredHingeApprox[n] +
      (1 - UseApprox[n])*CenteredHingeExact[n]

    Mu[n] <- Alpha0[Id[n]] + Beta1[Id[n]]*Month[n] +
      Beta2[Id[n]]*CenteredHinge[n]

    Z[n] <- (Y[n] - Mu[n])/Sigma_resid
    Zup[n] <- (Upper[n] - Mu[n])/Sigma_resid
    Ztilde[n] <- Bconst*Z[n] + Aconst
    ZtildeUp[n] <- Bconst*Zup[n] + Aconst

    Left[n] <- step(Knot - Z[n])
    Right[n] <- 1 - Left[n]
    BaseL[n] <- 1 + pow(Ztilde[n]/(1 - Lambda), 2)/(Nu - 2)
    BaseR[n] <- 1 + pow(Ztilde[n]/(1 + Lambda), 2)/(Nu - 2)
    PartL[n] <- Left[n]*(Bconst*Cc)*pow(BaseL[n], -(Nu + 1)/2)
    PartR[n] <- Right[n]*(Bconst*Cc)*pow(BaseR[n], -(Nu + 1)/2)
    PDFstd[n] <- PartL[n] + PartR[n]
    PDF[n] <- PDFstd[n]/Sigma_resid

    LeftUp[n] <- step(Knot - Zup[n])
    RightUp[n] <- 1 - LeftUp[n]
    tL_up[n] <- (ZtildeUp[n]/(1 - Lambda))*Kfac
    tR_up[n] <- (-ZtildeUp[n]/(1 + Lambda))*Kfac
    CDF_left_up[n] <- (1 - Lambda)*pt(tL_up[n], 0, 1, Nu)
    CDF_right_up[n] <- 1 - (1 + Lambda)*pt(tR_up[n], 0, 1, Nu)
    CDFstd_up[n] <- LeftUp[n]*CDF_left_up[n] + RightUp[n]*CDF_right_up[n]
    CDF[n] <- CDFstd_up[n]

    Like[n] <- pow(PDF[n], Censor[n])*pow(CDF[n], 1 - Censor[n])
    p[n] <- Like[n]/C
    Ones[n] ~ dbern(p[n])
  }
}

The two curve options considered

We consider two equations for \(\mu\):

  • The classic bi-linear model where two lines meet at a point
    • Note that this can also end up looking like a smooth transition if we allow for parameter uncertainty, subject uncertainty, or both (the main point of the paper)
  • A new model form which introduces a transition from one line to the other with flexible smoothness

Two-line (bi-phasic) curve illustration

  • A lot of real-world problems involve two straight-line trajectories meeting at a changepoint
    • HIV, some kidney issues, and TB are cited in the paper
    • I have seen it in farming (herbicide, insecticide, etc.)
  • The big question is, “How kinky is the transition from one curve to another?”

Bilinear curve (standard model)

  • Early examples of the model include Bacon and Watts (1971), and especially Laird and Ware (1982) - who include random effects, i.e. each subject has their own curve
  • The model is two lines that meet at an intersection point (for each subject)
  • This type of model is also called a bi-phasic regression
  • Given an unknown changepoint \(\gamma\), the simplest form of the model can be formally defined as

\[y_i\sim N(\beta_0+\beta_1(x_i-\gamma)I(x_i>\gamma)+\beta_2(x_i-\gamma)I(x_i<\gamma),\ \sigma^2)\]

New Soft Hinge curve

  • We model individual trajectories, optionally stratified by treatment, using a smooth transition around a subject-specific changepoint.
  • Let \(j = 1,\ldots,J\) index treatment groups, and let \(i = 1,\ldots,N_j\) index subjects assigned to treatment \(j\), where \(N_j\) is the number of subjects in group \(j\).
  • For each subject \(i\), let \(k = 1,\ldots,n_{ij}\) index repeated measurements, where \(n_{ij}\) is the number of observations available for that subject, and let the corresponding observation times be \(t_{ijk}\).
  • We write \(N = \sum_{j=1}^{J} \sum_{i=1}^{N_j} n_{ij}\) for the total number of observations across all treatment groups and subjects.
  • We take \(t = 0\) as baseline for all subjects, so that changepoint times are measured relative to baseline.
  • We use the term ‘hinge’ for the standard broken-stick basis function \(x^{+} = \max\left(x,0\right)\), which is zero for \(x \le 0\) and linear for \(x > 0\), so that including \(\beta_{2ij} x^{+}\) in a linear predictor creates a kink at \(x = 0\).

To obtain a differentiable version, we work with the smooth hinge

\[ h_{\kappa}\left(x\right) = x \Phi\left(\kappa x\right) + \frac{\phi\left(\kappa x\right)}{\kappa}, \qquad \kappa > 0, \]

which is continuous and differentiable, with derivative \(h_{\kappa}'\left(x\right) = \Phi\left(\kappa x\right)\), where \(\Phi\) and \(\phi\) denote the standard normal cumulative distribution function (CDF) and density, respectively. For large \(\kappa\), \(h_{\kappa}\left(x\right)\) is close to the hard hinge \(x^{+} = \max\left(x,0\right)\); for small \(\kappa\), the bend is gradual and spans a wider time interval around \(x = 0\).

The observed outcome is \(Y_{ijk}\), with conditional mean for \(k = 1,\ldots,n_{ij}\):

\[ \mu_{ijk} = \alpha_{0ij} + \beta_{1ij} t_{ijk} + \beta_{2ij} \left\{h_{\kappa_{ij}}\left(t_{ijk} - \tau_{ij}\right) - h_{\kappa_{ij}}\left(-\tau_{ij}\right)\right\}, \] where \(\tau_{ij}\) is a subject-specific changepoint and \(\kappa_{ij} > 0\) controls the smoothness of the transition.

The instantaneous slope with respect to \(t\) is

\[ \frac{d \mu_{ij}\left(t\right)}{d t} = \beta_{1ij} + \beta_{2ij} \Phi\left(\kappa_{ij}\left(t - \tau_{ij}\right)\right), \] which transitions smoothly from \(\beta_{1ij}\) for \(t \ll \tau_{ij}\) to \(\beta_{1ij} + \beta_{2ij}\) for \(t \gg \tau_{ij}\).

The parameter \(\beta_{1ij}\) is the initial slope before the gate, and \(\beta_{1ij} + \beta_{2ij}\) is the long-term slope after the gate.

The subject-specific smoothness parameter \(\kappa_{ij}\) controls how quickly the slope changes; a larger \(\kappa_{ij}\) means a sharper bend, and a smaller \(\kappa_{ij}\) means a slower bend over time.

Conditional on the mean trajectory, the outcome follows an error distribution:

\[ Y_{ijk} \left| \mu_{ijk}, \boldsymbol{\theta}_j \right. \sim f_{\varepsilon}\left(y; \mu_{ijk}, \boldsymbol{\theta}_j\right), \]

where \(\boldsymbol{\theta}_j\) collects dispersion, skewness, and tail-thickness parameters for treatment group \(j\), and the parameterization of \(f_{\varepsilon}\) ensures \(E\left(\varepsilon_{ijk}\right) = 0\) for \(\varepsilon_{ijk} = Y_{ijk} - \mu_{ijk}\). Examples include Gaussian, Student-\(t\), and skew-\(t\) families, parameterized so that \(E\left(\varepsilon_{ijk}\right) = 0\); allowing heavy-tailed residuals in nonlinear mixed-effects models is standard practice, e.g. Davidian and Giltinan (1995) Davidian and Giltinan (2003). In the Bayesian implementation for ACTG~315 we take \(f_{\varepsilon}\) to be a skew-\(t\) density.

Curve comparison

  • For model comparison between the two curve options, we use the leave-one-out cross-validation information criterion (LOOIC) (Vehtari et al. 2024)
    • It approximates how well the model might have predicted each observation had it been fitted without ever having seen that observation
    • Not true out-of-sample but a close approximation that is useful for smaller datasets and relatively fast

Posterior predictive distribution

To enable accurate calculation of predictive quantities we must be able to simulate from the posterior predictive distribution.

A simulation function for the residual distribution is created, which can then be applied for each combination of posterior vectors and new data matrices to arrive as predictions for any quantity of interest.

In this study the residual distribution used comes from Hansen (1994) (https://www.jstor.org/stable/2527081).

Non-symmetric conditional distributions

  • For symmetric distributions such as the normal, logistic, and t distributions:
    • \(\mu\ \equiv\) mean \(\equiv\) median \(\equiv\) mode
  • For the rest we have to explicitly state what we are trying to explain

When building models, ask yourself whether you are trying to model or predict, then ask yourself whether you are trying to model or predict

  • The expected value (mean)
  • The middle value (median)
  • The most likely value (mode)

Consider an ordinary normal regression on the log scale. What happens when you transform back? Are you still modelling the mean?

Single subject exploration

  • Here we consider a single subject and illustrate their fit with uncertainty.
  • This could easily be transformed into a web app on a company portal to allow a practitioner to interact with the model fit in future cases.
  • First code is created to illustrate a generic subject. This function calculates and plots the expected trajectory, the 95% credibility intervals of the expected trajectory, and 95% prediction intervals for a given subject.

For illustration we consider a few subjects.

Individual curve illustration

Subject 1 had few points and then dropped out so produced wide intervals.

Subject 2 did not follow the model pattern so produced even wider intervals which barely accommodate the observed pattern.

Subject 8 fitted the expected data pattern nearly perfectly, showing neat alignment and smaller intervals.

Subject 9 is extra, pun intended 😉

Conclusion

The real data curves are smooth, not kinky, so we smooth out the kink in 3 ways …

  • Bayesian parameter simulation
    • Even with the simplest two-line model, you get a smooth transition if you ‘integrate’ over the uncertainty in the turning point
  • Allowing the turning point to vary between subjects
    • Thus making the predictions valid for a random future subject, not just the middle subject from the observed data
  • Incorporating a smooth transition directly into the model curve
    • Allowing each subject to have their own level of curviness 😜

Technicalities

This presentation was created using the Reveal.js format in Quarto, using the RStudio IDE. Font and line colours according to UFS branding, and background image using image editor GIMP by compositing images from CoPilot.

  • Interactive plots were created by feeding ggplot2 plots to the plotly package, via the ggplotly function.

Parameter comparison

Model ParameterName Median Lower95 Upper95 Rhat
Abrupt \(\beta_{1,1}\) -3.764 -4.104 -3.448 1.004
Smooth \(\beta_{1,1}\) -4.659 -5.656 -4.085 1.040
Abrupt \(\beta_{2,1}\) 3.498 3.166 3.841 1.004
Smooth \(\beta_{2,1}\) 4.522 3.908 5.542 1.039
Abrupt \(\Lambda_{1}\) 0.252 0.000 0.485 1.000
Smooth \(\Lambda_{1}\) 0.212 -0.033 0.451 1.000
Abrupt \(\nu_{1}\) 2.320 2.049 3.138 1.000
Smooth \(\nu_{1}\) 2.552 2.078 3.863 1.000
Abrupt \(\sigma_{1}\) 0.624 0.410 1.448 1.001
Smooth \(\sigma_{1}\) 0.461 0.330 1.048 1.003
Abrupt \(L_{1,11}\) 2.316 1.601 10.667 1.019
Smooth \(L_{1,11}\) 2.841 1.708 12.595 1.028
Abrupt \(L_{1,21}\) -0.241 -30.211 22.713 1.011
Smooth \(L_{1,21}\) 1.948 -27.815 35.707 1.009
Abrupt \(L_{1,22}\) 4.682 2.602 23.693 1.066
Smooth \(L_{1,22}\) 5.306 2.355 29.102 1.040
Abrupt \(L_{1,31}\) 0.625 -16.788 28.243 1.013
Smooth \(L_{1,31}\) -0.441 -30.432 19.287 1.015
Abrupt \(L_{1,32}\) 1.607 -14.172 6.671 1.009
Smooth \(L_{1,32}\) 0.444 -17.068 7.551 1.010
Abrupt \(L_{1,33}\) 2.537 1.053 11.765 1.033
Smooth \(L_{1,33}\) 2.559 1.098 11.706 1.009
Abrupt Beta_RE_bar[1] 5.094 4.912 5.275 1.000
Smooth Beta_RE_bar[1] 5.098 4.922 5.274 1.001
Smooth Mu_kappa 3.636 -12.426 7.010 1.018
Abrupt Mu_tau 0.529 0.467 0.592 1.002
Smooth Mu_tau 0.440 0.329 0.519 1.029
Smooth Sigma_kappa 3.440 1.222 13.398 1.013
Abrupt Sigma_tau 0.099 0.034 0.154 1.000
Smooth Sigma_tau 0.092 0.042 0.148 1.001

References

Bacon, D. W., and D. G. Watts. 1971. “Estimating the transition between two intersecting straight lines.” Biometrika 58 (3): 525–34. https://doi.org/10.1093/biomet/58.3.525.
Burger, D. A., and R. Schall. 2015. “A Bayesian nonlinear mixed-effects regression model for the characterization of early bactericidal activity of tuberculosis drugs.” Journal of Biopharmaceutical Statistics 25 (6): 1247–71. https://doi.org/10.1080/10543406.2014.971170.
Davidian, M., and D. M. Giltinan. 1995. Nonlinear Models for Repeated Measurement Data. Chapman & Hall/CRC.
Davidian, M., and D. M. Giltinan. 2003. “Nonlinear models for repeated measurement data: an overview and update.” Journal of Agricultural, Biological, and Environmental Statistics 8 (4): 387. https://doi.org/10.1198/1085711032697.
Denwood, M. J. 2016. “runjags: An R package providing interface utilities, model templates, parallel computing methods and additional distributions for MCMC models in JAGS.” Journal of Statistical Software 71 (9): 1–25.
Donald, P. R., and A. H. Diacon. 2008. “The early bactericidal activity of anti-tuberculosis drugs: a literature review.” Tuberculosis 88: S75–83. https://doi.org/10.1016/S1472-9792(08)70038-6.
Ghosh, P., and F. Vaida. 2007. “Random changepoint modelling of HIV immunologic responses.” Statistics in Medicine 26 (9): 2074–87. https://doi.org/10.1002/sim.2671.
Hall, C. B., R. B. Lipton, M. Sliwinski, and W. F. Stewart. 2000. “A change point model for estimating the onset of cognitive decline in preclinical Alzheimer’s disease.” Statistics in Medicine 19 (11-12): 1555–66.
Hansen, B. E. 1994. “Autoregressive conditional density estimation.” International Economic Review, 705–30. https://doi.org/10.2307/2527081.
Hartig, F. 2021. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. Https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html.
Heerspink, H. J. L., N. Jongs, G. M. Chertow, et al. 2021. “Effect of dapagliflozin on the rate of decline in kidney function in patients with chronic kidney disease with and without type 2 diabetes: a prespecified analysis from the DAPA-CKD trial.” The Lancet Diabetes & Endocrinology 9 (11): 743–54. https://doi.org/10.1016/S2213-8587(21)00242-4.
Lachos, V. H., D. Bandyopadhyay, and D. K. Dey. 2011. “Linear and nonlinear mixed-effects models for censored HIV viral loads using normal/independent distributions.” Biometrics 67 (4): 1594–604. https://doi.org/10.1111/j.1541-0420.2011.01586.x.
Laird, N. M., and J. H. Ware. 1982. “Random-effects models for longitudinal data.” Biometrics, 963–74. https://doi.org/10.2307/2529876.
Morris, S. E., L. Dziobek-Garrett, and A. Yates. 2020. ushr: understanding suppression of HIV. Https://cran.r-project.org/package=ushr.
Perelson, A. S., A. U. Neumann, M. Markowitz, J. M. Leonard, and D. D. Ho. 1996. “HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time.” Science 271 (5255): 1582–86. https://doi.org/10.1038/387188a0.
van den Hout, A., G. Muniz-Terrera, and F. E. Matthews. 2011. “Smooth random change point models.” Statistics in Medicine 30 (6): 599–610. https://doi.org/10.1002/sim.4127.
Vehtari, A., J. Gabry, M. Magnusson, et al. 2024. loo: Efficient leave-one-out cross-validation and WAIC for Bayesian models. https://mc-stan.org/loo/.