Stepwise covariate modeling with nlmixr2scm

By Justin Wilkins in nlmixr2 nlmixr2scm

October 23, 2026

Watch on YouTube

nlmixr2scm is on CRAN!

Happily, nlmixr2scm, which does stepwise covariate modeling (SCM) for nlmixr2 models, is now on CRAN. It took a while - seems they’re very busy now that everyone is churning out packages with AI - but it’s finally there.

install.packages("nlmixr2scm")

If you’ve built covariate models in NONMEM, you’ve probably used scm in Perl-speaks-NONMEM. nlmixr2scm follows the same approach (Jonsson & Karlsson, 1998), and it works directly on an nlmixr2 fit. You don’t need to write control streams or config files or any of that kind of thing, although the command will probably be quite lengthy…

How SCM works

SCM searches in two phases, and both use likelihood-ratio tests (you probably know this, but a bit of underlining is always useful):

  1. Forward inclusion. Starting from the base model, it tests every candidate parameter-covariate relationship. The one that gives the largest significant drop in OFV is added, and the search repeats until no remaining candidate passes the forward threshold (p < 0.05 by default).
  2. Backward elimination. Starting from the final forward model, it removes each included relationship one at a time. Any relationship whose removal does not raise the OFV enough to pass a stricter threshold (p < 0.01 by default) is dropped.

runSCM() writes the covariate terms into the model body for you. Continuous covariates are centered (at the median, or at a reference value you choose), and categorical covariates are expanded into indicator columns. It takes the original, untransformed data.

A quick example

We’ll start from a one-compartment warfarin PK model:

library(nlmixr2)
library(nlmixr2scm)

pkdata <- nlmixr2data::warfarin[nlmixr2data::warfarin$dvid == "cp", ]

warf_pk <- function() {
  ini({
    tka <- log(1.15)
    tcl <- log(0.135)
    tv  <- log(7.0)
    eta.ka ~ 0.40
    eta.cl ~ 0.25
    eta.v  ~ 0.10
    prop.err <- 0.10
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv + eta.v)
    linCmt() ~ prop(prop.err)
  })
}

fit_base <- nlmixr2(warf_pk, pkdata, est = "focei",
                    control = foceiControl(print = 0))
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====] 0:00:00 
## done

Let’s test body weight (as a power or linear relationship) and sex on clearance and volume…

scm <- runSCM(
  fit        = fit_base,
  data       = pkdata,             # needed: the base model doesn't use wt or sex
  varsVec    = c("cl", "v"),       # parameters to test
  covarsVec  = "wt",               # continuous covariates
  catvarsVec = "sex",              # categorical covariates
  shapes     = c("power", "lin"),
  pVal       = list(fwd = 0.05, bck = 0.01),
  saveModels = FALSE,              # quick exploratory run: nothing written to disk
  print      = 0,
  confirm    = FALSE
)
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
##                step    covar var shape    objf  deltObjf      AIC      BIC
## cov_wt_power_v    1 wt_power   v power 449.113 -25.97642 926.4202 954.6238
##                numParams  qchisqr      pchisqr included searchType
## cov_wt_power_v         8 3.841459 3.456124e-07      yes    forward
##                      covNames covarEffect bsvReduction
## cov_wt_power_v cov_wt_power_v   0.8814415     61.74603
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
##                 step    covar var shape     objf  deltObjf      AIC     BIC
## cov_wt_power_cl    2 wt_power  cl power 442.6788 -6.434287 921.9859 953.715
##                 numParams  qchisqr    pchisqr included searchType
## cov_wt_power_cl         9 3.841459 0.01119381      yes    forward
##                        covNames covarEffect bsvReduction
## cov_wt_power_cl cov_wt_power_cl   0.5844377     17.27476
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
##                 step    covar var shape     objf deltObjf     AIC      BIC
## cov_wt_power_cl    1 wt_power  cl power 448.2978 5.619047 925.605 953.8086
##                 numParams  qchisqr    pchisqr included searchType
## cov_wt_power_cl         8 6.634897 0.01776631  dropped   backward
##                        covNames covarEffect bsvReduction
## cov_wt_power_cl cov_wt_power_cl   0.5844377     20.15977
## calculating covariance matrix
## [====|====|====|====|====|====|====|====|====|====done
## Direction  Step  Relation                Ref OFV         OFV      dOFV    p-value  Decision 
## ------------------------------------------------------------------------------------------- 
## Forward    1     wt_power~v [power]      475.089     449.113   -25.976     0.0000  Added
## Forward    2     wt_power~cl [power]     449.113     442.679    -6.434     0.0112  Added
## Forward    3     sex_female~v [cat]      442.679     440.495    -2.184     0.1395  Not selected
## Backward   1     wt_power~cl [power]     442.679     448.298     5.619     0.0178  Removed
## Backward   2     wt_power~v [power]      448.298     474.521    26.223     0.0000  Retained
## Direction  Step  Relation                Ref OFV         OFV      dOFV    p-value  Decision 
## ------------------------------------------------------------------------------------------- 
## Forward    1     wt_power~v [power]      475.089     449.113   -25.976     0.0000  Added
## Forward    1     wt_lin~v [lin]          475.089     450.010   -25.080     0.0000  Not selected
## Forward    1     sex_female~v [cat]      475.089     458.423   -16.666     0.0000  Not selected
## Forward    1     wt_power~cl [power]     475.089     470.203    -4.887     0.0271  Not selected
## Forward    1     wt_lin~cl [lin]         475.089     470.712    -4.378     0.0364  Not selected
## Forward    1     sex_female~cl [cat]     475.089     474.173    -0.917     0.3383  Not selected
## 
## Forward    2     wt_power~cl [power]     449.113     442.679    -6.434     0.0112  Added
## Forward    2     wt_lin~cl [lin]         449.113     444.150    -4.963     0.0259  Not selected
## Forward    2     sex_female~v [cat]      449.113     448.361    -0.752     0.3858  Not selected
## Forward    2     sex_female~cl [cat]     449.113     448.984    -0.129     0.7190  Not selected
## 
## Forward    3     sex_female~v [cat]      442.679     440.495    -2.184     0.1395  Not selected
## Forward    3     sex_female~cl [cat]     442.679     441.515    -1.164     0.2807  Not selected
## 
## Backward   1     wt_power~v [power]      442.679     470.173    27.494     0.0000  Retained
## Backward   1     wt_power~cl [power]     442.679     448.298     5.619     0.0178  Removed
## 
## Backward   2     wt_power~v [power]      448.298     474.521    26.223     0.0000  Retained

Printing the result gives a short overview:

scm
## ── nlmixr2scm result ───────────────────────────────────────────────────────────
## Search      : scm (focei), 6 candidate relations
## OFV         : base 475.089 -> final 448.298 (dOFV -26.792)
## Covariates  : wt_power~v [power]
## Use summary() for options, model comparison and step tables.

For the full picture, use summary(). It lists the options the search ran with, compares the base model with the final forward and final backward models (OFV, ΔOFV, AIC, BIC, parameter count and covariates), prints the step-by-step and all-candidate tables, and tells you where the models and report files were saved on disk, if they were:

summary(scm)
## ── nlmixr2scm summary ──────────────────────────────────────────────────────────
## ── Options ─────────────────────────────────────────────────────────────────────
## Search type        : scm
## p-value fwd / bck  : 0.05 / 0.01
## Estimation         : focei; control: foceiControl (inherited from fit)
## Candidates         : 6: wt_power~cl [power], wt_lin~cl [lin], wt_power~v
##                      [power], wt_lin~v [lin], sex_female~cl [cat], sex_female~v
##                      [cat]
## Included relations : none
## Shapes             : power, lin
## Centers            : median of each continuous covariate
## Categorical cutoff : 0.05
## Profile warm-start : stalled candidates only (stallTol = 0)
## Retries            : up to 3; maxDeltaOFV = Inf; OFV tolerance = 0; p-value
##                      underflow retried; on exhaustion keep best attempt
## Parallel           : workers = sequential; rxThreads per worker = 16
## Save models        : no
## ── Model comparison ────────────────────────────────────────────────────────────
##  Model          OFV     dOFV    AIC     BIC     Params
##  Base           475.089 0.000   950.397 975.075 7     
##  Forward final  442.679 -32.411 921.986 953.715 9     
##  Backward final 448.298 -26.792 925.605 953.809 8     
## 
## SCM covariates:
##   Base           : (none)
##   Forward final  : wt_power~v [power], wt_power~cl [power]
##   Backward final : wt_power~v [power]
## ── SCM Step Summary ────────────────────────────────────────────────────────────
## Direction  Step  Relation                Ref OFV         OFV      dOFV    p-value  Decision 
## ------------------------------------------------------------------------------------------- 
## Forward    1     wt_power~v [power]      475.089     449.113   -25.976     0.0000  Added
## Forward    2     wt_power~cl [power]     449.113     442.679    -6.434     0.0112  Added
## Forward    3     sex_female~v [cat]      442.679     440.495    -2.184     0.1395  Not selected
## Backward   1     wt_power~cl [power]     442.679     448.298     5.619     0.0178  Removed
## Backward   2     wt_power~v [power]      448.298     474.521    26.223     0.0000  Retained
## ── SCM All Candidates ──────────────────────────────────────────────────────────
## Direction  Step  Relation                Ref OFV         OFV      dOFV    p-value  Decision 
## ------------------------------------------------------------------------------------------- 
## Forward    1     wt_power~v [power]      475.089     449.113   -25.976     0.0000  Added
## Forward    1     wt_lin~v [lin]          475.089     450.010   -25.080     0.0000  Not selected
## Forward    1     sex_female~v [cat]      475.089     458.423   -16.666     0.0000  Not selected
## Forward    1     wt_power~cl [power]     475.089     470.203    -4.887     0.0271  Not selected
## Forward    1     wt_lin~cl [lin]         475.089     470.712    -4.378     0.0364  Not selected
## Forward    1     sex_female~cl [cat]     475.089     474.173    -0.917     0.3383  Not selected
## 
## Forward    2     wt_power~cl [power]     449.113     442.679    -6.434     0.0112  Added
## Forward    2     wt_lin~cl [lin]         449.113     444.150    -4.963     0.0259  Not selected
## Forward    2     sex_female~v [cat]      449.113     448.361    -0.752     0.3858  Not selected
## Forward    2     sex_female~cl [cat]     449.113     448.984    -0.129     0.7190  Not selected
## 
## Forward    3     sex_female~v [cat]      442.679     440.495    -2.184     0.1395  Not selected
## Forward    3     sex_female~cl [cat]     442.679     441.515    -1.164     0.2807  Not selected
## 
## Backward   1     wt_power~v [power]      442.679     470.173    27.494     0.0000  Retained
## Backward   1     wt_power~cl [power]     442.679     448.298     5.619     0.0178  Removed
## 
## Backward   2     wt_power~v [power]      448.298     474.521    26.223     0.0000  Retained
## ── Final model ─────────────────────────────────────────────────────────────────
## Retained:
## wt_power~v [power]
## ℹ Final model OFV: 448.298
## ────────────────────────────────────────────────────────────────────────────────
## ── Files ───────────────────────────────────────────────────────────────────────
## None: models and reports were not saved (saveModels = FALSE).

Underneath, the result is a list you can dig into:

  • summaryTable has one row per candidate tested in both phases, with the OFV, ΔOFV, p-value, AIC/BIC, and whether the relationship was included.
  • resFwd[[1]] and resBck[[1]] are the final forward and backward fits. These are ordinary nlmixr2 fit objects, so you can use your usual diagnostics and plots on them.
scm$resBck[[1]]$parFixedDf
##                  Estimate         SE      %RSE Back-transformed  CI Lower
## tka            -0.4775266 0.25661205 53.737750        0.6203158 0.3751321
## tcl            -1.9991278 0.05323381  2.662852        0.1354534 0.1220330
## tv              2.1282396 0.03044324  1.430442        8.4000660 7.9135138
## prop.err        0.2233200 0.02903750 13.002646        0.2233200 0.1664075
## cov_wt_power_v  0.8686918 0.13432749 15.463193        0.8686918 0.6054147
##                 CI Upper  BSV(CV%) Shrink(SD)%
## tka            1.0257497 71.587391   49.660122
## tcl            0.1503497 27.228417    3.203639
## tv             8.9165332  9.061567   51.577626
## prop.err       0.2802324        NA          NA
## cov_wt_power_v 1.1319688        NA          NA

If you want to control exactly which relationships are tested, pass pairsVec instead of the full varsVec × covarsVec grid. You can also use includedRelations to force known relationships (allometric weight on clearance, for example) into the backward-elimination model.

Real-world searches

Real-world SCM runs often fit hundreds of candidate models, and nlmixr2scm has several features for that:

  • Caching and resume. With saveModels = TRUE (the default), every candidate fit and step table is saved to outputDir. If a run is interrupted, call runSCM() again with the same outputDir and it picks up where it stopped. restart = TRUE starts from scratch instead.

  • Parallel fitting. workers fits the candidates within each step in parallel through future. rxThreads sets how many rxode2 threads each worker uses. runSCM() checks workers * rxThreads against your physical core count before it starts, so that it doesn’t quietly murder your computer by maxing out its capacity.

    runSCM(fit_base, data = pkdata, varsVec = c("cl", "v"), covarsVec = "wt",
       workers = 4, rxThreads = 1)

  • Protection against bad fits. Some candidate fits stall at the zero-effect initial estimate. This happens sometimes with ODE models. runSCM() detects these stalled candidates and refits them from a starting value obtained with a quick 1-D profile. Candidates that land on an implausible OFV are retried from perturbed initial estimates, and the best attempt is kept. Both behaviours are switched on by default. By default, a p-value that underflows to zero also counts as implausible, which means very strong effects (an OFV drop of more than about 70 for one parameter) get retried too. If you trust those, retryOnUnderflow = FALSE saves the extra fits.

  • Fixed centering. centers = c(wt = 70) centers weight at 70 kg instead of the per-dataset median. The coefficients then keep the same reference across datasets.

Learn more

  • The vignette gives the full argument reference and a longer worked example: vignette("runSCM", package = "nlmixr2scm").
  • Source, bug reports, and feature requests are on GitHub.
  • The NEWS file has release notes.

Acknowledgements

nlmixr2scm was written together with Matthew Fidler and Yaping Liu and contains contributions from Bill Denney, Vipul Mann, Vishal Sarsani, and Christian Bartels. The design owes a great deal to the original PsN scm implementation by Lars Lindbom, Niclas Jonsson, Pontus Pihlgren, Mats Karlsson, Andrew Hooker, Kajsa Harling, Rikard Nordgren and Svetlana Freiberga. Thank you all!

Posted on:
October 23, 2026
Length:
10 minute read, 1974 words
Categories:
nlmixr2 nlmixr2scm
See Also: