```
using Pumas
using PharmaDatasets
```

# Covariate Selection Methods - Forward Selection

As stated in Introduction to Covariate Selection Methods, broadly, there are two main types of **S**tepwise **C**ovariate **M**odel (SCM):

**Forward Selection**(FS)**Backward Elimination**(BE)

In this tutorial we’ll explore **Forward Selection** (FS) in detail.

## 1 Forward Selection (FS)

Forward selection stepwise covariate modeling procedures begin with a base model, \(M_0\), and for every possible covariate effect \(c \in C\), where \(C\) is the set of all covariate effects, it adds all available \(c\) individually to selected parameters of \(M_0\) and checks for model fit improvements using the selected metric. Once the best \(c\)-covariate parameter pair is found, SCM adds it to the model, defining a new reference model \(M_k\), where \(k\) is the iteration number, and removes \(c_k\) (the best \(c\) in iteration \(k\)) from the set \(C_{k+1}\) of covariate effects for iteration \(k+1\). SCM proceeds with additional iterations until no additional significant \(c_k\) is found, or there aren’t \(c_k\) left to be included in the model.

Here’s a breakdown of the FS algorithm:

- Base model without covariate effects: \(M_0\)
- Set of covariate effects: \(C\)
- Data: \(D\)
- Criterion, a function that takes a model along with data, and outputs a real number (with the assumption that a lower value is better): \(f: (M, D) \to \mathbb{R}\)

- Best model with a subset of covariate effects \(C^*\) according to criterion \(f\) evaluated on data \(D\)

`define`

\(M = M_0\)`define`

\(f_m = f(M, D)\)`for`

every \(c \in C\):`define`

\(M_c = M + c\)`define`

\(f_c = f(M_c, D)\)`if`

\(f_c < f_m\)`define`

\(M_{\text{best}} = M_c\)`define`

\(c_{\text{best}} = c\)`else`

`unset`

\(c_{\text{best}}\)

`if`

\(c_{\text{best}}\) is`defined`

`define`

\(M = M_{\text{best}}\)`define`

\(C = C \setminus \{ c_{\text{best}} \}\)`goto`

step 2

define \(C^* = C\)

`return`

\(M\) and \(C^*\)

In the procedures above, we start by defining the initial model \(M\) as the base model without covariate effects \(M_0\). Then we proceed by defining the initial criterion \(f_m\) as the criterion \(f\) of the initial model \(M\) and data \(D\). Next, we do a for loop on every covariate effect \(c\) in the set of all covariate effects \(C\), where we:

define the current model \(M_c\) as the model \(M\) with the covariate effect \(c\), that is \(M + c\).

define the criterion of the current model \(f_c\) as the criterion \(f\) applied to the current model \(M_c\) and data \(D\).

if the criterion improved, i.e \(f_c < f_m\), then:

- define the new best model \(M_{\text{best}}\) as the current model \(M_c\).
- define the new best covariate effect \(c_{\text{best}}\) as the covariate effect \(c\).

if we have a best covariate effect \(c_{\text{best}}\), then:

- define the new “initial model” \(M\) as the best model \(M_{\text{best}}\).
- define the new set of covariate effects as the set \(C\) minus best covariate effect \(c_{\text{best}}\)
- go to the step 2 and resume from there.

define the new set of covariate effects \(C^*\) as the set complement \(C^\complement\) of the set \(C\)

return the final “initial model” \(M\) and the new set of covariate effects \(C^*\)

Some models are very sensitive to initial parameter values, resulting in failed fitting procedures and values of \(-\infty\) in the log-likelihood estimation. In such cases, Forward Selection (FS) might be a better alternative than Backward Elimination (BS), since FS starts with fewer model parameters than BE and finding good initial estimates may be comparatively easier.

## 2 The `covariate_select`

function

All stepwise covariate selection technique in Pumas uses the `covariate_select`

function. It has the following function signature:

```
covariate_select(
model,
population,
param,
approx;
control_param,= aic,
criterion = CovariateSelection.Forward,
method
fit_options, )
```

There are positional arguments and keyword arguments.

These should sound familiar to you:

`model`

: a Pumas model`population`

: a`Population`

that is commonly obtained from`read_pumas`

`param`

: a set of initial parameter estimates`approx`

: a likelihood approximation method, e.g.`FOCE()`

,`LaplaceI()`

, etc.

These are exactly the same positional arguments that we pass to the `fit`

function to perform a model fit procedure.

The keyword arguments we have some new concepts:

`control_param`

: this is a tuple of parameters that encapsulate unique covariate effects such that, if we restrain their values to \(0\), we effectively eliminate any influence of the covariate effects in our model.`control_param`

can be thought as a list of covariate effects that we want to screen in our covariate selection routine.`criterion`

: Criterion for covariate selection. This needs to be a function that accepts a Pumas model and returns a real value. If you don’t provide a function Pumas will, by default, use the`aic`

function. Note that`criterion`

has the assumption that**lower is better**. Hence, you’ll need to take that into account when supplying a different function or defining your own user-defined function.`method`

: the desired method to be used for the stepwise covariate selection procedure. The default is`CovariateSelection.Forward`

which stands for Forward Selection (FS). You can also opt for`CovariateSelection.Backward`

if you want to use Backward Elimination (BS)`fit_options`

: a`NamedTuple`

of options that will be used during the model fit procedures that`covariate_select`

will perform. These will be the keyword arguments available for the`fit`

function, e.g.`(; constantcoef = (; myparam = 1.0))`

.

## 3 Example in Pumas

Let’s showcase an example that builds from the model we used in Covariate Selection Methods - Introduction.

First, let’s import the following packages:

We are going to use the `po_sad_1`

dataset from `PharmaDatasets`

:

```
= dataset("po_sad_1")
df first(df, 5)
```

Row | id | time | dv | amt | evid | cmt | rate | age | wt | doselevel | isPM | isfed | sex | route |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Int64 | Float64 | Float64? | Float64? | Int64 | Int64? | Float64 | Int64 | Int64 | Int64 | String3 | String3 | String7 | String3 | |

1 | 1 | 0.0 | missing | 30.0 | 1 | 1 | 0.0 | 51 | 74 | 30 | no | yes | male | ev |

2 | 1 | 0.25 | 35.7636 | missing | 0 | missing | 0.0 | 51 | 74 | 30 | no | yes | male | ev |

3 | 1 | 0.5 | 71.9551 | missing | 0 | missing | 0.0 | 51 | 74 | 30 | no | yes | male | ev |

4 | 1 | 0.75 | 97.3356 | missing | 0 | missing | 0.0 | 51 | 74 | 30 | no | yes | male | ev |

5 | 1 | 1.0 | 128.919 | missing | 0 | missing | 0.0 | 51 | 74 | 30 | no | yes | male | ev |

This is an oral dosing (`route = "ev"`

) NMTRAN-formatted dataset. It has 18 subjects, each with 1 dosing event (`evid = 1`

) and 18 measurement events (`evid = 0`

); and the following covariates:

`age`

: age in years (continuous)`wt`

: weight in kg (continuous)`doselevel`

: dosing amount, either`30`

,`60`

or`90`

milligrams (categorical)`isPM`

: subject is a**p**oor**m**etabolizer (binary)`isfed`

: subject is**fed**(binary)`sex`

: subject sex (binary)

Let’s parse `df`

into a `Population`

with `read_pumas`

:

```
=
population read_pumas(df; observations = [:dv], covariates = [:wt, :isPM, :isfed], route = :route)
```

```
Population
Subjects: 18
Covariates: wt, isPM, isfed
Observations: dv
```

We’ll use the same 2-compartment oral absorption model as in Covariate Selection Methods - Introduction, but now adding additional covariates and covariate effects:

```
= @model begin
model @metadata begin
= "full covariate model"
desc = u"hr"
timeu end
@param begin
"""
Clearance (L/hr)
"""
∈ RealDomain(; lower = 0)
tvcl """
Central Volume (L)
"""
∈ RealDomain(; lower = 0)
tvvc """
Peripheral Volume (L)
"""
∈ RealDomain(; lower = 0)
tvvp """
Distributional Clearance (L/hr)
"""
∈ RealDomain(; lower = 0)
tvq """
Fed - Absorption rate constant (h-1)
"""
∈ RealDomain(; lower = 0)
tvkafed """
Fasted - Absorption rate constant (h-1)
"""
∈ RealDomain(; lower = 0)
tvkafasted """
Power exponent on weight for Clearance
"""
∈ RealDomain()
dwtcl """
Proportional change in CL due to PM
"""
∈ RealDomain(; lower = -1, upper = 1)
ispmoncl """
- ΩCL
- ΩVc
- ΩKa
- ΩVp
- ΩQ
"""
∈ PDiagDomain(5)
Ω """
Proportional RUV (SD scale)
"""
∈ RealDomain(; lower = 0)
σₚ end
@random begin
~ MvNormal(Ω)
η end
@covariates begin
"""
Poor Metabolizer
"""
isPM"""
Fed status
"""
isfed"""
Weight (kg)
"""
wtend
@pre begin
= tvcl * (wt / 70)^dwtcl * (1 + ispmoncl * (isPM == "yes")) * exp(η[1])
CL = tvvc * (wt / 70) * exp(η[2])
Vc = (tvkafed + (isfed == "no") * (tvkafasted)) * exp(η[3])
Ka = tvq * (wt / 70)^0.75 * exp(η[4])
Q = tvvp * (wt / 70) * exp(η[5])
Vp end
@dynamics Depots1Central1Periph1
@derived begin
:= @. 1_000 * (Central / Vc)
cp """
DrugY Concentration (ng/mL)
"""
~ @. Normal(cp, cp * σₚ)
dv end
end
```

```
PumasModel
Parameters: tvcl, tvvc, tvvp, tvq, tvkafed, tvkafasted, dwtcl, ispmoncl, Ω, σₚ
Random effects: η
Covariates: isPM, isfed, wt
Dynamical system variables: Depot, Central, Peripheral
Dynamical system type: Closed form
Derived: dv
Observed: dv
```

The underlying Pumas model that `covariate_select`

will use needs to include **all** of the covariate effects to be tested. This is different than just specifying particular covariates in individual models for a more manual process. Note that one covariate could have more than one effect.

During each iteration of FS, `covariate_select`

will set all of the covariate effects to zero except the one that’s being evaluated for inclusion in the model. Each of the covariate effects should be totally controlled by a respective parameter defined in the `@param`

block. The model should be coded in such way that, if the corresponding parameter value is set to zero, then the covariate effect is completely removed from the model.

In the model above, we are defining the following parameters in the `@param`

block with respect to covariate effects:

`tvkafasted`

: typical value of the magnitude of the`isfed`

covariate effect on`tvka`

`dwtcl`

: typical value of the magnitude of the`wt`

covariate efffect on`tvcl`

`ispmoncl`

: typical value of the magnitude of the`isPM`

covariate efffect on`tvcl`

These will be the parameters that we will be testing for the covariate selection.

Further down in the `@pre`

block we see that if we set any of these parameters value to zero we will completely remove the covariate effect from the model. For the case of `tvkafasted`

, this expression in the `@pre`

block:

`= (tvkafed + (isfed == "no") * (tvkafasted)) * exp(η[3]) Ka `

shows us that if `tvkafasted`

is zero, then the whole expression `(isfed == "no")*(tvkafasted)`

evaluates to zero, and we are only left with:

`= tvkafed * exp(η[3]) Ka `

which removes the effect of `tvkafasted`

from the value of Ka.

For the case of `dwtcl`

and `ispmoncl`

, we have the expression:

`= tvcl * (wt / 70)^dwtcl * (1 + ispmoncl * (isPM == "yes")) * exp(η[1]) CL `

in the `@pre`

block.

- if
`dwtcl`

is set to zero, then`(wt / 70)^dwtcl`

evaluates to one. This would remove the`wt`

effect on`tvcl`

. - if
`ispmoncl`

is set to zero, then`ispmoncl * (isPM == "yes")`

evaluates to zero. This would remove the`isPM`

effect on`tvcl`

.

Note that it is the effect size parameter (i.e., `dwtcl`

or `ispmoncl`

), that is set to zero and effectively removes the covariate effect, whether continuous or categorical.

With the model done, now we define our initial parameters:

```
= (;
iparams = 0.4,
tvkafed = 1.5,
tvkafasted = 4.0,
tvcl = 70.0,
tvvc = 4.0,
tvq = 50.0,
tvvp = 0.75,
dwtcl = -0.7,
ispmoncl = Diagonal(fill(0.04, 5)),
Ω = 0.1,
σₚ )
```

```
(tvkafed = 0.4,
tvkafasted = 1.5,
tvcl = 4.0,
tvvc = 70.0,
tvq = 4.0,
tvvp = 50.0,
dwtcl = 0.75,
ispmoncl = -0.7,
Ω = [0.04 0.0 … 0.0 0.0; 0.0 0.04 … 0.0 0.0; … ; 0.0 0.0 … 0.04 0.0; 0.0 0.0 … 0.0 0.04],
σₚ = 0.1,)
```

### 3.1 FS with Different Metrics

Finally, we can run the covariate selection procedure.

We’ll use `FOCE()`

as our estimation method. For the `control_param`

tuple, we include all of the covariate effects we mentioned above: `dwtcl`

, `tvkafasted`

and `ispmoncl`

. These will be the parameters that we’ll test for covariate selection.

All others keyword arguments will be left to their default values, which implies `aic`

as the criterion and the FS method (`CovariateSelection.Forward`

):

```
= covariate_select(
covar_result
model,
population,
iparams,FOCE();
= (:dwtcl, :tvkafasted, :ispmoncl), # this is a Tuple
control_param )
```

```
[ Info: fitting model with zero restrictions on: dwtcl, tvkafasted, ispmoncl
[ Info: criterion: 2714.959963420873
[ Info: fitting model with zero restrictions on: tvkafasted, ispmoncl
[ Info: criterion: 2713.252922396735
[ Info: fitting model with zero restrictions on: dwtcl, ispmoncl
[ Info: criterion: 2680.0726047957905
[ Info: fitting model with zero restrictions on: dwtcl, tvkafasted
[ Info: criterion: 2702.591259599917
[ Info: current best model has zero restrictions on dwtcl, ispmoncl
[ Info: criterion improved, continuing!
[ Info: fitting model with zero restrictions on: ispmoncl
[ Info: criterion: 2678.298097810124
[ Info: fitting model with zero restrictions on: dwtcl
[ Info: criterion: 2667.720451400392
[ Info: current best model has zero restrictions on dwtcl
[ Info: criterion improved, continuing!
[ Info: fitting model with zero restrictions on:
[ Info: criterion: 2641.585946318844
[ Info: current best model has zero restrictions on
[ Info: criterion improved, continuing!
┌ Info: final best model summary
│ best_criterion = 2641.585946318844
└ best_zero_restrictions = ()
```

```
Pumas.CovariateSelectionResults
Criterion: aic
Method: Forward
Best model:
Unrestricted parameters: dwtcl, tvkafasted, ispmoncl
Zero restricted parameters:
Criterion value: 2641.59
Fitted model summary:
7×4 DataFrame
Row │ free_parameters zero_parameters fit ⋯
│ Tuple… Tuple… Fit ⋯
─────┼──────────────────────────────────────────────────────────────────────────
1 │ () (:dwtcl, :tvkafasted, :ispmoncl) Fit ⋯
2 │ (:dwtcl,) (:tvkafasted, :ispmoncl) Fit
3 │ (:tvkafasted,) (:dwtcl, :ispmoncl) Fit
4 │ (:ispmoncl,) (:dwtcl, :tvkafasted) Fit
5 │ (:dwtcl, :tvkafasted) (:ispmoncl,) Fit ⋯
6 │ (:tvkafasted, :ispmoncl) (:dwtcl,) Fit
7 │ (:dwtcl, :tvkafasted, :ispmoncl) () Fit
2 columns omitted
```

As you can see, `covariate_select`

, upon completion, prints its result.

Now, let’s showcase how you can use a custom `criterion`

.

Here’s an example of a function that takes a fitted Pumas model, `m`

, and outputs the -2LL (minus 2 times log-likelihood) without the constant. This is a model fit measure commonly used by NONMEM users and is is known as OFV: **O**bjective **F**unction **V**alue.

Hence, we will name the function `ofv`

:

`ofv(m) = (-2 * loglikelihood(m)) - (nobs(m) * log(2π))`

`ofv (generic function with 1 method)`

Now let’s run `covariate_select`

again but with `ofv`

as a criterion:

```
= covariate_select(
covar_result_ofv
model,
population,
iparams,FOCE();
= (:dwtcl, :tvkafasted, :ispmoncl),
control_param = ofv,
criterion = CovariateSelection.Forward,
method )
```

```
[ Info: fitting model with zero restrictions on: dwtcl, tvkafasted, ispmoncl
[ Info: criterion: 2196.73315549035
[ Info: fitting model with zero restrictions on: tvkafasted, ispmoncl
[ Info: criterion: 2193.026114466212
[ Info: fitting model with zero restrictions on: dwtcl, ispmoncl
[ Info: criterion: 2159.8457968652674
[ Info: fitting model with zero restrictions on: dwtcl, tvkafasted
[ Info: criterion: 2182.364451669394
[ Info: current best model has zero restrictions on dwtcl, ispmoncl
[ Info: criterion improved, continuing!
[ Info: fitting model with zero restrictions on: ispmoncl
[ Info: criterion: 2156.071289879601
[ Info: fitting model with zero restrictions on: dwtcl
[ Info: criterion: 2145.4936434698693
[ Info: current best model has zero restrictions on dwtcl
[ Info: criterion improved, continuing!
[ Info: fitting model with zero restrictions on:
[ Info: criterion: 2117.3591383883204
[ Info: current best model has zero restrictions on
[ Info: criterion improved, continuing!
┌ Info: final best model summary
│ best_criterion = 2117.3591383883204
└ best_zero_restrictions = ()
```

```
Pumas.CovariateSelectionResults
Criterion: ofv
Method: Forward
Best model:
Unrestricted parameters: dwtcl, tvkafasted, ispmoncl
Zero restricted parameters:
Criterion value: 2117.36
Fitted model summary:
7×4 DataFrame
Row │ free_parameters zero_parameters fit ⋯
│ Tuple… Tuple… Fit ⋯
─────┼──────────────────────────────────────────────────────────────────────────
1 │ () (:dwtcl, :tvkafasted, :ispmoncl) Fit ⋯
2 │ (:dwtcl,) (:tvkafasted, :ispmoncl) Fit
3 │ (:tvkafasted,) (:dwtcl, :ispmoncl) Fit
4 │ (:ispmoncl,) (:dwtcl, :tvkafasted) Fit
5 │ (:dwtcl, :tvkafasted) (:ispmoncl,) Fit ⋯
6 │ (:tvkafasted, :ispmoncl) (:dwtcl,) Fit
7 │ (:dwtcl, :tvkafasted, :ispmoncl) () Fit
2 columns omitted
```

First, we have the best model that was found with its parameters and criterion value (showcasing the `aic`

case).

Next, we have a `DataFrame`

-rendered table with the iterations that FS tested as the rows, and the following columns:

`free_parameters`

: covariate effects that were allowed to influence the model, i.e. were*not*forced to zero.`zero_parameters`

: covariate effects that were*not*allowed to influence the model, i.e. were forced to zero.`fitted_model`

: the underlying Pumas’ model fit`criterion`

: criterion value

We can access the results of the covariate selection procedure by accessing the `fits`

field of the `covar_result`

object, and convert it to a `DataFrame`

in order to perform any data wrangling or plotting.

Here’s an example of how to sort the results by the `criterion`

value:

`sort(DataFrame(covar_result.fits), :criterion)`

Row | free_parameters | zero_parameters | fitted_model | criterion |
---|---|---|---|---|

Tuple… | Tuple… | FittedPu… | Float64 | |

1 | (:dwtcl, :tvkafasted, :ispmoncl) | () | FittedPumasModel(...) | 2641.59 |

2 | (:tvkafasted, :ispmoncl) | (:dwtcl,) | FittedPumasModel(...) | 2667.72 |

3 | (:dwtcl, :tvkafasted) | (:ispmoncl,) | FittedPumasModel(...) | 2678.3 |

4 | (:tvkafasted,) | (:dwtcl, :ispmoncl) | FittedPumasModel(...) | 2680.07 |

5 | (:ispmoncl,) | (:dwtcl, :tvkafasted) | FittedPumasModel(...) | 2702.59 |

6 | (:dwtcl,) | (:tvkafasted, :ispmoncl) | FittedPumasModel(...) | 2713.25 |

7 | () | (:dwtcl, :tvkafasted, :ispmoncl) | FittedPumasModel(...) | 2714.96 |

As you can see this puts the “best” model (according to the `method`

and `criterion`

) in the first row.

After running the covariate selection procedure, we get a summary of the best model found and the fitted model for each step of the procedure.

You can recover the best model fit by accessing the `best_model`

field of the `covar_result`

object:

` covar_result.best_model`

```
FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -1306.793
Number of subjects: 18
Number of parameters: Fixed Optimized
0 14
Observation records: Active Missing
dv: 270 0
Total: 270 0
-------------------------
Estimate
-------------------------
tvcl 3.6888
tvvc 70.044
tvvp 48.706
tvq 4.3506
tvkafed 0.41114
tvkafasted 1.0894
dwtcl 0.81856
ispmoncl -0.62382
Ω₁,₁ 0.022271
Ω₂,₂ 0.050192
Ω₃,₃ 0.031194
Ω₄,₄ 0.10595
Ω₅,₅ 0.035072
σₚ 0.098043
-------------------------
```

This works as expected and returns the same best model as the AIC approach above.