Chapter 3: Coding Population Nonlinear Mixed Effects Models

MI-210 Essentials of Population PKPD M&S — Pumas Edition

Author

Adapted from Marc R. Gastonguay, Ph.D. (Metrum Institute)

Published

April 18, 2026

1 Overview

This chapter covers coding pharmacometric models in Pumas — the Julia equivalent of coding NONMEM control streams via NMTRAN. Topics include:

  • Coding pharmaco-statistical models through the Pumas @model macro
  • Interpreting modeling output
  • The Pumas analytical model library (equivalent to PREDPP)
  • Practice problems: enzyme kinetics, QTc, AST, and PopPK

2 Coding Pharmaco-Statistical Models in Pumas

2.1 NONMEM vs. Pumas Model Structure

In NONMEM, models are specified through control stream records. In Pumas, the @model macro provides a structured, block-based equivalent:

NONMEM Record Pumas Block Purpose
$PROBLEM @metadata Problem description
$INPUT / $DATA read_pumas() Data specification
$THETA @param Fixed effect parameters
$OMEGA @param + @random Random effect variance
$SIGMA @param Residual error parameters
$PK @pre Individual parameter definitions
$PRED @pre + @derived Algebraic model (no dynamics)
$SUBROUTINES @dynamics PK model library
$DES @dynamics begin ... end User-defined ODEs
$ERROR @derived Observation model
$ESTIMATION fit() Estimation method
$COVARIANCE infer() Standard errors
$TABLE inspect() Predictions and residuals

2.2 PRED Model Example (No Dynamics)

The book shows an Emax model coded with $PRED. Here is the NONMEM vs. Pumas side-by-side:

$PROB EMAX MODEL
$INPUT ID DV CONC
$DATA PDDATA.CSV IGNORE=C

$PRED
EMAX=THETA(1)*EXP(ETA(1))
EC50=THETA(2)*EXP(ETA(2))
E=EMAX*CONC/(EC50+CONC)
Y=E + EPS(1)

$THETA (0, 100) (0, 20)   ; EMAX and EC50
$OMEGA 0.4  0.4
$SIGMA 9
$ESTIMATION MAXEVAL=9999 PRINT=10
emax_model = @model begin
    @param begin
        tvEmax ∈ RealDomain(lower = 0.0, init = 100.0)
        tvEC50 ∈ RealDomain(lower = 0.0, init = 20.0)
        Ω      ∈ PDiagDomain(2)
        σ_add  ∈ RealDomain(lower = 0.0, init = 3.0)  # sqrt(9)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @covariates CONC
    @pre begin
        Emax = tvEmax * exp(η[1])
        EC50 = tvEC50 * exp(η[2])
    end
    @derived begin
        mu := @. Emax * CONC / (EC50 + CONC)
        dv ~ @. Normal(mu, σ_add)
    end
end

2.3 PREDPP Model Example (Compartmental Dynamics)

Physiological Meaning of PK Parameters

The one-compartment oral PK model describes a drug that enters a depot compartment (representing the gastrointestinal tract), undergoes first-order absorption at rate \(K_a\) into the systemic circulation (central compartment), and is eliminated by first-order processes governed by clearance (\(CL\)). The volume of distribution (\(V\)) is not a real physiological volume — it is a proportionality constant relating the total amount of drug in the body to the measured plasma concentration: \(V = X / C\). For a 70-kg human, \(V\) can far exceed total body water (~42 L) because drug partitions into tissues according to \(V = V_P + V_T \cdot f_U / f_{UT}\), where \(f_U\) and \(f_{UT}\) are the fractions unbound in plasma and tissue, respectively. Clearance (\(CL\)) is the volume of blood irreversibly cleared of drug per unit time (L/hr). It is not itself a rate; rather, it converts an observed concentration into an elimination rate via \(\text{Rate} = CL \times C\). For a linear system, \(CL = \text{Dose} / \text{AUC}\) and is independent of dose.

The book shows a 1-compartment oral PK model with ADVAN2 TRANS2. In Pumas, Depots1Central1 is the analytical equivalent:

$PROB Pop PK model
$INPUT C ID DV AMT II ADDL TIME
$DATA ../data.CSV IGNORE=C
$SUB ADVAN2 TRANS2

$PK
CL=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
KA=THETA(3)*EXP(ETA(3))
S2=V

$ERROR
Y=F*(1+ERR(1))
IPRED=F

$THETA (0, 8) (0, 50) (0, 0.45)
$OMEGA 0.04 0.04 0.04
$SIGMA 0.04
$ESTIMATION MAXEVAL=9999 METHOD=0 POSTHOC
poppk_model = @model begin
    @param begin
        tvCL ∈ RealDomain(lower = 0.0, init = 8.0)
        tvV  ∈ RealDomain(lower = 0.0, init = 50.0)
        tvKa ∈ RealDomain(lower = 0.0, init = 0.45)
        Ω    ∈ PDiagDomain(3)
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.2)  # sqrt(0.04)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @pre begin
        CL = tvCL * exp(η[1])
        Vc = tvV  * exp(η[2])
        Ka = tvKa * exp(η[3])
    end
    @dynamics Depots1Central1
    @derived begin
        cp := @. Central / Vc
        dv ~ @. Normal(cp, abs(cp) * σ_prop)
    end
end

2.4 The Pumas Analytical Model Library

Pumas provides closed-form analytical solutions equivalent to NONMEM’s PREDPP ADVAN subroutines:

NONMEM ADVAN Pumas Dynamics Description Required @pre
ADVAN1 Central1 1-cpt IV CL, Vc
ADVAN2 Depots1Central1 1-cpt oral CL, Vc, Ka
ADVAN3 Central1Periph1 2-cpt IV CL, Vc, Q, Vp
ADVAN4 Depots1Central1Periph1 2-cpt oral CL, Vc, Q, Vp, Ka
ADVAN11 Central1Periph2 3-cpt IV CL, Vc, Q2, Vp2, Q3, Vp3
ADVAN12 Depots1Central1Periph2 3-cpt oral CL, Vc, Q2, Vp2, Q3, Vp3, Ka
ADVAN6/13 @dynamics begin ... end User ODEs User-defined
Choosing the Number of Compartments

The choice of compartmental model is determined by the shape of the concentration-time curve on a semi-logarithmic plot. A one-compartment model assumes instantaneous distribution throughout the body, producing a monoexponential decline (\(C(t) = C_0 \cdot e^{-Kt}\)) that appears as a single straight line on a log-linear plot. A two-compartment model produces a biexponential decline with a rapid initial \(\alpha\) (distribution) phase followed by a slower \(\beta\) (elimination) phase — two distinct slopes on the semi-log plot. Adding an oral depot compartment (e.g., Depots1Central1, Depots1Central1Periph1) introduces a rising absorption phase governed by \(K_a\); the depot represents the GI tract where drug undergoes first-order absorption. Zero-order absorption (constant-rate input) is used for sustained-release formulations or IV infusions. Importantly, the “number of compartments” refers only to the systemic distribution compartments — the depot is not counted. Model selection should be guided by the data (graphical assessment of the log-concentration profile), not by physiological plausibility alone.

Note

Key differences from NONMEM:

  • No S2=V scaling needed — divide compartment by volume explicitly in @derived
  • Parameter names are case-sensitive and must match exactly: CL, Vc, Ka, Q, Vp
  • TRANS2 parameterization (CL, V, KA) maps directly to Pumas analytical models

2.4.1 Four Ways to Specify Dynamics in Pumas

Beyond the analytical library, Pumas supports multiple equivalent approaches for specifying model dynamics. All produce identical results for the same model:

Fastest. Use when the model matches a library structure.

@dynamics Depots1Central1

Most flexible. Required for nonlinear elimination, indirect response, or any custom dynamics.

@dynamics begin
    Depot'   = -Ka * Depot
    Central' = Ka * Depot - (CL / Vc) * Central
end

Pumas automatically detects linear ODE systems and solves via matrix exponential (faster than numerical integration). Use @options begin checklinear = false end to force numerical integration for nonlinear systems.

Define the ODE symbolically outside the model, enabling structural identifiability analysis and automatic simplification.

using ModelingToolkit

@mtkmodel MTK1Cpt begin
    @parameters begin
        Ka
        CL
        Vc
    end
    @variables begin
        Depot(t) = 0.0
        Central(t) = 0.0
    end
    @equations begin
        D(Depot)   ~ -Ka * Depot
        D(Central) ~ Ka * Depot - (CL / Vc) * Central
    end
end

@mtkbuild sys = MTK1Cpt()

# Then in your model:
@dynamics sys

Define dynamics as chemical/pharmacological reactions — intuitive for PBPK, QSP, and TMDD models.

using Catalyst

rn = @reaction_network begin
    Ka,    Depot --> Central
    CL/Vc, Central --> 0
end

# Then in your model:
@dynamics rn
Which Dynamics Approach to Use?
  • Analytical (closed-form): Default choice for standard 1/2/3-compartment models — fastest execution
  • Explicit ODEs: Nonlinear PK (Michaelis-Menten), indirect response, effect compartment, custom models
  • ModelingToolkit: When you need symbolic manipulation, identifiability analysis, or model composition
  • Catalyst reactions: PBPK, QSP, TMDD — large reaction networks where stoichiometry is natural

2.4.2 Absorption Model Patterns

Beyond first-order absorption (Ka with Depots1Central1), Pumas supports a range of absorption models via @dynamics and @dosecontrol:

Absorption Type Implementation Key Parameters
First-order Depots1Central1 or Depot' = -Ka * Depot Ka
Zero-order @dosecontrol begin duration = (Central = dur,) end dur (duration)
Parallel ZO + FO Dual virtual doses with bioav splitting Ka, dur, Frac
Two parallel FO Two depot compartments with different Ka Ka1, Ka2, Frac
Transit compartments Chain: T1' = -Ktr*T1; T2' = Ktr*T1 - Ktr*T2; … Ktr, n
Erlang @delay(Erlang(N, MAT/N), Central) N (integer shape), MAT
Gamma @delay(Gamma(N, MAT/N), Central) N (continuous shape), MAT
Weibull @delay(Weibull(k, λ), Central) k (shape), λ (scale)

See the Absorption Models tutorial for complete working examples of each.

2.5 Interpreting Modeling Output

In Pumas, model output is accessed programmatically rather than parsed from text files:

NONMEM Output Pumas Equivalent
Minimum OFV loglikelihood(fit) or -2 * loglikelihood(fit)
Parameter estimates coef(fit)
Standard errors infer(fit)
ETA-shrinkage Printed with fit summary
PRED, IPRED, WRES inspect(fit) → DataFrame
AIC aic(fit)
BIC bic(fit)

Useful variance terms:

  • For exponential IIV: \(\%CV \approx \sqrt{\omega^2} \times 100\)
  • For additive IIV: \(SD = \sqrt{\omega^2}\)
  • %RSE: \(\%RSE = \frac{SE}{PE} \times 100\)

3 Practice Problems

Michaelis-Menten Kinetics: The Foundation of Emax and Nonlinear PK

The Michaelis-Menten equation describes the rate of an enzymatic reaction as a function of substrate concentration:

\[V = \frac{V_{\max} \cdot C}{K_M + C}\]

\(V_{\max}\) reflects the enzyme capacity — the maximum rate achievable when all enzyme active sites are saturated. \(K_M\) is the substrate concentration at which the reaction rate is 50% of \(V_{\max}\) and serves as a surrogate for affinity: a low \(K_M\) indicates high affinity of substrate for enzyme. This equation is structurally identical to the \(E_{\max}\) model in pharmacodynamics, where \(E_{\max} \leftrightarrow V_{\max}\) and \(EC_{50} \leftrightarrow K_M\). At low concentrations (\(C \ll K_M\)), the equation simplifies to \(V \approx (V_{\max}/K_M) \cdot C\), yielding approximately linear (first-order) kinetics. At high concentrations (\(C \gg K_M\)), the rate saturates at \(V_{\max}\), producing zero-order kinetics. This saturation is the mechanistic basis for nonlinear (capacity-limited) pharmacokinetics — drugs such as phenytoin and ethanol exhibit dose-disproportional increases in exposure because their metabolizing enzymes become saturated at therapeutic concentrations.

3.1 Problem 1: Enzyme Kinetics (Michaelis-Menten)

Population analysis of in-vitro enzyme activity data in 20 liver slice samples. Concentrations: 0.1, 0.3, 1, 3, 10, 30, 100, 300, 1000 uM.

raw_enz = CSV.read(joinpath(data_dir, "enz3.csv"), DataFrame)

# Fix data entry error (concentration > 1000 → 1000, as in R script)
enz_df = @rtransform(raw_enz,
    :Concentration = :Concentration > 1000 ? 1000.0 : Float64(:Concentration),
)
rename!(enz_df, :ID => :id, :Concentration => :CONC, :ObsActivity => :dv)

fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1], xlabel = "Concentration (uM)", ylabel = "Enzyme Activity")
scatter!(ax, enz_df.CONC, enz_df.dv, color = :steelblue, markersize = 6)
fig
# Use concentration as time axis (non-temporal PRED-style model)
enz_df.time = Float64.(enz_df.CONC)

pop_enz = read_pumas(enz_df; observations = [:dv], covariates = [:CONC], event_data = false)
Population
  Subjects: 20
  Covariates: CONC
  Observations: dv
enz_model = @model begin
    @metadata begin
        desc = "Michaelis-Menten Enzyme Kinetics"
    end
    @param begin
        tvKM   ∈ RealDomain(lower = 0.0, init = 20.0)
        tvVMAX ∈ RealDomain(lower = 0.0, init = 120.0)
        Ω      ∈ PDiagDomain(2)
        σ_add  ∈ RealDomain(lower = 0.0, init = 3.16)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @covariates CONC
    @pre begin
        KM   = tvKM * exp(η[1])
        VMAX = tvVMAX * exp(η[2])
    end
    @derived begin
        ACT := @. VMAX * CONC / (KM + CONC)
        dv ~ @. Normal(ACT, σ_add)
    end
end
PumasModel
  Parameters: tvKM, tvVMAX, Ω, σ_add
  Random effects: η
  Covariates: CONC
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived: dv
  Observed: dv
enz_params = (tvKM = 20.0, tvVMAX = 120.0, Ω = Diagonal([0.04, 0.04]), σ_add = sqrt(10.0))

# FO estimation (METHOD=0)
enz_fo = fit(enz_model, pop_enz, enz_params, Pumas.FO())
enz_fo
[ 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.612324e+03     2.021640e+03
 * time: 0.012627124786376953
     1     7.765242e+02     1.665303e+02
 * time: 1.1877291202545166
     2     7.560971e+02     1.129869e+02
 * time: 1.1880099773406982
     3     7.234432e+02     4.035156e+01
 * time: 1.1882009506225586
     4     7.213581e+02     3.072872e+01
 * time: 1.1884150505065918
     5     7.203405e+02     2.585542e+01
 * time: 1.18858003616333
     6     7.179375e+02     2.887201e+01
 * time: 1.1887269020080566
     7     7.155956e+02     4.163822e+01
 * time: 1.1888940334320068
     8     7.140300e+02     3.219939e+01
 * time: 1.1890709400177002
     9     7.134381e+02     1.624644e+01
 * time: 1.1892240047454834
    10     7.130683e+02     1.093506e+01
 * time: 1.1893811225891113
    11     7.127838e+02     5.855108e+00
 * time: 1.1895339488983154
    12     7.126722e+02     5.831986e+00
 * time: 1.1896910667419434
    13     7.126297e+02     5.750130e+00
 * time: 1.1898210048675537
    14     7.125869e+02     5.652110e+00
 * time: 1.1899559497833252
    15     7.124757e+02     5.432402e+00
 * time: 1.190100908279419
    16     7.122104e+02     9.704714e+00
 * time: 1.1902430057525635
    17     7.115703e+02     1.873484e+01
 * time: 1.1903901100158691
    18     7.101726e+02     2.645589e+01
 * time: 1.190553903579712
    19     7.082110e+02     1.564613e+01
 * time: 1.1907129287719727
    20     7.072692e+02     6.606017e+00
 * time: 1.2259149551391602
    21     7.071223e+02     1.096616e+00
 * time: 1.2261879444122314
    22     7.071019e+02     2.603118e+00
 * time: 1.226362943649292
    23     7.070939e+02     1.204422e+00
 * time: 1.22658109664917
    24     7.070927e+02     1.089222e+00
 * time: 1.2267711162567139
    25     7.070915e+02     1.085128e+00
 * time: 1.226943016052246
    26     7.070860e+02     2.227252e+00
 * time: 1.2270920276641846
    27     7.070741e+02     4.471181e+00
 * time: 1.2272429466247559
    28     7.070422e+02     8.125073e+00
 * time: 1.227381944656372
    29     7.069685e+02     1.286940e+01
 * time: 1.227531909942627
    30     7.068075e+02     1.746722e+01
 * time: 1.227682113647461
    31     7.065105e+02     1.768163e+01
 * time: 1.2278358936309814
    32     7.061779e+02     9.495926e+00
 * time: 1.2279911041259766
    33     7.060588e+02     4.058715e-01
 * time: 1.2281451225280762
    34     7.060361e+02     2.759697e+00
 * time: 1.2283120155334473
    35     7.060203e+02     3.801312e+00
 * time: 1.2284860610961914
    36     7.060017e+02     3.009431e+00
 * time: 1.2286388874053955
    37     7.059896e+02     1.014903e+00
 * time: 1.2287828922271729
    38     7.059856e+02     3.506766e-01
 * time: 1.2289791107177734
    39     7.059839e+02     7.042265e-01
 * time: 1.2291929721832275
    40     7.059826e+02     6.669825e-01
 * time: 1.2294039726257324
    41     7.059818e+02     2.130024e-01
 * time: 1.2296209335327148
    42     7.059815e+02     2.150589e-02
 * time: 1.2298469543457031
    43     7.059813e+02     1.832773e-01
 * time: 1.2300539016723633
    44     7.059812e+02     1.405467e-01
 * time: 1.2302141189575195
    45     7.059812e+02     6.205356e-02
 * time: 1.2303850650787354
    46     7.059811e+02     1.204098e-02
 * time: 1.2305710315704346
    47     7.059811e+02     3.513147e-02
 * time: 1.2307910919189453
    48     7.059811e+02     3.008400e-02
 * time: 1.2309439182281494
    49     7.059811e+02     1.100063e-02
 * time: 1.2310960292816162
    50     7.059811e+02     3.082722e-03
 * time: 1.2312428951263428
    51     7.059811e+02     7.852107e-03
 * time: 1.2313830852508545
    52     7.059811e+02     6.154334e-03
 * time: 1.2315270900726318
    53     7.059811e+02     2.019669e-03
 * time: 1.231679916381836
    54     7.059811e+02     8.889316e-04
 * time: 1.231842041015625

FittedPumasModel

Dynamical system type:          No dynamical model

Number of subjects:                             20

Observation records:         Active        Missing
    dv:                         180              0
    Total:                      180              0

Number of parameters:      Constant      Optimized
                                  0              5

Likelihood approximation:                       FO
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -705.98112

---------------------
         Estimate
---------------------
tvKM      46.341
tvVMAX   104.94
Ω₁,₁       2.9107e-9
Ω₂,₂       0.0059394
σ_add     11.727
---------------------
enz_coefs = coef(enz_fo)
conc_range = range(0.1, 1000, length = 200)
pred_act = enz_coefs.tvVMAX .* conc_range ./ (enz_coefs.tvKM .+ conc_range)

fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1], xlabel = "Concentration (uM)", ylabel = "Enzyme Activity",
    xscale = log10)
scatter!(ax, enz_df.CONC, enz_df.dv, color = :steelblue, markersize = 5, alpha = 0.6)
lines!(ax, collect(conc_range), pred_act, color = :red, linewidth = 2, label = "Population Pred")
axislegend(ax, position = :lt)
fig
Figure 1: Problem 1: Michaelis-Menten Fit to Enzyme Kinetics Data

3.2 Problem 2: QTc Population Models

Population analysis of QTc exposure response data from Phase 1. Three models are compared:

  1. Naive pool (linear) — no IIV, OMEGA BLOCK fixed at 0
  2. Population linear — with IIV on intercept and slope
  3. Population Emax — with IIV on E0, Emax, EC50
dqtc_raw = CSV.read(joinpath(data_dir, "dqtc.csv"), DataFrame)
rename!(dqtc_raw, "SID" => :id, "CONC" => :conc, "dqtc" => :dv)

# Monotonic time per subject (row index)
dqtc_raw.time = Vector{Float64}(undef, nrow(dqtc_raw))
for gdf in groupby(dqtc_raw, :id)
    gdf.time .= Float64.(1:nrow(gdf))
end

pop_dqtc = read_pumas(dqtc_raw; observations = [:dv], covariates = [:conc], event_data = false)
println("Subjects: ", length(pop_dqtc))
Subjects: 64

3.2.1 Model 1: Naive Pool Linear (no IIV)

Equivalent to NONMEM $OMEGA BLOCK(2) 0 / 0 0 FIXED:

qtc_naive = @model begin
    @param begin
        tvINT ∈ RealDomain(init = 1.0)
        tvSLP ∈ RealDomain(init = 0.02)
        σ_add ∈ RealDomain(lower = 0.0, init = 5.0)
    end
    @covariates conc
    @pre begin
        INT = tvINT
        SLP = tvSLP
    end
    @derived begin
        mu := @. INT + SLP * conc
        dv ~ @. Normal(mu, σ_add)
    end
end
PumasModel
  Parameters: tvINT, tvSLP, σ_add
  Random effects:
  Covariates: conc
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived: dv
  Observed: dv
qtc_naive_fit = fit(qtc_naive, pop_dqtc,
    (tvINT = 1.0, tvSLP = 0.02, σ_add = 5.0),
    Pumas.NaivePooled())
qtc_naive_fit
[ 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     2.177739e+03     1.418514e+04
 * time: 5.793571472167969e-5
     1     2.145965e+03     6.258659e+02
 * time: 0.36159205436706543
     2     2.118681e+03     7.104756e+02
 * time: 0.3618330955505371
     3     1.915684e+03     2.055385e+03
 * time: 0.362030029296875
     4     1.856796e+03     2.183562e+03
 * time: 0.36221909523010254
     5     1.856667e+03     8.088581e+03
 * time: 0.362368106842041
     6     1.853771e+03     2.659683e+03
 * time: 0.36250805854797363
     7     1.851798e+03     1.320097e+03
 * time: 0.3626539707183838
     8     1.847591e+03     1.163686e+02
 * time: 0.3627970218658447
     9     1.814822e+03     6.855837e+03
 * time: 0.3629319667816162
    10     1.791657e+03     8.789379e+03
 * time: 0.3630690574645996
    11     1.785043e+03     3.339229e+03
 * time: 0.3632030487060547
    12     1.783504e+03     1.702782e+02
 * time: 0.3633451461791992
    13     1.783244e+03     1.914953e+02
 * time: 0.3634779453277588
    14     1.783242e+03     2.541537e+00
 * time: 0.3636140823364258
    15     1.783242e+03     1.408389e-01
 * time: 0.363753080368042
    16     1.783242e+03     3.340979e-03
 * time: 0.36388707160949707
    17     1.783242e+03     6.584688e-05
 * time: 0.3640329837799072

FittedPumasModel

Dynamical system type:          No dynamical model

Number of subjects:                             64

Observation records:         Active        Missing
    dv:                         812              0
    Total:                      812              0

Number of parameters:      Constant      Optimized
                                  0              3

Likelihood approximation:              NaivePooled
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -1783.2424

------------------
        Estimate
------------------
tvINT   -0.088942
tvSLP    0.016926
σ_add    2.1753
------------------

3.2.2 Model 2: Population Linear (diagonal IIV)

qtc_pop_lin = @model begin
    @param begin
        tvINT ∈ RealDomain(init = 1.0)
        tvSLP ∈ RealDomain(init = 0.2)
        Ω     ∈ PDiagDomain(2)
        σ_add ∈ RealDomain(lower = 0.0, init = 3.16)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @covariates conc
    @pre begin
        INT = tvINT + η[1]
        SLP = tvSLP + η[2]
    end
    @derived begin
        mu := @. INT + SLP * conc
        dv ~ @. Normal(mu, σ_add)
    end
end
PumasModel
  Parameters: tvINT, tvSLP, Ω, σ_add
  Random effects: η
  Covariates: conc
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived: dv
  Observed: dv
qtc_pop_fit = fit(qtc_pop_lin, pop_dqtc,
    (tvINT = 1.0, tvSLP = 0.2, Ω = Diagonal([0.1, 0.2]), σ_add = sqrt(10.0)),
    Pumas.FOCE())
qtc_pop_fit
[ 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.946113e+03     5.631218e+02
 * time: 1.7881393432617188e-5
     1     1.725899e+03     4.212633e+02
 * time: 0.4035990238189697
     2     1.677271e+03     1.943342e+02
 * time: 0.40445590019226074
     3     1.650182e+03     6.334516e+01
 * time: 0.40516090393066406
     4     1.637005e+03     7.223578e+01
 * time: 0.40584707260131836
     5     1.632158e+03     4.568235e+01
 * time: 0.40653204917907715
     6     1.624934e+03     2.013463e+01
 * time: 0.40718603134155273
     7     1.623045e+03     2.024793e+01
 * time: 0.4078540802001953
     8     1.617261e+03     3.906291e+01
 * time: 0.4085230827331543
     9     1.606577e+03     1.060489e+02
 * time: 0.4092130661010742
    10     1.565896e+03     3.676768e+02
 * time: 0.4099729061126709
    11     1.543994e+03     2.209361e+02
 * time: 0.4116358757019043
    12     1.541115e+03     1.818206e+02
 * time: 0.4123380184173584
    13     1.536554e+03     1.049645e+02
 * time: 0.41300487518310547
    14     1.532072e+03     1.480411e+02
 * time: 0.4136829376220703
    15     1.522174e+03     5.041024e+01
 * time: 0.4145689010620117
    16     1.514611e+03     8.654016e+01
 * time: 0.41552090644836426
    17     1.512978e+03     9.238392e+01
 * time: 0.41629791259765625
    18     1.510279e+03     1.252041e+02
 * time: 0.4170680046081543
    19     1.506438e+03     1.084174e+02
 * time: 0.44490694999694824
    20     1.499992e+03     9.634847e+01
 * time: 0.4457409381866455
    21     1.488966e+03     2.485298e+02
 * time: 0.4464840888977051
    22     1.483359e+03     4.427660e+02
 * time: 0.44748592376708984
    23     1.480256e+03     3.254343e+02
 * time: 0.4484560489654541
    24     1.479730e+03     1.634952e+02
 * time: 0.4493720531463623
    25     1.479427e+03     1.382633e+02
 * time: 0.45003294944763184
    26     1.479392e+03     1.275796e+01
 * time: 0.4507009983062744
    27     1.479392e+03     7.709314e-01
 * time: 0.45136499404907227
    28     1.479392e+03     5.331967e-02
 * time: 0.4520750045776367
    29     1.479392e+03     6.878027e-03
 * time: 0.4527130126953125
    30     1.479392e+03     6.363288e-04
 * time: 0.45334696769714355

FittedPumasModel

Dynamical system type:          No dynamical model

Number of subjects:                             64

Observation records:         Active        Missing
    dv:                         812              0
    Total:                      812              0

Number of parameters:      Constant      Optimized
                                  0              5

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -1479.3915

------------------
        Estimate
------------------
tvINT   0.15898
tvSLP   0.013788
Ω₁,₁    1.0842
Ω₂,₂    2.2605e-5
σ_add   1.3379
------------------

3.2.3 Model 3: Population Emax

qtc_emax = @model begin
    @param begin
        tvE0   ∈ RealDomain(init = 0.1)
        tvEmax ∈ RealDomain(init = 40.0)
        tvEC50 ∈ RealDomain(lower = 0.0, init = 1000.0)
        Ω      ∈ PSDDomain(3)
        σ_add  ∈ RealDomain(lower = 0.0, init = 3.16)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @covariates conc
    @pre begin
        E0   = tvE0 + η[1]
        EC50 = tvEC50 * exp(η[2])
        Emax = tvEmax + η[3]
    end
    @derived begin
        mu := @. E0 + Emax * conc / (EC50 + conc)
        dv ~ @. Normal(mu, σ_add)
    end
end
PumasModel
  Parameters: tvE0, tvEmax, tvEC50, Ω, σ_add
  Random effects: η
  Covariates: conc
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived: dv
  Observed: dv
qtc_emax_fit = fit(qtc_emax, pop_dqtc,
    (tvE0 = 0.1, tvEmax = 40.0, tvEC50 = 1000.0,
     Ω = [0.1 0.01 0.01; 0.01 0.1 0.01; 0.01 0.01 0.1],
     σ_add = sqrt(10.0)),
    Pumas.FOCE())
qtc_emax_fit
[ 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.871557e+03     5.256135e+02
 * time: 3.314018249511719e-5
     1     1.726619e+03     6.335882e+02
 * time: 0.3728041648864746
     2     1.632397e+03     2.206707e+02
 * time: 0.37501001358032227
     3     1.601012e+03     7.437475e+01
 * time: 0.37687206268310547
     4     1.590285e+03     4.988657e+01
 * time: 0.7839939594268799
     5     1.584782e+03     5.314367e+01
 * time: 0.7858569622039795
     6     1.570803e+03     2.829331e+01
 * time: 0.7875661849975586
     7     1.569767e+03     1.343165e+01
 * time: 0.7891290187835693
     8     1.567879e+03     3.358515e+01
 * time: 0.7907440662384033
     9     1.563679e+03     7.950324e+01
 * time: 0.7925851345062256
    10     1.558588e+03     1.150816e+02
 * time: 0.7945101261138916
    11     1.554793e+03     7.332243e+01
 * time: 0.7968120574951172
    12     1.550747e+03     1.623739e+01
 * time: 0.7988619804382324
    13     1.548074e+03     3.992431e+01
 * time: 0.801185131072998
    14     1.545245e+03     6.347403e+01
 * time: 0.8033521175384521
    15     1.538866e+03     8.631144e+01
 * time: 0.8057639598846436
    16     1.533309e+03     6.252670e+01
 * time: 0.80873703956604
    17     1.525244e+03     1.922217e+01
 * time: 0.8114571571350098
    18     1.520555e+03     3.148577e+01
 * time: 0.8141369819641113
    19     1.512474e+03     5.730347e+01
 * time: 0.8166530132293701
    20     1.497581e+03     7.174670e+01
 * time: 0.8190209865570068
    21     1.478386e+03     3.842529e+01
 * time: 0.8215301036834717
    22     1.472186e+03     1.127034e+01
 * time: 0.8509399890899658
    23     1.468135e+03     3.955126e+01
 * time: 0.8530151844024658
    24     1.465953e+03     3.408929e+01
 * time: 0.8551170825958252
    25     1.465131e+03     1.046805e+01
 * time: 0.8575210571289062
    26     1.464138e+03     1.625806e+01
 * time: 0.8594839572906494
    27     1.462860e+03     3.585974e+01
 * time: 0.8613979816436768
    28     1.461969e+03     3.464115e+01
 * time: 0.8633530139923096
    29     1.461397e+03     1.734594e+01
 * time: 0.8652641773223877
    30     1.461141e+03     3.733233e+00
 * time: 0.8672211170196533
    31     1.460974e+03     5.899231e+00
 * time: 0.8691320419311523
    32     1.460872e+03     6.726826e+00
 * time: 0.8712961673736572
    33     1.460772e+03     4.665350e+00
 * time: 0.8731729984283447
    34     1.460560e+03     4.449482e+00
 * time: 0.8822321891784668
    35     1.460140e+03     7.954170e+00
 * time: 0.8842029571533203
    36     1.459419e+03     1.550701e+01
 * time: 0.8861100673675537
    37     1.458427e+03     1.652405e+01
 * time: 0.8880491256713867
    38     1.457459e+03     5.999460e+00
 * time: 0.8899681568145752
    39     1.457069e+03     2.192073e+00
 * time: 0.8922080993652344
    40     1.456888e+03     7.414185e+00
 * time: 0.8941121101379395
    41     1.456707e+03     5.304777e+00
 * time: 0.895967960357666
    42     1.456511e+03     2.282375e+00
 * time: 0.89786696434021
    43     1.456292e+03     3.791058e+00
 * time: 0.8998689651489258
    44     1.455902e+03     6.517744e+00
 * time: 0.9083371162414551
    45     1.455467e+03     7.044861e+00
 * time: 0.910163164138794
    46     1.455376e+03     3.528953e+00
 * time: 0.9119560718536377
    47     1.455327e+03     3.173182e-01
 * time: 0.9137330055236816
    48     1.455286e+03     2.731534e+00
 * time: 0.9154629707336426
    49     1.455187e+03     5.604254e+00
 * time: 0.9174880981445312
    50     1.454975e+03     7.157259e+00
 * time: 0.9193241596221924
    51     1.454244e+03     7.270402e+00
 * time: 0.9212431907653809
    52     1.453604e+03     3.790794e+00
 * time: 0.9287431240081787
    53     1.453379e+03     7.986662e-01
 * time: 0.9309720993041992
    54     1.453354e+03     1.990846e-01
 * time: 0.9329590797424316
    55     1.453353e+03     1.949871e-01
 * time: 0.9348721504211426
    56     1.453353e+03     9.176580e-02
 * time: 0.9365851879119873
    57     1.453353e+03     2.141122e-02
 * time: 0.9380781650543213
    58     1.453353e+03     5.188243e-03
 * time: 0.9446201324462891
    59     1.453353e+03     5.188243e-03
 * time: 0.9479260444641113
    60     1.453353e+03     5.188243e-03
 * time: 0.9514420032501221

FittedPumasModel

Dynamical system type:          No dynamical model

Number of subjects:                             64

Observation records:         Active        Missing
    dv:                         812              0
    Total:                      812              0

Number of parameters:      Constant      Optimized
                                  0             10

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:              NoObjectiveChange
Log-likelihood value:                   -1453.3531

---------------------
         Estimate
---------------------
tvE0       -0.038189
tvEmax     75.865
tvEC50   4092.9
Ω₁,₁        1.0971
Ω₂,₁        0.77456
Ω₃,₁       31.495
Ω₂,₂        0.54683
Ω₃,₂       22.235
Ω₃,₃     1499.0
σ_add       1.2882
---------------------

3.2.4 Model Comparison

println("Model Comparison (QTc)")
println("=" ^ 50)
println("Naive Linear  AIC: ", round(aic(qtc_naive_fit), digits = 1))
println("Pop Linear    AIC: ", round(aic(qtc_pop_fit), digits = 1))
println("Pop Emax      AIC: ", round(aic(qtc_emax_fit), digits = 1))
Model Comparison (QTc)
==================================================
Naive Linear  AIC: 3572.5
Pop Linear    AIC: 2968.8
Pop Emax      AIC: 2926.7

3.3 Problem 3: AST Emax Population Model

Population analysis of AST exposure response data from Phase 1.

ast_raw = CSV.read(joinpath(data_dir, "ast.csv"), DataFrame)
rename!(ast_raw, "SID" => :id, "CONC" => :conc, "AST" => :dv)

ast_raw.time = Vector{Float64}(undef, nrow(ast_raw))
for gdf in groupby(ast_raw, :id)
    gdf.time .= Float64.(1:nrow(gdf))
end

pop_ast = read_pumas(ast_raw; observations = [:dv], covariates = [:conc], event_data = false)
Population
  Subjects: 64
  Covariates: conc (heterogenous)
  Observations: dv (heterogenous)
ast_emax = @model begin
    @metadata begin
        desc = "AST Emax Population Model"
    end
    @param begin
        tvE0   ∈ RealDomain(init = 20.0)
        tvEC50 ∈ RealDomain(lower = 0.0, init = 100.0)
        tvEmax ∈ RealDomain(init = 50.0)
        Ω      ∈ PSDDomain(3)
        σ_add  ∈ RealDomain(lower = 0.0, init = 3.16)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @covariates conc
    @pre begin
        E0   = tvE0 + η[1]
        EC50 = tvEC50 * exp(η[2])
        Emax = tvEmax + η[3]
    end
    @derived begin
        mu := @. E0 + Emax * conc / (EC50 + conc)
        dv ~ @. Normal(mu, σ_add)
    end
end
PumasModel
  Parameters: tvE0, tvEC50, tvEmax, Ω, σ_add
  Random effects: η
  Covariates: conc
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived: dv
  Observed: dv
ast_fit = fit(ast_emax, pop_ast,
    (tvE0 = 20.0, tvEC50 = 100.0, tvEmax = 50.0,
     Ω = [0.1 0.01 0.01; 0.01 0.1 0.01; 0.01 0.01 0.1],
     σ_add = sqrt(10.0)),
    Pumas.FOCE())
ast_fit
[ 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.694226e+03     1.577918e+03
 * time: 2.193450927734375e-5
     1     9.191714e+02     9.604917e+01
 * time: 0.5070838928222656
     2     9.067138e+02     8.195403e+01
 * time: 0.5085499286651611
     3     8.858533e+02     5.126383e+01
 * time: 0.5097339153289795
     4     8.783999e+02     6.345648e+01
 * time: 0.5108730792999268
     5     8.671584e+02     8.291959e+01
 * time: 0.5119888782501221
     6     8.538607e+02     7.804764e+01
 * time: 0.5131618976593018
     7     8.429917e+02     4.142100e+01
 * time: 0.5144391059875488
     8     8.406836e+02     9.453546e+00
 * time: 0.5156209468841553
     9     8.405201e+02     9.658065e+00
 * time: 0.5166380405426025
    10     8.399121e+02     1.046413e+01
 * time: 0.5177910327911377
    11     8.394418e+02     1.642671e+01
 * time: 0.5189030170440674
    12     8.390421e+02     1.503612e+01
 * time: 0.5199248790740967
    13     8.381302e+02     1.089462e+01
 * time: 0.5209779739379883
    14     8.358632e+02     1.002887e+01
 * time: 0.5485200881958008
    15     8.323047e+02     1.806133e+01
 * time: 0.5498509407043457
    16     8.310067e+02     1.099528e+01
 * time: 0.5511000156402588
    17     8.307111e+02     9.698122e+00
 * time: 0.5521519184112549
    18     8.301625e+02     9.316964e+00
 * time: 0.5532269477844238
    19     8.297585e+02     1.311702e+01
 * time: 0.5542259216308594
    20     8.277769e+02     2.275459e+01
 * time: 0.5553889274597168
    21     8.248505e+02     2.747646e+01
 * time: 0.556617021560669
    22     8.214115e+02     2.097432e+01
 * time: 0.5578880310058594
    23     8.181504e+02     8.449226e+00
 * time: 0.5591731071472168
    24     8.154305e+02     8.420683e+00
 * time: 0.5605309009552002
    25     8.153312e+02     3.275010e+00
 * time: 0.5616230964660645
    26     8.153082e+02     2.604366e-01
 * time: 0.5627260208129883
    27     8.153079e+02     2.975894e-01
 * time: 0.5637369155883789
    28     8.153076e+02     3.199933e-01
 * time: 0.5647048950195312
    29     8.153042e+02     4.871807e-01
 * time: 0.5657429695129395
    30     8.152983e+02     7.268925e-01
 * time: 0.5668799877166748
    31     8.152819e+02     1.340362e+00
 * time: 0.5681259632110596
    32     8.152525e+02     2.005302e+00
 * time: 0.5694329738616943
    33     8.152044e+02     2.425633e+00
 * time: 0.5875449180603027
    34     8.151088e+02     2.321137e+00
 * time: 0.5887150764465332
    35     8.149446e+02     1.163517e+00
 * time: 0.5898840427398682
    36     8.148452e+02     6.361037e-01
 * time: 0.5909380912780762
    37     8.148094e+02     1.206059e+00
 * time: 0.5920679569244385
    38     8.147652e+02     3.265117e-01
 * time: 0.5932090282440186
    39     8.147596e+02     2.561389e-01
 * time: 0.5949239730834961
    40     8.147580e+02     3.000458e-01
 * time: 0.5960021018981934
    41     8.147568e+02     2.611461e-01
 * time: 0.5970098972320557
    42     8.147558e+02     2.553203e-01
 * time: 0.598085880279541
    43     8.147554e+02     2.575697e-01
 * time: 0.5990650653839111
    44     8.147550e+02     2.591353e-01
 * time: 0.6000850200653076
    45     8.147544e+02     2.607006e-01
 * time: 0.6010849475860596
    46     8.147531e+02     2.629020e-01
 * time: 0.6021809577941895
    47     8.147507e+02     2.643156e-01
 * time: 0.6032769680023193
    48     8.147442e+02     2.665898e-01
 * time: 0.6043589115142822
    49     8.147253e+02     2.685605e-01
 * time: 0.6055519580841064
    50     8.146989e+02     6.027405e-01
 * time: 0.6132829189300537
    51     8.146523e+02     1.217843e+00
 * time: 0.6156408786773682
    52     8.144474e+02     3.016641e+00
 * time: 0.6168379783630371
    53     8.127958e+02     5.703653e+00
 * time: 0.6182301044464111
    54     8.121134e+02     8.840596e+00
 * time: 0.6201930046081543
    55     8.115071e+02     1.341296e+01
 * time: 0.6230208873748779
    56     8.113122e+02     9.704525e+00
 * time: 0.6256310939788818
    57     8.111868e+02     3.354920e+00
 * time: 0.6269519329071045
    58     8.111649e+02     3.315285e-01
 * time: 0.6283490657806396
    59     8.111621e+02     3.176063e-01
 * time: 0.629478931427002
    60     8.110946e+02     2.552949e+00
 * time: 0.6308720111846924
    61     8.109807e+02     5.009345e+00
 * time: 0.6323039531707764
    62     8.107321e+02     7.445230e+00
 * time: 0.633652925491333
    63     8.104088e+02     6.970537e+00
 * time: 0.6401770114898682
    64     8.101827e+02     3.136169e+00
 * time: 0.6415081024169922
    65     8.101494e+02     6.649904e-01
 * time: 0.6427919864654541
    66     8.101477e+02     2.623425e-01
 * time: 0.6440389156341553
    67     8.101468e+02     1.525673e-01
 * time: 0.6451449394226074
    68     8.101467e+02     1.528493e-01
 * time: 0.6461400985717773
    69     8.101467e+02     1.530786e-01
 * time: 0.6471230983734131
    70     8.101466e+02     1.534902e-01
 * time: 0.6482040882110596
    71     8.101463e+02     1.540950e-01
 * time: 0.6493480205535889
    72     8.101455e+02     1.584052e-01
 * time: 0.6505169868469238
    73     8.101435e+02     1.904824e-01
 * time: 0.6516819000244141
    74     8.101383e+02     3.452639e-01
 * time: 0.6529090404510498
    75     8.101244e+02     6.017616e-01
 * time: 0.6577548980712891
    76     8.100874e+02     1.034478e+00
 * time: 0.6590180397033691
    77     8.099885e+02     1.778811e+00
 * time: 0.6602818965911865
    78     8.098057e+02     3.325461e+00
 * time: 0.6617410182952881
    79     8.096503e+02     2.639750e+00
 * time: 0.6632349491119385
    80     8.093911e+02     1.669475e+00
 * time: 0.6647319793701172
    81     8.093681e+02     5.739834e-01
 * time: 0.6661860942840576
    82     8.093600e+02     9.433986e-01
 * time: 0.6675410270690918
    83     8.093472e+02     2.721046e-01
 * time: 0.6689310073852539
    84     8.093450e+02     2.378452e-01
 * time: 0.6701970100402832
    85     8.093449e+02     2.127704e-01
 * time: 0.6751348972320557
    86     8.093446e+02     2.220855e-01
 * time: 0.6762149333953857
    87     8.093438e+02     2.456031e-01
 * time: 0.6773738861083984
    88     8.093416e+02     3.418575e-01
 * time: 0.6786160469055176
    89     8.093361e+02     5.958788e-01
 * time: 0.679847002029419
    90     8.093216e+02     1.003145e+00
 * time: 0.6810839176177979
    91     8.092834e+02     1.664700e+00
 * time: 0.6824159622192383
    92     8.091831e+02     2.768165e+00
 * time: 0.6839759349822998
    93     8.089580e+02     4.820284e+00
 * time: 0.6855800151824951
    94     8.086509e+02     5.701941e+00
 * time: 0.6871850490570068
    95     8.080609e+02     3.807745e+00
 * time: 0.692781925201416
    96     8.075956e+02     9.191869e-01
 * time: 0.694633960723877
    97     8.075041e+02     8.324795e-01
 * time: 0.6961989402770996
    98     8.074673e+02     6.186203e-01
 * time: 0.6977629661560059
    99     8.074230e+02     5.781272e-01
 * time: 0.6996059417724609
   100     8.074084e+02     4.101525e-01
 * time: 0.7029290199279785
   101     8.074004e+02     1.810009e-01
 * time: 0.7061610221862793
   102     8.073931e+02     4.032068e-01
 * time: 0.7079830169677734
   103     8.073880e+02     4.295068e-01
 * time: 0.7096009254455566
   104     8.073503e+02     7.177194e-01
 * time: 0.7158739566802979
   105     8.072764e+02     1.078033e+00
 * time: 0.7176098823547363
   106     8.070835e+02     1.676936e+00
 * time: 0.7193310260772705
   107     8.067051e+02     2.970590e+00
 * time: 0.7209830284118652
   108     8.063589e+02     2.672784e+00
 * time: 0.7222399711608887
   109     8.057565e+02     3.656213e+00
 * time: 0.7238130569458008
   110     8.054303e+02     3.110965e+00
 * time: 0.7252740859985352
   111     8.054172e+02     2.375133e+00
 * time: 0.7269999980926514
   112     8.053551e+02     5.285517e-01
 * time: 0.7284688949584961
   113     8.053523e+02     1.364149e-01
 * time: 0.7299840450286865
   114     8.053515e+02     1.365453e-02
 * time: 0.7312290668487549
   115     8.053515e+02     8.137844e-03
 * time: 0.7357978820800781
   116     8.053515e+02     8.152911e-03
 * time: 0.7368988990783691
   117     8.053515e+02     8.156157e-03
 * time: 0.737847089767456
   118     8.053515e+02     8.156075e-03
 * time: 0.7390139102935791
   119     8.053515e+02     8.155839e-03
 * time: 0.7402489185333252
   120     8.053515e+02     8.155389e-03
 * time: 0.7413949966430664
   121     8.053515e+02     8.152666e-03
 * time: 0.7423350811004639
   122     8.053515e+02     8.152662e-03
 * time: 0.7434461116790771
   123     8.053515e+02     8.139791e-03
 * time: 0.7445979118347168
   124     8.053515e+02     1.135269e-02
 * time: 0.7457559108734131
   125     8.053515e+02     2.385565e-02
 * time: 0.7506210803985596
   126     8.053515e+02     3.979734e-02
 * time: 0.7517349720001221
   127     8.053514e+02     6.611758e-02
 * time: 0.7528579235076904
   128     8.053512e+02     1.019142e-01
 * time: 0.7539749145507812
   129     8.053509e+02     1.433596e-01
 * time: 0.7550520896911621
   130     8.053501e+02     1.684077e-01
 * time: 0.7561230659484863
   131     8.053490e+02     1.413035e-01
 * time: 0.7571990489959717
   132     8.053480e+02     5.888520e-02
 * time: 0.7584519386291504
   133     8.053477e+02     1.028002e-02
 * time: 0.7597129344940186
   134     8.053476e+02     2.437124e-02
 * time: 0.7608299255371094
   135     8.053476e+02     2.886175e-02
 * time: 0.7618279457092285
   136     8.053475e+02     1.936667e-02
 * time: 0.7663478851318359
   137     8.053475e+02     4.619635e-03
 * time: 0.7673399448394775
   138     8.053474e+02     3.401024e-03
 * time: 0.7683470249176025
   139     8.053474e+02     5.460919e-03
 * time: 0.7692348957061768
   140     8.053474e+02     4.049276e-03
 * time: 0.7701120376586914
   141     8.053474e+02     1.227562e-03
 * time: 0.7709879875183105
   142     8.053474e+02     5.661175e-04
 * time: 0.7719039916992188

FittedPumasModel

Dynamical system type:          No dynamical model

Number of subjects:                             64

Observation records:         Active        Missing
    dv:                         235              0
    Total:                      235              0

Number of parameters:      Constant      Optimized
                                  0             10

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -805.34742

-------------------
         Estimate
-------------------
tvE0       15.276
tvEC50    905.66
tvEmax     94.772
Ω₁,₁        9.6107
Ω₂,₁        4.3304
Ω₃,₁      106.66
Ω₂,₂        1.9512
Ω₃,₂       48.059
Ω₃,₃     1183.7
σ_add       6.6229
-------------------
ast_infer = infer(ast_fit)
ast_infer
[ Info: Calculating: variance-covariance matrix.
┌ Error: Failed.
└ @ Pumas ~/.julia/packages/Pumas/GZeMg/src/estimation/inference.jl:601

Asymptotic inference results

Dynamical system type:          No dynamical model

Number of subjects:                             64

Observation records:         Active        Missing
    dv:                         235              0
    Total:                      235              0

Number of parameters:      Constant      Optimized
                                  0             10

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -805.34742

---------------------------------------
         Estimate    SE    95.0% C.I.
---------------------------------------
tvE0       15.276    NaN   [ NaN; NaN]
tvEC50    905.66     NaN   [ NaN; NaN]
tvEmax     94.772    NaN   [ NaN; NaN]
Ω₁,₁        9.6107   NaN   [ NaN; NaN]
Ω₂,₁        4.3304   NaN   [ NaN; NaN]
Ω₃,₁      106.66     NaN   [ NaN; NaN]
Ω₂,₂        1.9512   NaN   [ NaN; NaN]
Ω₃,₂       48.059    NaN   [ NaN; NaN]
Ω₃,₃     1183.7      NaN   [ NaN; NaN]
σ_add       6.6229   NaN   [ NaN; NaN]
---------------------------------------


Variance-covariance matrix could not be
evaluated. The random effects may be over-
parameterized. Check the coefficients for
variance estimates near zero.
ast_insp = inspect(ast_fit)
ast_insp_df = DataFrame(ast_insp)
obs_rows = @rsubset(ast_insp_df, !ismissing(:dv))

fig = Figure(size = (800, 600))

ax1 = Axis(fig[1, 1], xlabel = "PRED", ylabel = "DV", title = "DV vs PRED")
scatter!(ax1, obs_rows.dv_pred, obs_rows.dv, color = :navy, markersize = 5)
ablines!(ax1, 0, 1, color = :black)

ax2 = Axis(fig[1, 2], xlabel = "IPRED", ylabel = "DV", title = "DV vs IPRED")
scatter!(ax2, obs_rows.dv_ipred, obs_rows.dv, color = :navy, markersize = 5)
ablines!(ax2, 0, 1, color = :black)

ax3 = Axis(fig[2, 1], xlabel = "PRED", ylabel = "WRES", title = "WRES vs PRED")
scatter!(ax3, obs_rows.dv_pred, obs_rows.dv_wres, color = :navy, markersize = 5)
hlines!(ax3, 0, color = :black)

ax4 = Axis(fig[2, 2], xlabel = "Time", ylabel = "WRES", title = "WRES vs TIME")
scatter!(ax4, obs_rows.time, obs_rows.dv_wres, color = :navy, markersize = 5)
hlines!(ax4, 0, color = :black)

fig
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating dose control parameters.
[ Info: Evaluating individual parameters.
[ Info: Done.
Figure 2: Problem 3: AST Emax GOF Diagnostics

3.4 Problem 4: Population PK (1-Compartment Oral)

Population PK analysis of a typical Phase 1 dataset. The data has known errors that must be fixed (as discovered in the R script).

pk_raw = CSV.read(joinpath(data_dir, "prob_2.csv"), DataFrame)
rename!(pk_raw,
    :ID => :id, :DV => :dv, :AMT => :amt,
    :II => :ii, :ADDL => :addl, :TIME => :time,
)

# Data cleaning (from R script: IndPopModels.Problem2.R)
# Fix 1: ID=13, AMT=20 should be 200
@rtransform!(pk_raw, :amt = (:id == 13 && :amt == 20) ? 200.0 : Float64(:amt))

# Fix 2: ID=22, II=6 should be 8
@rtransform!(pk_raw, :ii = (:id == 22 && :ii == 6) ? 8 : :ii)

# Fix 3: ID=2, dose TIME=24 should be 0
@rtransform!(pk_raw, :time = (:id == 2 && :amt > 0 && :time == 24) ? 0.0 : Float64(:time))

# Fix 4: ID=15, dose AMT=0 should be 200
@rtransform!(pk_raw, :amt = (:id == 15 && :time == 0 && :amt == 0) ? 200.0 : Float64(:amt))

# Drop comment column, add required columns
select!(pk_raw, Not(:C))
pk_raw.evid = ifelse.(pk_raw.amt .> 0, 1, 0)
pk_raw.cmt = ifelse.(pk_raw.amt .> 0, 1, 2)  # depot=1, central=2
pk_raw.dv = ifelse.(pk_raw.evid .== 1, missing, pk_raw.dv)

pop_pk = read_pumas(pk_raw)
println("Subjects: ", length(pop_pk))
Subjects: 25
# 1-cpt oral model: ADVAN2 TRANS2 equivalent
poppk_1cpt = @model begin
    @metadata begin
        desc = "1-Compartment Oral PopPK (ADVAN2 TRANS2)"
    end
    @param begin
        tvCL   ∈ RealDomain(lower = 0.0, init = 8.0)
        tvV    ∈ RealDomain(lower = 0.0, init = 50.0)
        tvKa   ∈ RealDomain(lower = 0.0, init = 0.45)
        Ω      ∈ PDiagDomain(3)
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.2)
    end
    @random begin
        η ~ MvNormal(Ω)
    end
    @pre begin
        CL = tvCL * exp(η[1])
        Vc = tvV  * exp(η[2])
        Ka = tvKa * exp(η[3])
    end
    @dynamics Depots1Central1
    @derived begin
        cp := @. Central / Vc
        dv ~ @. Normal(cp, abs(cp) * σ_prop)
    end
end
PumasModel
  Parameters: tvCL, tvV, tvKa, Ω, σ_prop
  Random effects: η
  Covariates:
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: dv
  Observed: dv
pk_params = (
    tvCL = 8.0, tvV = 50.0, tvKa = 0.45,
    Ω = Diagonal([0.04, 0.04, 0.04]),
    σ_prop = 0.2,
)

# FO estimation first (matching NONMEM METHOD=0)
pk_fo = fit(poppk_1cpt, pop_pk, pk_params, Pumas.FO())
pk_fo
[ 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     8.264980e+01     2.753078e+01
 * time: 2.2172927856445312e-5
     1     8.141061e+01     5.944208e+01
 * time: 1.3636481761932373
     2     8.077844e+01     5.122005e+01
 * time: 1.364319086074829
     3     7.993858e+01     2.362374e+01
 * time: 1.364886999130249
     4     7.771907e+01     8.524163e+00
 * time: 1.3654201030731201
     5     7.744921e+01     5.509531e+00
 * time: 1.3659591674804688
     6     7.713191e+01     3.205463e+00
 * time: 1.3664031028747559
     7     7.703286e+01     2.240696e+00
 * time: 1.366866111755371
     8     7.698240e+01     3.473049e-01
 * time: 1.3673031330108643
     9     7.697257e+01     4.504869e-01
 * time: 1.3677361011505127
    10     7.695368e+01     6.067292e-01
 * time: 1.3681612014770508
    11     7.694136e+01     4.894391e-01
 * time: 1.3686342239379883
    12     7.692262e+01     4.692972e-01
 * time: 1.3691060543060303
    13     7.689627e+01     8.668093e-01
 * time: 1.3695740699768066
    14     7.686013e+01     1.080698e+00
 * time: 1.3700392246246338
    15     7.683395e+01     6.307220e-01
 * time: 1.370530128479004
    16     7.682714e+01     9.958439e-02
 * time: 1.3709781169891357
    17     7.682624e+01     8.441247e-02
 * time: 1.371408224105835
    18     7.682581e+01     1.294959e-01
 * time: 1.3718361854553223
    19     7.682555e+01     8.027178e-02
 * time: 1.3722660541534424
    20     7.682551e+01     2.023900e-02
 * time: 1.372709035873413
    21     7.682550e+01     3.212026e-03
 * time: 1.373138189315796
    22     7.682550e+01     2.372556e-04
 * time: 1.373603105545044

FittedPumasModel

Dynamical system type:                 Closed form

Number of subjects:                             25

Observation records:         Active        Missing
    dv:                         125              0
    Total:                      125              0

Number of parameters:      Constant      Optimized
                                  0              7

Likelihood approximation:                       FO
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -76.825504

--------------------
         Estimate
--------------------
tvCL      8.0439
tvV      50.66
tvKa      0.49759
Ω₁,₁      0.058496
Ω₂,₂      0.057974
Ω₃,₃      0.0083115
σ_prop    0.15247
--------------------
# Re-fit with FOCEI for better estimates
pk_focei = fit(poppk_1cpt, pop_pk, pk_params, Pumas.FOCEI())
pk_focei
[ 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     7.571519e+01     3.057811e+01
 * time: 2.4080276489257812e-5
     1     7.272019e+01     2.670023e+01
 * time: 0.380601167678833
     2     7.044261e+01     1.658803e+01
 * time: 0.3841361999511719
     3     6.973885e+01     1.215980e+01
 * time: 0.3862471580505371
     4     6.956608e+01     5.783517e+00
 * time: 0.38887810707092285
     5     6.914330e+01     4.683046e+00
 * time: 0.39066410064697266
     6     6.892282e+01     1.076909e+00
 * time: 0.39240503311157227
     7     6.886545e+01     7.776024e-01
 * time: 0.3941681385040283
     8     6.879593e+01     7.129593e-01
 * time: 0.39599609375
     9     6.874807e+01     7.434706e-01
 * time: 0.3976931571960449
    10     6.863010e+01     1.021736e+00
 * time: 0.3994481563568115
    11     6.851866e+01     9.436458e-01
 * time: 0.4391801357269287
    12     6.841345e+01     9.920031e-01
 * time: 0.4410359859466553
    13     6.834219e+01     7.754581e-01
 * time: 0.4427511692047119
    14     6.827168e+01     7.149361e-01
 * time: 0.44454216957092285
    15     6.821152e+01     6.196741e-01
 * time: 0.44622302055358887
    16     6.818100e+01     2.081450e-01
 * time: 0.4478580951690674
    17     6.816906e+01     1.109124e-01
 * time: 0.4495580196380615
    18     6.816207e+01     1.665597e-01
 * time: 0.4512310028076172
    19     6.815634e+01     1.542026e-01
 * time: 0.45285820960998535
    20     6.815261e+01     7.146703e-02
 * time: 0.4544861316680908
    21     6.815096e+01     3.142187e-02
 * time: 0.4560542106628418
    22     6.815023e+01     3.787429e-02
 * time: 0.4574730396270752
    23     6.814975e+01     3.439474e-02
 * time: 0.45899009704589844
    24     6.814946e+01     1.448103e-02
 * time: 0.46054911613464355
    25     6.814933e+01     4.476179e-03
 * time: 0.46212220191955566
    26     6.814927e+01     8.170391e-03
 * time: 0.46361804008483887
    27     6.814923e+01     7.013814e-03
 * time: 0.4652111530303955
    28     6.814921e+01     2.612858e-03
 * time: 0.46680712699890137
    29     6.814920e+01     9.269186e-04
 * time: 0.4684431552886963

FittedPumasModel

Dynamical system type:                 Closed form

Number of subjects:                             25

Observation records:         Active        Missing
    dv:                         125              0
    Total:                      125              0

Number of parameters:      Constant      Optimized
                                  0              7

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -68.149202

--------------------
         Estimate
--------------------
tvCL      8.3055
tvV      45.422
tvKa      0.40781
Ω₁,₁      0.057279
Ω₂,₂      0.052748
Ω₃,₃      5.6956e-7
σ_prop    0.15245
--------------------
pk_infer = infer(pk_focei)
pk_infer
[ Info: Calculating: variance-covariance matrix.
┌ Error: Failed.
└ @ Pumas ~/.julia/packages/Pumas/GZeMg/src/estimation/inference.jl:601

Asymptotic inference results

Dynamical system type:                 Closed form

Number of subjects:                             25

Observation records:         Active        Missing
    dv:                         125              0
    Total:                      125              0

Number of parameters:      Constant      Optimized
                                  0              7

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -68.149202

----------------------------------------
         Estimate     SE    95.0% C.I.
----------------------------------------
tvCL      8.3055      NaN   [ NaN; NaN]
tvV      45.422       NaN   [ NaN; NaN]
tvKa      0.40781     NaN   [ NaN; NaN]
Ω₁,₁      0.057279    NaN   [ NaN; NaN]
Ω₂,₂      0.052748    NaN   [ NaN; NaN]
Ω₃,₃      5.6956e-7   NaN   [ NaN; NaN]
σ_prop    0.15245     NaN   [ NaN; NaN]
----------------------------------------


Variance-covariance matrix could not be
evaluated. The random effects may be over-
parameterized. Check the coefficients for
variance estimates near zero.
pk_insp = inspect(pk_focei)
pk_df = DataFrame(pk_insp)
pk_obs = @rsubset(pk_df, !ismissing(:dv))

fig = Figure(size = (800, 600))

ax1 = Axis(fig[1, 1], xlabel = "PRED", ylabel = "DV", title = "DV vs PRED")
scatter!(ax1, pk_obs.dv_pred, pk_obs.dv, color = :navy, markersize = 5, alpha = 0.6)
ablines!(ax1, 0, 1, color = :black)

ax2 = Axis(fig[1, 2], xlabel = "IPRED", ylabel = "DV", title = "DV vs IPRED")
scatter!(ax2, pk_obs.dv_ipred, pk_obs.dv, color = :navy, markersize = 5, alpha = 0.6)
ablines!(ax2, 0, 1, color = :black)

ax3 = Axis(fig[2, 1], xlabel = "PRED", ylabel = "WRES", title = "WRES vs PRED")
scatter!(ax3, pk_obs.dv_pred, pk_obs.dv_wres, color = :navy, markersize = 5, alpha = 0.6)
hlines!(ax3, 0, color = :black)

ax4 = Axis(fig[2, 2], xlabel = "Time (hr)", ylabel = "WRES", title = "WRES vs TIME")
scatter!(ax4, pk_obs.time, pk_obs.dv_wres, color = :navy, markersize = 5, alpha = 0.6)
hlines!(ax4, 0, color = :black)

fig
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating dose control parameters.
[ Info: Evaluating individual parameters.
[ Info: Done.
Figure 3: Problem 4: PopPK GOF Diagnostics 
fig = Figure(size = (900, 700))
for (i, subj_id) in enumerate(unique(pk_obs.id)[1:9])
    row, col = divrem(i - 1, 3) .+ (1, 1)
    ax = Axis(fig[row, col], xlabel = "Time", ylabel = "Conc",
        title = "ID $subj_id")
    sdf = @rsubset(pk_obs, :id == subj_id)
    scatter!(ax, sdf.time, sdf.dv, color = :blue, markersize = 6)
    lines!(ax, sdf.time, sdf.dv_ipred, color = :red, linewidth = 1.5)
    lines!(ax, sdf.time, sdf.dv_pred, color = :gray, linewidth = 1, linestyle = :dash)
end
fig
Figure 4: Problem 4: Individual Fits (first 9 subjects)
Comparing Pumas Objective Function to NONMEM

When validating a Pumas model against NONMEM results:

  • NONMEM OFV (with constant) = -2 * loglikelihood(fit) in Pumas — these should match exactly
  • NONMEM OFV (without constant) = -2 * loglikelihood(fit) - N * log(2π) where N = total observations
  • To compute the traditional OFV: ofv = -2 * loglikelihood(fit) - nobs(fit) * log(2π)

Common reasons for mismatch:

  1. Data differences — column mapping, missing value handling, or filtering
  2. Parameterization — different IIV or error model form
  3. Covariate interpolation — NONMEM defaults to NOCB; Pumas defaults to LOCF (covariates_direction = :left)
  4. Initial estimates — compare EBEs via empirical_bayes(fit) against NONMEM’s .phi file
Pumas Error Model Convenience Functions

Pumas provides named error distributions that simplify @derived blocks:

Error Model Manual Form Convenience Form
Additive dv ~ Normal(cp, σ_add) dv ~ Normal(cp, σ_add)
Proportional dv ~ Normal(cp, abs(cp) * σ_prop) dv ~ ProportionalNormal(cp, σ_prop)
Combined dv ~ Normal(cp, √(σ_add² + (cp*σ_prop)²)) dv ~ CombinedNormal(cp, σ_add, σ_prop)

ProportionalNormal and CombinedNormal handle the abs/sqrt internally.

Connecting Error Models to Bioanalytical Assay Characteristics

The choice of residual error model should be informed by the characteristics of the bioanalytical assay, not solely by statistical fit. Additive error implies constant measurement uncertainty across the concentration range, which is typical of assays operating near the lower limit of quantification (LLOQ) where absolute noise dominates. Proportional error implies uncertainty that scales with concentration, reflecting the behavior of most validated bioanalytical assays where the coefficient of variation (%CV) is approximately constant across the calibration range. Combined error captures both regimes simultaneously — additive noise dominates at low concentrations and proportional noise dominates at higher concentrations. In practice, the combined model is often the most realistic because assay precision degrades near the LLOQ (additive component) while maintaining a roughly constant %CV at mid-to-high concentrations (proportional component). Understanding the assay’s validated range and precision profile provides a principled basis for selecting the error model prior to estimation.

4 Study Guide Questions

  1. Explain how you might interpret the standard error of the estimates.
  2. List three things that you should look for in the Monitoring of Search output (optimizer convergence in Pumas).
  3. List three diagnostic plots that you should create for a population NLME model.
  4. How can you make sense of the off-diagonal elements of a covariance matrix?
  5. Provide an initial estimate for the first element of OMEGA, assuming a proportional inter-individual %CV equal to 30%. (Answer: 0.09, since for exponential IIV, omega = (CV/100)^2 = 0.3^2 = 0.09)
  6. Provide an initial estimate for SIGMA, assuming a constant residual error with standard deviation equal to 5. (Answer: In NONMEM sigma=25 (variance); in Pumas sigma_add=5 (SD))

5 Supplementary Material

Topic Tutorial Key Additions
Lecture: 1-Cpt IV Bolus PK W2 The physiology behind the monoexponential model coded here
Lecture: Volume of Distribution PK W4 Why V is not a real volume — tissue partitioning, protein binding
Lecture: Clearance & Hepatic CL PK W6 CL as volume cleared per unit time, extraction ratio, intrinsic clearance
Lecture: Intrinsic CL & Enzyme Kinetics PK W7a Vmax/KM foundation of Michaelis-Menten and Emax models
Lecture: Oral Absorption PK W8 First-order absorption, bioavailability, flip-flop kinetics
Lecture: Error Model Specification PKPD W7 Choosing additive vs proportional vs combined error from assay characteristics
Complete Pumas fitting workflow Fitting Domain types, @param/@random/@pre/@dynamics/@derived deep dive, inference and GOF
9 ways to specify dynamics Dynamical Systems Closed-form, explicit ODE, ModelingToolkit, Catalyst — with linearity detection comparison
NONMEM ↔︎ Pumas validation Comparing NONMEM & Pumas OFV matching, parameter extraction from .lst, EBE comparison, troubleshooting
Absorption model catalog Absorption Models First-order, zero-order, parallel, transit, Erlang, Gamma, Weibull — all with working code
VEM estimation VEM Introduction Variational EM as alternative to FOCE, options reference, comparison study
MCEM estimation MCEM Monte Carlo EM for complex likelihoods