{pmrm}: an R package for fitting progression models for repeated measures

Will Landau, Lars Lau Raket, Kasper Kristensen

Progression models for repeated measures (PMRMs)

Framing progressive diseases

In progressive diseases, patients are expected to deteriorate over time.

When deterioration cannot be reversed or cured, interventions on the disease aim to mitigate decline or slow progression.

Traditional statistical estimates may not capture the impact of a treatment on a patient’s life.

Modeling progressive diseases

Mixed Models for Repeated Measures (MMRMs)

• Popular longitudinal linear mixed model for continuous outcomes in clinical trials (Mallinckrodt et al. (2008), Mallinckrodt and Lipkovich (2017)).
• For certain estimands, can account for responses missing at random (MAR) without multiple imputation (Holzhauer et al. (2024)).
• Flexible but not powerful.

Variations on MMRMs

• Constrained longitudinal data analysis (cDLA) pools treatment and control at baseline (Dysken et al. (2014), Wang et al. (2022)).
• cLDA model with natural cubic splines (Donohue et al. (2023)).

Progression Models for Repeated Measures (PMRMs)

• Nonlinear longitudinal models with interpretable parameters for disease progression (Raket (2022)).

MMRMs use discrete time (visits).

The effect depends on the visit.

Proportional decline PMRM

Proportional slowing PMRM

PMRMs are nonlinear.

\(\alpha_\text{placebo}(t)\): natural cubic spline fitted to placebo response evaluated at time \(t\).

Proportional decline PMRM:

\[ \begin{aligned} \alpha_\text{treatment}(t) - \alpha_\text{treatment}(0) = (\alpha_\text{placebo}(t) - \alpha_\text{placebo}(0)) \cdot (1 - \beta) \end{aligned} \]

\(\beta\) = reduction in decline since baseline (relative to placebo)

Proportional slowing PMRM:

\[ \begin{aligned} \alpha_\text{treatment}(t) = \alpha_\text{placebo}\left( t \cdot (1 - \beta) \right) \end{aligned} \]

\(\beta\) = slowing of disease progression since baseline (relative to placebo)

Non-proportional PMRMs use visit-specific \(\beta_1, \ldots, \beta_J\).

PMRMs are difficult to fit.

Slow execution:

library(nlme)
system.time({
  fitted_pmrm <- gnls(
    model = custom_pmrm_formula,
    data = clinical_data,
    params = list(
      v0 + v1 + v2 + v3 + v4 + v5 + v6 ~ 1,
      as.formula(paste0('k ~ ', covariates, ' + 0')),
      coefficient_formula
    ),
    correlation = corSymm(form = ~ as.numeric(visit_sequence) | USUBJID),
    weights = varIdent(form = ~ 1 | visit_sequence),
    start = start_coefficients,
    control = gnlsControl(nlsTol = tolerance)
  )
})
#>    user  system elapsed
#> 513.851  99.188 491.208

Frequent errors:

#> Error in gnls(model = custom_pmrm_formula, data = clinical_data, params = :
#>   step halving factor reduced below minimum in NLS step

{pmrm} R package

{pmrm} R package

https://openpharma.github.io/pmrm

Fits the decline and slowing PMRMs (proportional and non-proportional).

Runs fast and reliably thanks to the {RTMB} package by Kasper Kristensen.

• 200x speedup on real clinical data.
• Divergences are rare.

Analyst-friendly interface with familiar S3 methods.

Part of CRAN and the OpenPharma GitHub organization (like packages mmrm, DoseFinding, and RobinCar2).

Core interface

Estimates and model comparison

Visualization and prediction

Features on the roadmap

Likelihood ratio tests for model comparison (e.g. tests of non-proportionality).

Significance tests of individual decline/slowing parameters.

Allow arbitrary parameterization the decline/slowing parameters (not just proportional vs non-proportional).

Simulating data

library(pmrm)
data <- pmrm_simulate_slowing_proportional(
  patients = 100,

  # Visit times inform covariance matrix and default spline knots:
  visit_times = c(0, 12, 24, 36, 52, 64, 76),

  # Placebo mean response at each time point:
  alpha = c(50, 44, 38, 32, 30, 28, 25),

  # Slowing in placebo (always 0) and treatment arms, respectively:
  beta = c(0, 0.3),

  # Visit-specific residual standard deviations:
  sigma = rep(3, 7),

  # Latent parameters for residual correlations among visits:
  rho = rep(1, 7 * (7 - 1) / 2),

  # Standard deviation of observed times around scheduled visits:
  tau = 3
)

Simulating data

data

#> # A tibble: 700 × 11
#>    patient   visit   arm       i     j     k     y     t  beta    mu     e
#>    <chr>     <ord>   <ord> <int> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
#>  1 patient_1 visit_1 arm_1     1     1     1  51.8  0      0    50   1.78 
#>  2 patient_1 visit_2 arm_1     1     2     1  42.4 17.4    0    41.4 0.973
#>  3 patient_1 visit_3 arm_1     1     3     1  42.7 19.4    0    40.4 2.33 
#>  4 patient_1 visit_4 arm_1     1     4     1  33.0 38.2    0    31.4 1.58 
#>  5 patient_1 visit_5 arm_1     1     5     1  32.3 50.3    0    30.1 2.18 
#>  6 patient_1 visit_6 arm_1     1     6     1  27.6 70.1    0    26.5 1.10 
#>  7 patient_1 visit_7 arm_1     1     7     1  28.8 73.3    0    25.7 3.10 
#>  8 patient_2 visit_1 arm_2     2     1     2  52.2  0      0.3  50   2.18 
#>  9 patient_2 visit_2 arm_2     2     2     2  49.8  7.93   0.3  47.2 2.63 
#> 10 patient_2 visit_3 arm_2     2     3     2  43.0 27.7    0.3  40.4 2.61 
#> # ℹ 690 more rows

Fitting a model

fit <- pmrm_model_slowing_proportional(
  data = data,
  outcome = "y",
  time = "t",
  patient = "patient",
  visit = "visit",
  arm = "arm",
  visit_times = c(0, 12, 24, 36, 52, 64, 76)
)

Model output

print(fit)

#> Model:
#> 
#>   PMRM type:        slowing
#>   Parameterization: proportional
#> 
#> Fit:
#> 
#>   Convergence:    converged
#>   Observations:   700
#>   Parameters:     36
#>   Log likelihood: -1301.719
#>   Deviance:       2603.438
#>   AIC:            2675.438
#>   BIC:            2839.277
#> 
#> Treatment effects:
#> 
#>         estimate  std.error
#>   arm_2 0.300327 0.01234619

names(fit)

#>  [1] "data"           
#>  [2] "constants"      
#>  [3] "options"        
#>  [4] "objective"      
#>  [5] "model"          
#>  [6] "optimization"   
#>  [7] "report"         
#>  [8] "initial"        
#>  [9] "final"          
#> [10] "estimates"      
#> [11] "standard_errors"
#> [12] "metrics"        
#> [13] "spline"

Estimates

pmrm_estimates(fit, parameter = "beta")

#> # A tibble: 2 × 6
#>   parameter arm   estimate standard_error  lower  upper
#>   <chr>     <ord>    <dbl>          <dbl>  <dbl>  <dbl>
#> 1 beta      arm_1    0            NA      NA     NA    
#> 2 beta      arm_2    0.300         0.0123  0.276  0.325

pmrm_marginals(fit) |> head(n = 4)

#> # A tibble: 4 × 7
#>   arm   visit    time estimate standard_error lower upper
#>   <ord> <ord>   <dbl>    <dbl>          <dbl> <dbl> <dbl>
#> 1 arm_1 visit_1     0     50.4          0.318  49.8  51.1
#> 2 arm_1 visit_2    12     44.3          0.308  43.7  44.9
#> 3 arm_1 visit_3    24     38.2          0.330  37.6  38.9
#> 4 arm_1 visit_4    36     32.1          0.320  31.5  32.8

Visualization

fit_plot <- plot(
  fit,
  show_data = FALSE,
  show_predictions = TRUE,
  show_marginals = TRUE,
  facet = FALSE
)

Visualization

fit_plot

Validation

✅ {pmrm} recovers true parameters across a range of scenarios: https://openpharma.github.io/pmrm/articles/validation.html.

✅ Full coverage unit testing with testthat to ensure:

1. Simulations and post-processing are correct.
2. Model output has the correct formatting.
3. Regressions are caught early.

{RTMB}

About {RTMB}

An efficient toolkit for fitting complex statistical models in R (Kristensen (2026)).

Developed by Kasper Kristensen at Danish Technical University.

Key features:

✅ Fast, exact, high-order automatic differentiation

✅ Fast flexible Laplace approximation of marginal likelihoods for models with many random effects.

✅ Uses R’s native syntax for model specification.

Tool stack highlights

%%{init: {"theme": "default", "themeVariables": {"fontSize": "40px"}, "flowchart": {"nodeSpacing": 30, "rankSpacing": 150, "padding": 20}} }%%
graph BT
  A["CppAD"] --> C
  B["Matrix"] --> C
  C["TMB"] --> D
  D["RTMB"]

{RTMB} = Native R syntax for {TMB} models (Kristensen (2026)).

{TMB} = Template Model Builder (Kristensen et al. (2016)).

{CppAD} = automatic differentiation in C++, runs exact high-order derivatives (Bell (2005)).

{Matrix} = numerical linear algebra in R, including sparse matrices (Bates et al. (2025)).

How to use {RTMB}

Step 1: write minus log likelihood in R.

\[ \begin{aligned} y_1, \ldots, y_n \stackrel{\text{iid}}{\sim} \text{Normal}(\mu, \ \sigma^2) \end{aligned} \]

library(RTMB)

data <- ChickWeight

objective <- function(parameters) {
  # Optional: attach list elements to objective function environment:
  getAll(parameters, data)

  # Optional: decalre the response variable to enable extra RTMB features:
  y <- OBS(weight)

  # Calculate the minus log likelihood using R syntax:
  -sum(dnorm(y, alpha + beta * Time, 1, log = TRUE))
}

How to use {RTMB}

Step 2: Register the objective function with RTMB

model <- MakeADFun(
  func = objective,
  parameters = list(alpha = 0, beta = 0),
  silent = TRUE
)

Step 3: fit the model with any optimizer in R

fit <- stats::nlminb(model$par, model$fn, model$gr)

How to use {RTMB}

Step 4: get convergence and parameter estimates

fit$convergence # 0 indicates convergence

#> [1] 0

report <- sdreport(model)
unlist(as.list(report, "Est")) # parameter estimates

#>     alpha      beta 
#> 27.467425  8.803039

unlist(as.list(report, "Std")) # standard errors

#>       alpha        beta 
#> 0.078031199 0.006159823

Fit an equivalent Bayesian model

Optional: set priors in the objective function

objective <- function(parameters) {
  getAll(parameters, data)
  y <- OBS(weight)

  -sum(dnorm(y, alpha + beta * Time, 1, log = TRUE)) -
    dnorm(alpha, 121, 100, log = TRUE) - # Optional prior on alpha
    dnorm(beta, 0, 100, log = TRUE) # Optional prior on beta
}

Use {tmbstan} to fit the model.

mcmc <- tmbstan::tmbstan(model, chains = 4, refresh = 0)
posterior::summarize_draws(mcmc)

#> # A tibble: 3 × 10
#>   variable       mean     median      sd     mad       q5     q95  rhat ess_bulk
#>   <chr>         <dbl>      <dbl>   <dbl>   <dbl>    <dbl>   <dbl> <dbl>    <dbl>
#> 1 alpha         27.5       27.5  0.0802  0.0819    2.73e1  2.76e1 1.00     1053.
#> 2 beta           8.80       8.80 0.00627 0.00624   8.79e0  8.81e0 1.00     1105.
#> 3 lp__     -436638.   -436638.   1.06    0.727    -4.37e5 -4.37e5 1.000    1363.
#> # ℹ 1 more variable: ess_tail <dbl>

Proportional slowing PMRM with {RTMB}

objective <- function(constants, parameters) {
  getAll(parameters, constants)
  y <- OBS(y)

  # Spline fitted to placebo respones:
  spline <- splinefun(x = spline_knots, y = alpha, method = "natural")

  # Covariate adjustment with sparse matrices (Matrix R package):
  adjustment <- (X %*% gamma) - sum(X_column_means * gamma)

  # Vector of fitted responses:
  mu <- spline((1 - beta[arm]) * time) + adjustment

  # Unstructured covariance among visits:
  sigma <- diag(exp(phi))
  Sigma <- sigma %*% unstructured(n_visits)$corr(rho) %*% sigma

  # Multivariate normal likelihood (on the minus log scale):
  target <- 0
  for (patient in seq_len(n_patients)) {
    visits <- seq(to = patient * n_visits, length.out = n_visits)
    target <- target - RTMB::dmvnorm(y[visits], mu[visits], Sigma, log = TRUE)
  }
  target
}

Full runnable example

https://wlandau.github.io/mbsw2026/example.R

103 lines of code

Recap

Recap

Progression models for repeated measures (PMRMs) provide meaningful estimates for patients.

The {pmrm} package fits PMRMs robustly and efficiently.

The backend package {RTMB} runs fast reliable models in native R syntax.

Slides: wlandau.github.io/mbsw2026

References

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2025. Matrix: Sparse and Dense Matrix Classes and Methods. https://doi.org/10.32614/CRAN.package.Matrix.

Bell, Bradley M. 2005. CppAD: A Package for c++ Algorithmic Differentiation. http://www.coin-or.org/CppAD.

Donohue, Michael C., Oliver Langford, Philip S. Insel, et al. 2023. “Natural Cubic Splines for the Analysis of Alzheimer’s Clinical Trials.” Pharmaceutical Statistics 22 (3): 508–19. https://doi.org/https://doi.org/10.1002/pst.2285.

Dysken, Maurice W., Mary Sano, Sanjay Asthana, et al. 2014. “Effect of Vitamin e and Memantine on Functional Decline in Alzheimer Disease: The TEAM-AD VA Cooperative Randomized Trial.” JAMA 311 (1): 33–44. https://doi.org/10.1001/jama.2013.282834.

Holzhauer, Björn, Sebastian Weber, Lukas Widmer, and Andrew Bean. 2024. “Bayesian Mixed Effects Model for Repeated Measures.” In Applied Modeling in Drug Development. Novartis AG. https://opensource.nibr.com/bamdd/src/02h_mmrm.html.

Kristensen, Kasper. 2026. RTMB: ’R’ Bindings for ’TMB’. https://doi.org/10.32614/CRAN.package.RTMB.

Kristensen, Kasper, Anders Nielsen, Casper W. Berg, Hans Skaug, and Bradley M. Bell. 2016. “TMB: Automatic Differentiation and Laplace Approximation.” Journal of Statistical Software 70 (5): 1–21. https://doi.org/10.18637/jss.v070.i05.

Mallinckrodt, C. H., P. W. Lane, D. Schnell, et al. 2008. “Recommendations for the Primary Analysis of Continuous Endpoints in Longitudinal Clinical Trials.” Therapeutic Innovation and Regulatory Science 42: 303–19. https://doi.org/https://doi.org/10.1177/009286150804200402.

Mallinckrodt, Craig H., and Ilya Lipkovich. 2017. Analyzing Longitudinal Clinical Trial Data: A Practical Guide. CRC Press, Taylor; Francis Group.

Raket, Lars Lau. 2022. “Progression Models for Repeated Measures: Estimating Novel Treatment Effects in Progressive Diseases.” Statistics in Medicine 10410 (28): 5537–57. https://doi.org/https://doi.org/10.1002/sim.9581.

Wang, Guoqiao, Lei Liu, Yan Li, et al. 2022. “Proportional Constrained Longitudinal Data Analysis Models for Clinical Trials in Sporadic Alzheimer’s Disease.” Alzheimer’s & Dementia: Translational Research & Clinical Interventions 8 (1): e12286. https://doi.org/https://doi.org/10.1002/trc2.12286.