using Pumas
using AlgebraOfGraphics
using CairoMakie
using StableRNGs
Statistical Models Without Differential Equations
1 Introduction
This tutorial represents a slight deviation from most other models presented and performed in Pumas that often includes pharmacokinetic modeling. Here, we will focus on simple Exposure-Response (ER) models. ER models can be models of efficacy, safety, toxicity, or any other response. The response can be continuous in nature, it can be binary, or even multinomial such as ordinal pain score models.
The main point of the tutorial is to introduce the situation where there is data that does not come with an inherent dynamical system component. Instead, we want to study the effect of an exposure measure such as plasma concentration (potentially as a proxy for biophase concentration) or AUC on some kind of response. The response could be almost any response you can think of from the total area of the body affected by eczema, a pain-score, pain relief indicator, body weight, presence of vomiting, or any other efficacy, safety, or toxicity measure. To keep things simple and focus on higher level ideas we go with a sigmoidal Emax model (with a Hill exponent) for some hypothetical response given some level of plasma concentrations.
# Data generating parameters
Emax = 100
C50 = 40
h = 2
# Function to evaluate the Emax model
hill_model(cp, emax, c50, hill) = emax * cp^hill / (c50^hill + cp^hill)hill_model (generic function with 1 method)The Hill function is well-known in pharmacodynamics when we have an exposure-response relationship that has an upper effect limit. The Hill parameter allows us to change where in the exposure domain the change in effect strongest. It includes two additional parameters: C50 and Emax. First, is the half maximal effective concentration C50 (sometimes EC50) parameter that is the exposure (here concentration, hence C) that leads to half of the response in excess of a baseline value. Here, the baseline is zero. Second, is the parameter Emax that is simply the maximum effect or response that can occur. If we sample with a technique that truly has additive error the observed maximal effect could in principle be above Emax.
Let us draw such a function with the data generating parameters.
Show plotting code
hill_plot =
data((exposure = 0:200, effect = hill_model.(0:200, Emax, C50, h))) *
mapping(:exposure, :effect) *
visual(Lines) +
data((exposure = [C50], effect = [0, hill_model(C50, Emax, C50, h)])) *
mapping(:exposure, :effect) *
visual(Lines; linestyle = :dot) +
data((exposure = [0, C50], effect = [hill_model(C50, Emax, C50, h)])) *
mapping(:exposure, :effect) *
visual(Lines; linestyle = :dot) +
data((exposure = [0, 200], effect = [Emax])) *
mapping(:exposure, :effect) *
visual(Lines; linestyle = :dot)
axis_spec = (;
axis = (
limits = (0, 200, 0, Emax + 10),
xticks = ([0, C50, 100, 200], ["0", "C50", "100", "200"]),
yticks = ([0, 25, Emax / 2, 75, Emax], ["0", "25", "Emax/2", "75", "Emax"]),
)
)
draw(hill_plot; axis_spec...)
Since we have no model to explain the exposure but it is rather just a measured quantity we will simply sample exposures for this example. To ensure proper behavior of the estimator, we will make sure to sample well above C50.
2 Data Without Any Events
One basic property of the data for this kind of analysis we need to consider is that there are no events in the provided data. By events we think of the usual dose events associated with infusion rates and durations, bolus dose amounts, and other information about the specific mechanisms at play. We normally parse this information and use it to explicitly build a model for the trajectories of pharmacokinetic variables and how they might influence pharmacodynamics - but not here! We still require a time column, but often it may only be used to separate individual observations within subjects instead of being used for dynamical modeling.
Instead of using a predefined dataset we will construct a simple one as described in the previous section using sampled concentrations, the Hill model, and an additive error model. We need at least id, time, and observations to define an eventless dataset, but to drive the Emax model we need to include cp_i that are the measured or predicted exposures. The observations will be called resp here for response.
# Define the number of concentrations to sample
N = 40
# Define the random number generator
rng = StableRNG(983)
# Sample concentrations from a log-normal distribution
cp_i = rand(rng, LogNormal(log(C50 + 5), 0.6), N)
# Generate response variables given the exposure, cp_i and parameters for the Hill model
resp = @. rand(rng, Normal(hill_model(cp_i, Emax, C50, h), 3.2))
# Combine results into a DataFrame
response_df = DataFrame(id = 1:N, time = 1, cp_i = cp_i, resp = resp)40×4 DataFrame
15 rows omitted
| Row | id | time | cp_i | resp |
|---|---|---|---|---|
| Int64 | Int64 | Float64 | Float64 | |
| 1 | 1 | 1 | 29.9617 | 32.829 |
| 2 | 2 | 1 | 67.714 | 78.3184 |
| 3 | 3 | 1 | 58.7976 | 69.3724 |
| 4 | 4 | 1 | 37.54 | 46.1287 |
| 5 | 5 | 1 | 58.5205 | 67.6474 |
| 6 | 6 | 1 | 51.4522 | 61.8933 |
| 7 | 7 | 1 | 106.938 | 84.1643 |
| 8 | 8 | 1 | 17.5912 | 21.9294 |
| 9 | 9 | 1 | 14.974 | 11.6478 |
| 10 | 10 | 1 | 40.0817 | 47.5374 |
| 11 | 11 | 1 | 67.7002 | 77.037 |
| 12 | 12 | 1 | 38.6564 | 43.8903 |
| 13 | 13 | 1 | 31.8099 | 40.5693 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 29 | 29 | 1 | 109.159 | 82.9605 |
| 30 | 30 | 1 | 21.5369 | 22.6964 |
| 31 | 31 | 1 | 79.237 | 82.8221 |
| 32 | 32 | 1 | 94.6278 | 82.5703 |
| 33 | 33 | 1 | 53.9975 | 59.4201 |
| 34 | 34 | 1 | 44.5625 | 52.5194 |
| 35 | 35 | 1 | 130.494 | 92.3875 |
| 36 | 36 | 1 | 66.5144 | 75.9562 |
| 37 | 37 | 1 | 18.7343 | 22.4305 |
| 38 | 38 | 1 | 62.8681 | 71.8646 |
| 39 | 39 | 1 | 63.7703 | 73.3971 |
| 40 | 40 | 1 | 48.1412 | 64.6768 |
Show plotting code
draw(hill_plot + data(response_df) * mapping(:cp_i, :resp) * visual(Scatter); axis_spec...)
2.1 Defining Pumas Population Without Events
To map from tabular data in response_df to a Population in response_pop we use read_pumas just as we did in the case with event data. The important part is to turn off event_data to disable checks that are not relevant to this eventless example. If event_data is not set to false we would get errors about missing event columns for example.
response_pop = read_pumas(
response_df,
id = :id,
time = :time,
covariates = [:cp_i],
observations = [:resp],
event_data = false,
)Population
Subjects: 40
Covariates: cp_i
Observations: resp2.2 A Model Without Dynamics
Above we simply defined the input exposure and output response as simple arrays and function evaluations to defer the PumasModel definition until after seeing how to load eventless datasets. However, to perform a maximum likelihood fit we need to define a proper PumasModel.
response_model = @model begin
@param begin
θemax ∈ RealDomain(lower = 0, init = 90)
θc50 ∈ RealDomain(lower = 0, init = 30)
θhill ∈ RealDomain(lower = 0, init = 3)
σ ∈ RealDomain(lower = 1e-5, init = 0.1)
end
@covariates cp_i
@pre begin
emax_i = hill_model(cp_i, θemax, θc50, θhill)
end
@derived begin
resp ~ @. Normal(emax_i, σ)
end
endPumasModel
Parameters: θemax, θc50, θhill, σ
Random effects:
Covariates: cp_i
Dynamical system variables:
Dynamical system type: No dynamical model
Derived: resp
Observed: resp2.3 Fitting
To fit the model, we simply invoke fit with the model, population, parameters, and likelihood approximation method. Since there are no random effects in this model there is no integral to approximate so we use NaivePooled(). This will perform a maximum likelihood estimation according to the distribution used in @derived. The parameters in this case will match a non-linear least squares fit since we have an additive gaussian error model, but inference may deviate as usual.
emax_fit = fit(response_model, response_pop, init_params(response_model), NaivePooled())[ 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.191392e+05 2.289229e+06 * time: 0.04492497444152832 1 8.932842e+04 8.822257e+05 * time: 2.4604389667510986 2 4.110550e+04 1.497057e+05 * time: 2.4607698917388916 3 2.803549e+04 7.704834e+04 * time: 2.4609320163726807 4 1.746711e+04 7.957100e+04 * time: 2.461082935333252 5 1.107308e+04 7.925575e+04 * time: 2.4612419605255127 6 7.560761e+03 6.331888e+04 * time: 2.4613919258117676 7 6.125273e+03 4.654326e+04 * time: 2.4615519046783447 8 5.512225e+03 3.127545e+04 * time: 2.461697816848755 9 5.111728e+03 1.651115e+04 * time: 2.461840867996216 10 4.701021e+03 1.206406e+04 * time: 2.461987018585205 11 4.161755e+03 1.326976e+04 * time: 2.462139844894409 12 3.261693e+03 1.201441e+04 * time: 2.46228289604187 13 2.356246e+03 1.385234e+04 * time: 2.4624269008636475 14 1.928814e+03 1.207445e+04 * time: 2.462568998336792 15 1.058086e+03 6.260232e+03 * time: 2.4627180099487305 16 6.221586e+02 2.995796e+03 * time: 2.462877035140991 17 3.669135e+02 1.316280e+03 * time: 2.4630229473114014 18 2.436398e+02 6.822477e+02 * time: 2.4631829261779785 19 1.819199e+02 3.472574e+02 * time: 2.4633309841156006 20 1.540523e+02 1.827051e+02 * time: 2.4634718894958496 21 1.433511e+02 1.040047e+02 * time: 2.463618040084839 22 1.402804e+02 6.746720e+01 * time: 2.4637598991394043 23 1.396952e+02 5.123420e+01 * time: 2.46390700340271 24 1.395818e+02 4.380188e+01 * time: 2.4640488624572754 25 1.394928e+02 4.050505e+01 * time: 2.4642138481140137 26 1.393822e+02 5.505660e+01 * time: 2.464357852935791 27 1.393262e+02 6.733245e+01 * time: 2.4645028114318848 28 1.393130e+02 7.071214e+01 * time: 2.4646449089050293 29 1.393111e+02 7.015237e+01 * time: 2.4647929668426514 30 1.393106e+02 6.972962e+01 * time: 2.4649360179901123 31 1.393087e+02 6.863999e+01 * time: 2.4650819301605225 32 1.393046e+02 6.711256e+01 * time: 2.465235948562622 33 1.392930e+02 6.443478e+01 * time: 2.4653828144073486 34 1.392631e+02 6.009199e+01 * time: 2.4655399322509766 35 1.391837e+02 5.266633e+01 * time: 2.465693950653076 36 1.389749e+02 4.852068e+01 * time: 2.4658420085906982 37 1.384224e+02 6.914404e+01 * time: 2.465986967086792 38 1.369787e+02 1.095747e+02 * time: 2.4661450386047363 39 1.337496e+02 1.819193e+02 * time: 2.466295003890991 40 1.292445e+02 1.453304e+02 * time: 2.466439962387085 41 1.266898e+02 3.437051e+01 * time: 2.466596841812134 42 1.265079e+02 1.411626e+01 * time: 2.4667420387268066 43 1.264825e+02 5.597764e+00 * time: 2.466887950897217 44 1.264816e+02 5.740614e+00 * time: 2.467033863067627 45 1.264815e+02 5.793538e+00 * time: 2.4671878814697266 46 1.264814e+02 5.581702e+00 * time: 2.467341899871826 47 1.264814e+02 5.530522e+00 * time: 2.4674909114837646 48 1.264813e+02 5.275939e+00 * time: 2.4676380157470703 49 1.264810e+02 4.970802e+00 * time: 2.4677820205688477 50 1.264803e+02 4.413515e+00 * time: 2.467923879623413 51 1.264785e+02 3.555867e+00 * time: 2.4680678844451904 52 1.264737e+02 4.535225e+00 * time: 2.468219041824341 53 1.264612e+02 8.589969e+00 * time: 2.468362808227539 54 1.264287e+02 1.520662e+01 * time: 2.468505859375 55 1.263464e+02 2.630965e+01 * time: 2.468648910522461 56 1.261502e+02 4.613011e+01 * time: 2.4687929153442383 57 1.257527e+02 7.955082e+01 * time: 2.4689409732818604 58 1.249560e+02 1.123712e+02 * time: 2.469085931777954 59 1.235645e+02 9.395865e+01 * time: 2.4692349433898926 60 1.223284e+02 5.036318e+01 * time: 2.469374895095825 61 1.200148e+02 2.219219e+01 * time: 2.46952486038208 62 1.184537e+02 4.474061e+01 * time: 2.469728946685791 63 1.173665e+02 6.035637e+01 * time: 2.4701650142669678 64 1.161599e+02 7.012795e+01 * time: 2.4703688621520996 65 1.144975e+02 1.939249e+02 * time: 2.470515012741089 66 1.134709e+02 4.912669e+01 * time: 2.4706578254699707 67 1.117234e+02 5.079030e+01 * time: 2.470798969268799 68 1.082554e+02 8.642695e+01 * time: 2.4709978103637695 69 1.064381e+02 1.023070e+02 * time: 2.4712278842926025 70 1.051470e+02 1.054707e+02 * time: 2.4714338779449463 71 1.046726e+02 9.971096e+01 * time: 2.471637010574341 72 1.032863e+02 8.590532e+01 * time: 2.4717929363250732 73 1.028726e+02 1.319924e+02 * time: 2.471944808959961 74 1.025334e+02 8.632347e+01 * time: 2.472092866897583 75 1.019406e+02 1.075698e+01 * time: 2.47225284576416 76 1.019208e+02 1.331307e+00 * time: 2.472398042678833 77 1.019191e+02 3.467414e-01 * time: 2.4725449085235596 78 1.019190e+02 6.481356e-02 * time: 2.4726908206939697 79 1.019189e+02 5.769418e-03 * time: 2.472853899002075 80 1.019189e+02 1.779522e-04 * time: 2.473001003265381
FittedPumasModel
Dynamical system type: No dynamical model
Number of subjects: 40
Observation records: Active Missing
resp: 40 0
Total: 40 0
Number of parameters: Constant Optimized
0 4
Likelihood approximation: NaivePooled
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -101.91895
-----------------
Estimate
-----------------
θemax 104.13
θc50 41.558
θhill 1.8147
σ 3.0927
-----------------and we may use the usual workflow to get estimates of parameter uncertainty
infer(emax_fit)[ Info: Calculating: variance-covariance matrix. [ Info: Done.
Asymptotic inference results using sandwich estimator
Dynamical system type: No dynamical model
Number of subjects: 40
Observation records: Active Missing
resp: 40 0
Total: 40 0
Number of parameters: Constant Optimized
0 4
Likelihood approximation: NaivePooled
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -101.91895
--------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------
θemax 104.13 4.1798 [ 95.942 ; 112.33 ]
θc50 41.558 2.1123 [ 37.418 ; 45.698 ]
θhill 1.8147 0.1155 [ 1.5883; 2.0411]
σ 3.0927 0.27055 [ 2.5624; 3.6229]
--------------------------------------------------as well as inspect. Diagnostics for these models are typically simple enough that they can be constructed from the DataFrame constructed from inspect output.
2.4 Extensions
Since we used a normal PumasModel we can extend the response analysis with:
- covariate effects including time, dose level, etc
- random effects if there are multiple observations per subject
- more complicated response models such as binary response and ordinal response
3 Concluding Remarks
This tutorial introduced a simple Exposure-Response (ER) modeling approach using the Emax model, focusing on cases where no explicit time component or dynamic system is involved. By working with eventless data, we demonstrated how to structure the dataset, define a model without pharmacokinetics, and fit it using maximum likelihood estimation in Pumas.
The methodology presented provides a foundation for analyzing ER relationships in various contexts, from efficacy assessments to safety evaluations. While we used a basic model, the approach can be extended to incorporate covariate effects, random effects, and more complex response types, making it highly adaptable to different modeling needs. By leveraging these capabilities, users can perform robust ER analyses even in the absence of traditional pharmacokinetic modeling constraints.