Covariate Selection Methods - Forward Selection

Authors

Jose Storopoli

Andreas Noack

Joel Owen

As stated in Introduction to Covariate Selection Methods, broadly, there are two main types of Stepwise Covariate Model (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\):
1. define \(M_c = M + c\)
2. define \(f_c = f(M_c, D)\)
3. 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^*\)

Note

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,
    criterion = aic,
    method = CovariateSelection.Forward,
    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:

using Pumas
using PharmaDatasets

We are going to use the po_sad_1 dataset from PharmaDatasets:

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

5×14 DataFrame

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 poor metabolizer (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 = @model begin
    @metadata begin
        desc = "full covariate model"
        timeu = u"hr"
    end

    @param begin
        """
        Clearance (L/hr)
        """
        tvcl ∈ RealDomain(; lower = 0)
        """
        Central Volume (L)
        """
        tvvc ∈ RealDomain(; lower = 0)
        """
        Peripheral Volume (L)
        """
        tvvp ∈ RealDomain(; lower = 0)
        """
        Distributional Clearance (L/hr)
        """
        tvq ∈ RealDomain(; lower = 0)
        """
        Fed - Absorption rate constant (h-1)
        """
        tvkafed ∈ RealDomain(; lower = 0)
        """
        Fasted - Absorption rate constant (h-1)
        """
        tvkafasted ∈ RealDomain(; lower = 0)
        """
        Power exponent on weight for Clearance
        """
        dwtcl ∈ RealDomain()
        """
        Proportional change in CL due to PM
        """
        ispmoncl ∈ RealDomain(; lower = -1, upper = 1)
        """
          - Ω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)
        """
        wt
    end

    @pre begin
        CL = tvcl * (wt / 70)^dwtcl * (1 + ispmoncl * (isPM == "yes")) * exp(η[1])
        Vc = tvvc * (wt / 70) * exp(η[2])
        Ka = (tvkafed + (isfed == "no") * (tvkafasted)) * exp(η[3])
        Q = tvq * (wt / 70)^0.75 * exp(η[4])
        Vp = tvvp * (wt / 70) * exp(η[5])
    end

    @dynamics Depots1Central1Periph1

    @derived begin
        cp := @. 1000 * (Central / Vc)
        """
        DrugY Concentration (ng/mL)
        """
        dv ~ @. Normal(cp, cp * σₚ)
    end
end

PumasModel
  Parameters: tvcl, tvvc, tvvp, tvq, tvkafed, tvkafasted, dwtcl, ispmoncl, Ω, σₚ
  Random effects: η
  Covariates: isPM, isfed, wt
  Dynamical variables: Depot, Central, Peripheral
  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:

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

shows us that if tvkafasted is zero, then the whole expression (isfed == "no")*(tvkafasted) evaluates to zero, and we are only left with:

Ka = tvkafed * exp(η[3])

which removes the effect of tvkafasted from the value of Ka.

For the case of dwtcl and ispmoncl, we have the expression:

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

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 = (;
    tvkafed = 0.4,
    tvkafasted = 1.5,
    tvcl = 4.0,
    tvvc = 70.0,
    tvq = 4.0,
    tvvp = 50.0,
    dwtcl = 0.75,
    ispmoncl = -0.7,
    Ω = 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):

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

[ 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.2529223967354
[ 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.7204514003934
[ Info: current best model has zero restrictions on dwtcl
[ Info: criterion improved, continuing!
[ Info: fitting model with zero restrictions on: 
[ Info: criterion: 2641.58594655915
[ Info: current best model has zero restrictions on 
[ Info: criterion improved, continuing!
┌ Info: final best model summary
│   best_criterion = 2641.58594655915
└   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: Objective Function Value.

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:

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

[ 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.0261144662118
[ Info: fitting model with zero restrictions on: dwtcl, ispmoncl
[ Info: criterion: 2159.845796865268
[ 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.4936434698707
[ Info: current best model has zero restrictions on dwtcl
[ Info: criterion improved, continuing!
[ Info: fitting model with zero restrictions on: 
[ Info: criterion: 2117.359138628627
[ Info: current best model has zero restrictions on 
[ Info: criterion improved, continuing!
┌ Info: final best model summary
│   best_criterion = 2117.359138628627
└   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)

7×4 DataFrame

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
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.

Reuse

CC BY-SA 4.0