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)

1 rows × 10 columns

timecmtamtevidiiaddlratedurationssroute
Float64Int64Float64Int8Float64Int64Float64Float64Int8Route…
10.01100.014.020.00.00NullRoute

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.