Stepwise Covariate Modelling (SCM) is the standard automated covariate selection procedure in population PK/PD analysis. It uses a likelihood-ratio test to evaluate whether adding a covariate-parameter relationship significantly improves model fit, and proceeds in two phases:
runSCM() implements both phases for nlmixr2
fits. It automatically generates the covariate expressions in the model
body — centering continuous covariates at their observed median,
creating indicator columns for categorical covariates — so no data
pre-processing is required.
We use nlmixr2data::warfarin, filtered to PK
observations only (dvid == "cp"), with covariates body
weight (wt, continuous) and sex (sex,
categorical).
Pass the original, untransformed dataset to
runSCM(). All covariate centering and other shape
transformations are generated inside the model body by the SCM
machinery; applying them to the data first would lead to
double-transformation.
pkdata <- warfarin[warfarin$dvid == "cp", ]
d_subj <- pkdata[!duplicated(pkdata$id), ]
cat(sprintf("Subjects: %d Rows: %d\n", length(unique(pkdata$id)), nrow(pkdata)))
#> Subjects: 32 Rows: 283
cat(sprintf(
"wt median: %.1f range: %.1f – %.1f\n",
median(d_subj$wt), min(d_subj$wt), max(d_subj$wt)
))
#> wt median: 71.7 range: 40.0 – 102.0
cat(sprintf(
"sex levels: %s\n",
paste(names(table(d_subj$sex)), table(d_subj$sex), sep = "=", collapse = ", ")
))
#> sex levels: female=5, male=27The base model contains no covariate terms — runSCM()
adds them.
warf_pk <- function() {
ini({
tka <- log(1.15) # log absorption rate constant (h^-1)
tcl <- log(0.135) # log clearance (L/h)
tv <- log(7.0) # log volume of distribution (L)
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 = nlmixr2est::foceiControl(print = 0),
table = nlmixr2est::tableControl(cwres = TRUE)
)fit_base$parFixedDf
#> Estimate SE %RSE Back-transformed CI Lower CI Upper
#> tka -0.5945217 0.34277068 57.654862 0.5518264 0.2818614 1.0803623
#> tcl -1.9935985 0.06516081 3.268502 0.1362044 0.1198743 0.1547591
#> tv 2.0897598 0.04288462 2.052131 8.0829737 7.4313499 8.7917357
#> prop.err 0.2239179 0.03281636 14.655529 0.2239179 0.1595991 0.2882368
#> BSV(CV%) Shrink(SD)%
#> tka 77.12023 51.705421
#> tcl 27.30318 4.952248
#> tv 17.53332 18.821830
#> prop.err NA NA
cat(sprintf("Base OFV: %.3f\n", fit_base$objf))
#> Base OFV: 475.089runSCM() supports two main ways to define the search
space.
The most explicit approach is pairsVec: a list of
list(var=, covar=, shapes=) items. Each item names the PK
parameter (var), the covariate (covar), and
optionally the functional shapes to test. This is useful when you want
exact control over which relationships are evaluated.
When the search space is a full Cartesian product of parameters and
covariates, a shorter option is varsVec plus
covarsVec, with catvarsVec used for
categorical covariates. That is the approach used in the worked example
below.
The data argument must be supplied explicitly whenever
the base model does not reference a covariate column — otherwise
nlme::getData(fit) may drop it.
runSCM() argument reference| Argument | Default | Meaning |
|---|---|---|
varsVec |
NULL |
PK parameters to test (used with covarsVec; ignored
when pairsVec supplied) |
covarsVec |
NULL |
Continuous covariates (ignored when pairsVec
supplied) |
catvarsVec |
NULL |
Categorical covariates; indicator columns created automatically |
pairsVec |
NULL |
Explicit parameter-covariate pairs; overrides
varsVec/covarsVec |
shapes |
"power" |
Shape(s) for continuous covariates: "power",
"lin", "log", "identity" |
customShapes |
NULL |
Named list of additional shape builder functions |
centers |
NULL |
Named vector fixing the reference (centering) value for continuous
covariates, e.g. c(BW = 70), instead of the per-dataset
median |
pVal |
list(fwd=0.05, bck=0.01) |
Forward and backward p-value thresholds |
searchType |
"scm" |
"scm" (both phases), "forward", or
"backward" |
data |
NULL |
Dataset; required when base model omits covariate columns |
inits |
list() |
Initial theta estimates for covariate parameters, by shape |
missingToken |
NA |
Sentinel value for missing covariate observations |
catCutoff |
0.05 |
Minimum non-reference level prevalence to test |
includedRelations |
NULL |
Relations forced into the backward-start model |
outputDir |
NULL |
Subdirectory for saved model .rds files; auto-named
when NULL |
saveModels |
TRUE |
Write fitted models to outputDir |
verbose |
FALSE |
Print full candidate tables at each step |
control |
NULL |
Control object (e.g. saemControl()) passed to every
nlmixr2() call |
print |
100 |
Optimiser progress-print interval, applied as
control$print |
restart |
FALSE |
Discard existing cache and start fresh |
workers |
NULL |
Sequential; set > 1 for parallel candidate fitting |
confirm |
TRUE |
Prompt for confirmation in interactive sessions |
profileInit |
FALSE |
Warm-start every forward candidate’s covariate coefficient with a frozen 1-D profile (Brent method) |
profileInitOnStall |
TRUE |
Warm-start only candidates whose ordinary fit stalls (OFV
improvement <= stallTol) |
stallTol |
0 |
OFV-improvement threshold below which a forward candidate is considered stalled |
maxRetries |
3L |
Maximum retry attempts per candidate when the OFV is deemed unrealistic |
maxDeltaOFV |
Inf |
Absolute ceiling on plausible \|dOFV\| before a
candidate is retried |
retryPerturbSD |
0.5 |
SD of the log-normal perturbation applied to cov_init
on odd-numbered retries |
retrySmallInit |
0.01 |
Near-zero covariate theta init used on even-numbered retries |
retryOFVTolerance |
NULL |
Margin by which candidate OFV must exceed the parent before
retrying; auto-detected (10 for SAEM, 0
otherwise) when NULL |
retryFailOnExhaustion |
FALSE |
Exclude a candidate as failed (instead of accepting the best retry) when all retries are exhausted |
retryOnUnderflow |
TRUE |
Treat a p-value underflow (OFV drop above ~70 at
df = 1) as unrealistic and retry; set FALSE to
skip the extra fits for genuinely strong effects |
runSCM() supports four built-in shapes for continuous
covariates. "power" and "lin" use the
observed median as the centering value, computed
automatically from the data.
| Shape | Expression in model body | Interpretation |
|---|---|---|
"power" |
log(cov / median(cov)) |
Allometric / power-law relationship |
"lin" |
cov - median(cov) |
Linear deviation from the median |
"log" |
log(cov) |
Log-transformed covariate (no centering) |
"identity" |
cov |
Raw covariate (no transformation) |
Categorical covariates always use the "cat" shape: a
binary indicator column cov_<level> multiplied by a
theta parameter. The reference level is the most frequent category.
Levels representing fewer than catCutoff (default 5 %) of
subjects are merged with the reference and not tested.
Custom shapes can be added via the customShapes argument
— a named list of functions function(col, center, level)
that return a character expression string.
By default "power" and "lin" center on the
per-dataset median, which can differ between datasets fit with the same
model. Pass centers (a named numeric vector, e.g.
centers = c(wt = 70)) to fix the reference value for one or
more continuous covariates instead, so covariate coefficients stay on a
consistent reference across datasets. Covariates not named in
centers continue to use the median.
When the search space is a full Cartesian product of parameters ×
covariates, use varsVec and covarsVec instead
of enumerating every pair.
The example below tests wt~cl, wt~v,
sex~cl, and sex~v, with both power and linear
shapes for the continuous covariate — six candidates in total.
scm <- runSCM(
fit = fit_base,
data = pkdata,
varsVec = c("cl", "v"),
covarsVec = "wt",
catvarsVec = "sex",
shapes = c("power", "lin"),
pVal = list(fwd = 0.05, bck = 0.01),
searchType = "scm",
saveModels = TRUE,
verbose = FALSE,
restart = TRUE,
print = 0,
workers = 1L,
rxThreads = 2L,
confirm = FALSE
)
#> 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
#> Warning: ! sex_female ~ v: unrealistic OFV on attempt 1/4: OFV increased vs parent
#> (452.129 > 449.113).
#> ℹ Retrying with perturbed init.
#> 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
#> Warning: ! sex_female ~ cl: unrealistic OFV on attempt 1/4: OFV increased vs parent
#> (445.315 > 442.679).
#> ℹ Retrying with perturbed init.
#> Warning: ! sex_female ~ v: unrealistic OFV on attempt 1/4: OFV increased vs parent
#> (443.426 > 442.679).
#> ℹ Retrying with perturbed init.
#> Warning: ! sex_female ~ v: unrealistic OFV on attempt 2/4: OFV increased vs parent
#> (443.212 > 442.679).
#> ℹ Retrying with small init.
#> Warning: ! sex_female ~ v: unrealistic OFV on attempt 3/4: OFV increased vs parent
#> (442.88 > 442.679).
#> ℹ Retrying with perturbed init.
#> 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
#> 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.398 -2.281 0.1309 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 445.844 -3.269 0.0706 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.398 -2.281 0.1309 Not selected
#> Forward 3 sex_female~cl [cat] 442.679 441.316 -1.362 0.2431 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 RetainedrunSCM() returns an object of class
"nlmixr2scm". Printing it 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]
#> Saved to : /path/to/project/fit_base_scm_1
#> Use summary() for options, model comparison and step tables.summary() gives the full human-readable report: the
options used for the search, a comparison of the base, final forward and
final backward models, the step and all-candidate tables, and where
models and report files were saved on disk:
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 = 1; rxThreads per worker = 2
#> Save models : yes (restarted)
#> ── 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.398 -2.281 0.1309 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 445.844 -3.269 0.0706 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.398 -2.281 0.1309 Not selected
#> Forward 3 sex_female~cl [cat] 442.679 441.316 -1.362 0.2431 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 ───────────────────────────────────────────────────────────────────────
#> Output directory: /path/to/project/fit_base_scm_1
#> Report files : scm_log.txt, scm_step_summary.csv, scm_all_candidates.csv
#> Saved models : 9 files (.rds)Underneath, the result is a named list:
summaryTable contains one row per candidate tested
across all steps and both phases, with the columns most useful for
review:
cols <- intersect(
c(
"searchType", "step", "covar", "var", "shape",
"deltObjf", "pchisqr", "included"
),
colnames(scm$summaryTable)
)
print(scm$summaryTable[, cols])
#> searchType step covar var shape deltObjf
#> cov_wt_power_cl forward 1 wt_power cl power -4.8866928
#> cov_wt_lin_cl forward 1 wt_lin cl lin -4.3775934
#> cov_wt_power_v forward 1 wt_power v power -25.9764236
#> cov_wt_lin_v forward 1 wt_lin v lin -25.0795749
#> cov_sex_female_cl forward 1 sex_female cl cat -0.9167855
#> cov_sex_female_v forward 1 sex_female v cat -16.6663280
#> cov_wt_power_cl1 forward 2 wt_power cl power -6.4342868
#> cov_wt_lin_cl1 forward 2 wt_lin cl lin -4.9629385
#> cov_sex_female_cl1 forward 2 sex_female cl cat -0.1294283
#> cov_sex_female_v1 forward 2 sex_female v cat -3.2688731
#> cov_sex_female_cl2 forward 3 sex_female cl cat -1.3624840
#> cov_sex_female_v2 forward 3 sex_female v cat -2.2812182
#> cov_wt_power_cl2 backward 1 wt_power cl power 5.6190466
#> cov_wt_power_v2 backward 1 wt_power v power 27.4944709
#> cov_wt_power_v1 backward 2 wt_power v power 26.2227916
#> pchisqr included
#> cov_wt_power_cl 2.706448e-02 no
#> cov_wt_lin_cl 3.641438e-02 no
#> cov_wt_power_v 3.456124e-07 yes
#> cov_wt_lin_v 5.501249e-07 no
#> cov_sex_female_cl 3.383204e-01 no
#> cov_sex_female_v 4.456505e-05 no
#> cov_wt_power_cl1 1.119381e-02 yes
#> cov_wt_lin_cl1 2.589617e-02 no
#> cov_sex_female_cl1 7.190256e-01 no
#> cov_sex_female_v1 7.060611e-02 no
#> cov_sex_female_cl2 2.431075e-01 no
#> cov_sex_female_v2 1.309490e-01 no
#> cov_wt_power_cl2 1.776631e-02 dropped
#> cov_wt_power_v2 1.575443e-07 retained
#> cov_wt_power_v1 3.042108e-07 retainedKey columns:
| Column | Meaning |
|---|---|
searchType |
"forward" or "backward" |
step |
Step number within the phase |
covar |
Covariate name |
var |
PK parameter |
shape |
Functional shape tested |
deltObjf |
ΔOFV (candidate − reference: negative = improvement in forward, positive = deterioration in backward) |
pchisqr |
p-value from chi-squared test (1 df) |
bsvReduction |
Reduction in BSV for var (%) |
covarEffect |
Estimated effect magnitude at the covariate extremes |
included |
"yes" (added in forward), "no" (rejected
in forward), "dropped" (removed in backward),
"retained" (tested in backward but kept in model) |
fwd_tbl <- scm$resFwd[[2]]
if (!is.null(fwd_tbl) && nrow(fwd_tbl) > 0) {
print(fwd_tbl[
fwd_tbl$included == "yes",
intersect(c(
"step", "covar", "var", "shape", "deltObjf",
"pchisqr", "bsvReduction"
), colnames(fwd_tbl))
])
} else {
cat("No covariates accepted in forward search.\n")
}fwd_tbl
#> step covar var shape objf deltObjf AIC
#> cov_wt_power_cl 1 wt_power cl power 470.2028 -4.8866928 947.5099
#> cov_wt_lin_cl 1 wt_lin cl lin 470.7119 -4.3775934 948.0190
#> cov_wt_power_v 1 wt_power v power 449.1130 -25.9764236 926.4202
#> cov_wt_lin_v 1 wt_lin v lin 450.0099 -25.0795749 927.3170
#> cov_sex_female_cl 1 sex_female cl cat 474.1727 -0.9167855 951.4798
#> cov_sex_female_v 1 sex_female v cat 458.4231 -16.6663280 935.7303
#> cov_wt_power_cl1 2 wt_power cl power 442.6788 -6.4342868 921.9859
#> cov_wt_lin_cl1 2 wt_lin cl lin 444.1501 -4.9629385 923.4573
#> cov_sex_female_cl1 2 sex_female cl cat 448.9836 -0.1294283 928.2908
#> cov_sex_female_v1 2 sex_female v cat 445.8442 -3.2688731 925.1513
#> cov_sex_female_cl2 3 sex_female cl cat 441.3163 -1.3624840 922.6234
#> cov_sex_female_v2 3 sex_female v cat 440.3975 -2.2812182 921.7047
#> BIC numParams qchisqr pchisqr included searchType
#> cov_wt_power_cl 975.7135 8 3.841459 2.706448e-02 no forward
#> cov_wt_lin_cl 976.2226 8 3.841459 3.641438e-02 no forward
#> cov_wt_power_v 954.6238 8 3.841459 3.456124e-07 yes forward
#> cov_wt_lin_v 955.5207 8 3.841459 5.501249e-07 no forward
#> cov_sex_female_cl 979.6835 8 3.841459 3.383204e-01 no forward
#> cov_sex_female_v 963.9339 8 3.841459 4.456505e-05 no forward
#> cov_wt_power_cl1 953.7150 9 3.841459 1.119381e-02 yes forward
#> cov_wt_lin_cl1 955.1863 9 3.841459 2.589617e-02 no forward
#> cov_sex_female_cl1 960.0198 9 3.841459 7.190256e-01 no forward
#> cov_sex_female_v1 956.8804 9 3.841459 7.060611e-02 no forward
#> cov_sex_female_cl2 957.8780 10 3.841459 2.431075e-01 no forward
#> cov_sex_female_v2 956.9592 10 3.841459 1.309490e-01 no forward
#> covNames covarEffect bsvReduction
#> cov_wt_power_cl cov_wt_power_cl 0.527569356 13.4824367
#> cov_wt_lin_cl cov_wt_lin_cl 0.007405764 13.8607798
#> cov_wt_power_v cov_wt_power_v 0.881441524 61.7460287
#> cov_wt_lin_v cov_wt_lin_v 0.012441116 66.9813363
#> cov_sex_female_cl cov_sex_female_cl -0.107372644 0.1758785
#> cov_sex_female_v cov_sex_female_v -0.400863506 18.3482567
#> cov_wt_power_cl1 cov_wt_power_cl 0.584437702 17.2747592
#> cov_wt_lin_cl1 cov_wt_lin_cl 0.004636182 12.8092199
#> cov_sex_female_cl1 cov_sex_female_cl -0.104984340 -5.5881667
#> cov_sex_female_v1 cov_sex_female_v -0.128347020 55.8389186
#> cov_sex_female_cl2 cov_sex_female_cl 0.144149022 0.4603984
#> cov_sex_female_v2 cov_sex_female_v -0.164344744 28.0074901bck_tbl <- scm$resBck[[2]]
if (!is.null(bck_tbl) && nrow(bck_tbl) > 0) {
print(bck_tbl[, intersect(c(
"step", "covar", "var", "shape", "deltObjf",
"pchisqr", "included"
), colnames(bck_tbl))])
} else {
cat("No backward elimination steps.\n")
}bck_tbl
#> step covar var shape objf deltObjf AIC BIC
#> cov_wt_power_cl 1 wt_power cl power 448.2978 5.619047 925.6050 953.8086
#> cov_wt_power_v 1 wt_power v power 470.1732 27.494471 947.4804 975.6840
#> cov_wt_power_v1 2 wt_power v power 474.5206 26.222792 949.8277 974.5059
#> numParams qchisqr pchisqr included searchType
#> cov_wt_power_cl 8 6.634897 1.776631e-02 dropped backward
#> cov_wt_power_v 8 6.634897 1.575443e-07 retained backward
#> cov_wt_power_v1 7 6.634897 3.042108e-07 retained backward
#> covNames covarEffect bsvReduction
#> cov_wt_power_cl cov_wt_power_cl 0.5844377 20.15977
#> cov_wt_power_v cov_wt_power_v 0.8895159 83.20203
#> cov_wt_power_v1 cov_wt_power_v 0.8686918 76.37747The final model fit is in resBck[[1]] after a full SCM
(or resFwd[[1]] after a forward-only search):
fit_final <- scm$resBck[[1]]
fit_final$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
cat(sprintf(
"Final OFV: %.3f ΔOFV vs base: %.3f\n",
fit_final$objf, fit_base$objf - fit_final$objf
))
#> Final OFV: 448.298 ΔOFV vs base: 26.792includedRelations forces specific covariate
relationships into the model at the start of backward elimination — even
if they were not accepted during forward inclusion. This mirrors PsN’s
[included_relations] block:
The wt~cl power relationship will be present at the
start of backward elimination regardless of whether it was significant
in the forward phase. It remains eligible for removal during backward
elimination like any other covariate.
Set searchType to run a single phase:
With searchType = "forward", the run stops after the
inclusion phase and the forward-final model is returned in
resFwd[[1]]. With searchType = "backward",
runSCM() applies elimination to a model that already
contains covariate terms.
Each candidate model at a given step is an independent fit, so
parallelisation gives a near-linear speed-up up to the number of
candidates. Install the optional future,
future.apply, and progressr packages and set
workers to the desired number of parallel processes.
Sequential execution (workers = NULL or
workers = 1) requires no additional packages.
Note: parallel SCM workers load from the installed package. If you are developing with
devtools::load_all(), install the package first withdevtools::install()before usingworkers > 1.
When a covariate has missing observations, set
missingToken to the sentinel used in the data
(e.g. -99 or ".") in addition to
NA. runSCM() wraps the covariate expression in
an ifelse() guard that imputes the population-typical
value:
cov = median(cov). This gives
0 for "power"
(log(median(cov)/median(cov)) = 0) and "lin"
(median(cov)-median(cov)= 0), log(median(cov))
for "log", and the median(cov) itself for
"identity".When saveModels = TRUE (the default),
runSCM() writes the following files to the output directory
(auto-named <fitName>_scm_<N> in the working
directory).
Two CSV summaries and one log file are written at the end of the search:
| File | Contents |
|---|---|
scm_step_summary.csv |
Best candidate per step (one row per step) |
scm_all_candidates.csv |
Every candidate model tested, across all steps and both phases |
scm_log.txt |
Human-readable run log with header metadata |
Three .rds files are written per accepted
step (not per candidate). The <key> segment
is <covar>_<var> for the covariate that was
added (forward) or removed (backward) at that step:
| File | Contents |
|---|---|
forward_step_<N>_fit_<key>.rds |
Fit object for the model accepted at forward step N |
forward_step_<N>_table_<key>.rds |
One-row data frame for the accepted best candidate |
forward_step_<N>_completetable_<key>.rds |
All candidates tested across forward steps 1..N (cumulative) |
backward_step_<N>_fit_<key>.rds |
Fit object after the covariate is removed at backward step N |
backward_step_<N>_table_<key>.rds |
One-row data frame for the removed covariate |
backward_step_<N>_completetable_<key>.rds |
All candidates tested across backward steps 1..N (cumulative) |
Only the accepted (best-at-step) candidate is persisted to
.rds; per-candidate fit results are not written
individually but are available in scm_all_candidates.csv
and in the cumulative _completetable_ snapshots.
The
_table_and_completetable_.rdsfiles duplicate information in the CSVs but at full numeric precision and as per-step checkpoints (useful for resume/inspection mid-run). Routine post-hoc analysis can rely on the CSVs alone.
Set saveModels = FALSE to run without writing any files,
for example when exploring the search space in a scratch session.
Sample size: the SCM uses a chi-squared approximation with 1 degree of freedom. The rule of thumb is at least 50–100 subjects for a reliable test; sparse data leads to inflated type I error in forward inclusion.
p-value thresholds: the conventional PsN defaults
are fwd = 0.05 and bck = 0.01. In exploratory
analyses, a more liberal forward threshold (e.g. 0.10) avoids missing
important relationships at the cost of higher false-positive rates.
Shape choice: the power shape
(log(cov/median)) is appropriate for weight and other
allometric predictors. Use shapes = c("power", "lin") to
test both and let the data decide.
Warm-starting: when a covariate is accepted, the estimated theta from that step is automatically used as the starting estimate for all subsequent steps that involve the same covariate shape, which speeds convergence and improves numerical stability.
Stalled candidates: with ODE models, solver noise
can flatten the outer objective enough that a candidate’s fit never
leaves its zero-effect initial estimate. By default
(profileInitOnStall = TRUE), a candidate whose OFV
improvement over its parent is <= stallTol is
automatically rescued with a one-shot frozen 1-D profile (Brent method)
that supplies a gradient-informative starting value, then refit. Set
profileInit = TRUE to warm-start every forward candidate
this way rather than only stalled ones.
Unrealistic OFVs: runSCM() retries a
candidate (up to maxRetries, default 3) when its fit
produces an implausible OFV — using a perturbed or near-zero covariate
init on alternate attempts — before falling back to the best attempt
seen (or, with retryFailOnExhaustion = TRUE, excluding the
candidate as failed). Stochastic estimators (SAEM) get a wider tolerance
automatically to avoid spurious retries from Monte Carlo noise. A
p-value that underflows to zero also counts as implausible by default,
which means genuinely strong effects (an OFV drop above about 70 for one
parameter) are retried too; set retryOnUnderflow = FALSE to
skip those extra fits.
Categorical covariates: pass catvarsVec
rather than pre-creating dummy columns. runSCM() handles
level detection, reference-level selection, and indicator column
creation automatically.
Jonsson, E.N. & Karlsson, M.O. (1998). Automated covariate model building within NONMEM. Pharmaceutical Research, 15(9), 1463–1468.
Lindbom, L., Ribbing, J., & Jonsson, E.N. (2004). Perl-speaks-NONMEM (PsN) — a Perl module for NONMEM related programming. Computer Methods and Programs in Biomedicine, 75, 85–94.
Ribbing, J. & Jonsson, E.N. (2004). Power, selection bias and predictive performance of the population pharmacokinetic covariate model. Journal of Pharmacokinetics and Pharmacodynamics, 31(2), 109–134.