Hi all,
is there an easy way to extract individual parameter estimates based on the EBE’s? I imagine all parameters appearing in the @pre block could be added to the FittedPumasModelInspection object. As they can be time-dependent (via t or via covariates) this seems the correct place for them?
Thanks in advance.
Hi, welcome to the community! I would like to take this opportunity to showcase the different post-processing functions available through a complete example.
For this lets use the Theoph dataset that comes with installation.
theopp = read_pumas(example_data("event_data/THEOPP"), covariates = [:SEX,:WT])
Here is a simple model and the corresponding parameters
theopmodel = @model begin
@param begin
θ₁ ∈ RealDomain(lower=0.1, upper=5.0, init=2.77)
θ₂ ∈ RealDomain(lower=0.0008, upper=0.5, init=0.0781)
θ₃ ∈ RealDomain(lower=0.004, upper=0.9, init=0.0363)
θ₄ ∈ RealDomain(lower=0.1, upper=5.0, init=1.5)
Ω ∈ PSDDomain(3)
σ²_add ∈ RealDomain(lower=0.001, init=sqrt(0.388))
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
Ka = SEX == 0 ? θ₁ + η[1] : θ₄ + η[1]
K = θ₂ + η[2]
CL = θ₃*WT + η[3]
Vc = CL/K
SC = CL/K/WT
end
@covariates SEX WT
@vars begin
conc = Central / SC
end
@dynamics Depots1Central1
@derived begin
dv ~ @. Normal(conc, sqrt(σ²_add))
end
end
param = (
θ₁ = 2.77, #Ka MEAN ABSORPTION RATE CONSTANT for SEX = 1(1/HR)
θ₂ = 0.0781, #K MEAN ELIMINATION RATE CONSTANT (1/HR)
θ₃ = 0.0363, #SLP SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
θ₄ = 1.5, #Ka MEAN ABSORPTION RATE CONSTANT for SEX=0 (1/HR)
Ω = diagm(0 => [5.55, 0.0024, 0.515]), # update to block diagonal
σ²_add = 0.388
)
You can fit the model using the syntax below
myfit = fit(theopmodel,
theopp,
param,
Pumas.FOCEI(),
ensemblealg = EnsembleThreads())
There are multiple post-processing steps that one would like to perform after a fit
The result of the model fit is an object of type FittedPumasModel, in this case, myfit.
Typing myfit in the Julia REPL/console will result in the printing of fitted parameter estimates.
julia> myfit
FittedPumasModel
Successful minimization: true
Likelihood approximation: Pumas.FOCEI
Deviance: 114.55288
Total number of observation records: 132
Number of active observation records: 132
Number of subjects: 12
-----------------------
Estimate
-----------------------
θ₁ 1.747
θ₂ 0.086384
θ₃ 0.04019
θ₄ 1.9824
Ω₁,₁ 1.5428
Ω₂,₁ 0.003224
Ω₃,₁ -0.090697
Ω₂,₂ 0.00013038
Ω₃,₂ 0.0074753
Ω₃,₃ 0.48048
σ²_add 0.4757
-----------------------
While this is a convenient way to look at the model fit parameters, one may want to see the results
in a dataframe which can be either rendered as a html or latex table if you are working in a jmd file
or you may want to export it to CSV output.
You can use the coeftable(fit) to get a table of parameter estimates
julia> coeftable(myfit)
11×2 DataFrame
│ Row │ parameter │ estimate │
│ │ String │ Float64 │
├─────┼───────────┼─────────────┤
│ 1 │ θ₁ │ 1.74696 │
│ 2 │ θ₂ │ 0.0863839 │
│ 3 │ θ₃ │ 0.0401902 │
│ 4 │ θ₄ │ 1.98236 │
│ 5 │ Ω₁,₁ │ 1.54278 │
│ 6 │ Ω₂,₁ │ 0.00322403 │
│ 7 │ Ω₃,₁ │ -0.0906974 │
│ 8 │ Ω₂,₂ │ 0.000130385 │
│ 9 │ Ω₃,₂ │ 0.00747534 │
│ 10 │ Ω₃,₃ │ 0.480476 │
│ 11 │ σ²_add │ 0.475701 │
If working in a html or latex rendering engine like Weave.jl, these tables will be rendered directly. If
you would rather save the files, you can write it out to disk, e.g.
CSV.write("./reports/myfit_params.csv", myfit)
The parameters of the fit can also be accessed using the coef function that prints the result
as a NamedTuple().
julia> myfit_coef = coef(myfit)
(θ₁ = 1.746957083389781, θ₂ = 0.08638392390943932, θ₃ = 0.04019020329335648, θ₄ = 1.982355904736838, Ω = PDMats.PDMat{Float64,Array{Float64,2}}(3, [1.5427849261948847 0.0032240317997140745 -0.09069740250273521; 0.0032240317997140745 0.0001303848008102982 0.007475341834314581; -0.09069740250273521 0.007475341834314581 0.4804764539570251], LinearAlgebra.Cholesky{Float64,Array{Float64,2}}([1.2420889365077223 0.00259565293994069 -0.0730200550354643; 0.0 0.011119684645963456 0.6893070081465335; 0.0 0.0 0.0006115878659088313], 'U', 0)), σ²_add = 0.47570143542677834)
The advantage of using coef is that you can use the output as initial estimates for a new model, e.g.
if you want to restart your run with these minimized set of estimates. You just pass this in to the fit as
coef(myfit) instead of the param.
However, if you would like to update your model by adding a new covariate effect
(e.g. you want to add SEX on Cl in the model above), and want to the use the
final estimates from the myfit model to begin you work, then you can use Julia syntax to extend
the NamedTuple like this
julia> inits_for_new_model = (myfit_coef..., sexoncl = 0.5)
(θ₁ = 1.746957083389781, θ₂ = 0.08638392390943932, θ₃ = 0.04019020329335648, θ₄ = 1.982355904736838, Ω = PDMats.PDMat{Float64,Array{Float64,2}}(3, [1.5427849261948847 0.0032240317997140745 -0.09069740250273521; 0.0032240317997140745 0.0001303848008102982 0.007475341834314581; -0.09069740250273521 0.007475341834314581 0.4804764539570251], LinearAlgebra.Cholesky{Float64,Array{Float64,2}}([1.2420889365077223 0.00259565293994069 -0.0730200550354643; 0.0 0.011119684645963456 0.6893070081465335; 0.0 0.0 0.0006115878659088313], 'U', 0)), σ²_add = 0.47570143542677834, sexoncl = 0.5)
The covariance step to compute the standard errors and confidence intervals is not part of the fit by default.
You can make an inference on your myfit using the infer function as below.
myfit_infer = infer(myfit)
The coeftable function that was introduced above can also be used to generate a dataframe on
the infer output so that you get a table of parameter estimates/confidence intervals together.
After making an inference, we can inspect the model performance by looking at the diagnostics.
julia> myfit_inspect = DataFrame(inspect(myfit))
Calculating: predictions, weighted residuals, empirical bayes. Done.
132×14 DataFrame
│ Row │ id │ time │ evid │ dv │ dv_pred │ dv_ipred │ dv_wres │ dv_iwres │ wres_approx │ η_1 │ η_2 │ η_3 │ SEX │ WT │
│ │ String │ Float64 │ Int64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Pumas.FOCEI │ Float64 │ Float64 │ Float64 │ Int64 │ Float64 │
├─────┼────────┼─────────┼───────┼─────────┼─────────┼──────────┼───────────┼────────────┼─────────────┼────────────┼─────────────┼──────────┼───────┼─────────┤
│ 1 │ 1 │ 0.0 │ 0 │ 0.74 │ 0.0 │ 0.0 │ 1.07291 │ 1.07291 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
│ 2 │ 1 │ 0.25 │ 0 │ 2.84 │ 3.02238 │ 3.81171 │ -1.42921 │ -1.40886 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
│ 3 │ 1 │ 0.57 │ 0 │ 6.57 │ 5.29503 │ 6.78682 │ -0.220898 │ -0.314366 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
│ 4 │ 1 │ 1.12 │ 0 │ 10.5 │ 6.96696 │ 9.14446 │ 2.31839 │ 1.96537 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
│ 5 │ 1 │ 2.02 │ 0 │ 9.66 │ 7.36781 │ 9.96997 │ 0.173548 │ -0.449416 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
│ 6 │ 1 │ 3.82 │ 0 │ 8.58 │ 6.5236 │ 9.24889 │ -0.222183 │ -0.969808 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
│ 7 │ 1 │ 5.1 │ 0 │ 8.36 │ 5.84979 │ 8.54664 │ 0.475421 │ -0.270606 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
│ 8 │ 1 │ 7.03 │ 0 │ 7.47 │ 4.95246 │ 7.56659 │ 0.583903 │ -0.14004 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
│ 9 │ 1 │ 9.05 │ 0 │ 6.89 │ 4.15951 │ 6.65907 │ 1.02907 │ 0.334819 │ FOCEI() │ -0.0933388 │ -0.0231288 │ -1.41617 │ 0 │ 79.6 │
⋮
│ 123 │ 12 │ 0.25 │ 0 │ 1.25 │ 4.4001 │ 2.62294 │ -2.23931 │ -1.9906 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 124 │ 12 │ 0.5 │ 0 │ 3.96 │ 6.98669 │ 4.59662 │ -1.29335 │ -0.923026 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 125 │ 12 │ 1.0 │ 0 │ 7.82 │ 9.28438 │ 7.15797 │ 0.561269 │ 0.959873 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 126 │ 12 │ 2.0 │ 0 │ 9.72 │ 9.7949 │ 9.15466 │ 0.644427 │ 0.819682 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 127 │ 12 │ 3.52 │ 0 │ 9.75 │ 8.77673 │ 9.19084 │ 0.912403 │ 0.810715 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 128 │ 12 │ 5.07 │ 0 │ 8.57 │ 7.68604 │ 8.3484 │ 0.524142 │ 0.321295 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 129 │ 12 │ 7.07 │ 0 │ 6.59 │ 6.46694 │ 7.16639 │ -0.610685 │ -0.835697 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 130 │ 12 │ 9.03 │ 0 │ 6.11 │ 5.45969 │ 6.1358 │ 0.178745 │ -0.0374063 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 131 │ 12 │ 12.05 │ 0 │ 4.57 │ 4.206 │ 4.82374 │ -0.174835 │ -0.367894 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
│ 132 │ 12 │ 24.15 │ 0 │ 1.17 │ 1.47885 │ 1.83844 │ -0.861641 │ -0.969165 │ FOCEI() │ -0.948334 │ -0.00666158 │ -0.23435 │ 1 │ 60.5 │
Note that inspect does 3 things and each of these tasks can be performed individually outside an inspect call
julia> myfit_predict = DataFrame(predict(myfit))
132×8 DataFrame
│ Row │ id │ time │ evid │ dv │ dv_pred │ dv_ipred │ SEX │ WT │
│ │ String │ Float64 │ Int64 │ Float64 │ Float64 │ Float64 │ Int64 │ Float64 │
├─────┼────────┼─────────┼───────┼─────────┼─────────┼──────────┼───────┼─────────┤
│ 1 │ 1 │ 0.0 │ 0 │ 0.74 │ 0.0 │ 0.0 │ 0 │ 79.6 │
│ 2 │ 1 │ 0.25 │ 0 │ 2.84 │ 3.02238 │ 3.81171 │ 0 │ 79.6 │
│ 3 │ 1 │ 0.57 │ 0 │ 6.57 │ 5.29503 │ 6.78682 │ 0 │ 79.6 │
│ 4 │ 1 │ 1.12 │ 0 │ 10.5 │ 6.96696 │ 9.14446 │ 0 │ 79.6 │
│ 5 │ 1 │ 2.02 │ 0 │ 9.66 │ 7.36781 │ 9.96997 │ 0 │ 79.6 │
│ 6 │ 1 │ 3.82 │ 0 │ 8.58 │ 6.5236 │ 9.24889 │ 0 │ 79.6 │
│ 7 │ 1 │ 5.1 │ 0 │ 8.36 │ 5.84979 │ 8.54664 │ 0 │ 79.6 │
│ 8 │ 1 │ 7.03 │ 0 │ 7.47 │ 4.95246 │ 7.56659 │ 0 │ 79.6 │
│ 9 │ 1 │ 9.05 │ 0 │ 6.89 │ 4.15951 │ 6.65907 │ 0 │ 79.6 │
⋮
│ 123 │ 12 │ 0.25 │ 0 │ 1.25 │ 4.4001 │ 2.62294 │ 1 │ 60.5 │
│ 124 │ 12 │ 0.5 │ 0 │ 3.96 │ 6.98669 │ 4.59662 │ 1 │ 60.5 │
│ 125 │ 12 │ 1.0 │ 0 │ 7.82 │ 9.28438 │ 7.15797 │ 1 │ 60.5 │
│ 126 │ 12 │ 2.0 │ 0 │ 9.72 │ 9.7949 │ 9.15466 │ 1 │ 60.5 │
│ 127 │ 12 │ 3.52 │ 0 │ 9.75 │ 8.77673 │ 9.19084 │ 1 │ 60.5 │
│ 128 │ 12 │ 5.07 │ 0 │ 8.57 │ 7.68604 │ 8.3484 │ 1 │ 60.5 │
│ 129 │ 12 │ 7.07 │ 0 │ 6.59 │ 6.46694 │ 7.16639 │ 1 │ 60.5 │
│ 130 │ 12 │ 9.03 │ 0 │ 6.11 │ 5.45969 │ 6.1358 │ 1 │ 60.5 │
│ 131 │ 12 │ 12.05 │ 0 │ 4.57 │ 4.206 │ 4.82374 │ 1 │ 60.5 │
│ 132 │ 12 │ 24.15 │ 0 │ 1.17 │ 1.47885 │ 1.83844 │ 1 │ 60.5 │
julia> myfit_resids = DataFrame(wresiduals(myfit))
132×7 DataFrame
│ Row │ id │ time │ dv_wres │ dv_iwres │ wres_approx │ SEX │ WT │
│ │ String │ Float64 │ Float64 │ Float64 │ Pumas.FOCEI │ Int64 │ Float64 │
├─────┼────────┼─────────┼───────────┼────────────┼─────────────┼───────┼─────────┤
│ 1 │ 1 │ 0.0 │ 1.07291 │ 1.07291 │ FOCEI() │ 0 │ 79.6 │
│ 2 │ 1 │ 0.25 │ -1.42921 │ -1.40886 │ FOCEI() │ 0 │ 79.6 │
│ 3 │ 1 │ 0.57 │ -0.220898 │ -0.314366 │ FOCEI() │ 0 │ 79.6 │
│ 4 │ 1 │ 1.12 │ 2.31839 │ 1.96537 │ FOCEI() │ 0 │ 79.6 │
│ 5 │ 1 │ 2.02 │ 0.173548 │ -0.449416 │ FOCEI() │ 0 │ 79.6 │
│ 6 │ 1 │ 3.82 │ -0.222183 │ -0.969808 │ FOCEI() │ 0 │ 79.6 │
│ 7 │ 1 │ 5.1 │ 0.475421 │ -0.270606 │ FOCEI() │ 0 │ 79.6 │
│ 8 │ 1 │ 7.03 │ 0.583903 │ -0.14004 │ FOCEI() │ 0 │ 79.6 │
│ 9 │ 1 │ 9.05 │ 1.02907 │ 0.334819 │ FOCEI() │ 0 │ 79.6 │
⋮
│ 123 │ 12 │ 0.25 │ -2.23931 │ -1.9906 │ FOCEI() │ 1 │ 60.5 │
│ 124 │ 12 │ 0.5 │ -1.29335 │ -0.923026 │ FOCEI() │ 1 │ 60.5 │
│ 125 │ 12 │ 1.0 │ 0.561269 │ 0.959873 │ FOCEI() │ 1 │ 60.5 │
│ 126 │ 12 │ 2.0 │ 0.644427 │ 0.819682 │ FOCEI() │ 1 │ 60.5 │
│ 127 │ 12 │ 3.52 │ 0.912403 │ 0.810715 │ FOCEI() │ 1 │ 60.5 │
│ 128 │ 12 │ 5.07 │ 0.524142 │ 0.321295 │ FOCEI() │ 1 │ 60.5 │
│ 129 │ 12 │ 7.07 │ -0.610685 │ -0.835697 │ FOCEI() │ 1 │ 60.5 │
│ 130 │ 12 │ 9.03 │ 0.178745 │ -0.0374063 │ FOCEI() │ 1 │ 60.5 │
│ 131 │ 12 │ 12.05 │ -0.174835 │ -0.367894 │ FOCEI() │ 1 │ 60.5 │
│ 132 │ 12 │ 24.15 │ -0.861641 │ -0.969165 │ FOCEI() │ 1 │ 60.5 │
julia> myfit_ebes = empirical_bayes(myfit)
12-element Array{NamedTuple{(:η,),Tuple{Array{Float64,1}}},1}:
(η = [-0.0933387988324176, -0.0231287596895002, -1.4161692310894067],)
(η = [0.37821971388420084, 0.002521075924176358, 0.08505001149756336],)
(η = [0.4549020027332671, 0.0030150769810539484, 0.10123223719036228],)
(η = [-0.478079575517628, -0.006781308473883922, -0.33033537897568765],)
(η = [-0.2764232608639846, 0.002144685346230285, 0.18500792432344426],)
(η = [-0.590292471666874, 0.009147817587369522, 0.6782426031954858],)
(η = [-1.243482576652008, 0.004582979739288362, 0.5182848785296695],)
(η = [-0.6028975993389734, 0.004558260278199378, 0.39611047362483376],)
(η = [3.2336939515362926, -0.001873477997516174, -0.7251432116140162],)
(η = [-1.2897096767721485, -0.010639040016305933, -0.4166198053464913],)
(η = [1.0147889114680149, 0.02207025314656471, 1.1770158207301915],)
(η = [-0.9483341234405305, -0.006661575179981912, -0.23435013137913568],)
All the individual parameters can be accessed using the icoef function. This function
currently outputs all the parameters specified in the @pre block
julia> myfit_indpar = reduce(vcat, DataFrame.(icoef(myfit)))
12×6 DataFrame
│ Row │ id │ Ka │ K │ CL │ Vc │ SC │
│ │ String │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼────────┼──────────┼───────────┼─────────┼─────────┼──────────┤
│ 1 │ 1 │ 1.65362 │ 0.0632552 │ 1.78297 │ 28.187 │ 0.354108 │
│ 2 │ 2 │ 2.12518 │ 0.088905 │ 2.99482 │ 33.6856 │ 0.465271 │
│ 3 │ 3 │ 2.20186 │ 0.089399 │ 2.93464 │ 32.8263 │ 0.465622 │
│ 4 │ 4 │ 1.26888 │ 0.0796026 │ 2.59149 │ 32.5554 │ 0.447804 │
│ 5 │ 5 │ 1.47053 │ 0.0885286 │ 2.37939 │ 26.8771 │ 0.492255 │
│ 6 │ 6 │ 1.15666 │ 0.0955317 │ 3.89346 │ 40.7557 │ 0.509446 │
│ 7 │ 7 │ 0.738873 │ 0.0909669 │ 3.11457 │ 34.2385 │ 0.530008 │
│ 8 │ 8 │ 1.37946 │ 0.0909422 │ 3.22952 │ 35.5118 │ 0.503713 │
│ 9 │ 9 │ 5.21605 │ 0.0845104 │ 2.74729 │ 32.5083 │ 0.376253 │
│ 10 │ 10 │ 0.692646 │ 0.0757449 │ 1.92245 │ 25.3806 │ 0.436093 │
│ 11 │ 11 │ 2.99714 │ 0.108454 │ 3.78938 │ 34.9399 │ 0.537537 │
│ 12 │ 12 │ 1.03402 │ 0.0797223 │ 2.19716 │ 27.5601 │ 0.455539 │
We initially thought of adding individual parameters to the inspect output, but went with the option of
having a separate table output that can be easily merged with the inspect output.
myfit_inspect_infer = leftjoin(myfit_inspect, myfit_indpar, on = :id)
Note that icoef handles both constant and time varying covariates. In the example above, we
have only constant ones, but if there were time varying parameters, the table would extend to the covariate
value at each time.
You can know more about these functions using the ?fx syntax in the REPL.
Hope this helps.
Hi Vijay,
thanks for your detailed answer. Unfortunately, icoef is not available in my version v0.13.0. Also the covariates option for read_pumas in my installation is cvs. Am I using an outdated version and should I do an update?
Thanks again, KR
Hi KR,
Did you look at the latest announcement our beta release. Everything that I showcased above can be done in that release. You can update your package by simply doing ] up on the REPL
Vijay
Hi Vijay,
I had the beta release installed and I thought that came with the latest Pumas version, but apparently not. Now on v0.14.0 it works as expected, thanks!
v 0.14.0 is the latest release, and is the beta version.