# Parametric bootstrap with Pumas?

Hello,

I know Pumas defaults to paired bootstrapped, but I’m curious whether it is it capable of doing parametric bootstrap of random effects and residuals, as noted in 2.3.3 here:

https://www.hal.inserm.fr/inserm-00939284/document

Pumas defaults to non parameteric bootstrap, not sure what you mean by paired.

Didn’t see the attachment you provided, but parametric bootstrap is essentially

1. Stimulate from model, `n` times
2. Re-estimate the simulated `n` datasets

Thank you. For the nonparametric bootstrap, is there a way to get the medians for the parameter estimates?

when making an `infer`ence we are mostly interested in the uncertainty around the estimated values and not an alternate estimator (`median` in your case) for that parameter. Hence, in Pumas, when the results are reported for a bootstrap inference, we only show the bootstrap-based percentile confidence intervals along with the maximum likelihood estimate based on the original fit.

However, the generated bootstrap object stores the individual parameters for each fit which you could use to derive a metric of your choice, as below.

``````bs = infer(fit1, Pumas.Bootstrap(samples=200))

# to get the individual estimates across all fits
DataFrame(bs.vcov)

# compute your median from there
``````

Hope that helps.

Great! Thank you.

I’m attempting to simulate trials (e.g., 1000x) with the uncertainty derived from the bootstrap estimates, similar to the ‘vcov’ method in ‘simobs’. Am I correct in that it’s more appropriate to simulate using the median parameter estimates from the nonparametric bootstrapping as the mean of the sampling distribution, over using the estimates from the original model fit?

Aha! If you want to simulate from the uncertainty distribution, you derived the `vcov` from the bootstrap results. In the next release of Pumas in 2 weeks (if I am not wrong) you can do this directly, just like you do `simobs` on a standard `infer`ence.

``````simobs(inferredbootstrapobject, nsamples=200)
``````
1 Like

Hello,

I see that the docs show this feature now for simulating with uncertainty from bootstrapping. It looks like I can also pass my own pumas population that is different from the inspect object. Is this correct? For example, I want to change the number of subjects for the simulation, so I need different populations than from the fitting.

Currently, you can also use `simobs` in at least four different ways:

• `simobs(fpm::FittedPumasModel)`: this simulates directly from the underlying population that you specified in `fpm = fit(my_model, my_pop, ...)`. It recovers the parameters from the estimated parameters from the fitted model.
• `simobs(fpm::FittedPumasModel, population::Population)`: this simulates using the fitted model `fpm` but for a different `population` that the model was initially fitted. It also recovers the parameters from the estimated parameters from the fitted model.
• `simobs(model::AbstractPumasModel, subject, params)`: this will simulate from a `@model` without being fitted against a certain `subject` with certain `params` as a `NamedTuple` of parameters values. This does not require a fitted Pumas model since you are passing the parameters.
• `simobs(model::AbstractPumasModel, population, params)`: same as the above, but now with a Population `population` instead of a single Subject.

Here is the documentation section link for a deep dive into simulation: Simulation of Pumas Models · Pumas.

@storopoli The documentation is full of such notation and I am trying to understanding how to use it. Are we supposed to use the following syntax for eg see below. I am sure this is all very simple to you.

1. simobs(fpm::myfit)? or simobs(fpm) and define fpm = myfit?
2. simobs(model::AbstractPumasModel). I get errors when I use this syntax. Perhaps I am making typographical errors of A, P or M. My syntax is: simobs(myfit::AbstractPumasModel). Is this correct?
B

@bobbrown these are `methods` of the function `simobs()`.
These are ways to call `simobs()` using different argument types.

For example:

``````simobs(fpm::FittedPumasModel)
``````

take as a single argument something called `fpm` that is of type `FittedPumasModel`.

So if I have a model like this:

``````my_model = @model begin
...
end
``````

And I estimate it using `fit()`:

``````my_fit = fit(
my_model, # model
my_pop,   # population,
init_params(my_model), # initial parameters
Pumas.FOCEI() # estimation method
)
``````

I can then use `my_fit` which is of type `FittedPumasModel` inside `simobs()`:

``````simobs(my_fit)
``````

Is that more clear to you?

@storopoli I see, ok I think so. But what are those two colons for then? “::” thats confusing syntax to me. Same with ‘AbstractPumasModel’ - my model is a simple compartmental PK model - what is abstract about it? Thank you.
Bob

The `::` is to denote the type of the argument.

For example:

``````round(x::Float64) = Int64(x)
round(x::Int64) = x
``````

The round function behaves differently if you provide the argument `x` as either `Float64` (floating point number, \mathbb{R}) or an `Int64` (integer number, \mathbb{Z}).

Thank you. @storopoli I tried your recommended syntax:
round(pkdata.TIME :: Float64)

obtained this error.
ERROR: LoadError: TypeError: in typeassert, expected Float64, got a value of type SentinelArrays.ChainedVector{Float64, Vector{Float64}}
Stacktrace:
[1] top-level scope

@storopoli I am getting an error while trying to simulate from an `infer` object.

`test = simobs(abbrev_pl_BE_sims_ML_reest_infers[1], full_pl_BE_pops[1], obstimes=obstimes_SD, ensemblealg = EnsembleThreads(), samples=100)`

ERROR: MethodError: no method matching simobs(::Pumas.FittedPumasModelInference{Pumas.FittedPumasModel{PumasModel{ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_add, :σ_prop), NTuple{8, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#41#57", var"#42#58", var"#44#60", var"#46#62", Pumas.LinearODE, var"#47#63", var"#52#68", ModelingToolkit.ODESystem}, Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:ab_rep, :params_num, :ab_NSubj, :occ, :period, :sequence, :seq_n, :formulation, :isT, :input_tr_ratio), Tuple{Int64, Int64, Int64, Int64, Int64, String, Int64, String, Int64, Float64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Optim.MultivariateOptimizationResults{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, Float64, Vector{Float64}, Float64, Float64, Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, Bool, NamedTuple{(:f_limit_reached, :g_limit_reached, :h_limit_reached, :time_limit, :callback, :f_increased), NTuple{6, Bool}}}, Pumas.FOCE, Vector{Vector{Float64}}, NamedTuple{(:optimize_fn, :constantcoef, :omegas, :ensemblealg, :diffeq_options), Tuple{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, NamedTuple{(), Tuple{}}, Tuple{}, EnsembleThreads, NamedTuple{(), Tuple{}}}}, ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_add, :σ_prop), NTuple{8, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}}, Pumas.Bootstraps, Float64}, ::Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:params_num, :fl_NSubj, :occ, :period, :sequence, :seq_n, :formulation, :isT), Tuple{Int64, Int64, Int64, Int64, String, Int64, String, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}; obstimes=[0.001, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0, 12.0, 14.0 … 85.0, 90.0, 104.0, 118.0, 132.0, 146.0, 160.0, 190.0, 220.0, 250.0], ensemblealg=EnsembleThreads(), samples=100)
Closest candidates are:
simobs(::PumasModel, ::AbstractVector{T} where T<:Subject, ::NamedTuple, ::NamedTuple) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/HDuXQ/src/estimation/inference.jl:635 got unsupported keyword arguments “obstimes”, “ensemblealg”, “samples”
simobs(::PumasModel, ::AbstractVector{T} where T<:Subject, ::NamedTuple, ::NamedTuple, ::Union{Nothing, AbstractVector{var"#s652"} where var"#s652"<:NamedTuple}; samples, rng, diffeq_options) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/HDuXQ/src/estimation/inference.jl:635 got unsupported keyword arguments “obstimes”, “ensemblealg”
simobs(::PumasModel, ::AbstractVector{T} where T<:Subject, ::MvNormal) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/HDuXQ/src/estimation/inference.jl:649 got unsupported keyword arguments “obstimes”, “ensemblealg”, “samples”

Stacktrace:
[1] top-level scope
@ REPL[113]:1

Correct. You can do:

``````pmi = infer(res, Pumas.Bootstrap(samples = 20))
simobs(pmi, new_population, samples = 10)
``````

where `res` is the result from `fit` and `new_population` is a vector of new subjects.

Hi Mohamed,

I am still getting an error with that method. Here is the error even without a `population`:

``````julia> test = simobs(abbrev_pl_BE_sims_ML_reest_infers[1], samples=100)
ERROR: MethodError: no method matching simobs(::Pumas.FittedPumasModelInference{Pumas.FittedPumasModel{PumasModel{ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_add, :σ_prop), NTuple{8, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#41#57", var"#42#58", var"#44#60", var"#46#62", Pumas.LinearODE, var"#47#63", var"#52#68", ModelingToolkit.ODESystem}, Vector{Subject{NamedTuple{(:dv,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:ab_rep, :params_num, :ab_NSubj, :occ, :period, :sequence, :seq_n, :formulation, :isT, :input_tr_ratio), Tuple{Int64, Int64, Int64, Int64, Int64, String, Int64, String, Int64, Float64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, Optim.MultivariateOptimizationResults{Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}, Float64, Vector{Float64}, Float64, Float64, Vector{Optim.OptimizationState{Float64, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Nothing, Float64, Optim.Flat}}}, Bool, NamedTuple{(:f_limit_reached, :g_limit_reached, :h_limit_reached, :time_limit, :callback, :f_increased), NTuple{6, Bool}}}, Pumas.FOCE, Vector{Vector{Float64}}, NamedTuple{(:optimize_fn, :constantcoef, :omegas, :ensemblealg, :diffeq_options), Tuple{Pumas.DefaultOptimizeFN{Nothing, NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases), Tuple{Bool, Bool, Bool, Float64, Bool}}}, NamedTuple{(), Tuple{}}, Tuple{}, EnsembleThreads, NamedTuple{(), Tuple{}}}}, ParamSet{NamedTuple{(:tvtr_ratio, :tvCL, :tvVc, :tvDuration, :ωCL_BSV, :ωVc_BSV, :σ_add, :σ_prop), NTuple{8, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}}, Pumas.Bootstraps, Float64}; samples=100)
Closest candidates are:
simobs(::PumasModel, ::DataFrame, ::Any...; kwargs...) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/HDuXQ/src/models/model_api.jl:1143
simobs(::PumasModel, ::AbstractVector{T} where T<:Subject, ::NamedTuple, ::NamedTuple) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/HDuXQ/src/estimation/inference.jl:635 got unsupported keyword argument "samples"
simobs(::PumasModel, ::AbstractVector{T} where T<:Subject, ::NamedTuple, ::NamedTuple, ::Union{Nothing, AbstractVector{var"#s652"} where var"#s652"<:NamedTuple}; samples, rng, diffeq_options) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/HDuXQ/src/estimation/inference.jl:635
...
Stacktrace:
[1] top-level scope
@ REPL[121]:1
``````

which version of Pumas are you using?

I’m using what’s available on umb.juliahub.com.

This feature is only available in Pumas 2.1 which I think you don’t have access to yet. @vijay can comment on when it will be available to you.