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.
Hi, Vijay,
itโs a plausible way to derive individual parameters.
But I have a question related to deriving individual parameters from the simulated dataset.
Iโm trying to fit pumas model with constantcoef for all parameters and then inspecting my fitted result, but it has some error as you can see.
Thank you very much for helping.
Itโs an example from Fitting models in Pumas
Welcome to the community @NamTienNguyen
Can you tell me what are your trying to do by fixing all parameters in a fit?
Thank you for your response,
Iโm needing individual parameters for the simulated dataset, so my strategy is that: I fixed all parameters (including fixed and random effect) and then fit to create result = fitted pumas model. After that, I inspect this result to derive ETAs value and individual params
If you familar with NONMEM, it equivalent to
$SIMULATION ONLYSIM
$TABLE ID TIME CL Vd ETA1 ETA2
Thank you!
Oh sorry, I have confusion with Pumas.
Is there any way for me to obtain ETAs from the simulation process? Because simobs function doesnโt have ETAs in output.
In NONME, itโs equivalent to
$SIMULATION ONLYSIM
$TABLE ID TIME CL Vd ETA1 ETA2
Thank you!
Unfortunately thereโs no easy way to grab that information after you have simulated the object in Pumas v1.1. However, if you modify your workflow slightly, it is possible. When you simulate a subject, you write
sobs = simobs(model, subject, param)
However, you also have the possibility of inputting the random effects yourself as a fourth argument, this means that you can sample them yourself and then input them. Something like
vrandeffs = [sample_randeffs(model, param) for subject in population]
simulated_obs = [simobs(model, subject, param, randeffs) for (subject, randeffs) in zip(population, vrandeffs)]
The ith element of vrandeffs will be the random effects for the ith element of simulated_obs.
I hope this is useful.
I just read it carefully in Simulation of Pumas Models and got the same idea familiar with you, but itโs time-consuming to code.
Thank you for your solution! It works successfully.
Dear @pkofod
I have a question after executing your solution.
Relating to circumstances having many random effect such as IOV for instance, then I changed ฮท in @random block to ฮท_IIV ~ MvNormal(ฮฉ_IIV) and ฮท_IOV ~ MvNormal(ฮฉ_IOV) to avoid confusion.
But ERROR is that UndefVarError: ฮท not defined
The same ERROR when I remove ฮท_IOV. I think there was a named variable in function which you used, but I donโt find any tutorials concerning this.
Itโs clear that I just merge ฮท_IIV and ฮท_IOV to ฮท and the ERROR will be passed.
But is there a solution for my interest, Thank you!
Since I donโt have your model it is a bit hard to help, but I would suspect you forgot to adjust the names in your model somewhere? Probably somewhere in the pre-block youโre making a reference to ฮท which is no longer in your random block. If you can share your code it will be easier to help you.
Best,
Yeah. Itโs my code, I just copy from Pumas tutorials and modify it.
Thank you for your consideration!
# Load library and creating dosage, population
using Pumas, DataFrames, Random, Tables
repeated_dose_regimen = DosageRegimen(100, time=0, ii=4, addl=2)
choose_covariates() = (Wt = rand(55:80),)
pop = Population(map(i -> Subject(id = i,
events = repeated_dose_regimen,
observations =(dv=Float64[],),
covariates = choose_covariates()),
1:4))
#
model = @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)
ฮฉ_IIV โ PDiagDomain([0.9, 0.07, 0.05])
ฮฉ_IOV โ PDiagDomain([0.9, 0.07, 0.05])
ฯ_prop โ RealDomain(lower = 0,init = 0.03)
end
@random begin
ฮท_IIV ~ MvNormal(ฮฉ_IIV)
ฮท_IOV ~ MvNormal(ฮฉ_IOV)
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
param = init_param(model)
vrandeffs = [sample_randeffs(model, param) for subject in pop]
# df1 = DataFrame(vrandeffs)
simulated_obs = [simobs(model, subject, param, randeffs, obstimes = 1:1:12) for (subject, randeffs) in zip(pop, vrandeffs)]
Right, so as I suspected you also need to adjust your names in the pre block. Please compare your posted snippet with the one below. Notice, to actually do the IOV you need to include a covariate to control which index to use.
# Load library and creating dosage, population
using Pumas, DataFrames, Random, Tables
repeated_dose_regimen = DosageRegimen(100, time=0, ii=4, addl=2)
choose_covariates() = (Wt = rand(55:80),)
pop = Population(map(i -> Subject(id = i,
events = repeated_dose_regimen,
observations =(dv=Float64[],),
covariates = choose_covariates()),
1:4))
#
model = @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)
ฮฉ_IIV โ RealDomain(lower=0.0)
ฮฉ_IOV โ PDiagDomain([0.07, 0.05])
ฯ_prop โ RealDomain(lower = 0,init = 0.03)
end
@random begin
ฮท_IIV ~ Normal(sqrt(ฮฉ_IIV))
ฮท_IOV ~ MvNormal(ฮฉ_IOV)
end
@covariates Wt
@pre begin
CL = cl * (Wt/70)^0.75 * exp(ฮท_IIV[1])
Vc = tv * (Wt/70) * exp(ฮท_IOV[1])
Ka = ka * exp(ฮท_IOV[2])
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
param = init_param(model)
vrandeffs = [sample_randeffs(model, param) for subject in pop]
# df1 = DataFrame(vrandeffs)
simulated_obs = [simobs(model, subject, param, randeffs, obstimes = 1:1:12) for (subject, randeffs) in zip(pop, vrandeffs)]My mistake for forgetting of changing ฮท in @pre block.
Thank you again for your help!