The fit Subscript
Having simulated some data to which we can fit our model, the Fit Script is run in order to estimate the “best fit” fixed effects (and random effects) for the simulated dataset.
fit Inputs
As in the gen script, a METHOD_OPTIONS block contains a py_module
entry that defines the action to take (fit in this case).
A new FILE_PATHS section includes an input_data_file entry that points
to the data file (my_pkpd_model_gen.pyml_output/generated_data/cx_obs_params.csv)
that defines the study and the individuals taking part in it.
FILE_PATHS:
input_data_file: builtin_tut_example_gen.pyml_output/generated_data/cx_obs_params.csv
output_file_ext: ['svg', 'pdf']
Because the properties of the population and study are entirely
contained in the input data, a fit script need only consider what
fixed effects and random effects to estimate, so the parent tut script’s
FIT_EFFECTS that is copied to the fit script’s EFFECTS block
does not include covariates (c[X]) or times (t[X]):
EFFECTS:
POP: |
f[KA] ~ P1.0
f[CL] ~ P1.0
f[V1] ~ P20
f[Q] ~ P0.5
f[V2] ~ P100
f[KA_isv,CL_isv,V1_isv,Q_isv,V2_isv] ~ spd_matrix() [
[0.05],
[0.01, 0.05],
[0.01, 0.01, 0.05],
[0.01, 0.01, 0.01, 0.05],
[0.01, 0.01, 0.01, 0.01, 0.05],
]
f[PNOISE] ~ P0.1
ID: |
r[KA, CL, V1, Q, V2] ~ mnorm([0,0,0,0,0], f[KA_isv,CL_isv,V1_isv,Q_isv,V2_isv])
Furthermore, because fixed effects are unknown but allowed to vary they are given a distribution of values (usually to limit their range) and an initial value to begin the estimation process:
f[X] ~ <distribution> <initial value>
For example,
f[X] ~ unif(0.0, +inf) 1.0
defines a uniformly distributed variable that must take a value in the range
(0, +infinity) - i.e. a positive value - and that has an initial value of 1.0.
Because this range is a common one in PK/PD modelling,
PoPy provides a convenient shortcut, ~ P:
f[X] ~ P 1.0
where the ‘P’ is short for ‘positive’.
Random effects definitions are identical to those in the gen script.
In total there are 271 parameters to estimate:
six univariate parameters
one symmetric covariance matrix with 15 degrees of freedom
five random effects for each of the 50 individuals in the population
MODEL_PARAMS, DERIVATIVES and PREDICTIONS define the mathematical model as before, and the ODE_SOLVER section sets the parameters of the solver that integrates derivatives into amounts.
OUTPUT_SCRIPTS defines the subscripts we want to run after completing
the model fitting, which in this case are sim, msim and fitsum.
As with gen, sim simulates smooth curves at dense time points for the
model with estimated (rather than true) values of the parameters,
and fitsum creates a HTML report of the fitting process.
msim is a different subscript that, given population fixed effects, will
generate many new populations whose statistics can be used to evaluate the
goodness of fit via a Visual Predictive Check (VPC).
The other section - also new to the fit script - is the FIT_METHODS
section that defines the fitting method (e.g. FOCE, JOE or ND) and its
parameters (e.g. tolerances or number of iterations). These parameters control
how well and how quickly a fitting algorithm will converge on a final estimate
and will be discussed in greater depth later in this book.
FIT_METHODS: [ND: {max_n_main_iterations: 30}]
fit Outputs
Running this fit script (either manually or as part of this tut
tutorial), will attempt to minimize iteratively an error function
(known as the Objective Function, OBJV) over a number of steps such that
the model parameters generate a better prediction of the data at every step.
The Objective Function is defined in such a way that it is minimized when
the likelihood (the probability of the data given the model) is maximized,
and when OBJV stops decreasing the fit is said to have converged
(Table 3).
Iteration |
Time |
OBJV |
0.00 |
4224.091990462789 |
|
11.27 |
888.5810225001037 |
|
11.57 |
888.5810225001037 |
|
1.1 |
11.64 |
-624.5908516989929 |
1.2 |
14.19 |
-688.3203750485884 |
1.3 |
15.71 |
-783.4113895200995 |
1.4 |
17.30 |
-817.544655433727 |
1.5 |
18.82 |
-819.0278364926105 |
1.6 |
20.24 |
-820.6396717219868 |
1.7 |
21.97 |
-821.8064086681877 |
1.8 |
23.39 |
-822.4392081874006 |
1.9 |
25.01 |
-824.1441967766439 |
1.10 |
27.08 |
-825.5293291066278 |
1.11 |
28.70 |
-828.130755672234 |
1.12 |
31.17 |
-829.0629711533805 |
1.13 |
32.51 |
-831.9135180999128 |
1.14 |
34.16 |
-833.5065464878089 |
1.15 |
35.87 |
-835.658611770827 |
1.16 |
37.67 |
-838.15981287146 |
1.17 |
39.96 |
-841.6970964235949 |
1.18 |
41.69 |
-844.3049259627668 |
1.19 |
43.31 |
-848.650027481997 |
1.20 |
45.48 |
-850.9383284732306 |
1.21 |
47.61 |
-853.0830534238515 |
1.22 |
48.91 |
-854.1773655364336 |
1.23 |
50.52 |
-854.4969371123143 |
1.24 |
51.88 |
-855.8025865592084 |
1.25 |
53.21 |
-858.9790044800343 |
1.26 |
55.11 |
-860.7660976740185 |
1.27 |
56.43 |
-861.8112092897645 |
1.28 |
57.77 |
-862.283489303459 |
1.29 |
59.34 |
-862.8287881866131 |
1.30 |
60.89 |
-863.1496093288893 |
Because the data to which we are fitting with our model was created by simulating observations with the same model, we should expect a good model fit to the data. (Some differences in estimated parameter values should, however, be expected because the observations are noisy and imperfect [Sheiner1980].)
When OBJV has converged to a locally minimal value (a “local minimum”) then the corresponding values of the fixed effects and random effects are returned as the “best” estimates of those variables:
f[KA] = 0.2148
f[CL] = 1.7896
f[V1] = 55.4580
f[Q] = 1.0564
f[V2] = 486.5705
f[KA_isv,CL_isv,V1_isv,Q_isv,V2_isv] = [
[ 0.1393, -0.0553, 0.0414, -0.0256, 0.1217 ],
[ -0.0553, 0.0808, -0.0070, -0.0061, -0.1784 ],
[ 0.0414, -0.0070, 0.1098, -0.0827, 0.2139 ],
[ -0.0256, -0.0061, -0.0827, 0.3173, -0.3239 ],
[ 0.1217, -0.1784, 0.2139, -0.3239, 0.9443 ],
]
f[PNOISE] = 0.1415
After running fit followed by sim followed by grph
we can see a smooth predicted time series of drug concentration
for every individual, overlaid with the input observations
(Table 4)
|
|
|
where a dashed blue line represents the noiseless model predictions
given the initial f[X] parameters and sampled r[X] values for each individual,
a solid blue line represents the predicted curve for the final values, and
solid blue dots represent the simulated, noisy observations
(c[DV_CENTRAL] in synthetic_data.csv).
These graphs show that PoPy has adjusted the f[X] and r[X] parameters
such that the PK curves more closely match the input data
and therefore maximise the likelihood of the data being generated from this model.