Parametric bootstrap with Pumas?

How-to
1

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

2

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
  3. Get your parameter’s distribution
3

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

4

when making an inference 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.

5

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?

6

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 inference.

simobs(inferredbootstrapobject, nsamples=200)

cc @mohamed82008

1 Like
7

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.

8

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.

9

@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
10

@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?

11

@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

12

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}).

13

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

14

@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

15

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.

16

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
17

which version of Pumas are you using?

18

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

19

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.