Fitting models in Pumas
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.
using Pumas, Plots, Random Random.seed!(0);
MersenneTwister(UInt32[0x00000000], Random.DSFMT.DSFMT_state(Int32[74839879 7, 1073523691, -1738140313, 1073664641, -1492392947, 1073490074, -162528183 9, 1073254801, 1875112882, 1073717145 … 943540191, 1073626624, 1091647724 , 1073372234, -1273625233, -823628301, 835224507, 991807863, 382, 0]), [0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x000 00000000000000000000000000000, 0x00000000000000000000000000000000, 0x000000 00000000000000000000000000, 0x00000000000000000000000000000000, 0x000000000 00000000000000000000000, 0x00000000000000000000000000000000, 0x000000000000 00000000000000000000, 0x00000000000000000000000000000000, 0x000000000000000 00000000000000000 … 0x00000000000000000000000000000000, 0x000000000000000 00000000000000000, 0x00000000000000000000000000000000, 0x000000000000000000 00000000000000, 0x00000000000000000000000000000000, 0x000000000000000000000 00000000000, 0x00000000000000000000000000000000, 0x000000000000000000000000 00000000, 0x00000000000000000000000000000000, 0x000000000000000000000000000 00000], 1002, 0)
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 | |
|---|---|---|---|---|---|---|---|---|---|
| Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | |
| 1 | 0.0 | 1 | 100.0 | 1 | 4.0 | 2 | 0.0 | 0.0 | 0 |
As ususal, 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 = Population(map(i -> Subject(id = i, events = repeated_dose_regimen, observations =(dv=Float64[],), covariates = choose_covariates()), 1:24))
Population Subjects: 24 Covariates: Wt Observables: 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 indiviual 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_param(mymodel) obs = simobs(mymodel, pop, param, obstimes=1:1:72) plot(obs)
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.(obs), param, Pumas.FOCEI())
Iter Function value Gradient norm
0 9.509662e+03 4.974134e+02
* time: 0.00012087821960449219
1 9.509087e+03 5.235690e+01
* time: 0.6147558689117432
2 9.508605e+03 1.141497e+01
* time: 0.7336869239807129
3 9.508415e+03 3.304623e+00
* time: 1.092930793762207
4 9.508362e+03 3.055397e+00
* time: 1.1886138916015625
5 9.508092e+03 7.897585e+00
* time: 1.2845778465270996
6 9.507916e+03 7.176151e+00
* time: 1.3791728019714355
7 9.507822e+03 1.852197e+00
* time: 1.4786708354949951
8 9.507810e+03 1.416720e+00
* time: 1.5947978496551514
9 9.507794e+03 1.847344e+00
* time: 1.6789729595184326
10 9.507765e+03 3.646199e+00
* time: 1.7948508262634277
11 9.507727e+03 4.127054e+00
* time: 1.887901782989502
12 9.507701e+03 2.259724e+00
* time: 1.9814648628234863
13 9.507695e+03 4.267957e-01
* time: 2.0880179405212402
14 9.507695e+03 1.529258e-01
* time: 2.1761229038238525
15 9.507695e+03 1.458034e-01
* time: 2.256622791290283
16 9.507695e+03 1.374069e-01
* time: 2.357555866241455
17 9.507694e+03 1.983480e-01
* time: 2.446509838104248
18 9.507694e+03 1.832846e-01
* time: 2.5351829528808594
19 9.507694e+03 8.495916e-02
* time: 2.638772964477539
20 9.507694e+03 1.434653e-02
* time: 2.728404998779297
21 9.507694e+03 1.142185e-03
* time: 2.8144118785858154
22 9.507694e+03 1.139690e-03
* time: 2.9424610137939453
23 9.507694e+03 4.554493e-02
* time: 3.034693956375122
FittedPumasModel
Successful minimization: false
Likelihood approximation: Pumas.FOCEI
Log-likelihood value: -9507.6936
Number of subjects: 24
Number of parameters: Fixed Optimized
0 6
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
---------------------
Estimate
---------------------
cl 0.87552
tv 9.9608
ka 0.97115
q 0.4989
Ω₁,₁ 0.99073
Ω₂,₂ 0.061093
Ω₃,₃ 0.037914
σ_prop 0.029516
---------------------
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.(obs), alternative_param, Pumas.FOCEI())
Iter Function value Gradient norm
0 2.171648e+04 4.024499e+04
* time: 4.220008850097656e-5
1 1.374758e+04 1.846694e+04
* time: 0.2631080150604248
2 1.062678e+04 4.399382e+03
* time: 0.3806910514831543
3 1.035016e+04 1.328733e+03
* time: 0.49639415740966797
4 1.028054e+04 1.295589e+03
* time: 0.5840921401977539
5 1.021065e+04 1.841083e+03
* time: 0.6791222095489502
6 9.949460e+03 6.144939e+03
* time: 0.7906382083892822
7 9.595062e+03 3.239578e+03
* time: 0.8814511299133301
8 9.569951e+03 4.301083e+02
* time: 0.9913311004638672
9 9.568758e+03 1.318693e+02
* time: 1.0890111923217773
10 9.567892e+03 2.141681e+02
* time: 1.1796600818634033
11 9.561440e+03 7.908701e+02
* time: 1.277250051498413
12 9.552319e+03 1.086506e+03
* time: 1.396630048751831
13 9.545273e+03 6.937697e+02
* time: 1.4925100803375244
14 9.543505e+03 1.872858e+02
* time: 1.6035311222076416
15 9.543171e+03 5.062953e+01
* time: 1.689317226409912
16 9.542870e+03 1.524584e+02
* time: 1.774001121520996
17 9.542118e+03 3.573733e+02
* time: 1.8821971416473389
18 9.540396e+03 6.307445e+02
* time: 1.97446608543396
19 9.536735e+03 9.350193e+02
* time: 2.071498155593872
20 9.530476e+03 1.095649e+03
* time: 2.1951160430908203
21 9.522366e+03 8.881636e+02
* time: 2.3025741577148438
22 9.516800e+03 3.390603e+02
* time: 2.4035191535949707
23 9.515844e+03 3.149418e+01
* time: 2.523315191268921
24 9.515776e+03 3.246223e+01
* time: 2.622610092163086
25 9.515691e+03 4.957381e+01
* time: 2.743077039718628
26 9.515502e+03 5.679289e+01
* time: 2.8480241298675537
27 9.515158e+03 3.770801e+01
* time: 2.9513301849365234
28 9.514635e+03 1.984292e+01
* time: 3.0782461166381836
29 9.514115e+03 5.875926e+01
* time: 3.181462049484253
30 9.513726e+03 7.526339e+01
* time: 3.2869460582733154
31 9.513454e+03 5.726155e+01
* time: 3.4090402126312256
32 9.513370e+03 2.345084e+01
* time: 3.5236990451812744
33 9.513346e+03 1.495845e+01
* time: 3.6385672092437744
34 9.513305e+03 3.106059e+01
* time: 3.740854024887085
35 9.513241e+03 6.757386e+01
* time: 3.8473870754241943
36 9.513059e+03 1.342898e+02
* time: 3.96631121635437
37 9.512653e+03 2.234116e+02
* time: 4.072812080383301
38 9.511771e+03 3.222847e+02
* time: 4.179495096206665
39 9.510273e+03 3.589024e+02
* time: 4.303644180297852
40 9.508599e+03 2.512825e+02
* time: 4.411496162414551
41 9.507861e+03 8.450873e+01
* time: 4.516286134719849
42 9.507785e+03 1.614494e+01
* time: 4.631969213485718
43 9.507780e+03 3.390312e+00
* time: 4.7293541431427
44 9.507779e+03 1.391292e+00
* time: 4.840573072433472
45 9.507779e+03 1.398817e+00
* time: 4.935310125350952
46 9.507779e+03 1.402850e+00
* time: 5.0290892124176025
47 9.507779e+03 1.407965e+00
* time: 5.138683080673218
48 9.507778e+03 1.411078e+00
* time: 5.230583190917969
49 9.507777e+03 1.405185e+00
* time: 5.325125217437744
50 9.507773e+03 1.368683e+00
* time: 5.443305015563965
51 9.507763e+03 1.726587e+00
* time: 5.546467065811157
52 9.507745e+03 2.360911e+00
* time: 5.645001173019409
53 9.507719e+03 2.370152e+00
* time: 5.757469177246094
54 9.507700e+03 1.436938e+00
* time: 5.854623079299927
55 9.507695e+03 4.604145e-01
* time: 5.956088066101074
56 9.507694e+03 2.545239e-01
* time: 6.069289207458496
57 9.507694e+03 8.748370e-02
* time: 6.1631481647491455
58 9.507694e+03 8.781354e-02
* time: 6.251407146453857
59 9.507694e+03 8.779529e-02
* time: 6.369810104370117
60 9.507694e+03 8.780567e-02
* time: 6.501018047332764
FittedPumasModel
Successful minimization: false
Likelihood approximation: Pumas.FOCEI
Log-likelihood value: -9507.694
Number of subjects: 24
Number of parameters: Fixed Optimized
0 6
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
---------------------
Estimate
---------------------
cl 0.87564
tv 9.9607
ka 0.9711
q 0.49891
Ω₁,₁ 0.99095
Ω₂,₂ 0.061185
Ω₃,₃ 0.037614
σ_prop 0.029516
---------------------
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)
Calculating: variance-covariance matrix. Done.
Asymptotic inference results
Successful minimization: false
Likelihood approximation: Pumas.FOCEI
Log-likelihood value: -9507.6936
Number of subjects: 24
Number of parameters: Fixed Optimized
0 6
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
--------------------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------------------
cl 0.87552 0.17807 [0.52651 ; 1.2245 ]
tv 9.9608 0.50239 [8.9761 ; 10.945 ]
ka 0.97115 0.038954 [0.8948 ; 1.0475 ]
q 0.4989 0.00103 [0.49689 ; 0.50092 ]
Ω₁,₁ 0.99073 0.27442 [0.45288 ; 1.5286 ]
Ω₂,₂ 0.061093 0.014388 [0.032894; 0.089292]
Ω₃,₃ 0.037914 0.011946 [0.014501; 0.061328]
σ_prop 0.029516 0.00051089 [0.028515; 0.030517]
--------------------------------------------------------------------
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 = Population(map(i -> Subject(id = i, events = repeated_dose_regimen, observations =(dv=Float64[],), covariates = choose_covariates()), 1:100)) obs_big = simobs(mymodel, pop_big, param, obstimes=1:1:72) result_big = fit(mymodel, Subject.(obs_big), param, Pumas.FOCEI())
Iter Function value Gradient norm
0 3.923167e+04 1.061584e+03
* time: 4.291534423828125e-5
1 3.923104e+04 5.439570e+01
* time: 1.1168508529663086
2 3.923062e+04 5.493562e+01
* time: 1.6585218906402588
3 3.922981e+04 7.061316e+00
* time: 2.275825023651123
4 3.922971e+04 6.257481e+00
* time: 2.724810838699341
5 3.922917e+04 8.285576e+00
* time: 3.2103898525238037
6 3.922915e+04 2.215932e+00
* time: 3.639216899871826
7 3.922913e+04 2.501836e+00
* time: 4.043552875518799
8 3.922910e+04 5.866463e+00
* time: 4.478799819946289
9 3.922909e+04 4.914102e+00
* time: 4.9104249477386475
10 3.922908e+04 1.679291e+00
* time: 5.345788955688477
11 3.922908e+04 2.096004e-01
* time: 5.757008790969849
12 3.922908e+04 3.402351e-02
* time: 6.146726846694946
13 3.922908e+04 2.982270e-02
* time: 6.5279388427734375
14 3.922908e+04 2.982270e-02
* time: 7.275569915771484
15 3.922908e+04 2.982270e-02
* time: 8.016044855117798
16 3.922908e+04 2.982270e-02
* time: 8.782685995101929
infer(result_big)
Calculating: variance-covariance matrix. Done.
Asymptotic inference results
Successful minimization: true
Likelihood approximation: Pumas.FOCEI
Log-likelihood value: -39229.078
Number of subjects: 100
Number of parameters: Fixed Optimized
0 6
Observation records: Active Missing
dv: 7200 0
Total: 7200 0
--------------------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------------------
cl 0.93233 0.093072 [0.74991 ; 1.1147 ]
tv 10.39 0.28913 [9.8235 ; 10.957 ]
ka 1.0051 0.022993 [0.96008 ; 1.0502 ]
q 0.50059 0.00052754 [0.49955 ; 0.50162 ]
Ω₁,₁ 0.99768 0.15308 [0.69765 ; 1.2977 ]
Ω₂,₂ 0.07766 0.015107 [0.048051; 0.10727 ]
Ω₃,₃ 0.049499 0.010411 [0.029093; 0.069906]
σ_prop 0.029905 0.00023157 [0.029451; 0.030358]
--------------------------------------------------------------------
This time we see similar estimates, but much narrower confidence intervals across the board.
Estimating a misspecified model
To explore some of the diagnostics tools available in Julia, 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: cp, dv Observed: cp, 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.(obs), alternative_param_no_q, Pumas.FOCEI())
Iter Function value Gradient norm
0 1.130666e+05 2.034877e+05
* time: 4.9114227294921875e-5
1 2.587757e+04 2.631157e+04
* time: 0.23863816261291504
2 2.252342e+04 1.917207e+04
* time: 0.4143691062927246
3 1.742887e+04 7.806532e+03
* time: 0.6188540458679199
4 1.586746e+04 3.869815e+03
* time: 0.8084421157836914
5 1.514158e+04 1.626030e+03
* time: 0.9823591709136963
6 1.491812e+04 6.143775e+02
* time: 1.1739201545715332
7 1.486465e+04 4.559947e+02
* time: 1.3434960842132568
8 1.485667e+04 4.548559e+02
* time: 1.5338661670684814
9 1.485453e+04 4.529243e+02
* time: 1.7061591148376465
10 1.485029e+04 4.470285e+02
* time: 1.9110851287841797
11 1.484016e+04 4.302327e+02
* time: 2.1178431510925293
12 1.481713e+04 3.885750e+02
* time: 2.3135690689086914
13 1.477276e+04 4.540622e+02
* time: 2.534414052963257
14 1.470854e+04 4.794530e+02
* time: 2.7332310676574707
15 1.464361e+04 3.845968e+02
* time: 2.8944289684295654
16 1.460286e+04 1.716652e+02
* time: 3.077684164047241
17 1.459619e+04 8.137585e+01
* time: 3.228778123855591
18 1.459564e+04 8.232241e+01
* time: 3.4018490314483643
19 1.459533e+04 8.310168e+01
* time: 3.5583319664001465
20 1.459522e+04 8.314237e+01
* time: 3.720219135284424
21 1.459360e+04 8.253116e+01
* time: 3.8814001083374023
22 1.459079e+04 7.990532e+01
* time: 4.0680601596832275
23 1.458301e+04 7.043052e+01
* time: 4.241152048110962
24 1.456850e+04 6.365034e+01
* time: 4.425739049911499
25 1.454808e+04 5.721254e+01
* time: 4.592349052429199
26 1.453536e+04 2.934679e+01
* time: 4.770431041717529
27 1.453260e+04 1.249869e+01
* time: 4.959182977676392
28 1.453239e+04 1.445128e+01
* time: 5.1073760986328125
29 1.453238e+04 1.444180e+01
* time: 5.256748199462891
30 1.453238e+04 1.416957e+01
* time: 5.4203760623931885
31 1.453237e+04 1.384491e+01
* time: 5.571076154708862
32 1.453237e+04 1.345241e+01
* time: 5.7346251010894775
33 1.453235e+04 1.278907e+01
* time: 5.884273052215576
34 1.453231e+04 1.168377e+01
* time: 6.0557541847229
35 1.453220e+04 9.925444e+00
* time: 6.207094192504883
36 1.453193e+04 1.017005e+01
* time: 6.375234127044678
37 1.453132e+04 9.562097e+00
* time: 6.531067132949829
38 1.453020e+04 7.132143e+00
* time: 6.712423086166382
39 1.452882e+04 6.237787e+00
* time: 6.866991996765137
40 1.452773e+04 4.886911e+00
* time: 7.041015148162842
41 1.452721e+04 2.309835e+00
* time: 7.204063177108765
42 1.452709e+04 8.448814e-01
* time: 7.387569189071655
43 1.452708e+04 1.028132e+00
* time: 7.5365471839904785
44 1.452708e+04 1.062585e+00
* time: 7.699488162994385
45 1.452708e+04 1.067032e+00
* time: 7.828843116760254
46 1.452708e+04 1.072589e+00
* time: 7.981573104858398
47 1.452708e+04 1.075766e+00
* time: 8.11764407157898
48 1.452708e+04 1.081692e+00
* time: 8.277626037597656
49 1.452708e+04 1.089409e+00
* time: 8.420660018920898
50 1.452708e+04 1.100283e+00
* time: 8.588013172149658
51 1.452708e+04 1.111684e+00
* time: 8.736998081207275
52 1.452707e+04 1.114455e+00
* time: 8.90843915939331
53 1.452707e+04 1.080389e+00
* time: 9.06112813949585
54 1.452705e+04 1.087978e+00
* time: 9.23018503189087
55 1.452703e+04 1.115299e+00
* time: 9.387166023254395
56 1.452701e+04 6.293728e-01
* time: 9.561125993728638
57 1.452700e+04 4.329917e-01
* time: 9.715710163116455
58 1.452700e+04 3.560297e-01
* time: 9.887900114059448
59 1.452700e+04 3.347408e-01
* time: 10.018774032592773
60 1.452700e+04 3.332461e-01
* time: 10.141225099563599
61 1.452700e+04 3.309391e-01
* time: 10.29244613647461
62 1.452700e+04 3.289458e-01
* time: 10.422115087509155
63 1.452700e+04 3.262512e-01
* time: 10.574211120605469
64 1.452700e+04 3.228271e-01
* time: 10.716219186782837
65 1.452700e+04 3.177641e-01
* time: 10.87358808517456
66 1.452700e+04 3.097423e-01
* time: 11.01814317703247
67 1.452700e+04 5.004012e-01
* time: 11.15893816947937
68 1.452700e+04 8.731792e-01
* time: 11.32102108001709
69 1.452699e+04 1.396165e+00
* time: 11.469442129135132
70 1.452698e+04 1.994401e+00
* time: 11.637778043746948
71 1.452696e+04 2.426608e+00
* time: 11.786231994628906
72 1.452691e+04 2.240644e+00
* time: 11.943849086761475
73 1.452686e+04 8.847438e-01
* time: 12.090716123580933
74 1.452683e+04 2.401597e-01
* time: 12.253134965896606
75 1.452683e+04 6.924838e-01
* time: 12.397193193435669
76 1.452682e+04 7.426297e-01
* time: 12.559774160385132
77 1.452682e+04 4.989830e-01
* time: 12.69515609741211
78 1.452682e+04 1.387091e-01
* time: 12.86349105834961
79 1.452681e+04 7.360207e-02
* time: 13.012005090713501
80 1.452681e+04 1.312669e-01
* time: 13.171314001083374
81 1.452681e+04 1.092161e-01
* time: 13.295031070709229
82 1.452681e+04 4.168175e-02
* time: 13.425817966461182
83 1.452681e+04 9.372777e-03
* time: 13.574540138244629
84 1.452681e+04 2.615442e-02
* time: 13.70778512954712
85 1.452681e+04 2.321991e-02
* time: 13.849776983261108
86 1.452681e+04 1.010919e-02
* time: 13.971124172210693
87 1.452681e+04 9.626878e-04
* time: 14.091981172561646
FittedPumasModel
Successful minimization: true
Likelihood approximation: Pumas.FOCEI
Log-likelihood value: -14526.813
Number of subjects: 24
Number of parameters: Fixed Optimized
0 5
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
--------------------
Estimate
--------------------
cl 1.2266
tv 24.739
ka 279250.0
Ω₁,₁ 0.42901
Ω₂,₂ 0.12159
Ω₃,₃ 0.04714
σ_prop 0.49129
--------------------
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 fit: the weighted residuals. We get these by using the wresiduals function
wres = wresiduals(result) wres_misspec = wresiduals(result_misspec) p1 = plot([w.wres.dv for w in wres], title="Correctly specified", legend=false) p2 = plot([w.wres.dv for w in wres_misspec], title = "Misspecified", legend=false) plot(p1, p2)
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 misspecified 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.