using CSV, DataFrames, DataFramesMeta, LinearAlgebra
using Random, CategoricalArrays, CairoMakie, Pumas
using Pumas: cv_populations
using StatsBase: zscore
Random.seed!(378)
LASSO for Covariate Modeling in NLME Models
1 Introduction
LASSO, which stands for Least Absolute Shrinkage and Selection Operator (Ribbing et al. 2007), is a powerful covariate selection method in statistics. It enhances the base model by adding a penalty term to the objective function value, aka loss function, which is proportional to the absolute value of the coefficients. This penalty term enables the shrinking of some coefficients to be close to zero. As a result, LASSO helps in reducing model complexity and selecting a subset of relevant covariates from a potentially large set of covariates. Note that the LASSO solution may not be unique but it may still offer some insights into the data. The covariates with a coefficient close to 0 can be interpreted as not adding any more predictive power to the model, assuming the remaining covariates with non-zero coefficients are already present in the model.
Assuming the LASSO solution included a non-zero coefficient for covariate \(x_1\), LASSO can also be used to answer questions such as: “does including covariate \(x_2\) replace the need for including covariate \(x_1\) in the model?” One can answer this question by not penalizing the coefficient of \(x_2\), performing LASSO, and then observing the effect of that on the coefficient of \(x_1\) which was previously selected when all coefficients were equally penalized. Not penalizing the coefficient of \(x_2\) makes it more favorable to be selected by LASSO if it adds a comparable amount of predictive power to the model as covariate \(x_1\). If covariate \(x_1\) is still significantly more predictive of the response, LASSO will still assign a non-zero coefficient to \(x_1\) despite the coefficient of \(x_2\) not having a penalty.
The rest of this tutorial is organized as follows:
- Reading and pre-processing the data
- Presenting the LASSO PK model
- Cross-validation and hyper-parameter tuning
- Interpreting the results
In this tutorial, the model and dataset from “Pharmacokinetics of Quinidine in Male Patients” (Verme et al. 1992) will be used and adapted. In the original study, the authors analyzed data from 139 adult hospitalized men receiving oral therapy of quinidine, an anti-arrhythmic drug. In previous works, quinidine has shown considerable inter-subject variability in pharmacokinetic (PK) parameters, prompting the search by the authors after predictive covariates, specially for clearance. We will use the resulting data to train a LASSO covariate PK model.
2 Loading packages
First, some environment setup is necessary. The following code brings the packages we need into the environment. A seed for the random number generator is also set for reproducibility.
3 Data processing
3.1 Data reading
We begin by reading the data. The data can be downloaded from Kaggle, an online platform and community for sharing datasets and trained models. Alternatively, the data may be obtained from the RDatasets.jl package, which makes available to Julia many datasets already in the R programming language. Right from the start, we take the opportunity for a bit of processing by removing the first column, which just counts lines. Also, this dataset represents missing values with the string “NA”, so we indicate that in the call to CSV.read().
df = CSV.read("./Quinidine.csv", DataFrame; missingstring = ["NA"])[:, 2:end]
describe(df)[:, Not(7)]14×6 DataFrame
| Row | variable | mean | min | median | max | nmissing |
|---|---|---|---|---|---|---|
| Symbol | Union… | Any | Union… | Any | Int64 | |
| 1 | Subject | 73.9266 | 1 | 74.0 | 223 | 0 |
| 2 | time | 373.252 | 0.0 | 60.0 | 8095.5 | 0 |
| 3 | conc | 2.45429 | 0.4 | 2.3 | 9.4 | 1110 |
| 4 | dose | 225.302 | 83 | 201.0 | 603 | 443 |
| 5 | interval | 7.10843 | 4 | 6.0 | 12 | 1222 |
| 6 | Age | 66.7322 | 42 | 66.0 | 92 | 0 |
| 7 | Height | 69.1931 | 60 | 69.0 | 79 | 0 |
| 8 | Weight | 79.6608 | 41 | 78.0 | 119 | 0 |
| 9 | Race | Black | Latin | 0 | ||
| 10 | Smoke | no | yes | 0 | ||
| 11 | Ethanol | current | none | 0 | ||
| 12 | Heart | Moderate | Severe | 0 | ||
| 13 | Creatinine | < 50 | >= 50 | 0 | ||
| 14 | glyco | 1.28168 | 0.39 | 1.23 | 3.16 | 0 |
3.2 Categorical variable recoding
As can be seen by the dataset summary printed previously, some preprocessing is still necessary. First, some columns include categories as strings, which will be encoded as integers. Also, two new columns are necessary for future steps, namely the compartment (cmt) to which the drug is administered and evid indicating dosing events.
@chain df begin
# Encode string categories as integers
@transform! :Smoke = recode(:Smoke, "no" => 0, "yes" => 1)
@transform! :Heart = recode(:Heart, "No/Mild" => 1, "Moderate" => 2, "Severe" => 3)
@transform! :Creatinine = recode(:Creatinine, ">= 50" => 0, "< 50" => 1)
@transform! :Race = recode(:Race, "Caucasian" => 1, "Latin" => 2, "Black" => 3)
@transform! :Ethanol = recode(:Ethanol, "none" => 1, "former" => 2, "current" => 3)
# Fields for pumas_read()
@transform! :cmt = ifelse(ismissing(:dose), missing, 1)
@rtransform! :evid = ismissing(:dose) ? 0 : 1
end
describe(df)[:, Not(7)]16×6 DataFrame
| Row | variable | mean | min | median | max | nmissing |
|---|---|---|---|---|---|---|
| Symbol | Float64 | Real | Float64 | Real | Int64 | |
| 1 | Subject | 73.9266 | 1 | 74.0 | 223 | 0 |
| 2 | time | 373.252 | 0.0 | 60.0 | 8095.5 | 0 |
| 3 | conc | 2.45429 | 0.4 | 2.3 | 9.4 | 1110 |
| 4 | dose | 225.302 | 83 | 201.0 | 603 | 443 |
| 5 | interval | 7.10843 | 4 | 6.0 | 12 | 1222 |
| 6 | Age | 66.7322 | 42 | 66.0 | 92 | 0 |
| 7 | Height | 69.1931 | 60 | 69.0 | 79 | 0 |
| 8 | Weight | 79.6608 | 41 | 78.0 | 119 | 0 |
| 9 | Race | 1.42284 | 1 | 1.0 | 3 | 0 |
| 10 | Smoke | 0.303875 | 0 | 0.0 | 1 | 0 |
| 11 | Ethanol | 1.45615 | 1 | 1.0 | 3 | 0 |
| 12 | Heart | 1.93202 | 1 | 2.0 | 3 | 0 |
| 13 | Creatinine | 0.28416 | 0 | 0.0 | 1 | 0 |
| 14 | glyco | 1.28168 | 0.39 | 1.23 | 3.16 | 0 |
| 15 | cmt | 1.0 | 1 | 1.0 | 1 | 0 |
| 16 | evid | 0.698844 | 0 | 1.0 | 1 | 0 |
3.3 Data exploration
Since familiarizing with the data we are using can be useful for future troubleshooting, we’ll do a bit of exploration. In the code below, the dataset is first sorted by subject IDs. Then the number of unique values per column is determined and listed.
sort!(df, :Subject)
uniqs = Dict([Symbol(n) => unique(col) for (n, col) in zip(names(df), eachcol(df))])
for (name, uni) in uniqs
println("$(rpad(name, 12))$(length(uni))")
endSmoke 2
interval 5
Heart 3
time 728
Creatinine 2
evid 2
Age 38
Race 3
Height 20
glyco 145
cmt 1
Ethanol 3
Subject 136
dose 16
Weight 58
conc 57And to get a better idea of the distributions of variables, we can check the histograms of all columns.
CairoMakie.activate!()
fig = Figure()
for (title, cartInd) in zip(names(df), CartesianIndices((4, 4)))
ax = Axis(fig[cartInd.I...]; title)
hist!(ax, df[!, title] |> skipmissing |> collect)
end
fig
We can also check how each patient in the dataset looks. In the following block, columns associated with pharmacokinetics are set aside, and shown for a subject.
pkDF = select(df, ["Subject", "time", "dose", "conc", "evid"])
pkDFsubGroup = groupby(pkDF, :Subject) # Group by subject
pkDFsubGroup[1] # Show rows from subject16×5 SubDataFrame
| Row | Subject | time | dose | conc | evid |
|---|---|---|---|---|---|
| Int64 | Float64 | Int64? | Float64? | Int64 | |
| 1 | 1 | 0.0 | 249 | missing | 1 |
| 2 | 1 | 3.0 | 249 | missing | 1 |
| 3 | 1 | 6.0 | 249 | missing | 1 |
| 4 | 1 | 16.0 | 201 | missing | 1 |
| 5 | 1 | 24.0 | 201 | missing | 1 |
| 6 | 1 | 32.0 | 201 | missing | 1 |
| 7 | 1 | 40.0 | 201 | missing | 1 |
| 8 | 1 | 49.25 | missing | 1.4 | 0 |
| 9 | 1 | 138.0 | 201 | missing | 1 |
| 10 | 1 | 145.5 | 201 | missing | 1 |
| 11 | 1 | 150.0 | 201 | missing | 1 |
| 12 | 1 | 160.0 | 402 | missing | 1 |
| 13 | 1 | 168.0 | 402 | missing | 1 |
| 14 | 1 | 176.0 | 402 | missing | 1 |
| 15 | 1 | 184.0 | 402 | missing | 1 |
| 16 | 1 | 193.33 | missing | 3.0 | 0 |
The same is done to the covariate columns in the code below to check how they look.
# Covariates
covas = [:Smoke, :Heart, :Creatinine, :Age, :Race, :Height, :glyco, :Ethanol, :Weight]
covDF = select(df, [:Subject, covas...]) # New DataFrame with covariates and subject column
covSubGroup = groupby(covDF, :Subject) # Group by subjects
covSubGroup[1] # Show rows from subject16×10 SubDataFrame
| Row | Subject | Smoke | Heart | Creatinine | Age | Race | Height | glyco | Ethanol | Weight |
|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Int64 | Int64 | Int64 | Int64 | Int64 | Int64 | Float64 | Int64 | Int64 | |
| 1 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 0.41 | 3 | 106 |
| 2 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 0.41 | 3 | 106 |
| 3 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 0.41 | 3 | 106 |
| 4 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 0.41 | 3 | 106 |
| 5 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 0.41 | 3 | 106 |
| 6 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 0.41 | 3 | 106 |
| 7 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 0.41 | 3 | 106 |
| 8 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 0.41 | 3 | 106 |
| 9 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 1.11 | 3 | 106 |
| 10 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 1.11 | 3 | 106 |
| 11 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 1.11 | 3 | 106 |
| 12 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 1.11 | 3 | 106 |
| 13 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 1.11 | 3 | 106 |
| 14 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 1.11 | 3 | 106 |
| 15 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 1.11 | 3 | 106 |
| 16 | 1 | 0 | 2 | 0 | 60 | 1 | 69 | 1.11 | 3 | 106 |
3.4 Data encoding and standardization
Now, let’s actually prepare the data for the models. We begin by taking only the first row for each subject and dropping the :Subject column.
baselineInput = unique(covDF, :Subject)[:, Not(:Subject)]All covariates must be coerced into numerically continuous values. Continuous and ordinal discrete variables need not be modified. Nominal categorical variables with \(n\) values are one-hot encoded using \(n - 1\) binary variables. Consider :Race as an example which has 3 possible values. The one-hot encoding of :Race involving representing the race of a subject using 2 binary covariates:
:Race__1: 1 if and only if the race has the value 1 (Caucasian):Race__2: 1 if and only if the race has the value 2 (Latin)
The third race (Black) is considered the reference value. The covariate effect of :Race__1 is therefore interpreted as the effect of the Caucasian race relative to the reference Black race. Similarly, the covariate effect of :Race__2 is interpreted as the effect of the Latin race relative to the reference Black race.
After encoding all covariates using continuous values, standardization is performed to ensure all covariates are on the same scale. As examples from the data summaries earlier, ‘glyco’ (\(\alpha_1\)-acid glycoprotein) is in the [0.39, 3.16] mg/100 dL range, while ‘weight’ is in [41, 119] kg. In LASSO, it is essential for all covariates to have the same scale to ensure their covariate effects are fairly penalized. After standardization, all the covariates will have a mean of 0 and standard deviation of 1. If a covariate has the same value across all subjects (standard deviation of 0), it should be filtered out.
stdInputs = @chain baselineInput begin
# One-hot encode the nominal categorical :Race covariate. Its 3 levels are
# represented by 2 binary covariates; the last level (3, Black) is the
# reference and is dropped.
@rtransform begin
:Race__1 = :Race == 1
:Race__2 = :Race == 2
end
@select Not(:Race)
# Also standardize, since covariates have different units/sizes
mapcols(zscore, _)
end
# Check results
describe(stdInputs, :mean, :std, :min, :max)10×5 DataFrame
| Row | variable | mean | std | min | max |
|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Float64 | Float64 | |
| 1 | Smoke | 1.3388e-16 | 1.0 | -0.700621 | 1.41681 |
| 2 | Heart | -4.24497e-17 | 1.0 | -1.05675 | 1.33855 |
| 3 | Creatinine | 8.08177e-17 | 1.0 | -0.59779 | 1.66053 |
| 4 | Age | 2.12249e-17 | 1.0 | -2.79031 | 2.81836 |
| 5 | Height | 1.97228e-15 | 1.0 | -2.85786 | 2.77932 |
| 6 | glyco | -3.42047e-16 | 1.0 | -1.72042 | 4.27123 |
| 7 | Ethanol | 7.1838e-17 | 1.0 | -0.653312 | 2.21283 |
| 8 | Weight | 1.7633e-16 | 1.0 | -2.46618 | 2.51976 |
| 9 | Race__1 | 3.91843e-17 | 1.0 | -1.41681 | 0.700621 |
| 10 | Race__2 | 4.73477e-17 | 1.0 | -0.586504 | 1.69248 |
WarningStandardize within the cross-validation folds
The mean and standard deviation above are computed once over all subjects, before the cross-validation split performed later in this tutorial. The scaling of the training data is therefore influenced by the validation subjects, which is a form of data leakage and makes the validation log likelihood slightly optimistic.
We keep the simpler form here to stay focused on LASSO itself, but a stricter analysis would recompute the mean and standard deviation from each training fold alone and apply those values to the corresponding validation fold.
Finally, let’s combine all the standardized covariates together in a single column of vectors \(x\):
x_names = names(stdInputs)
stdInputs.x = Vector{Float64}.(eachrow(stdInputs))
stdInputs.Subject = unique(df.Subject)
cov_df = leftjoin(df, stdInputs[:, [:Subject, :x]], on = :Subject)4 Covariate-free model
4.1 Model definition
Based on (Verme et al. 1992), below is the definition for a Pumas PK model with one compartment, first-order absorption and elimination, 3 fixed effect parameters (clearance tvcl, central volume tvvc, and absorption constant tvka), associated random effects, and an additive residual error term (\(\sigma\)). Additionally, reference parameter values were extracted from (Pinheiro and Bates 2000) to be used as initial parameter values when fitting.
model = @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvka ∈ RealDomain(; lower = 0)
Ω ∈ PSDDomain(3)
σ ∈ RealDomain(; lower = 0)
end
@random η ~ MvNormal(Ω)
@pre begin
CL = tvcl * exp(η[1])
Vc = tvvc * exp(η[2])
Ka = tvka * exp(η[3])
end
@dynamics Depots1Central1
@derived conc ~ @. Normal(Central / Vc, σ)
end
# Parameters
params = (;
tvcl = ℯ^2.5, # 12.18
tvvc = ℯ^5.4, # 221
tvka = ℯ^-0.21, # 0.81
Ω = [0.056 0.0 0.0; 0.0 0.03 0.0; 0.0 0.0 0.01],
σ = 0.0123,
)The above model does not include any covariates yet. It is generally a good idea to fit such a model first and obtain values for the population parameters before introducing covariates. The optimal population parameter values when fitting the covariate-free model can then be used as initial parameter values when fitting the model with covariates included.
4.2 Pumas population
To fit the above Pumas model to the data, the data needs to be converted to a Pumas Population using the read_pumas function:
pop = read_pumas(
df;
id = :Subject,
time = :time,
amt = :dose,
cmt = :cmt,
evid = :evid,
observations = [:conc],
)Population
Subjects: 136
Observations: concNote that, we created the columns cmt and evid because they are required by read_pumas.
4.3 Model fit
Now, we fit the model using the FOCE algorithm:
fpm = fit(model, pop, params, FOCE(), optim_options = (; show_trace = false))[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed.
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 136
Observation records: Active Missing
conc: 361 82
Total: 361 82
Number of parameters: Constant Optimized
0 10
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -610.76828
-----------------
Estimate
-----------------
tvcl 9.5278
tvvc 35.393
tvka 0.42438
Ω₁,₁ 0.09541
Ω₂,₁ -0.067139
Ω₃,₁ -0.12856
Ω₂,₂ 0.45413
Ω₃,₂ 0.45256
Ω₃,₃ 0.64913
σ 1.102
-----------------5 LASSO model
5.1 Laplace prior
The LASSO objective function value to be minimized is:
\[ -\log p(y \mid x, \beta) + \lambda \cdot \left\lVert \beta \right\rVert_1 \]
where \(\beta\) is the set of covariate effects in the model and \(p(y \mid x, \beta)\) is the marginal likelihood of the model given the population data. An equivalent way to implement LASSO is using a Laplace prior distribution on \(\beta\) and using maximum-a-posteriori (MAP) estimation to optimize for \(\beta\).
\[ \begin{aligned} \beta_i & \sim \text{Laplace}(0, b) \\ -\log p(\beta) & = \sum_{i=1}^{n_\beta} \Bigg( \log (2 \cdot b) + \frac{|\beta_i|}{b} \Bigg) \\ & = \underbrace{n_\beta \cdot \log (2 \cdot b)}_{\text{constant}} + \underbrace{\frac{1}{b}}_{\lambda} \cdot \underbrace{\sum_{i=1}^{n_\beta} |\beta_i|}_{\left\lVert \beta \right\rVert_1} \\ & = \lambda \cdot \left\lVert \beta \right\rVert_1 + \text{constant} \end{aligned} \]
Therefore, by setting the parameter \(b\) of the Laplace prior, one can control the LASSO penalty parameter \(\lambda = \frac{1}{b}\).
5.2 Pumas model
To implement the above prior in a Pumas model, the following model can be used:
λ = 1.0 # LASSO penalty parameter
# multi-variate Laplace distribution
function mvlaplace(λ)
return product_distribution(fill(Distributions.Laplace(0.0, 1 / λ), 10))
end
model_map = @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvka ∈ RealDomain(; lower = 0)
Ω ∈ PSDDomain(3)
σ ∈ RealDomain(; lower = 0)
βcl ~ mvlaplace(λ)
βvc ~ mvlaplace(λ)
βka ~ mvlaplace(λ)
end
@covariates x
@random η ~ MvNormal(Ω)
@pre begin
CL = tvcl * exp(η[1] + dot(x, βcl))
Vc = tvvc * exp(η[2] + dot(x, βvc))
Ka = tvka * exp(η[3] + dot(x, βka))
end
@dynamics Depots1Central1
@derived conc ~ @. Normal(Central / Vc, σ)
endIn the above model, all the standardized covariates are assumed to be in a vector \(x\) for each subject. The vector \(x\) has length 10 for each subject.
5.3 Pumas population with covariates
cov_pop = read_pumas(
cov_df;
id = :Subject,
time = :time,
amt = :dose,
cmt = :cmt,
evid = :evid,
covariates = [:x],
observations = [:conc],
)Population
Subjects: 136
Covariates: x
Observations: conc5.4 Fitting the model
To fit the above model using MAP estimation, call:
# initial parameter values
map_params = merge(init_params(model_map), coef(fpm))
fpm_map =
fit(model_map, cov_pop, map_params, MAP(FOCE()), optim_options = (; show_trace = false))[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed.
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 136
Observation records: Active Missing
conc: 361 82
Total: 361 82
Number of parameters: Constant Optimized
0 40
Estimation method: MAP (FOCE)
-------------------
Estimate
-------------------
tvcl 9.9192
tvvc 24.727
tvka 0.29743
Ω₁,₁ 0.075284
Ω₂,₁ -0.071662
Ω₃,₁ -0.078523
Ω₂,₂ 0.095906
Ω₃,₂ 0.15178
Ω₃,₃ 0.29619
σ 1.0529
βcl₁ -0.0073902
βcl₂ -0.05126
βcl₃ -0.053049
βcl₄ 0.046375
βcl₅ 0.0018842
βcl₆ -0.12291
βcl₇ 0.0064426
βcl₈ 0.045659
βcl₉ 0.12492
βcl₁₀ 0.11307
βvc₁ -0.027897
βvc₂ 0.18895
βvc₃ 0.07224
βvc₄ 0.026123
βvc₅ 0.061366
βvc₆ -0.0067779
βvc₇ -0.25812
βvc₈ -0.0027002
βvc₉ -0.2959
βvc₁₀ 0.11518
βka₁ -7.3121e-7
βka₂ 0.052801
βka₃ 0.066978
βka₄ 0.11268
βka₅ -0.16531
βka₆ 0.15881
βka₇ -0.13743
βka₈ -0.063349
βka₉ -5.9993e-7
βka₁₀ 0.40054
-------------------The following function summarizes the covariate effects and sorts them by their absolute value for any input variable var, where var is either :CL, :Vc or :Ka. A higher covariate effect indicates a relative importance score.
function sort_by_abs(var::Symbol)
coef_map = coef(fpm_map)
effect_summary = DataFrame(
:covariate => x_names,
:CL => coef_map.βcl,
:Vc => coef_map.βvc,
:Ka => coef_map.βka,
)
sort!(effect_summary, [var], by = abs, rev = true)
return effect_summary[:, [:covariate, var]]
end
println(sort_by_abs(:CL), "\n")
println(sort_by_abs(:Vc), "\n")
println(sort_by_abs(:Ka), "\n");10×2 DataFrame Row │ covariate CL │ String Float64 ─────┼───────────────────────── 1 │ Race__1 0.124923 2 │ glyco -0.122905 3 │ Race__2 0.113073 4 │ Creatinine -0.0530487 5 │ Heart -0.05126 6 │ Age 0.0463754 7 │ Weight 0.0456588 8 │ Smoke -0.00739025 9 │ Ethanol 0.00644262 10 │ Height 0.00188418 10×2 DataFrame Row │ covariate Vc │ String Float64 ─────┼───────────────────────── 1 │ Race__1 -0.295897 2 │ Ethanol -0.258124 3 │ Heart 0.18895 4 │ Race__2 0.11518 5 │ Creatinine 0.0722403 6 │ Height 0.0613661 7 │ Smoke -0.0278974 8 │ Age 0.0261234 9 │ glyco -0.00677788 10 │ Weight -0.00270015 10×2 DataFrame Row │ covariate Ka │ String Float64 ─────┼───────────────────────── 1 │ Race__2 0.400536 2 │ Height -0.165307 3 │ glyco 0.158811 4 │ Ethanol -0.137434 5 │ Age 0.112676 6 │ Creatinine 0.0669777 7 │ Weight -0.0633491 8 │ Heart 0.0528015 9 │ Smoke -7.31211e-7 10 │ Race__1 -5.9993e-7
Covariates with an effect whose absolute value is less than a user-specified threshold, e.g. 1e-4, can be dropped from the model. A higher threshold can be used if a smaller model is required.
6 Cross-validation and hyper-parameter selection
6.1 Effect of \(\lambda\)
In the above fit, the value for \(\lambda\) was set to 1.0 upfront. However, that value is a hyper-parameter which needs to be objectively selected. A higher \(\lambda\) has the following effects:
- Covariate effects are more likely to be sparse, and
- The model is more likely to under-fit.
A lower \(\lambda\) has the opposite effects where:
- Covariate effects are less likely to be sparse, and
- The model is more likely to over-fit.
The optimal value of \(\lambda\) is one that maximizes the model’s ability to fit unseen data, not used to fit the model, thus avoiding both under-fitting and over-fitting. To quantify this ability, cross-validation is typically used.
Note that instead of using a single \(\lambda\), one may also use a different \(\lambda\) for each of the covariate effects on CL, Vc and Ka. This will likely give a better model at the end. However, for simplicity, we use a single \(\lambda\) in this tutorial.
6.2 Cross-validation
One common strategy to split the data in the context of hyperparameter search is cross-validation (CV), which Figure 1 illustrates. The hyperparameter optimization procedure can be summarized as follows:
Split the data into \(k\) ‘folds’ (partitions) of roughly equal size.
For each value of \(\lambda\):
For each fold of the \(k\) folds:
- Assign the selected fold as the validation data, and use the remaining folds for fitting.
- Fit/train the model using the training folds and the selected value of \(\lambda\).
- Evaluate the fitted model’s performance on the validation fold, e.g. using the log likelihood of the fitted model given the validation data.
Take the average performance metric across all folds.
Use an optimization or search strategy to choose the best \(\lambda\) that maximizes the average model performance on unseen data.
Re-fit the model using the optimal \(\lambda\) and the entire dataset.
K-fold cross-validation can also be used without changing hyperparameters between iterations to get a measure of a model’s performance on unseen data for model comparison.
While there is generally no easy way to select the number of folds K (Kohavi 2001), values from 5 to 10 are common. Higher values increase the computational cost of the procedure, because more fits need to be done, but reduce the bias when estimating the model’s performance. Lower K values are faster to run but more biased.
6.3 Select \(\lambda\)
The cross-validation procedure presented will be repeatedly applied to determine the best value of \(\lambda\). And the average log likelihood using the validation data will be used as the metric to compare the models.
First some setup is required. In the following code, the function createModel is defined for code organization. It returns a Pumas model using a provided \(\lambda\). Then, Pumas.cv_populations is used to split the population into K = 5 folds. In the split object returned, the fields train_data and validation_data are vectors of populations with one entry per fold: the former contains the 5 training populations; and the latter, the respective 5 validation populations. Lastly, the values of \(\lambda\) to be tried are defined in lambdas, spanning 1e-5 to 1e3 in log10 scale.
function createModel(λ)
return @model begin
@param begin
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvka ∈ RealDomain(; lower = 0)
Ω ∈ PSDDomain(3)
σ ∈ RealDomain(; lower = 0)
βcl ~ mvlaplace(λ)
βvc ~ mvlaplace(λ)
βka ~ mvlaplace(λ)
end
@covariates x
@random η ~ MvNormal(Ω)
@pre begin
CL = tvcl * exp(η[1] + dot(x, βcl))
Vc = tvvc * exp(η[2] + dot(x, βvc))
Ka = tvka * exp(η[3] + dot(x, βka))
end
@dynamics Depots1Central1
@derived conc ~ @. Normal(Central / Vc, σ)
end
end
# Split population
split = cv_populations(cov_pop, KFold(; K = 5), BySubject())
(; train_data, validation_data) = split
lambdas = 10 .^ (range(-5, 3, length = 20))20-element Vector{Float64}:
1.0e-5
2.6366508987303607e-5
6.951927961775606e-5
0.00018329807108324357
0.0004832930238571752
0.0012742749857031334
0.003359818286283781
0.008858667904100823
0.023357214690901226
0.06158482110660264
0.16237767391887217
0.4281332398719394
1.128837891684689
2.976351441631318
7.847599703514613
20.6913808111479
54.5559478116852
143.8449888287663
379.26901907322497
1000.0With these definitions, we are able to run the procedure. For optimization, we use a simple log-scale grid search technique, where we iterate over pre-specified \(\lambda\) values. For each \(\lambda\), k-fold cross-validation is done by looping over the train and test pairs of populations in train_data and validation_data. In each fold:
- The model is built with the current \(\lambda\).
- The initial parameters are defined.
- The model is fitted to training data.
- The validation log likelihood of the fitted model is stored.
For each \(\lambda\), the mean validation log likelihood is calculated, averaged across the \(K\) runs. Runs giving -Inf or NaN validation log likelihood values are assumed to be unlucky runs and ignored. If the number of unlucky runs is significant, the values of \(\lambda\) giving such non-finite validation log likelihood values should be penalized and not selected.
cv_result = map(lambdas) do λ
# cv with current λ
@info "Running cross-validation for λ = $λ."
map(zip(train_data, validation_data)) do (train, test)
# model with current lambda
fold_model = createModel(λ)
fold_params = merge(init_params(fold_model), coef(fpm))
# fit to current train fold
fold_fpm = fit(
fold_model,
train,
fold_params,
MAP(FOCE()),
optim_options = (; show_trace = false),
verbose = false,
)
# evaluate average loglikelihood in current validation fold
return loglikelihood(fold_model, test, coef(fold_fpm), FOCE())
end
end
medFin = mean.(filter.(isfinite, cv_result))
cvDF = DataFrame(λ = lambdas, mean_val_ll = medFin)20×2 DataFrame
| Row | λ | mean_val_ll |
|---|---|---|
| Float64 | Float64 | |
| 1 | 1.0e-5 | -135.65 |
| 2 | 2.63665e-5 | -135.65 |
| 3 | 6.95193e-5 | -135.729 |
| 4 | 0.000183298 | -135.649 |
| 5 | 0.000483293 | -135.65 |
| 6 | 0.00127427 | -135.842 |
| 7 | 0.00335982 | -135.647 |
| 8 | 0.00885867 | -135.643 |
| 9 | 0.0233572 | -135.631 |
| 10 | 0.0615848 | -135.742 |
| 11 | 0.162378 | -135.865 |
| 12 | 0.428133 | -136.154 |
| 13 | 1.12884 | -134.754 |
| 14 | 2.97635 | -131.104 |
| 15 | 7.8476 | -126.51 |
| 16 | 20.6914 | -126.025 |
| 17 | 54.5559 | -125.286 |
| 18 | 143.845 | -120.479 |
| 19 | 379.269 | -120.479 |
| 20 | 1000.0 | -120.479 |
As shown in the final DataFrame, increasing \(\lambda\) improves the model’s performance, with the best mean validation log likelihood at the largest \(\lambda\) values in the grid. Recall that the \(\lambda\) is a penalty coefficient in the following objective function value being minimized when fitting the model:
\[ -\log p(y \mid x, \beta) + \lambda \cdot \left\lVert \beta \right\rVert_1 \]
A high value of \(\lambda = 1000\) penalizes all the covariates, leading to their respective coefficients to be all close to 0, i.e. none of the covariates is significant. This can happen in real analyses where none of the covariates have significant predictive power.
The dataset used above is a real dataset, so we have no control over the result. However, for didactic purposes, let’s simulate data from a presumed covariate model to showcase the next couple of steps that would be carried out in this type of analysis.
Below, simobs uses model_map to simulate observations for cov_pop. The covariate effects (\(\beta\)) of smoke, heart and creatinine are arbitrarily (no biological interpretation intended) and manually specified to impact clearance, volume and absorption, respectively. This way, LASSO will have relationships to identify. The remaining parameters are taken from fpm. The data splits and \(\lambda\) values are then defined.
sim_coef = merge(
(;
βcl = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
βvc = [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
βka = [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
),
coef(fpm),
)
sims = simobs(model_map, cov_pop, sim_coef)
sim_pop = Subject.(sims)
split = cv_populations(sim_pop, KFold(; K = 5), BySubject())
(; train_data, validation_data) = split;Now, the previous steps are repeated to optimize the hyperparameter \(\lambda\) using cross-validation, including taking the average of the validation log likelihood across folds and organizing the results in a DataFrame.
cv_result = map(lambdas) do λ
# cv with current λ
@info "Running cross-validation for λ = $λ."
map(zip(train_data, validation_data)) do (train, test)
# model with current lambda
fold_model = createModel(λ)
fold_params = merge(init_params(fold_model), coef(fpm))
# fit to current train fold
fold_fpm = fit(
fold_model,
train,
fold_params,
MAP(FOCE()),
optim_options = (; show_trace = false),
verbose = false,
)
# evaluate average loglikelihood in current validation fold
return loglikelihood(fold_model, test, coef(fold_fpm), FOCE())
end
end
medFin = mean.(filter.(isfinite, cv_result))
cvDF = DataFrame(λ = lambdas, mean_val_ll = medFin)20×2 DataFrame
| Row | λ | mean_val_ll |
|---|---|---|
| Float64 | Float64 | |
| 1 | 1.0e-5 | -177.504 |
| 2 | 2.63665e-5 | -177.102 |
| 3 | 6.95193e-5 | -177.294 |
| 4 | 0.000183298 | -178.938 |
| 5 | 0.000483293 | -179.884 |
| 6 | 0.00127427 | -178.894 |
| 7 | 0.00335982 | -178.76 |
| 8 | 0.00885867 | -178.754 |
| 9 | 0.0233572 | -178.334 |
| 10 | 0.0615848 | -177.653 |
| 11 | 0.162378 | -176.153 |
| 12 | 0.428133 | -170.076 |
| 13 | 1.12884 | -163.515 |
| 14 | 2.97635 | -169.293 |
| 15 | 7.8476 | -164.503 |
| 16 | 20.6914 | -155.568 |
| 17 | 54.5559 | -167.906 |
| 18 | 143.845 | -176.776 |
| 19 | 379.269 | -208.99 |
| 20 | 1000.0 | -208.99 |
This time, we can see the performance doesn’t improve monotonically with larger \(\lambda\) values. It reaches an optimum with \(\lambda = 20.7\). To finish the analysis, the best \(\lambda\) can be used to fit the model to the entire (simulated) dataset. Note that the non-zero covariate effects learned indeed correspond to the true covariate effects assumed above.
bestλ = sort(cvDF, :mean_val_ll; rev = true)[1, :λ]
model_bestλ = createModel(bestλ)
params_bestλ = merge(init_params(model_bestλ), coef(fpm))
fpm_bestλ = fit(
model_bestλ,
sim_pop,
params_bestλ,
MAP(FOCE()),
optim_options = (; show_trace = false),
verbose = false,
)FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 136
Observation records: Active Missing
conc: 443 0
Total: 443 0
Number of parameters: Constant Optimized
0 40
Estimation method: MAP (FOCE)
--------------------
Estimate
--------------------
tvcl 8.4401
tvvc 52.475
tvka 0.84883
Ω₁,₁ 0.021451
Ω₂,₁ -0.050189
Ω₃,₁ -0.010609
Ω₂,₂ 0.24353
Ω₃,₂ 0.36531
Ω₃,₃ 0.95685
σ 1.1901
βcl₁ 0.94788
βcl₂ 0.00050451
βcl₃ -1.0641e-5
βcl₄ -0.0083133
βcl₅ -0.014887
βcl₆ 0.023541
βcl₇ -0.056767
βcl₈ 5.5956e-5
βcl₉ -4.3096e-5
βcl₁₀ -5.5317e-6
βvc₁ 3.953e-5
βvc₂ 0.6767
βvc₃ -0.41082
βvc₄ -1.3586e-5
βvc₅ -1.2544e-5
βvc₆ -2.0619e-7
βvc₇ -0.0001348
βvc₈ 1.5661e-5
βvc₉ 8.3458e-6
βvc₁₀ 9.9509e-6
βka₁ 1.9806e-6
βka₂ -7.0058e-6
βka₃ 0.23693
βka₄ -1.7341e-5
βka₅ -1.611e-6
βka₆ -4.5222e-6
βka₇ 2.0769e-5
βka₈ 1.7516e-6
βka₉ -1.5247e-5
βka₁₀ 1.6161e-6
--------------------References
Kohavi, Ron. 2001. A Study of Cross-Validation and Bootstrap for Accuracy Estimation and Model Selection. 14 (March).
Pinheiro, José C., and Douglas M. Bates. 2000. Mixed-Effects Models in s and s-PLUS. Springer.
Ribbing, Jakob, Joakim Nyberg, Ola Caster, and E. Niclas Jonsson. 2007. “The Lasso—a Novel Method for Predictive Covariate Model Building in Nonlinear Mixed Effects Models.” Journal of Pharmacokinetics and Pharmacodynamics 34 (4): 485–517. https://doi.org/10.1007/s10928-007-9057-1.
Verme, Catherine N., Thomas M. Ludden, William A. Clementi, and Steven C. Harris. 1992. “Pharmacokinetics of Quinidine in Male Patients: A Population Analysis.” Clinical Pharmacokinetics 22 (6): 468–80. https://doi.org/10.2165/00003088-199222060-00005.