Fitting models in Pumas
2021-08-20
using Random using Pumas using PumasUtilities using CairoMakie
Important Note: Pumas.jl is a proprietary package. It is free to use for no n-commercial academic teaching and research purposes. For commercial users, license fees apply. Further, Pumas may not be used under any license agreement for patient mana gement services in a hospital or clinical settings. Please refer to End User Licen se Agreement (https://juliacomputing.com/eula) for details. Please contact sales@pumas.a i for purchase.
Fitting a PK model
In this tutorial we will go through the steps required to fit a model in Pumas.jl.
The dosage regimen
We start by simulating a population from a two compartment model with oral absorption, and then we show how to fit and do some model validation using the fit output.
The dosage regimen is an dose of 100 into Depot at time 0, followed by two additional (addl=2) doses every fourth hour
repeated_dose_regimen = DosageRegimen(100, time = 0, ii = 4, addl = 2)
| time | cmt | amt | evid | ii | addl | rate | duration | ss | route | |
|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | Route… | |
| 1 | 0.0 | 1 | 100.0 | 1 | 4.0 | 2 | 0.0 | 0.0 | 0 | NullRoute |
As usual, let's define a function to choose body weight randomly per subject
choose_covariates() = (Wt = rand(55:80),)
choose_covariates (generic function with 1 method)
and generate a population of subjects with a random weight generated from the covariate function above
pop = map(i -> Subject(id = i, events = repeated_dose_regimen, observations = (dv = nothing,), covariates = choose_covariates()), 1:24)
Population Subjects: 24 Covariates: Wt Observations: dv
We now have 24 subjects equipped with a weight and a dosage regimen.
The PK model of drug concentration and elimination
To simulate a data set and attempt to estimate the data generating parameters, we have to set up the actual pharmacokinetics (PK) model and simulate the data. We use the closed form model called Depots1Central1Periph1 which is a two compartment model with first order absorption. This requires CL, Vc, Ka, Vp, and Q to be defined in the @pre-block, since they define the rates of transfer between (and out of) the compartments
mymodel = @model begin @param begin cl ∈ RealDomain(lower = 0.0, init = 1.0) tv ∈ RealDomain(lower = 0.0, init = 10.0) ka ∈ RealDomain(lower = 0.0, init = 1.0) q ∈ RealDomain(lower = 0.0, init = 0.5) Ω ∈ PDiagDomain(init = [0.9,0.07, 0.05]) σ_prop ∈ RealDomain(lower = 0,init = 0.03) end @random begin η ~ MvNormal(Ω) end @covariates Wt @pre begin CL = cl * (Wt/70)^0.75 * exp(η[1]) Vc = tv * (Wt/70) * exp(η[2]) Ka = ka * exp(η[3]) Vp = 30.0 Q = q end @dynamics Depots1Central1Periph1 @derived begin cp := @. 1000*(Central / Vc) # We use := because we don't want simobs to store the variable dv ~ @. Normal(cp, abs(cp)*σ_prop) end end
PumasModel Parameters: cl, tv, ka, q, Ω, σ_prop Random effects: η Covariates: Wt Dynamical variables: Depot, Central, Peripheral Derived: dv Observed: dv
Some parameters are left free by giving them domains in the @param-block, and one PK parameter (the volume of distribution of the peripheral compartment) is fixed to 20.0.
Simulating the individual observations
The simobs function is used to simulate individual time series. We input the model, the population of Subjects that currently only have dosage regimens and covariates, the parameter vector and the times where we want to simulate. Since we have a proportional error model we avoid observations at time zero to avoid degenerate distributions of the dependent variable. The problem is, that if the concentration is zero the variance in distribution of the explained variable will also be zero. Let's use the default parameters, as set in the @param-block, and simulate the data
param = init_params(mymodel) Random.seed!(1234) sims = simobs(mymodel, pop, param, obstimes = 1:1:72)
sim_plot(mymodel, sims, observations =[:dv], figure = (fontsize = 18, ), axis = (xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"))
Fitting the model
To fit the model, we use the fit function. It requires a model, a population, a named tuple of parameters and a likelihood approximation method.
result = fit(mymodel, Subject.(sims), param, Pumas.FOCE())
Iter Function value Gradient norm
0 1.006124e+04 1.189287e+02
* time: 0.01732802391052246
1 1.006122e+04 1.008886e+02
* time: 1.0786399841308594
2 1.006097e+04 8.354758e+01
* time: 1.0924510955810547
3 1.006047e+04 4.952202e+01
* time: 1.1041531562805176
4 1.005951e+04 4.278742e+01
* time: 1.1143040657043457
5 1.005935e+04 1.730721e+02
* time: 1.1244909763336182
6 1.005847e+04 3.982153e+00
* time: 1.1340019702911377
7 1.005840e+04 1.824897e+00
* time: 1.1427371501922607
8 1.005826e+04 3.248011e+00
* time: 1.1520881652832031
9 1.005816e+04 3.388170e+00
* time: 1.161694049835205
10 1.005814e+04 7.354519e-01
* time: 1.1709630489349365
11 1.005814e+04 6.604098e-02
* time: 1.1806020736694336
12 1.005814e+04 6.279394e-02
* time: 1.1884961128234863
13 1.005814e+04 1.040618e-01
* time: 1.1970741748809814
14 1.005814e+04 7.067727e-02
* time: 1.205679178237915
15 1.005814e+04 1.983862e-02
* time: 1.2150611877441406
16 1.005814e+04 1.886800e-03
* time: 1.223703145980835
17 1.005814e+04 1.835534e-03
* time: 1.2315449714660645
18 1.005814e+04 3.418999e-03
* time: 1.2388830184936523
19 1.005814e+04 3.538633e-03
* time: 1.247145175933838
20 1.005814e+04 2.945302e-03
* time: 1.2556490898132324
21 1.005814e+04 2.440051e-04
* time: 1.2662761211395264
FittedPumasModel
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -10058.137
Number of subjects: 24
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
---------------------
Estimate
---------------------
cl 0.70364
tv 9.7689
ka 1.0287
q 0.50024
Ω₁,₁ 0.76205
Ω₂,₂ 0.048448
Ω₃,₃ 0.058201
σ_prop 0.030283
---------------------
Of course, we started the fitting at the true parameters, so let us define our own starting parameters, and fit based on those values
alternative_param = ( cl = 0.5, tv = 9.0, ka = 1.3, q = 0.3, Ω = Diagonal([0.18,0.04, 0.03]), σ_prop = 0.04) fit(mymodel, Subject.(sims), alternative_param, Pumas.FOCE())
Iter Function value Gradient norm
0 2.242244e+04 4.079826e+04
* time: 4.291534423828125e-5
1 1.596822e+04 1.997004e+04
* time: 0.032621145248413086
2 1.163331e+04 5.687883e+03
* time: 0.046514034271240234
3 1.111722e+04 2.014736e+03
* time: 0.05762505531311035
4 1.096977e+04 1.355036e+03
* time: 0.06728196144104004
5 1.091489e+04 1.328994e+03
* time: 0.07570409774780273
6 1.055951e+04 1.150869e+03
* time: 0.09117412567138672
7 1.034606e+04 9.558073e+02
* time: 0.10634207725524902
8 1.011361e+04 3.090090e+02
* time: 0.11777806282043457
9 1.010353e+04 8.935514e+02
* time: 0.12695908546447754
10 1.009992e+04 1.493625e+03
* time: 0.13565301895141602
11 1.009585e+04 3.463120e+02
* time: 0.14665603637695312
12 1.009432e+04 1.004895e+02
* time: 0.15567708015441895
13 1.009390e+04 9.835821e+01
* time: 0.16340994834899902
14 1.008198e+04 2.144859e+02
* time: 0.17381811141967773
15 1.008164e+04 5.575778e+01
* time: 0.18186211585998535
16 1.008153e+04 3.610147e+01
* time: 0.1905660629272461
17 1.008128e+04 1.414547e+02
* time: 0.19892311096191406
18 1.008065e+04 3.197100e+02
* time: 0.20764899253845215
19 1.007921e+04 5.537094e+02
* time: 0.2161259651184082
20 1.007632e+04 7.768123e+02
* time: 0.22536110877990723
21 1.007179e+04 8.145490e+02
* time: 0.2353370189666748
22 1.006734e+04 5.121966e+02
* time: 0.24671411514282227
23 1.006582e+04 1.261857e+02
* time: 0.2562999725341797
24 1.006571e+04 2.431418e+01
* time: 0.26478099822998047
25 1.006567e+04 2.665189e+01
* time: 0.2738659381866455
26 1.006551e+04 1.189259e+02
* time: 0.28220200538635254
27 1.006528e+04 1.934590e+02
* time: 0.29140400886535645
28 1.006485e+04 2.496918e+02
* time: 0.30128908157348633
29 1.006432e+04 2.247446e+02
* time: 0.30983400344848633
30 1.006387e+04 1.122881e+02
* time: 0.3190159797668457
31 1.006367e+04 2.218938e+01
* time: 0.6476240158081055
32 1.006359e+04 2.740558e+01
* time: 0.6557700634002686
33 1.006355e+04 3.344587e+01
* time: 0.6637070178985596
34 1.006349e+04 2.646741e+01
* time: 0.6717629432678223
35 1.006335e+04 1.748026e+01
* time: 0.6793560981750488
36 1.006305e+04 2.838470e+01
* time: 0.6878321170806885
37 1.006246e+04 7.364463e+01
* time: 0.6947779655456543
38 1.006146e+04 1.095040e+02
* time: 0.7036211490631104
39 1.006005e+04 1.108130e+02
* time: 0.7124779224395752
40 1.005865e+04 6.833121e+01
* time: 0.7221460342407227
41 1.005835e+04 2.150404e+01
* time: 0.7313971519470215
42 1.005833e+04 7.732434e+00
* time: 0.740915060043335
43 1.005833e+04 1.931051e+00
* time: 0.7491390705108643
44 1.005833e+04 1.938613e+00
* time: 0.7581050395965576
45 1.005833e+04 1.951543e+00
* time: 0.7669370174407959
46 1.005833e+04 1.958896e+00
* time: 0.7756450176239014
47 1.005833e+04 1.969511e+00
* time: 0.7832961082458496
48 1.005832e+04 1.968297e+00
* time: 0.7923779487609863
49 1.005831e+04 2.188541e+00
* time: 0.8019530773162842
50 1.005829e+04 3.269215e+00
* time: 0.8105249404907227
51 1.005825e+04 4.154125e+00
* time: 0.8191149234771729
52 1.005819e+04 3.788114e+00
* time: 0.8278179168701172
53 1.005815e+04 1.935035e+00
* time: 0.8364870548248291
54 1.005814e+04 4.360267e-01
* time: 0.846066951751709
55 1.005814e+04 2.151379e-01
* time: 0.855355978012085
56 1.005814e+04 5.672937e-02
* time: 0.863037109375
57 1.005814e+04 2.931482e-02
* time: 0.8716959953308105
58 1.005814e+04 2.931269e-02
* time: 0.8821239471435547
59 1.005814e+04 5.189560e-02
* time: 0.8902690410614014
60 1.005814e+04 9.916391e-02
* time: 0.8977150917053223
61 1.005814e+04 2.047834e-01
* time: 0.9066629409790039
62 1.005814e+04 3.428638e-01
* time: 0.9154460430145264
63 1.005814e+04 5.335779e-01
* time: 0.923875093460083
64 1.005814e+04 6.985796e-01
* time: 0.932560920715332
65 1.005814e+04 6.637835e-01
* time: 0.9416670799255371
66 1.005814e+04 3.519294e-01
* time: 0.9496650695800781
67 1.005814e+04 7.386846e-02
* time: 0.9593639373779297
68 1.005814e+04 1.169473e-03
* time: 0.9673380851745605
69 1.005814e+04 1.103316e-03
* time: 0.9751200675964355
70 1.005814e+04 1.892715e-04
* time: 0.982928991317749
FittedPumasModel
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -10058.137
Number of subjects: 24
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
---------------------
Estimate
---------------------
cl 0.70364
tv 9.7689
ka 1.0287
q 0.50024
Ω₁,₁ 0.76205
Ω₂,₂ 0.048448
Ω₃,₃ 0.058201
σ_prop 0.030283
---------------------
and we see that the estimates are essentially the same up to numerical noise.
To augment the basic information listed when we print the results, we can use infer to provide RSEs and confidence intervals
infer(result)
Asymptotic inference results using sandwich estimator
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -10058.137
Number of subjects: 24
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
--------------------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------------------
cl 0.70364 0.12548 [0.45771 ; 0.94957 ]
tv 9.7689 0.43957 [8.9074 ; 10.63 ]
ka 1.0287 0.052072 [0.92667 ; 1.1308 ]
q 0.50024 0.0010611 [0.49816 ; 0.50232 ]
Ω₁,₁ 0.76205 0.29007 [0.19353 ; 1.3306 ]
Ω₂,₂ 0.048448 0.013318 [0.022346; 0.07455 ]
Ω₃,₃ 0.058201 0.016218 [0.026415; 0.089987]
σ_prop 0.030283 0.00051879 [0.029266; 0.0313 ]
--------------------------------------------------------------------
So as we observed earlier, the parameters look like they have sensible values. The confidence intervals are a bit wide, and especially so for the random effect variability parameters. To see how we can use simulation to better understand the statistical properties of our model, we can simulate a much larger population and try again
pop_big = map(i -> Subject(id = i, events = repeated_dose_regimen, observations = (dv = nothing,), covariates = choose_covariates()), 1:100) sims_big = simobs(mymodel, pop_big, param, obstimes=1:1:72) result_big = fit(mymodel, Subject.(sims_big), param, Pumas.FOCE()) infer(result_big)
Iter Function value Gradient norm
0 3.679364e+04 7.630791e+02
* time: 0.0002110004425048828
1 3.679332e+04 1.483985e+02
* time: 0.23825597763061523
2 3.679308e+04 1.287341e+02
* time: 0.2739448547363281
3 3.679287e+04 5.245973e+01
* time: 0.30597782135009766
4 3.679173e+04 7.537700e+01
* time: 0.3377988338470459
5 3.679106e+04 4.713081e+01
* time: 0.37213993072509766
6 3.679072e+04 1.131401e+01
* time: 0.4063758850097656
7 3.679068e+04 4.231405e+00
* time: 0.4401698112487793
8 3.679067e+04 5.464074e-01
* time: 0.47487378120422363
9 3.679067e+04 5.086507e-01
* time: 0.5078549385070801
10 3.679066e+04 7.688545e-01
* time: 0.9162437915802002
11 3.679066e+04 5.785163e-01
* time: 0.9450597763061523
12 3.679066e+04 1.437409e-01
* time: 0.9723308086395264
13 3.679066e+04 1.577841e-02
* time: 0.9983208179473877
14 3.679066e+04 1.577773e-02
* time: 1.0528099536895752
15 3.679066e+04 1.577773e-02
* time: 1.1543779373168945
16 3.679066e+04 1.577772e-02
* time: 1.2510960102081299
17 3.679066e+04 1.577772e-02
* time: 1.6714198589324951
18 3.679066e+04 1.577772e-02
* time: 1.7575418949127197
Asymptotic inference results using sandwich estimator
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -36790.662
Number of subjects: 100
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
dv: 7200 0
Total: 7200 0
--------------------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------------------
cl 1.2148 0.11819 [0.9832 ; 1.4465 ]
tv 10.041 0.27034 [9.5116 ; 10.571 ]
ka 1.0159 0.023039 [0.97072 ; 1.061 ]
q 0.5004 0.00051675 [0.49939 ; 0.50141 ]
Ω₁,₁ 0.94738 0.13271 [0.68727 ; 1.2075 ]
Ω₂,₂ 0.071937 0.010667 [0.051031; 0.092843]
Ω₃,₃ 0.049599 0.0065982 [0.036667; 0.062531]
σ_prop 0.029833 0.00024886 [0.029345; 0.030321]
--------------------------------------------------------------------
This time we see similar estimates, but much narrower confidence intervals across the board.
Estimating a mis-specified model
To explore some of the diagnostics tools available in Pumas, we can try to set up a model that does not fit out data generating process. This time we propose a one compartent model. The problem with estimating a one compartment model when the data comes from a two compartment model, is that we cannot capture the change in slope on the concentration profile you get with a two compartment model. This means that even if we can capture the model fit someone well on average, we should expect to see systematic trends in the residual diagnostics post estimation.
mymodel_misspec = @model begin @param begin cl ∈ RealDomain(lower = 0.0, init = 1.0) tv ∈ RealDomain(lower = 0.0, init = 20.0) ka ∈ RealDomain(lower = 0.0, init = 1.0) Ω ∈ PDiagDomain(init = [0.12, 0.05, 0.08]) σ_prop ∈ RealDomain(lower = 0, init = 0.03) end @random begin η ~ MvNormal(Ω) end @pre begin CL = cl * (Wt/70)^0.75 * exp(η[1]) Vc = tv * (Wt/70) * exp(η[2]) Ka = ka * exp(η[3]) end @covariates Wt @dynamics Depots1Central1 @derived begin cp := @. 1000*(Central / Vc) dv ~ @. Normal(cp, abs(cp)*σ_prop) end end
PumasModel Parameters: cl, tv, ka, Ω, σ_prop Random effects: η Covariates: Wt Dynamical variables: Depot, Central Derived: dv Observed: dv
alternative_param_no_q = ( cl = 0.5, tv = 9.0, ka = 1.3, Ω = Diagonal([0.18,0.04, 0.03]), σ_prop = 0.04) result_misspec = fit(mymodel_misspec, Subject.(sims), alternative_param_no_q, Pumas.FOCE())
Iter Function value Gradient norm
0 9.859164e+04 1.740033e+05
* time: 7.104873657226562e-5
1 2.415137e+04 2.231621e+04
* time: 0.03387093544006348
2 2.133314e+04 1.625040e+04
* time: 0.05183601379394531
3 1.707315e+04 6.560504e+03
* time: 0.0690610408782959
4 1.579818e+04 3.211978e+03
* time: 0.08685612678527832
5 1.522797e+04 1.311890e+03
* time: 0.10706901550292969
6 1.506519e+04 4.674492e+02
* time: 0.1273820400238037
7 1.503074e+04 4.176990e+02
* time: 0.14860200881958008
8 1.502596e+04 4.166385e+02
* time: 0.1678941249847412
9 1.502395e+04 4.144028e+02
* time: 0.18882393836975098
10 1.501859e+04 4.058294e+02
* time: 0.2115769386291504
11 1.500678e+04 3.833996e+02
* time: 0.23414397239685059
12 1.498072e+04 3.656211e+02
* time: 0.2589280605316162
13 1.493625e+04 4.421138e+02
* time: 0.28327202796936035
14 1.488124e+04 4.164407e+02
* time: 0.3117051124572754
15 1.483632e+04 2.634889e+02
* time: 0.33542394638061523
16 1.482048e+04 7.965700e+01
* time: 0.7111349105834961
17 1.481931e+04 8.198374e+01
* time: 0.7279119491577148
18 1.481900e+04 8.307845e+01
* time: 0.7474479675292969
19 1.481888e+04 8.323456e+01
* time: 0.7646389007568359
20 1.481854e+04 8.319847e+01
* time: 0.7841010093688965
21 1.481770e+04 8.240822e+01
* time: 0.8038160800933838
22 1.481553e+04 7.930206e+01
* time: 0.8244450092315674
23 1.481048e+04 7.079104e+01
* time: 0.8451180458068848
24 1.480006e+04 6.733081e+01
* time: 0.8672780990600586
25 1.478431e+04 7.223456e+01
* time: 0.8882019519805908
26 1.477247e+04 4.015798e+01
* time: 0.9095799922943115
27 1.476993e+04 1.277528e+01
* time: 0.9284911155700684
28 1.476977e+04 1.326386e+01
* time: 0.9467189311981201
29 1.476976e+04 1.289432e+01
* time: 0.9658610820770264
30 1.476976e+04 1.236934e+01
* time: 0.9837970733642578
31 1.476975e+04 1.209211e+01
* time: 1.00142502784729
32 1.476974e+04 1.152062e+01
* time: 1.0193650722503662
33 1.476972e+04 1.095485e+01
* time: 1.0393180847167969
34 1.476967e+04 1.132512e+01
* time: 1.0587670803070068
35 1.476953e+04 1.164147e+01
* time: 1.0790770053863525
36 1.476922e+04 1.142438e+01
* time: 1.1005918979644775
37 1.476860e+04 9.716891e+00
* time: 1.450066089630127
38 1.476772e+04 6.916953e+00
* time: 1.469465970993042
39 1.476680e+04 6.698177e+00
* time: 1.4884650707244873
40 1.476619e+04 3.851861e+00
* time: 1.5096969604492188
41 1.476597e+04 3.738518e+00
* time: 1.5283639430999756
42 1.476586e+04 3.605082e+00
* time: 1.5476250648498535
43 1.476577e+04 3.396881e+00
* time: 1.565032958984375
44 1.476569e+04 3.040295e+00
* time: 1.5843710899353027
45 1.476565e+04 3.129574e+00
* time: 1.603437900543213
46 1.476565e+04 3.216652e+00
* time: 1.6205730438232422
47 1.476565e+04 3.228074e+00
* time: 1.636596918106079
48 1.476565e+04 3.262238e+00
* time: 1.6551320552825928
49 1.476565e+04 3.287127e+00
* time: 1.6746959686279297
50 1.476565e+04 3.322913e+00
* time: 1.6950440406799316
51 1.476564e+04 3.361729e+00
* time: 1.7162771224975586
52 1.476562e+04 3.399238e+00
* time: 1.738090991973877
53 1.476559e+04 3.401055e+00
* time: 1.7605550289154053
54 1.476550e+04 3.276514e+00
* time: 1.7825961112976074
55 1.476531e+04 4.939640e+00
* time: 1.8046050071716309
56 1.476502e+04 5.861639e+00
* time: 1.8261959552764893
57 1.476472e+04 3.690961e+00
* time: 1.8460190296173096
58 1.476464e+04 5.173962e-01
* time: 2.1728599071502686
59 1.476463e+04 2.828018e-01
* time: 2.188831090927124
60 1.476463e+04 2.398396e-01
* time: 2.2033820152282715
61 1.476463e+04 2.390774e-01
* time: 2.2169270515441895
62 1.476463e+04 2.375504e-01
* time: 2.2285749912261963
63 1.476463e+04 2.338593e-01
* time: 2.242110013961792
64 1.476463e+04 2.285027e-01
* time: 2.255863904953003
65 1.476463e+04 2.553134e-01
* time: 2.269639015197754
66 1.476463e+04 3.838928e-01
* time: 2.2832109928131104
67 1.476463e+04 5.895030e-01
* time: 2.297606945037842
68 1.476463e+04 9.050444e-01
* time: 2.3129100799560547
69 1.476463e+04 1.353121e+00
* time: 2.3279991149902344
70 1.476462e+04 1.861178e+00
* time: 2.344688892364502
71 1.476461e+04 2.131338e+00
* time: 2.3633439540863037
72 1.476460e+04 1.738425e+00
* time: 2.3817999362945557
73 1.476458e+04 7.231803e-01
* time: 2.401494026184082
74 1.476457e+04 5.235524e-02
* time: 2.420306921005249
75 1.476457e+04 3.138355e-01
* time: 2.4377810955047607
76 1.476457e+04 3.484829e-01
* time: 2.4535489082336426
77 1.476457e+04 2.206247e-01
* time: 2.4712069034576416
78 1.476457e+04 5.036072e-02
* time: 2.488598108291626
79 1.476457e+04 3.848967e-02
* time: 2.505247116088867
80 1.476457e+04 6.104728e-02
* time: 2.519792079925537
81 1.476457e+04 4.232548e-02
* time: 2.5337600708007812
82 1.476457e+04 9.751373e-03
* time: 2.5492119789123535
83 1.476457e+04 9.243068e-03
* time: 2.870316982269287
84 1.476457e+04 1.377599e-02
* time: 2.880934953689575
85 1.476457e+04 9.204716e-03
* time: 2.891993999481201
86 1.476457e+04 2.118159e-03
* time: 2.9047040939331055
87 1.476457e+04 2.049893e-03
* time: 2.9159131050109863
88 1.476457e+04 2.986712e-03
* time: 2.9262280464172363
89 1.476457e+04 1.913842e-03
* time: 2.935673952102661
90 1.476457e+04 3.799866e-04
* time: 2.9460389614105225
FittedPumasModel
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -14764.568
Number of subjects: 24
Number of parameters: Fixed Optimized
0 7
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
---------------------
Estimate
---------------------
cl 1.032
tv 21.928
ka 1.1969e7
Ω₁,₁ 0.20656
Ω₂,₂ 0.031515
Ω₃,₃ 0.059073
σ_prop 0.44644
---------------------
First off, the absorption flow parameters Ka is quite off the charts, so that would be a warning sign off the bat, but let us try to use a tool in the toolbox to asses the goodness of fit. We get these by using the inspect function
res_inspect_2cmp = inspect(result) res_inspect_1cmp = inspect(result_misspec)
FittedPumasModelInspection Fitting was successful: true Likehood approximations used for weighted residuals : Pumas.FOCE()
goodness_of_fit(res_inspect_1cmp, figure = (resolution = (1200, 800),), axis = (title = "Mis-specified model",))
goodness_of_fit(res_inspect_2cmp, figure = (resolution = (1200, 800),), axis = (title = "True model",))
The weighted residuals should be standard normally distributed with throughout the time domain. We see that this is the case for the correctly specified model, but certainly not for the mis-specified model. That latter has a very clear pattern in time. This comes from the fact that the one compartment model is not able to capture the change in slope as time progresses, so it can never accurately capture the curves generated by a two compartment model.
Conclusion
This tutorial showed how to use fitting in Pumas.jl based on a simulated data set. There are many more models and simulation experiments to explore. Please try out fit on your own data and model, and reach out if further questions or problems come up.