using Pumas
using PharmaDatasets
Covariate Selection Methods - Introduction
In pharmacometric workflows, we often have competing models to select from. In this tutorial we will review selection criteria and automated procedures to select the best model out of a set of competing candidate models.
First, we’ll review how to measure model fit, then we’ll cover model selection algorithms.
1 Model Fit Measures
Traditionally in Statistics, model comparison has been done based on a theoretical divergence metric that originates from information theory’s entropy:
\[H = - \operatorname{E}\log(p) = -\sum_i p_i \log(p_i)\]
where \(p_i\) is the probability of occurrence of the \(i\)-th possible value.
Note
We use the \(\log\) scale because it transforms a product of probabilities into a sum, which is both numerically faster and numerically more stable due to the robustness against floating point underflow.
Entropy was the inspiration behind Akaike’s Information Criterion (AIC) (Akaike, 1973):
\[\operatorname{AIC} = -2\log{\hat{\mathcal{L}}} + 2k\]
where \(\hat{\mathcal{L}}\) is the estimated value of the likelihood for a given model and data, and \(k\) is the number of parameters in the model. Generally the likelihood is estimated by maximizing the likelihood function, thus the name maximum likelihood estimation (MLE). The likelihood describes how well the model fits the data, and in certain conditions, can be treated similarly to a probability: higher values means higher plausibility. Hence, models with higher likelihood values demonstrate better fits to the data. Since we are multiplying by a negative value, this means that lower values are preferred.
Note
The \(-2\) was proposed in Akaike’s 1973 original paper to simplify some calculations involving \(\chi^2\) distributions and was kept around since then.
AIC was devised to “punish” model complexity, i.e models that have more parameters to fit to the data. This is why we add \(2\) to the loglikelihood value for every parameter that the model has. Due to the preference of lower AIC values this penalizes models by the number of parameters, while also making it possible to compare models with different complexities.
Building from the AIC, the Bayesian Information Criterion (BIC) (Schwarz, 1978) uses the same idea, but the penalty term is different:
\[\operatorname{BIC} = -2\log{\hat{\mathcal{L}}} + k\log(n)\]
where \(\hat{\mathcal{L}}\) is the estimated value of the likelihood for a given model and data, \(k\) is the model’s number of parameters, and \(n\) is the number of observations. It is called Bayesian because it uses a “Bayesian” argument to derive the punishment term \(k\log(n)\) in the original 1975 paper.
1.1 Example in Pumas
Let’s go over an example of model fit measures in Pumas.
First, let’s import the following packages:
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, either30,60or90milligrams (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])Population
Subjects: 18
Covariates: wt, isPM, isfed
Observations: dvLet’s create a 2-compartment oral absorption base model with no covariate effects:
base_model = @model begin
@metadata begin
desc = "base model"
timeu = "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)
"""
Absorption rate constant (1/h)
"""
tvka ∈ RealDomain(; lower = 0)
"""
- ΩCL
- ΩVc
- ΩKa
- ΩVp
- ΩQ
"""
Ω ∈ PDiagDomain(5)
"""
Proportional RUV (SD scale)
"""
σₚ ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
CL = tvcl * exp(η[1])
Vc = tvvc * exp(η[2])
Ka = tvka * exp(η[3])
Q = tvq * exp(η[4])
Vp = tvvp * exp(η[5])
end
@dynamics Depots1Central1Periph1
@derived begin
cp := @. 1_000 * (Central / Vc)
"""
Drug Concentration (ng/mL)
"""
dv ~ @. ProportionalNormal(cp, σₚ)
end
endPumasModel
Parameters: tvcl, tvvc, tvvp, tvq, tvka, Ω, σₚ
Random effects: η
Covariates:
Dynamical system variables: Depot, Central, Peripheral
Dynamical system type: Closed form
Derived: dv
Observed: dvLet’s go over the model.
In the @metadata block we are adding a model description and adding information regarding the time units (hours).
Next, we define the model’s parameters in @param while also prepending them with a string that serves as an annotation for the parameter description. This is helpful for post-processing, since Pumas can use the description instead of the parameter name in tables and figures.
Our ηs are defined in the @random block and are sampled from a multivariate normal distribution with mean 0 and a positive-diagonal covariance matrix Ω. We have 5 ηs, one for each PK typical value (also known as θs).
We proceed by defining the individual PK parameters in the @pre block. Each typical value is incremented by the subject’s ηs in a non-linear exponential transformation. This is done to enforce that all individual PK parameters are constrained to being positive. This also has a side effect that the individual PK parameters will be log-normally distributed.
We use the aliased short notation Depots1Central1Periph1 for the ODE system in the @dynamics. This is equivalent to having the following equations:
Depot' = -Ka * Depot
Central' = Ka * Depot - (CL + Q) / Vc * Central + Q / Vp * Peripheral
Peripheral' = Q / Vc * Central - Q / Vp * PeripheralNote that, in order for Depots1Central1Periph1 work correctly, we need to define Ka, CL, Q, Vc, and Vp in the @pre block.
Finally, in the @derived block we define our error model (or likelihood for the statistically-minded). Here we are using a proportional error model with the Gaussian/normal likelihood. Note that Normal is parameterized with mean and standard deviation, not with variance. That’s why we name our proportional error parameter as σₚ and not σ²ₚ.
Let’s now define a initial set of parameter estimates to fit our model:
iparams = (;
tvka = 0.4,
tvcl = 4.0,
tvvc = 70.0,
tvq = 4.0,
tvvp = 50.0,
Ω = Diagonal(fill(0.04, 5)),
σₚ = 0.1,
)(tvka = 0.4, tvcl = 4.0, tvvc = 70.0, tvq = 4.0, tvvp = 50.0, Ω = Diagonal([0.04, 0.04, 0.04, 0.04, 0.04]), σₚ = 0.1)We call the fit function to estimate the parameters of the model:
base_fit = fit(base_model, population, iparams, FOCE())[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 1.630402e+03 2.604358e+02 * time: 0.03934502601623535 1 1.499510e+03 9.365700e+01 * time: 12.627730131149292 2 1.447619e+03 4.714464e+01 * time: 12.644601106643677 3 1.427906e+03 4.439232e+01 * time: 12.66196608543396 4 1.414326e+03 2.726109e+01 * time: 12.677695989608765 5 1.387798e+03 1.159019e+01 * time: 12.694393157958984 6 1.382364e+03 7.060796e+00 * time: 12.710526943206787 7 1.380839e+03 4.839103e+00 * time: 12.72626805305481 8 1.380281e+03 4.075615e+00 * time: 12.741450071334839 9 1.379767e+03 3.303901e+00 * time: 12.756112098693848 10 1.379390e+03 2.856359e+00 * time: 12.771324157714844 11 1.379193e+03 2.650736e+00 * time: 12.960346937179565 12 1.379036e+03 2.523349e+00 * time: 12.974655151367188 13 1.378830e+03 2.638648e+00 * time: 12.989346027374268 14 1.378593e+03 3.463990e+00 * time: 13.00433897972107 15 1.378335e+03 3.471127e+00 * time: 13.01966404914856 16 1.378143e+03 2.756670e+00 * time: 13.035571098327637 17 1.378019e+03 2.541343e+00 * time: 13.051499128341675 18 1.377888e+03 2.163251e+00 * time: 13.066779136657715 19 1.377754e+03 2.571076e+00 * time: 13.082422971725464 20 1.377620e+03 3.370764e+00 * time: 13.097433090209961 21 1.377413e+03 3.938291e+00 * time: 13.11287808418274 22 1.377094e+03 4.458016e+00 * time: 13.199388980865479 23 1.376674e+03 5.713348e+00 * time: 13.215831995010376 24 1.375946e+03 5.417530e+00 * time: 13.233389139175415 25 1.375343e+03 5.862876e+00 * time: 13.25101900100708 26 1.374689e+03 5.717165e+00 * time: 13.267974138259888 27 1.374056e+03 4.400490e+00 * time: 13.285594940185547 28 1.373510e+03 2.191437e+00 * time: 13.304062128067017 29 1.373277e+03 1.203587e+00 * time: 13.322754144668579 30 1.373233e+03 1.157761e+00 * time: 13.340660095214844 31 1.373218e+03 8.770728e-01 * time: 13.35707712173462 32 1.373204e+03 8.021952e-01 * time: 13.391674041748047 33 1.373190e+03 6.613857e-01 * time: 13.409414052963257 34 1.373183e+03 7.602394e-01 * time: 13.426414966583252 35 1.373173e+03 8.552154e-01 * time: 13.443189144134521 36 1.373162e+03 6.961928e-01 * time: 13.460067987442017 37 1.373152e+03 3.162546e-01 * time: 13.477543115615845 38 1.373148e+03 1.747381e-01 * time: 13.494282960891724 39 1.373147e+03 1.258699e-01 * time: 13.525418043136597 40 1.373147e+03 1.074908e-01 * time: 13.541476011276245 41 1.373147e+03 6.799619e-02 * time: 13.557314157485962 42 1.373147e+03 1.819329e-02 * time: 13.573434114456177 43 1.373147e+03 1.338880e-02 * time: 13.58825397491455 44 1.373147e+03 1.370144e-02 * time: 13.602304935455322 45 1.373147e+03 1.315666e-02 * time: 13.630499124526978 46 1.373147e+03 1.065953e-02 * time: 13.64494013786316 47 1.373147e+03 1.069775e-02 * time: 13.659554958343506 48 1.373147e+03 6.234846e-03 * time: 13.673542976379395 49 1.373147e+03 1.512139e-03 * time: 13.687623977661133 50 1.373147e+03 1.115831e-04 * time: 13.700510025024414
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 18
Observation records: Active Missing
dv: 270 0
Total: 270 0
Number of parameters: Constant Optimized
0 11
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -1373.1468
-----------------
Estimate
-----------------
tvcl 2.8346
tvvc 77.8
tvvp 48.753
tvq 3.9784
tvka 1.0279
Ω₁,₁ 0.26382
Ω₂,₂ 0.22877
Ω₃,₃ 0.40034
Ω₄,₄ 0.3796
Ω₅,₅ 0.21493
σₚ 0.097806
-----------------Now we are ready to showcase model fit measures. All of these functions should take a result from fit and output a real number.
Let’s start with aic and bic which are included in Pumas:
aic(base_fit)2768.293584154597bic(base_fit)2807.876225703579We are also free to create our own functions if we want to use something different than aic or bic.
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)We can use it on our base_fit model fit result:
ofv(base_fit)2250.06677622407372 Likelihood Ratio Tests
A likelihood-ratio test (LRT) is a statistical hypothesis test used in the field of statistics and probability theory to compare two statistical models and determine which one provides a better fit to a given set of observed data. It is particularly useful in the context of maximum likelihood estimation (MLE) and is commonly used for hypothesis testing in parametric statistical modeling.
The basic idea behind the likelihood ratio test is to compare the likelihoods of two competing models:
Null Hypothesis (\(H_0\)): This is the model that you want to test against. It represents a specific set of parameter values or restrictions on the model.
Alternative Hypothesis (\(H_a\)): This is the alternative model, often a more complex one or the one you want to support.
The test statistic is calculated as the ratio of the likelihood under the alternative model (\(H_a\)) to the likelihood under the null model (\(H_0\)). Mathematically, it can be expressed as:
\[\operatorname{LRT} = - 2 \log \left( \frac{\mathcal{L}(H_0)}{\mathcal{L}(H_a)} \right)\]
where:
- \(\operatorname{LRT}\): likelihood ratio test statistic
- \(\mathcal{L}(H_0)\): likelihood under \(H_0\), the likelihood of the data under the null hypothesis
- \(\mathcal{L}(H_a)\): likelihood under \(H_a\), the likelihood of the data under the alternative hypothesis
The LRT statistic follows a \(\chi^2\) (chi-squared) distribution with degrees of freedom equal to the difference in the number of parameters between the two models (i.e., the degrees of freedom is the number of additional parameters in the alternative model). In practice, you compare the LRT statistic to \(\chi^2\) distribution to determine whether the alternative model is a significantly better fit to the data than the null model.
The key idea is that if the p-value derived from the LRT statistic is lower than your desired \(\alpha\) (the type-1 error rate, commonly set to \(0.05\)), you would reject the null hypothesis in favor of the alternative hypothesis, indicating that the alternative model provides a better fit to the data.
Note
The likelihood-ratio test requires that the models be nested, i.e. the more complex model can be transformed into the simpler model by imposing constraints on the former’s parameters.
This is generally the case when performing LRT in a covariate selection context. However, be mindful of not violating this assumption when performing LRT.
2.1 Example in Pumas
Pumas provides us with the lrtest function to perform LRT. It takes 2 positional arguments as competing models:
- Model under \(H_0\) (i.e. the model with less parameters)
- Model under \(H_a\) (i.e. the model with more parameters)
Let’s define a covariate model that takes wt into consideration for all the clearance and volume PK parameters:
covariate_model = @model begin
@metadata begin
desc = "covariate model that uses weight covariate information"
timeu = "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)
"""
Absorption rate constant (h-1)
"""
tvka ∈ RealDomain(; lower = 0)
"""
Power exponent on weight for Clearance # new
"""
dwtcl ∈ RealDomain() # new
"""
Power exponent on weight for Distributional Clearance # new
"""
dwtq ∈ RealDomain() # new
"""
- ΩCL
- ΩVc
- ΩKa
- ΩVp
- ΩQ
"""
Ω ∈ PDiagDomain(5)
"""
Proportional RUV (SD scale)
"""
σₚ ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates begin
"""
Weight (kg) # new
"""
wt # new
end
@pre begin
CL = tvcl * exp(η[1]) * (wt / 70)^dwtcl # new
Vc = tvvc * exp(η[2]) * (wt / 70) # new
Ka = tvka * exp(η[3])
Q = tvq * exp(η[4]) * (wt / 70)^dwtq # new
Vp = tvvp * exp(η[5]) * (wt / 70) # new
end
@dynamics Depots1Central1Periph1
@derived begin
cp := @. 1000 * (Central / Vc)
"""
Drug Concentration (ng/mL)
"""
dv ~ @. ProportionalNormal(cp, σₚ)
end
endPumasModel
Parameters: tvcl, tvvc, tvvp, tvq, tvka, dwtcl, dwtq, Ω, σₚ
Random effects: η
Covariates: wt
Dynamical system variables: Depot, Central, Peripheral
Dynamical system type: Closed form
Derived: dv
Observed: dvThis is almost the same model as before. However, we are adding a few tweaks (commented with # new):
wtin the new@covariatesblock- allometric scaling based on
wtfor the individual PK parametersCL,Q,VcandVp - new parameters in
@paramfor the exponent of the power function ofwton both individual clearance PK parametersCLandQ
Since covariate_model has two new parameters in the @param block, we need to add them to the initial set of parameter estimates. We can do this by creating a new NamedTuple that builds upon the last one iparams, while also adding initial values for dwtcl and dwtq:
iparams_covariate = (; iparams..., dwtcl = 0.75, dwtq = 0.75)(tvka = 0.4, tvcl = 4.0, tvvc = 70.0, tvq = 4.0, tvvp = 50.0, Ω = Diagonal([0.04, 0.04, 0.04, 0.04, 0.04]), σₚ = 0.1, dwtcl = 0.75, dwtq = 0.75)
Tip
We are using Julia’s splatting ... operator to expand inline the iparams NamedTuple.
Now we fit our covariate_model:
covariate_fit = fit(covariate_model, population, iparams_covariate, FOCE())[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 1.555051e+03 2.584685e+02 * time: 1.9073486328125e-5 1 1.436886e+03 9.959639e+01 * time: 1.320741891860962 2 1.383250e+03 3.318037e+01 * time: 1.3366260528564453 3 1.372961e+03 2.525341e+01 * time: 1.3535380363464355 4 1.365242e+03 2.081002e+01 * time: 1.4779400825500488 5 1.350200e+03 1.667386e+01 * time: 1.494136095046997 6 1.346374e+03 9.195785e+00 * time: 1.5099759101867676 7 1.344738e+03 8.614309e+00 * time: 1.5259878635406494 8 1.343902e+03 4.950745e+00 * time: 1.540497064590454 9 1.343662e+03 1.478699e+00 * time: 1.555027961730957 10 1.343626e+03 9.575005e-01 * time: 1.5687479972839355 11 1.343609e+03 8.509968e-01 * time: 1.583216905593872 12 1.343589e+03 7.964671e-01 * time: 1.5969979763031006 13 1.343567e+03 8.202459e-01 * time: 1.6557340621948242 14 1.343550e+03 8.133359e-01 * time: 1.6695630550384521 15 1.343542e+03 6.865506e-01 * time: 1.6831939220428467 16 1.343538e+03 3.869567e-01 * time: 1.697037935256958 17 1.343534e+03 2.805019e-01 * time: 1.711210012435913 18 1.343531e+03 3.271442e-01 * time: 1.7248759269714355 19 1.343529e+03 4.584302e-01 * time: 1.7383098602294922 20 1.343527e+03 3.951940e-01 * time: 1.766157865524292 21 1.343525e+03 1.928385e-01 * time: 1.7797300815582275 22 1.343524e+03 1.958575e-01 * time: 1.7934129238128662 23 1.343523e+03 2.008844e-01 * time: 1.806931972503662 24 1.343522e+03 1.636364e-01 * time: 1.8202509880065918 25 1.343522e+03 1.041929e-01 * time: 1.8335609436035156 26 1.343521e+03 7.417497e-02 * time: 1.8604719638824463 27 1.343521e+03 7.297961e-02 * time: 1.8748910427093506 28 1.343521e+03 8.109591e-02 * time: 1.8879430294036865 29 1.343520e+03 7.067080e-02 * time: 1.9011099338531494 30 1.343520e+03 5.088025e-02 * time: 1.9144439697265625 31 1.343520e+03 4.980085e-02 * time: 1.9407758712768555 32 1.343520e+03 4.778940e-02 * time: 1.953847885131836 33 1.343520e+03 5.667067e-02 * time: 1.96714186668396 34 1.343520e+03 5.825591e-02 * time: 1.9803788661956787 35 1.343519e+03 5.354660e-02 * time: 1.9933290481567383 36 1.343519e+03 5.300792e-02 * time: 2.0194129943847656 37 1.343519e+03 4.011720e-02 * time: 2.032459020614624 38 1.343519e+03 3.606197e-02 * time: 2.0455620288848877 39 1.343519e+03 3.546034e-02 * time: 2.0584628582000732 40 1.343519e+03 3.525307e-02 * time: 2.071537971496582 41 1.343519e+03 3.468091e-02 * time: 2.0976250171661377 42 1.343519e+03 3.313732e-02 * time: 2.110214948654175 43 1.343518e+03 4.524162e-02 * time: 2.1230239868164062 44 1.343518e+03 5.769309e-02 * time: 2.136017084121704 45 1.343518e+03 5.716613e-02 * time: 2.1490628719329834 46 1.343517e+03 4.600797e-02 * time: 2.1744120121002197 47 1.343517e+03 3.221948e-02 * time: 2.1870129108428955 48 1.343517e+03 2.610758e-02 * time: 2.1994380950927734 49 1.343517e+03 2.120270e-02 * time: 2.212049961090088 50 1.343517e+03 1.887916e-02 * time: 2.224752902984619 51 1.343517e+03 1.229271e-02 * time: 2.2502379417419434 52 1.343517e+03 4.778802e-03 * time: 2.2631869316101074 53 1.343517e+03 2.158460e-03 * time: 2.2757608890533447 54 1.343517e+03 3.861385e-04 * time: 2.288145065307617
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 18
Observation records: Active Missing
dv: 270 0
Total: 270 0
Number of parameters: Constant Optimized
0 13
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -1343.5173
------------------
Estimate
------------------
tvcl 2.7287
tvvc 70.681
tvvp 47.396
tvq 4.0573
tvka 0.98726
dwtcl 0.5835
dwtq 1.1761
Ω₁,₁ 0.21431
Ω₂,₂ 0.050414
Ω₃,₃ 0.42467
Ω₄,₄ 0.040325
Ω₅,₅ 0.045982
σₚ 0.097904
------------------Now we are ready to perform LRT with lrtest:
mytest = lrtest(base_fit, covariate_fit)Statistic: 59.3
Degrees of freedom: 2
P-value: 0.0The degrees of freedom of the underlying \(\chi^2\) distribution is \(2\), i.e. we have two additional parameters in the model under \(H_a\); and the test statistic is \(59.3\).
The \(p\)-value corresponding for the test statistic and degree of freedom is very close to \(0\). It prints as 0.0, but we can access the value with the pvalue function:
pvalue(mytest)1.3554723626395645e-13This indicates strong evidence against the base_model (i.e. model under \(H_0\)) and in favor of the covariate_model (i.e. model under \(H_a\)).
3 Model Selection Algorithms
There are several model selection techniques that take into account covariate selection. In the statistical literature, the reader can check Thayer (1990), and for the pharmacometric context, the reader can check Hutmacher & Kowalski (2015) and Jonsson & Karlsson (1998).
Pumas currently only implements the Stepwise Covariate Model (SCM). SCM, also known as stepwise procedures, is a model building strategy that is used to identify the best covariate model for a given dataset by a series of iterations (Hutmacher & Kowalski, 2015). Broadly, there are two main types of SCM:
- Forward Selection (FS)
- Backward Elimination (BE)
We will be covering these in detail in a new set of tutorials, please check tutorials.pumas.ai.
4 References
Akaike, H. (1973). Information theory and the extension of the maximum likelihood principle. Proceedings of the Second International Symposium on Information Theory.
Hutmacher, M. M., & Kowalski, K. G. (2015). Covariate selection in pharmacometric analyses: a review of methods. British journal of clinical pharmacology, 79(1), 132–147. https://doi.org/10.1111/bcp.12451
Jonsson, E. N., & Karlsson, M. O. (1998). Automated covariate model building within NONMEM. Pharmaceutical research, 15(9), 1463–1468. https://doi.org/10.1023/a:1011970125687
Schwarz, Gideon E. (1978). Estimating the dimension of a model. Annals of Statistics, 6 (2): 461–464, doi:10.1214/aos/1176344136.
Thayer, J. D. (1990). Implementing Variable Selection Techniques in Regression. ERIC.