Covariate Models

The purpose of this vignette is to demonstrate how to work with covariate models using the Certara.RsNLME package, for both built-in and textual PML models.

In this vignette you will:

• Build a base model and add covariates
• Add covariate effects and run Stepwise Covariate Modeling (SCM), then review results
• Convert to a textual model and edit the PML to include/exclude covariate effects
• Use the PML enable statement to toggle covariate effects
• Run an SCM search on the textual model
• Investigate covariate inflation on different structural parameters with a Shotgun covariate search

library(dplyr)
library(Certara.RsNLME)

Build base model

We’ll begin by creating a built-in model with the pkmodel() function to illustrate how covariates are added. Since this is the base model, we will introduce covariates without setting any covariate effects at this time. This approach is useful to include covariate columns in the model outputs and visually explore potential covariate–parameter relationships to inform subsequent covariate modeling.

baseModel <- pkmodel(columnMap = FALSE, numCompartments = 2, modelName = "baseModel") %>%  
  dataMapping(pkData) %>%
  colMapping(c(id = "Subject", time = "Act_Time", CObs = "Conc", A1 = "Amount"))

Now, we’ll add the covariates from the dataset to the base model.

head(pkData)
#>    Subject Nom_Time Act_Time Amount  Conc   Age BodyWeight Gender
#>      <num>    <num>    <num>  <num> <num> <num>      <num> <char>
#> 1:       1        0     0.00  25000  2010    22         73   male
#> 2:       1        1     0.26     NA  1330    22         73   male
#> 3:       1        1     1.10     NA   565    22         73   male
#> 4:       1        2     2.10     NA   216    22         73   male
#> 5:       1        4     4.13     NA   180    22         73   male
#> 6:       1        8     8.17     NA   120    22         73   male

We’ll use the addCovariate() function to add covariates to the base model. For continuous covariates, set type = "Continuous". For categorical covariates, set type = "Categorical" and provide the appropriate levels and labels based on the unique values in the dataset.

Note that the labels argument is only required if the categorical covariate is a character column in the dataset.

unique(pkData$Gender)
#> [1] "male"   "female"

Add covariates to the base model, no effects on structural parameters at this time.

baseModel <- baseModel %>%
  addCovariate(covariate = "Age", type = "Continuous", direction = "Forward") %>%
  addCovariate(covariate = "BodyWeight", type = "Continuous", direction = "Forward") %>% 
  addCovariate(covariate = "Gender",
               type = "Categorical",
               levels=c(0,1),
               labels=c("male","female"),
               direction = "Forward")

In the above covariate argument we are setting the covariate name in the model/PML to be the same as the column name in the dataset.

When this is done, you do not need to map the column names in the dataset to the model variables in the colMapping() function afterwards.

If instead you used addCovariate(covariate = "BWT", type = "Continuous") for example, you would need to additionally supply the function command colMapping(c(BWT = "BodyWeight")) to map the model variable BWT to the column name BodyWeight in the pkData dataset.

An additional comment on the direction argument in addCovariate(): it controls how missing covariate values are propagated in the data. By default, direction = "Forward" carries missing values forward to subsequent observations. Other options are "Backward" (carry values backward to previous observations) and "Interpolate" (interpolate values between observations; applicable only for continuous covariates).

• direction = "Forward" is equivalent to PML code fcovariate(CovName)

• direction = "Backward" is equivalent to PML code covariate(CovName)

• direction = "Interpolate" is equivalent to PML code interpolate(CovName).

Now view the model and observe the fcovariate() statements added to the PML for each covariate.

    fcovariate(Age)
    fcovariate(BodyWeight)
    fcovariate(Gender())

Note: The PML syntax to indicate a categorical covariate is to include () after the covariate name (e.g., fcovariate(Gender())).

print(baseModel)

Next, we’ll fit the base model and inspect the results.

baseResult <- fitmodel(baseModel)
#> Using MPI host with 4 cores
#> 
#> NLME Job
#> sharedDirectory is not given in the host class instance.
#> Valid NLME_ROOT_DIRECTORY is not given, using current working directory:
#> C:/Repos/R-RsNLME/vignettes
#> TDL5 version: 25.7.1.1
#> 
#> Status: OK
#> License expires: 2026-11-14
#> Refresh until: 2025-12-14 09:53:10
#> Current Date: 2025-11-14
#> The model compiled
#> 
#> Trying to generate job results...
#> Done generating job results.

baseResult$Overall
#>    Scenario RetCode    LogLik     -2LL      AIC      BIC nParm  nObs  nSub
#>      <char>   <int>     <num>    <num>    <num>    <num> <int> <int> <int>
#> 1: WorkFlow       1 -632.7953 1265.591 1283.591 1308.057     9   112    16
#>    EpsShrinkage Condition
#>           <num>     <num>
#> 1:      0.17312   4.86919
baseResult$theta
#>    Scenario Parameter  Estimate  Units     Stderr       CV%     2.5%CI
#>      <char>    <char>     <num> <lgcl>      <num>     <num>      <num>
#> 1: WorkFlow       tvV 15.393731     NA 1.13948176  7.402246 13.1338373
#> 2: WorkFlow      tvCl  6.599838     NA 0.72057416 10.918058  5.1707492
#> 3: WorkFlow      tvV2 41.190009     NA 1.06233995  2.579120 39.0831087
#> 4: WorkFlow     tvCl2 14.025770     NA 1.03509504  7.379952 11.9729030
#> 5: WorkFlow      CEps  0.161206     NA 0.01964795 12.188101  0.1222389
#>       97.5%CI Var.Inf.factor
#>         <num>         <lgcl>
#> 1: 17.6536240             NA
#> 2:  8.0289269             NA
#> 3: 43.2969102             NA
#> 4: 16.0786368             NA
#> 5:  0.2001731             NA
baseResult$omega
#>     Label         nV       nCl          nV2       nCl2
#>    <char>      <num>     <num>        <num>      <num>
#> 1:     nV 0.06914959 0.0000000 0.000000e+00 0.00000000
#> 2:    nCl 0.00000000 0.1806635 0.000000e+00 0.00000000
#> 3:    nV2 0.00000000 0.0000000 4.649898e-08 0.00000000
#> 4:   nCl2 0.00000000 0.0000000 0.000000e+00 0.04239193

Generate a few goodness-of-fit plots to assess the base model across the three covariates we introduced in the model.

library(Certara.Xpose.NLME)

baseModel_xpdb <- xposeNlmeModel(baseModel, baseResult)

eta_vs_cov(baseModel_xpdb, type = "ps", covariate = "Age")

eta_vs_cov(baseModel_xpdb, type = "ps", covariate = "BodyWeight") 

eta_vs_cov(baseModel_xpdb, type = "ps", covariate = "Gender") 

Note: You can also use the function xposeNlme() to generate the xpose_data object from the run output directory, which is useful if you want to regenerate or modify the model diagnostics without re-running the model.

There appears to be a clear effect of BodyWeight on Cl and V. Other covariates are less compelling; we will further explore potential effects with Stepwise Covariate Modeling.

Create a covariate model from the base model

Copy the base model, accepting the parameter estimates from the initial run.

scmModel <-  copyModel(baseModel, acceptAllEffects = TRUE, modelName = "scmModel")
 

Next, add the covariate effects that we want to search for. We’ll use the same addCovariate() function as we did for the base model, but this time specify the effect argument to include the names of the structural parameters we want to include in the search.

We’ll search for effects on the V and Cl parameters.

scmModel <-  scmModel %>% 
  addCovariate("BodyWeight", 
               effect=(c("V","Cl"))) %>%  
  addCovariate("Gender", 
               type = "Categorical",
               effect=c("V","Cl"),
               levels=c(0,1),
               labels=c("male","female")) %>% 
  addCovariate("Age", 
               effect=c("V","Cl")) 
 
 

What do “covariate effects” mean?

When you call addCovariate(..., effect = c("V","Cl")), you’re telling RsNLME to link a covariate (e.g., BodyWeight, Age, Gender) to one or more structural parameters (here V and Cl). Internally, RsNLME creates fixed-effect terms in the PML that modify each linked parameter:

• Continuous covariates (e.g., BodyWeight, Age) add one fixed-effect coefficient per parameter (e.g., dVdBodyWeight, dCldBodyWeight).
• Categorical covariates (e.g., Gender with two levels) add number of levels − 1 fixed-effect coefficients per parameter (e.g., dVdGender1, dCldGender1).

These effect terms are what the stepwise/shotgun search will switch on/off in different scenarios.

How RsNLME translates effects into the model (options & styles)

The exact link used depends on the parameter’s style (e.g., LogNormal vs Normal) and the option you choose in addCovariate():

• option = "Yes" (default): standard link
  – multiplicative for LogNormal style; additive for Normal style.
• option = "PlusOne": “1 + effect” formulation; only when the affected parameter has style = "LogNormal"; not for Occasion covariates.
• option = "No": remove the specified covariate→parameter link (the covariate declaration remains).

Direction & centering affect how covariate values are propagated and scaled, not which parameters they link to:

• direction = "Backward" | "Forward" | "Interpolate" map to PML covariate(...), fcovariate(...), and interpolate(...).

• center = "Mean" | "Median" | "Value" | "None" (continuous only) controls centering; centerValue is required when center = "Value".

Tip: Keep categorical levels/labels explicit so the reference group is stable and generated coefficients (e.g., ...Gender1) are interpretable.

Worked examples

A) Continuous covariate on LogNormal parameters (standard link)

Suppose V and Cl are LogNormal. With the default "Yes" option, the standard multiplicative link is used:

Model_BW_on_V <- baseModel %>%
  addCovariate("BodyWeight", effect = c("V"),
               type = "Continuous", center = "Median")

Conceptually (illustrative PML):

covariate(BodyWeight)
stparm(V = tvV * ((BodyWeight/median(BodyWeight))^dVdBodyWeight) * exp(nV))

In RsNLME, LogNormal parameters combine covariate effects multiplicatively and add the random effect as * exp(η).

Tip: Since covariate direction is not given, the same direction is used if the covariate is already in the model.

B) Continuous covariate with “PlusOne” (LogNormal-only)

Model_Age_on_Cl <- baseModel %>%
  addCovariate("Age", effect = c("Cl"),
               type = "Continuous", option = "PlusOne",
               direction = "Backward", center = "Mean")

covariate(Age)
stparm(Cl = tvCl * ( (1+(Age-mean(Age))*dCldAge))  * exp(nCl))

Note: option="PlusOne" is intended for LogNormal parameters (and not for Occasion covariates). For non-LogNormal styles, use the standard link.

C) Two-level categorical covariate (Gender)

ModelGender <- baseModel %>%
  addCovariate("Gender", effect = c("Cl"),
               type = "Categorical", levels = c(0,1), labels = c("male","female"))

For a LogNormal parameter with option="Yes", this behaves like:

fcovariate(Gender())
stparm(Cl = tvCl * exp(dCldGender1*(Gender==1)) * exp(nCl))

With option="PlusOne", the categorical multiplier becomes

fcovariate(Gender())
stparm(Cl = tvCl * ( (1+dCldGender1*(Gender==1)))  * exp(nCl))

D) Three-level categorical covariate (example)

ModelRace <- baseModel %>%
  addCovariate("Race", effect = "V",
               type = "Categorical", levels = c(0,1,2),
               labels = c("White","Black","Asian"))

This adds 2 fixed effects on V: dVdRace1, dVdRace2 (one per non-reference level).

    fcovariate(Race())
    fixef( dVdRace1(enable=c(0)) = c(,0,))
    fixef( dVdRace2(enable=c(0)) = c(,0,))
    stparm(V = tvV * exp(dVdRace1*(Race==1)) * exp(dVdRace2*(Race==2)) * exp(nV))

Note: The "Race" column does not exist in the built-in pkData; this example illustrates usage of a categorical covariate effect with three levels.

Non-LogNormal styles: how covariate links look

For non-LogNormal parameter styles, covariate effects are additive inside a sum, and the random effect is added inside that sum (or applied after the sum depending on style), not as * exp(η) unless specified by the style.

Tip: Change a parameter’s style with structuralParameter(), then add covariates. You can also change style after; RsNLME will rebuild the link.

E) Continuous covariate on a Normal parameter (additive link)

Model_V_Normal <- baseModel %>%
  structuralParameter("V", style = "Normal") %>% # switch V to Normal
  addCovariate("BodyWeight", effect = "V", direction = "Backward",
               type = "Continuous", center = "Mean")

Illustrative PML (Normal → additive):

covariate(BodyWeight)
stparm(V = tvV + ((BodyWeight-mean(BodyWeight))*dVdBodyWeight) + nV)

F) Categorical covariate on a Normal parameter

Model_V_Normal_Gender <- baseModel %>%
  structuralParameter("V", style = "Normal") %>%
  addCovariate("Gender", effect = "V",
               type = "Categorical", levels = c(0,1), labels = c("male","female"))

Illustrative PML:

fcovariate(Gender())
stparm(V = tvV + dVdGender1*(Gender==1) + nV)

G) Continuous covariate on Combination (a.k.a. “LogNormal1”)

Combination style forms a sum of fixed and covariate terms, then multiplies by exp(η).

Model_Cl_Combo <- baseModel %>%
  structuralParameter("Cl", style = "LogNormal1") %>% # Combination style
  addCovariate("BodyWeight", effect = "Cl", direction = "Forward",
               type = "Continuous", center = "Mean")

Illustrative PML:

fcovariate(BodyWeight)
stparm(Cl = (tvCl + ((BodyWeight-mean(BodyWeight))*dCldBodyWeight) )  * exp(nCl))

---

H) Continuous covariate on Log (a.k.a. “LogNormal2”)

Log style exponentiates a sum of fixed, covariate, and random terms: exp(sum + η).

Model_V_Log <- baseModel %>%
  structuralParameter("V", style = "LogNormal2") %>% # Log style
  addCovariate("Age", effect = "V", direction = "Backward",
               type = "Continuous", center = "Mean")

Illustrative PML:

covariate(Age)
stparm(V = exp(tvV + ((Age-mean(Age))*dVdAge) + nV))

---

I) Continuous covariate on Logit (a.k.a. “LogitNormal”) — bounded [0,1]

Use this when the parameter is naturally bounded (e.g., a fraction like Fe). The model builds ilogit(sum + η).

Model_Fe_Logit <- 
  pkmodel(columnMap = FALSE, 
          numCompartments = 1, 
          absorption = "FirstOrder", 
          modelName = "FModel", 
          isFractionExcreted = TRUE, 
          isClosedForm = FALSE) %>%
  structuralParameter("Fe", style = "LogitNormal") %>% # only if Fe exists in your PK setup
  addCovariate("Age", effect = "Fe",
               type = "Continuous", center = "Mean")

Illustrative PML:

fcovariate(Age)
stparm(Fe = ilogit(tvFe + ((Age-mean(Age))*dFedAge) + nFe))

J) Removing covariate/covariate effect

# remove covariate Gender completely:
ModelWOGender <-
  removeCovariate(ModelGender, covariate = "Gender")

# Keep Gender declared but remove its effect on Cl
ModelGenderNoEff <- ModelGender %>%
  addCovariate("Gender", effect = "Cl",
               type = "Categorical", option = "No",
               levels = c(0,1), labels = c("male","female"))
#> updating option argument for covariate effect of Gender in structural parameter Cl to 'No'

---

How many degrees of freedom does each effect add?

• Each continuous effect on a parameter contributes 1 df.
• Each categorical effect contributes (number of levels − 1) df on that parameter.

Example (this vignette):

• BodyWeight on V and Cl → 2 df;

• Age on V and Cl → 2 df;

• Gender (2 levels) on V and Cl → 2 df;

total = 6 df.

---

Inspect which effects will be searched

# Create the covariate-effect model object used by Stepwise/Shotgun
cm <- covariateModel(scmModel)

# Quick peek
cm@covariateList        
#> [1] "V-Age,V-BodyWeight,V-Gender,Cl-Age,Cl-BodyWeight,Cl-Gender"

cm@degreesOfFreedom
#> [1] "1,1,1,1,1,1"

Other quick utilities:

# List the covariate names currently declared in the model
covariateNames(scmModel)
#> [1] "Age"        "BodyWeight" "Gender"

Run stepwise search

Running SCM is performed with the stepwiseSearch() function.

Optionally, you can control the search thresholds and the model selection criterion via the stepwiseParams argument:

• addPValue: threshold to add a covariate effect (default 0.01). Interpretation depends on method:
  • If method = "-2LL": p-value cutoff for a likelihood-ratio test; add if p ≤ addPValue.
  • If method = "AIC" or "BIC": minimum required decrease in the chosen criterion; add if the criterion decreases by ≥ addPValue units.
• removePValue: threshold to remove a covariate effect (default 0.001). Interpretation depends on method:
  • If method = "-2LL": p-value cutoff for a likelihood-ratio test; remove if p > removePValue.
  • If method = "AIC" or "BIC": effect can be removed if the criterion does not worsen by more than removePValue units (i.e., ΔAIC/ΔBIC ≤ removePValue).
• method: metric to rank scenarios; one of “-2LL”, “AIC”, “BIC” (default “-2LL”)

If stepwiseParams is not provided, stepwiseSearch() uses these defaults and prints a message similar to:

stepwiseParams argument is not given. 

Using Default: 

Thresholds for adding/removing a covariate effect: 0.01 / 0.001

Method: -2LL

scmResult <- stepwiseSearch(scmModel)
#> 
#> stepwiseParams argument is not given. 
#> Using Default: 
#> Thresholds for adding/removing a covariate effect: 0.01 / 0.001
#> Method: -2LL
#> 
#> covariateModel argument is not given. Processing model object to make it.
#> Stepwise Job
#> 
#> NLME Job
#> sharedDirectory is not given in the host class instance.
#> Valid NLME_ROOT_DIRECTORY is not given, using current working directory:
#> C:/Repos/R-RsNLME/vignettes
#> 
#> Preparing files for Stepwise Covariate Search run
#> 
#> ----------------------------------------------------------------------
#> Processing scenarios:
#> "Base Model, no covariates", " V-Age", " V-BodyWeight", " V-Gender", " Cl-Age", " Cl-BodyWeight", " Cl-Gender"
#> ----------------------------------------------------------------------
#> 
#> Compiling 1 of 1 NLME models
#> TDL5 version: 25.7.1.1
#> 
#> Status: OK
#> License expires: 2026-11-14
#> Refresh until: 2025-12-14 09:53:10
#> Current Date: 2025-11-14
#> The model compiled
#> 
#> Num Jobs/Completed/Failed:7/0/0
#> 
#> Num Jobs/Completed/Failed:7/0/0
#> 
#> Num Jobs/Completed/Failed:7/0/0
#> 
#> Num Jobs/Completed/Failed:7/0/0
#> 
#> Num Jobs/Completed/Failed:7/0/0
#> 
#> Num Jobs/Completed/Failed:7/1/0
#> 
#> Num Jobs/Completed/Failed:7/3/0
#> 
#> Num Jobs/Completed/Failed:7/4/0
#> 
#> Num Jobs/Completed/Failed:7/4/0
#> 
#> Num Jobs/Completed/Failed:7/4/0
#> 
#> Num Jobs/Completed/Failed:7/5/0
#> 
#> Num Jobs/Completed/Failed:7/6/0
#> 
#> Num Jobs/Completed/Failed:7/6/0
#> 
#> Num Jobs/Completed/Failed:7/6/0
#> 
#> Num Jobs/Completed/Failed:7/7/0
#> 
#> Trying to generate job results...
#> Done generating job results.
#> Find effect to add that reduces -2LL the most:
#> cstep001  V-Age 100000   X 1272.143217 ( 1265.508320 +6.634897 ) > 1265.590400 )
#> cstep002  V-BodyWeight 010000     1255.631877 ( 1248.996980 +6.634897 ) < 1265.590400 )
#> cstep003  V-Gender 001000   X 1271.990837 ( 1265.355940 +6.634897 ) > 1265.590400 )
#> cstep004  Cl-Age 000100   X 1271.969037 ( 1265.334140 +6.634897 ) > 1265.590400 )
#> cstep005  Cl-BodyWeight 000010     1241.226597 ( 1234.591700 +6.634897 ) < 1265.590400 )
#> cstep006  Cl-Gender 000001   X 1271.292697 ( 1264.657800 +6.634897 ) > 1265.590400 )
#> cstep005  Cl-BodyWeight 000010 chosen, -2LL = 1234.591700
#> ----------------------------------------------------------------------
#> Processing scenarios:
#> " V-Age Cl-BodyWeight", " V-BodyWeight Cl-BodyWeight", " V-Gender Cl-BodyWeight", " Cl-Age Cl-BodyWeight", " Cl-BodyWeight Cl-Gender"
#> ----------------------------------------------------------------------
#> 
#> Compiling 1 of 1 NLME models
#> TDL5 version: 25.7.1.1
#> 
#> Status: OK
#> License expires: 2026-11-14
#> Refresh until: 2025-12-14 09:53:10
#> Current Date: 2025-11-14
#> The model compiled
#> 
#> Num Jobs/Completed/Failed:5/0/0
#> 
#> Num Jobs/Completed/Failed:5/0/0
#> 
#> Num Jobs/Completed/Failed:5/0/0
#> 
#> Num Jobs/Completed/Failed:5/0/0
#> 
#> Num Jobs/Completed/Failed:5/0/0
#> 
#> Num Jobs/Completed/Failed:5/0/0
#> 
#> Num Jobs/Completed/Failed:5/2/0
#> 
#> Num Jobs/Completed/Failed:5/3/0
#> 
#> Num Jobs/Completed/Failed:5/4/0
#> 
#> Num Jobs/Completed/Failed:5/4/0
#> 
#> Num Jobs/Completed/Failed:5/4/0
#> 
#> Num Jobs/Completed/Failed:5/4/0
#> 
#> Num Jobs/Completed/Failed:5/5/0
#> 
#> Trying to generate job results...
#> Done generating job results.
#> Find effect to add that reduces -2LL the most:
#> cstep007  V-Age Cl-BodyWeight 100010   X 1241.393497 ( 1234.758600 +6.634897 ) > 1234.591700 )
#> cstep008  V-BodyWeight Cl-BodyWeight 010010     1223.942897 ( 1217.308000 +6.634897 ) < 1234.591700 )
#> cstep009  V-Gender Cl-BodyWeight 001010   X 1241.003937 ( 1234.369040 +6.634897 ) > 1234.591700 )
#> cstep010  Cl-Age Cl-BodyWeight 000110   X 1238.548457 ( 1231.913560 +6.634897 ) > 1234.591700 )
#> cstep011  Cl-BodyWeight Cl-Gender 000011   X 1239.541237 ( 1232.906340 +6.634897 ) > 1234.591700 )
#> cstep008  V-BodyWeight Cl-BodyWeight 010010 chosen, -2LL = 1217.308000
#> ----------------------------------------------------------------------
#> Processing scenarios:
#> " V-Age V-BodyWeight Cl-BodyWeight", " V-BodyWeight V-Gender Cl-BodyWeight", " V-BodyWeight Cl-Age Cl-BodyWeight", " V-BodyWeight Cl-BodyWeight Cl-Gender"
#> ----------------------------------------------------------------------
#> 
#> Compiling 1 of 1 NLME models
#> TDL5 version: 25.7.1.1
#> 
#> Status: OK
#> License expires: 2026-11-14
#> Refresh until: 2025-12-14 09:53:10
#> Current Date: 2025-11-14
#> The model compiled
#> 
#> Num Jobs/Completed/Failed:4/0/0
#> 
#> Num Jobs/Completed/Failed:4/0/0
#> 
#> Num Jobs/Completed/Failed:4/0/0
#> 
#> Num Jobs/Completed/Failed:4/0/0
#> 
#> Num Jobs/Completed/Failed:4/0/0
#> 
#> Num Jobs/Completed/Failed:4/0/0
#> 
#> Num Jobs/Completed/Failed:4/0/0
#> 
#> Num Jobs/Completed/Failed:4/2/0
#> 
#> Num Jobs/Completed/Failed:4/4/0
#> 
#> Trying to generate job results...
#> Done generating job results.
#> Find effect to add that reduces -2LL the most:
#> cstep012  V-Age V-BodyWeight Cl-BodyWeight 110010   X 1225.493677 ( 1218.858780 +6.634897 ) > 1217.308000 )
#> cstep013  V-BodyWeight V-Gender Cl-BodyWeight 011010   X 1223.921317 ( 1217.286420 +6.634897 ) > 1217.308000 )
#> cstep014  V-BodyWeight Cl-Age Cl-BodyWeight 010110   X 1226.756817 ( 1220.121920 +6.634897 ) > 1217.308000 )
#> cstep015  V-BodyWeight Cl-BodyWeight Cl-Gender 010011   X 1222.227977 ( 1215.593080 +6.634897 ) > 1217.308000 )
#> 
#> No effect chosen to add
#> Find effect to subtract that increases -2LL the least:
#> cstep005  Cl-BodyWeight 000010   X 1223.764134 ( 1234.591700 -10.827566 ) > 1217.308000 )
#> cstep002  V-BodyWeight 010000   X 1238.169414 ( 1248.996980 -10.827566 ) > 1217.308000 )
#> 
#> No effect chosen to subtract
#> 
#> Scenario to use = cstep008  V-BodyWeight Cl-BodyWeight 010010
#> 
#> Summarizing stepwise covariate search results for 16 scenarios
#> 
#> Finished summarizing results. Transferring data and loading the results...

scmResult[scmResult$BestScenario == TRUE,]
#>                               Scenario RetCode   LogLik     -2LL      AIC
#> 9 cstep008  V-BodyWeight Cl-BodyWeight       2 -608.654 1217.308 1239.308
#>        BIC nParm nObs nSub EpsShrinkage  Condition BestScenario
#> 9 1269.211    11  112   16      0.13174 1370156093         TRUE

To customize the thresholds and/or selection method:

# Example: use AIC and adjust thresholds
customParams <- StepwiseParams(addPValue = 0.05,
                               removePValue = 0.01,
                               method = "AIC")
scmResult <- stepwiseSearch(scmModel, stepwiseParams = customParams)

The SCM search has also confirmed the covariate relationships we observed in our initial exploratory plots.

Scenario to use = cstep008  V-BodyWeight Cl-BodyWeight 010010

SCM with a textual model

A textual model is simply a model that is developed/edited directly using PML code, rather than the built-in Certara.RsNLME model functions that we used previously (e.g., pkmodel(), addCovariate()).

Any built-in model can be converted to a textual model and is automatically done so if you directly edit the PML using the editModel() function:

scmModel <- editModel(scmModel)

Note: While you can convert any built-in model to a textual model, you cannot convert a textual model back to a built-in model.

Tip: For interactive syntax checking use the modelTextualUI() function. See here for additional details.

Let’s directly edit the PML statements to always include BodyWeight in our covariate model (i.e., exclude it from the SCM). One way to do this is by removing the enable statement in the corresponding fixef(), and we’ll also update the initial estimate to 1.

Since we’ve removed the enable statement for this fixef(), we’ll also need to update the indices for the remaining enable statements in the other fixef() used for SCM.

Change from:

  fixef( dVdAge       (enable=c(0)) = c(,0,))
  fixef( dVdBodyWeight(enable=c(1)) = c(,0,))
  fixef( dVdGender1   (enable=c(2)) = c(,0,))
  fixef( dCldAge      (enable=c(3)) = c(,0,))
  fixef( dCldBodyWeight(enable=c(4)) = c(,0,))
  fixef( dCldGender1  (enable=c(5)) = c(,0,))

Change to:

  # Forced-in effect (not part of SCM search)
  fixef( dVdBodyWeight = c(,1,))
  
  # SCM candidates (contiguous indices 0..4)
  fixef( dVdAge       (enable=c(0)) = c(,0,))
  fixef( dVdGender1   (enable=c(1)) = c(,0,))
  fixef( dCldAge      (enable=c(2)) = c(,0,))
  fixef( dCldBodyWeight(enable=c(3)) = c(,0,))
  fixef( dCldGender1  (enable=c(4)) = c(,0,))

Note: You can simply comment out the original fixef() statement using #.

The full updated PML code is below:

test(){
    cfMicro(A1,Cl/V, Cl2/V, Cl2/V2)
    dosepoint(A1)
    C = A1 / V
    error(CEps=0.161205984597858)
    observe(CObs=C * ( 1 + CEps))
    stparm(V = tvV * (Age^dVdAge)   * (BodyWeight^dVdBodyWeight)   * exp(dVdGender1*(Gender==1)) * exp(nV))
    stparm(Cl = tvCl * (Age^dCldAge)   * (BodyWeight^dCldBodyWeight)   * exp(dCldGender1*(Gender==1)) * exp(nCl))
    stparm(V2 = tvV2 * exp(nV2))
    stparm(Cl2 = tvCl2 * exp(nCl2))
    fcovariate(Age)
    fcovariate(BodyWeight)
    fcovariate(Gender())
    fixef( tvV = c(,15.393730662165,))
    fixef( tvCl = c(,6.59983805751661,))
    fixef( tvV2 = c(,41.1900094327716,))
    fixef( tvCl2 = c(,14.0257699237193,))
    fixef( dVdAge(enable=c(0)) = c(,0,))
    #fixef( dVdBodyWeight(enable=c(1)) = c(,0,))
    fixef( dVdBodyWeight = c(,1,))
    fixef( dVdGender1(enable=c(1)) = c(,0,))
    fixef( dCldAge(enable=c(2)) = c(,0,))
    fixef( dCldBodyWeight(enable=c(3)) = c(,0,))
    fixef( dCldGender1(enable=c(4)) = c(,0,))
    ranef(diag(nV,nCl,nV2,nCl2) = c(0.0691495878371252,0.180663528781695,4.64989751058228e-08,0.0423919347837654))
}

Tip: PML supports in-line commenting within a statement.

For example, you could have instead written the following inline comment in the PML:

• fixef( dVdBodyWeight/*(enable=c(1)/**/ = c(,1,))

Or alternatively, simply remove a given enable statement entirely from the fixef() rather than commenting it out.

See here for additional details.

Start from a textual model (.mdl file)

There will likely be cases where you simply start with PML (from Phoenix or other model files shared with you) and will not start with a built-in model in Certara.RsNLME. As long as you have saved your PML code in a text file, you can import into R and run SCM.

In the code below, we will write our PML to a file named “test.mdl”, then re-import back into our R session for SCM.

Note: The filename “test.mdl” is arbitrary and can be named anything you’d like, so long as it is a text file. An alternative programmatic way to write the PML from a model to a file is writeLines(paste(scmModel@statements), "test.mdl").

pml <- "test(){
    cfMicro(A1,Cl/V, Cl2/V, Cl2/V2)
    dosepoint(A1)
    C = A1 / V
    error(CEps=0.161205984597858)
    observe(CObs=C * ( 1 + CEps))
    stparm(V = tvV * (Age^dVdAge)   * (BodyWeight^dVdBodyWeight)   * exp(dVdGender1*(Gender==1)) * exp(nV))
    stparm(Cl = tvCl * (Age^dCldAge)   * (BodyWeight^dCldBodyWeight)   * exp(dCldGender1*(Gender==1)) * exp(nCl))
    stparm(V2 = tvV2 * exp(nV2))
    stparm(Cl2 = tvCl2 * exp(nCl2))
    fcovariate(Age)
    fcovariate(BodyWeight)
    fcovariate(Gender())
    fixef( tvV = c(,15.393730662165,))
    fixef( tvCl = c(,6.59983805751661,))
    fixef( tvV2 = c(,41.1900094327716,))
    fixef( tvCl2 = c(,14.0257699237193,))
    fixef( dVdAge(enable=c(0)) = c(,0,))
    #fixef( dVdBodyWeight(enable=c(1)) = c(,0,))
    fixef( dVdBodyWeight = c(,1,))
    fixef( dVdGender1(enable=c(1)) = c(,0,))
    fixef( dCldAge(enable=c(2)) = c(,0,))
    fixef( dCldBodyWeight(enable=c(3)) = c(,0,))
    fixef( dCldGender1(enable=c(4)) = c(,0,))
    ranef(diag(nV,nCl,nV2,nCl2) = c(0.0691495878371252,0.180663528781695,4.64989751058228e-08,0.0423919347837654))
}"
  
writeLines(pml, "test.mdl")

Now that we have our PML code as a text file in our working directory, we can create a textual model in Certara.RsNLME using the textualmodel() function:

Tip: To open your PML file in RStudio for editing before creating a model with textualmodel(), run file.edit("test.mdl"). Alternatively, edit in your favorite text editor (e.g., Notepad++). You can make additional changes to PML after creating the model with textualmodel() using editModel().

scmModelTextual <- textualmodel(modelName = "scmModelTextual", data = pkData, mdl = "test.mdl")

print(scmModelTextual)

#> 
#>  Model Overview 
#>  ------------------------------------------- 
#> Model Name        :  scmModelTextual
#> Working Directory :  C:/Repos/R-RsNLME/vignettes/scmModelTextual
#> Model Type        :  Textual
#> 
#>  PML 
#>  ------------------------------------------- 
#> test(){
#>     cfMicro(A1,Cl/V, Cl2/V, Cl2/V2)
#>     dosepoint(A1)
#>     C = A1 / V
#>     error(CEps=0.161205984597858)
#>     observe(CObs=C * ( 1 + CEps))
#>     stparm(V = tvV * (Age^dVdAge)   * (BodyWeight^dVdBodyWeight)   * exp(dVdGender1*(Gender==1)) * exp(nV))
#>     stparm(Cl = tvCl * (Age^dCldAge)   * (BodyWeight^dCldBodyWeight)   * exp(dCldGender1*(Gender==1)) * exp(nCl))
#>     stparm(V2 = tvV2 * exp(nV2))
#>     stparm(Cl2 = tvCl2 * exp(nCl2))
#>     fcovariate(Age)
#>     fcovariate(BodyWeight)
#>     fcovariate(Gender())
#>     fixef( tvV = c(,15.393730662165,))
#>     fixef( tvCl = c(,6.59983805751661,))
#>     fixef( tvV2 = c(,41.1900094327716,))
#>     fixef( tvCl2 = c(,14.0257699237193,))
#>     fixef( dVdAge(enable=c(0)) = c(,0,))
#>     #fixef( dVdBodyWeight(enable=c(1)) = c(,0,))
#>     fixef( dVdBodyWeight = c(,1,))
#>     fixef( dVdGender1(enable=c(1)) = c(,0,))
#>     fixef( dCldAge(enable=c(2)) = c(,0,))
#>     fixef( dCldBodyWeight(enable=c(3)) = c(,0,))
#>     fixef( dCldGender1(enable=c(4)) = c(,0,))
#>     ranef(diag(nV,nCl,nV2,nCl2) = c(0.0691495878371252,0.180663528781695,4.64989751058228e-08,0.0423919347837654))
#> }
#> 
#>  Structural Parameters 
#>  ------------------------------------------- 
#>  V Cl V2 Cl2
#>  ------------------------------------------- 
#>  Column Mappings 
#>  ------------------------------------------- 
#> Model Variable Name : Data Column name
#> id                  : ?
#> time                : ?
#> A1                  : ?
#> Age                 : Age
#> BodyWeight          : BodyWeight
#> Gender              : Gender
#> CObs                : ?

From the Model Overview in the print output, we can see that Model Type: Textual.

From the Column Mappings information, since our covariate model variables are named the same as the columns in our pkData they get automatically mapped. For other model variables with ? we will need to map with colMapping().

 ------------------------------------------- 
 Column Mappings 
 ------------------------------------------- 
Model Variable Name : Data Column name
id                  : ?
time                : ?
A1                  : ?
Age                 : Age
BodyWeight          : BodyWeight
Gender              : Gender
CObs                : ?

scmModelTextual <- scmModelTextual %>%
    colMapping(c(id = "Subject", time = "Act_Time", CObs = "Conc", A1 = "Amount"))

For this scmModelTextual, since we created it directly from PML, and are using a categorical covariate in the model that has a character column in the data (Gender), we’ll need to additionally provide the correct levels/labels. The labels are not required if the categorical covariate is numeric in the data

Note: For the built-in model, we added covariates using addCovariate(), where we specified levels/labels for the Gender categorical covariate.

scmModelTextual <- scmModelTextual %>% 
  addLabel("Gender",
           levels=c(0,1),
           labels=c("male","female"))

Now, when running print(scmModelTextual), you will observe the column mapping for Gender has been updated with specified levels and labels.

 ------------------------------------------- 
 Column Mappings 
 ------------------------------------------- 
Model Variable Name : Data Column name
id                  : Subject
time                : Act_Time
A1                  : Amount
Age                 : Age
BodyWeight          : BodyWeight
Gender              : Gender( male=0 female=1 )
CObs                : Conc

Run stepwise search for textual model and examine best model

Just as we did previously for the built-in model, we’ll use the same stepwiseSearch() function on the textual model to run SCM, except this time the covariate effect of BodyWeight on V will be present across all scenarios, since we’ve removed the corresponding enable statement for the fixef().

scmResultTextual <- stepwiseSearch(scmModelTextual) 
#> 
#> stepwiseParams argument is not given. 
#> Using Default: 
#> Thresholds for adding/removing a covariate effect: 0.01 / 0.001
#> Method: -2LL
#> 
#> covariateModel argument is not given. Processing model object to make it.
#> Stepwise Job
#> 
#> NLME Job
#> sharedDirectory is not given in the host class instance.
#> Valid NLME_ROOT_DIRECTORY is not given, using current working directory:
#> C:/Repos/R-RsNLME/vignettes
#> TDL5 version: 25.7.1.1
#> 
#> Status: OK
#> License expires: 2026-11-14
#> Refresh until: 2025-12-14 09:53:10
#> Current Date: 2025-11-14
#> The model compiled
#> 
#> Trying to generate job results...
#> Done generating job results.
#> TDL5 version: 25.7.1.1
#> 
#> Status: OK
#> License expires: 2026-11-14
#> Refresh until: 2025-12-14 09:53:10
#> Current Date: 2025-11-14
#> The model compiled
#> 
#> Trying to generate job results...
#> Done generating job results.

print(scmResultTextual)

SCM results

scmResultTextual[scmResultTextual$BestScenario == TRUE,]
#>                  Scenario RetCode    LogLik     -2LL      AIC      BIC nParm
#> 5 cstep004  Cl-BodyWeight       2 -609.2603 1218.521 1240.521 1270.424    11
#>   nObs nSub EpsShrinkage Condition BestScenario
#> 5  112   16      0.11327        NA         TRUE

Shotgun search

Use shotgun search to assess the inflation of additional covariates included in the model. See here for additional examples.

shotgunResult <- shotgunSearch(scmModelTextual) 
#> 
#> covariateModel argument is not given. Processing model object to make it.
#> Shotgun Job
#> 
#> NLME Job
#> sharedDirectory is not given in the host class instance.
#> Valid NLME_ROOT_DIRECTORY is not given, using current working directory:
#> C:/Repos/R-RsNLME/vignettes
#> TDL5 version: 25.7.1.1
#> 
#> Status: OK
#> License expires: 2026-11-14
#> Refresh until: 2025-12-14 09:53:10
#> Current Date: 2025-11-14
#> The model compiled
#> 
#> Trying to generate job results...
#> Done generating job results.

shotgunResult %>%
  arrange(AIC) %>%
  head(10)
#>                                         Scenario RetCode    LogLik     -2LL
#> 1                  cshot012 Cl-Age Cl-BodyWeight       2 -606.7152 1213.430
#> 2        cshot028 Cl-Age Cl-BodyWeight Cl-Gender       2 -605.9782 1211.956
#> 3         cshot014 V-Gender Cl-Age Cl-BodyWeight       3 -605.9924 1211.985
#> 4         cshot025 V-Age Cl-BodyWeight Cl-Gender       2 -606.0559 1212.112
#> 5               cshot024 Cl-BodyWeight Cl-Gender       2 -607.5986 1215.197
#> 6  cshot029 V-Age Cl-Age Cl-BodyWeight Cl-Gender       2 -605.9738 1211.948
#> 7                         cshot008 Cl-BodyWeight       2 -609.2603 1218.521
#> 8          cshot011 V-Age V-Gender Cl-BodyWeight       2 -607.4807 1214.961
#> 9                   cshot009 V-Age Cl-BodyWeight       2 -608.6864 1217.373
#> 10               cshot010 V-Gender Cl-BodyWeight       2 -608.7207 1217.441
#>         AIC      BIC nParm nObs nSub EpsShrinkage    Condition
#> 1  1237.430 1270.052    12  112   16      0.12722 8.778448e+12
#> 2  1237.956 1273.297    13  112   16      0.12759 7.459216e+08
#> 3  1237.985 1273.325    13  112   16      0.13125 3.738741e+05
#> 4  1238.112 1273.452    13  112   16      0.11483           NA
#> 5  1239.197 1271.819    12  112   16      0.13135 1.387281e+17
#> 6  1239.948 1278.007    14  112   16      0.12974 3.165627e+05
#> 7  1240.521 1270.424    11  112   16      0.11327           NA
#> 8  1240.961 1276.302    13  112   16      0.11561           NA
#> 9  1241.373 1273.995    12  112   16      0.12856 3.986554e+04
#> 10 1241.441 1274.063    12  112   16      0.12459 6.324903e+04

Traceability with textual models

When using the built-in modeling functions, changes are naturally captured in your R script, which makes provenance and reproducibility straightforward. After converting to a textual model and editing PML directly, it’s important to explicitly track and compare model definitions so that:

• You can verify exactly which covariate effects or structural edits were introduced between model versions via the R code
• Reviewers can audit changes without ambiguity
• You can reconstruct the analysis later

For textual models, there may be different ways that you’ve edited the PML e.g., editModel(), modelTextualUI() or manually editing the PML file in a text editor. Given this, it’s good practice to preserve traceability by writing the PML to disk and diffing the files.

As we’ve demonstrate above, we can use:

writeLines(paste(scmModel@statements), "scmModel.mdl")
writeLines(paste(scmModelTextual@statements), "scmModelTextual.mdl")

to write the PML for the the two models to our working directory. With the diffobj package, we can then ‘diff’ (show the differences between) the two files:

library(diffobj)
diffChr(
    readLines("scmModel.mdl"), 
    readLines("scmModelTextual.mdl"), 
    mode = "sidebyside", 
    format = "ansi8"
    )
#> < readLines("scmModel.mdl")              > readLines("scmModelTextual.mdl")     
#> @@ 17,10 @@                              @@ 17,10 @@                            
#>       fixef( tvCl2 = c(,14.025769923719        fixef( tvCl2 = c(,14.025769923719
#>   3,))                                     3,))                                 
#>       fixef( dVdAge(enable=c(0)) = c(,0        fixef( dVdAge(enable=c(0)) = c(,0
#>   ,))                                      ,))                                  
#> <     fixef( dVdBodyWeight(enable=c(1))  >     #fixef( dVdBodyWeight(enable=c(1)
#> :  = c(,0,))                             : ) = c(,0,))                          
#> <     fixef( dVdGender1(enable=c(2)) =   >     fixef( dVdBodyWeight = c(,1,))   
#> : c(,0,))                                ~                                      
#> ~                                        >     fixef( dVdGender1(enable=c(1)) = 
#> ~                                        : c(,0,))                              
#> <     fixef( dCldAge(enable=c(3)) = c(,  >     fixef( dCldAge(enable=c(2)) = c(,
#> : 0,))                                   : 0,))                                 
#> <     fixef( dCldBodyWeight(enable=c(4)  >     fixef( dCldBodyWeight(enable=c(3)
#> : ) = c(,0,))                            : ) = c(,0,))                          
#> <     fixef( dCldGender1(enable=c(5)) =  >     fixef( dCldGender1(enable=c(4)) =
#> :  c(,0,))                               :  c(,0,))                             
#>       ranef(diag(nV,nCl,nV2,nCl2) = c(0        ranef(diag(nV,nCl,nV2,nCl2) = c(0
#>   .0691495878371252,0.180663528781695,4    .0691495878371252,0.180663528781695,4
#>   .64989751058228e-08,0.042391934783765    .64989751058228e-08,0.042391934783765
#>   4))                                      4))                                  
#> <                                        ~                                      
#>   }                                        }

#> [1] TRUE TRUE

You can also simply diff the PML statements directly from the model objects e.g., diffChr(scmModel@statements, scmModelTextual@statements, mode = "sidebyside", format = "ansi8").