# Individual parameter estimates

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?

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)
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
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
)
``````

## Fitting a model

You can fit the model using the syntax below

``````myfit = fit(theopmodel,
theopp,
param,
Pumas.FOCEI(),
``````

There are multiple post-processing steps that one would like to perform after a `fit`

## Parameter estimates

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
-----------------------
``````

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)
``````

## Inference

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.

## Inspect

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

1. Computes the predictions
``````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    │
``````
1. Computes the weighted residuals
``````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    │
``````
1. Computes the empirical Bayes estimates
``````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],)
``````

## Individual parameters

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.

1 Like

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.