--- title: "Stepwise Covariate Modelling with runSCM()" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Stepwise Covariate Modelling with runSCM()} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ## Introduction 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: 1. **Forward inclusion**: starting from the base model, the candidate relationship that produces the largest statistically significant ΔOFV (OFV decrease) is added at each step. The process repeats until no remaining candidate meets the threshold. 2. **Backward elimination**: starting from the forward-final model, each included relationship is tested for removal. A relationship is dropped if its removal does not produce a statistically significant ΔOFV increase (at a stricter threshold than forward). `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. --- ## Example: warfarin PK covariate search We use `nlmixr2data::warfarin`, filtered to PK observations only (`dvid == "cp"`), with covariates body weight (`wt`, continuous) and sex (`sex`, categorical). ``` r library(nlmixr2data) library(nlmixr2utils) library(nlmixr2scm) ``` ### 1. Prepare the data 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. ``` r 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=27 ``` ### 2. Fit the base model The base model contains no covariate terms — `runSCM()` adds them. ``` r 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) ) ``` ``` r 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.089 ``` --- ### 3. Candidate specification strategies `runSCM()` 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 | --- ### 4. Covariate shapes `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_` 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. --- ### 5. Run the SCM 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. ``` r 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 Retained ``` --- ### 6. Interpreting the results `runSCM()` returns an object of class `"nlmixr2scm"`. Printing it gives a short overview: ``` r 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: ``` r 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: ``` r names(scm) #> [1] "summaryTable" "resFwd" "resBck" "baseFit" "options" #> [6] "outputDir" ``` #### Summary table `summaryTable` contains one row per candidate tested across all steps and both phases, with the columns most useful for review: ``` r 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 retained ``` Key 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) | #### Forward search detail ``` r 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") } ``` ``` r 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.0074901 ``` #### Backward search detail ``` r bck_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") } ``` ``` r 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.37747 ``` #### Final model The final model fit is in `resBck[[1]]` after a full SCM (or `resFwd[[1]]` after a forward-only search): ``` r 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.792 ``` --- ### 7. Forced inclusions `includedRelations` 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. --- ### 8. Forward-only or backward-only search 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. --- ### 9. Parallel candidate fitting 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 with > `devtools::install()` before using `workers > 1`. --- ### 10. Missing covariate values 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: * **Continuous**: imputes the median by evaluating the shape expression at `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"`. * **Categorical**: imputes the mode (most frequent level). --- ## Output files When `saveModels = TRUE` (the default), `runSCM()` writes the following files to the output directory (auto-named `_scm_` 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 `` segment is `_` for the covariate that was added (forward) or removed (backward) at that step: | File | Contents | |---|---| | `forward_step__fit_.rds` | Fit object for the model accepted at forward step N | | `forward_step__table_.rds` | One-row data frame for the accepted best candidate | | `forward_step__completetable_.rds` | All candidates tested across forward steps 1..N (cumulative) | | `backward_step__fit_.rds` | Fit object after the covariate is removed at backward step N | | `backward_step__table_.rds` | One-row data frame for the removed covariate | | `backward_step__completetable_.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_` `.rds` files 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. --- ## Tips for a robust SCM **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. --- ## References 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.