Simulation of Pumas Models

The simobs Function

Simulation of a PumasModel are performed via the simobs function. The function is given by the values:

simobs(m, data, param, [randeffs]; kwargs...)

The terms in the function call are:

  • m: the PumasModel, either defined via the @model DSL or the function-based interface.
  • data: either a Subject or a Population.
  • param: a NamedTuple of parameters which conform to the ParamSet of the model.
  • randeffs: an optional argument for the random effects for the simulation. If the random effects are not given, they are sampled as described in the model.
  • kwargs: extra keyword arguments.

Additionally, the following keyword arguments can be used:

The result of simobs function is a SimulatedObservation if the data was Subject and a SimulatedPopulation if the data was a Population.

Handling Simulated Returns

When running

sim = simobs(m, data, param)

sim is a SimulatedObservation which can be accessed via its fields. These fields are:

  • subject: the Subject used to generate the observation
  • time: the times associated with the observations
  • observations: the resulting observations of the simulation

If the @model DSL is used, then observed is a NamedTuple where the names give the associated values. From the function-based interface, observed is the chosen return type of the observed function in the model specification.

A SimulatedPopulation is a collection of SimulatedObservations, and when indexed like sim[i] it returns the SimulatedObservation of the ith simulation subject.

Visualizing Simulated Returns

These objects have automatic plotting and dataframe visualization. To plot a simulation return, simply call plot on the output using Plots.jl. For example, the following will run a simulation and plot the observed variables:

obs = simobs(m, data, param)
using Plots
plot(obs)

By default this generates a plot for each derived variable. To choose which variables to plot, the obsnames argument can be given which declares indices or derived variable names to plot. For example, plot(obs,obsnames=[:dv1,:dv2]) would only plot the values dv1 and dv2. In addition, all of the Plots.jl attributes can be used in this plot command. For more information on using Plots.jl, please see the Plots.jl tutorial. Note that if the simulated return is a SimulatedPopulation, then the plots overlay the results of the various subjects.

To generate the DataFrame associated with the observed outputs, simply call DataFrame on the simulated return. For example, the following builds the tabular output from the returned object:

obs = simobs(m, data, param)
df = DataFrame(obs)

Different usage patterns of simobs

Below, we showcase the different usage patterns of simobs. Consider the following model, population and parameters that we will use to simulate from.

onecomp = @model begin
  @param begin
    tvcl ∈ RealDomain(lower=0, init=4)
    tvvc ∈ RealDomain(lower=0, init=70)
    Ω    ∈ PDiagDomain(init=[0.04,0.04])
    σ    ∈ RealDomain(lower=0.00001, init = 0.1)
  end
  @random begin
    η ~ MvNormal(Ω)
  end
  @pre begin
    CL = tvcl * exp(η[1])
    Vc = tvvc * exp(η[2])
  end
  @dynamics Central1
  @derived begin
    cp = @. Central/Vc
    dv ~ @. Normal(cp, cp*σ)
  end
end

model_params = init_param(onecomp)

ext_param = (tvcl = 3, tvvc = 60, Ω = Diagonal([0.04, 0.04]), σ = 0.1)

regimen = DosageRegimen(100)
pop = map(i -> Subject(id=i, events = regimen), 1:3)

NOTE: In all the examples below, the default simulation time is till 24 hours.

Subject level simulation

  1. Subject level simulation

In this pattern, you can use simobs to do a single subject simulation using the syntax below.

simobs(model::PumasModel, subject::Subject, param::NamedTuple)

Particularly, for this example, it translates to

sims2 = simobs(onecomp, pop[1], ext_param)

where, onecomp is the model, pop[1] represents one subject in the Population, pop. The parameters of the model are passed in as the NamedTuple defined earlier.

  1. Subject level simulation with passing random effects

In this pattern, you can use simobs to do a single subject simulation using the syntax below.

simobs(model::PumasModel, subject::Subject, param::NamedTuple, randeffs::Union{Nothing, NamedTuple})

Here individual subject random effects can be passed in as vector. The length of the random effect vector should be the same as the number of parameters for each subject.

ext_randeffs = (η = [0.7, -0.44],)

Particularly, for this example, it translates to

sims3 = simobs(onecomp, pop[1], ext_param, ext_randeffs)

where, onecomp is the model, pop[1] represents one subject in the Population, pop. The parameters of the model are passed in as the NamedTuple defined earlier. ext_randeffs are the specific values of the η's that are passed in, which means that simobs uses these values as opposed to those sampled from the distributions specified in the @random block.

  1. Subject level simulation with passing in an array of parameters

In this pattern, you can use simobs to do a single subject simulation using the syntax below.

simobs(model::PumasModel, subject::Subject, vparam::AbstractArray{<:NamedTuple,1})

Here, a subject can be simulated with an array of parameters that generates a unique simulation solution for each element of the parameter array. E.g. let's simulate pop[1] with three different combinations of tvcl and tvvc. In the example below, we change the value of tvcl to represent lo, med, hi of 3, 6, 9 but we don't change the other values (Note that are more efficient programmatic ways to generate such arrays)

ext_param_array = [(tvcl = 3, tvvc = 60, Ω = Diagonal([0.04, 0.04]), σ = 0.1),
                   (tvcl = 6, tvvc = 60, Ω = Diagonal([0.04, 0.04]), σ = 0.1),
                   (tvcl = 9, tvvc = 60, Ω = Diagonal([0.04, 0.04]), σ = 0.1)]

Particularly, for this example, it translates to

sims4 = simobs(onecomp, pop[1], ext_param_array)

where, onecomp is the model, pop[1] represents one subject in the Population, pop. The parameters of the model are passed in as the vector of NamedTuples ext_param_array.

Tip

A DataFrame of parameter values can be converted to a NamedTuple easily using Tables.rowtable syntax. See the example below

julia> myparams = DataFrame(
    tvcl = [1, 2, 3],
    tvvc =[20, 30, 40],
    Ω = [Diagonal([0.04, 0.04]), Diagonal([0.04, 0.04]), Diagonal([0.04, 0.04])],
    σ = 0.1)

3×4 DataFrame
│ Row │ tvcl  │ tvvc  │ Ω                    │ σ       │
│     │ Int64 │ Int64 │ Diagonal…            │ Float64 │
├─────┼───────┼───────┼──────────────────────┼─────────┤
│ 1   │ 1     │ 20    │ [0.04 0.0; 0.0 0.04] │ 0.1     │
│ 2   │ 2     │ 30    │ [0.04 0.0; 0.0 0.04] │ 0.1     │
│ 3   │ 3     │ 40    │ [0.04 0.0; 0.0 0.04] │ 0.1     │

This can be converted to a NamedTuple

julia> using Tables

julia> myparam_tuple = Tables.rowtable(myparams)
3-element Array{NamedTuple{(:tvcl, :tvvc, :Ω, :σ),Tuple{Int64,Int64,Diagonal{Float64,Array{Float64,1}},Float64}},1}:
 (tvcl = 1, tvvc = 20, Ω = [0.04 0.0; 0.0 0.04], σ = 0.1)
 (tvcl = 2, tvvc = 30, Ω = [0.04 0.0; 0.0 0.04], σ = 0.1)
 (tvcl = 3, tvvc = 40, Ω = [0.04 0.0; 0.0 0.04], σ = 0.1)

myparam_tuple is of the same form as the ext_param_array.

  1. Subject level simulation with passing in an array of parameters and random effects

The two usage patterns above can be combined such that one can pass in an array of parameters and corresponding random effects with the usage pattern below

simobs(model::PumasModel, subject::Subject, vparam::AbstractArray{<:NamedTuple,1}, vrandeffs::Union{Nothing, AbstractArray{<:NamedTuple,1}})
ext_randeffs = [(η = [0.7, -0.44],),
                (η = [0.6, -0.9],),
                (η = [-0.7, 0.6],)]
sims5 = simobs(onecomp, pop[1], ext_param_array, ext_randeffs)

Population level simulation

  1. Population level simulation with passing in parameters

In this pattern, you can use simobs to do a population simulation using the syntax below.

simobs(model::PumasModel, population::AbstractArray{<:Subject,1}, param::NamedTuple)

Particularly, for this example, it translates to

sims7 = simobs(onecomp, pop, ext_param)

where, onecomp is the model, pop represents the Population, pop. The parameters of the model are passed in as the NamedTuple defined earlier.

  1. Population level simulation passing random effects

In this pattern, you can use simobs to do a population simulation using the syntax below.

simobs(model::PumasModel, population::AbstractArray{<:Subject,1}, param::NamedTuple, vrandeffs::Union{Nothing, AbstractArray{<:NamedTuple,1}})

Here individual subject random effects can be passed in as vector. The length of the random effect vector should be the same as the number of parameters for each subject.

Tip

You can use this method to pass in the empirical_bayes() of a FittedPumasModel

ext_randeffs =  [(η = rand(2),),
                 (η = rand(2),),
                 (η = rand(2),)]

Particularly, for this example, it translates to

sims8 = simobs(onecomp, pop, ext_param, ext_randeffs)

where, onecomp is the model, pop represents the Population, pop. The parameters of the model are explicitly passed in as the NamedTuple defined earlier. ext_randeffs are the specific values of the η's that are passed in, which means that simobs uses these values as opposed to those sampled from the distributions specified in the @random block.

  1. Population level simulation with passing in an array of parameters

In this pattern, you can use simobs to do a population simulation using the syntax below.

simobs(model::PumasModel, population::AbstractArray{<:Subject,1}, vparam::AbstractArray{<:NamedTuple,1})

Here, a population can be simulated with an array of parameters that generates a unique simulation solution for each element of the parameter array. E.g. let's simulate pop with three different combinations of tvcl and tvvc. In the example below, we change the value of tvcl to represent lo, med, hi of 3, 6, 9 but we don't change the other values (Note that are more efficient programmatic ways to generate such arrays)

ext_param_array = [(tvcl = 3, tvvc = 60, Ω = Diagonal([0.04, 0.04]), σ = 0.1),
                   (tvcl = 6, tvvc = 60, Ω = Diagonal([0.04, 0.04]), σ = 0.1),
                   (tvcl = 9, tvvc = 60, Ω = Diagonal([0.04, 0.04]), σ = 0.1)]

Particularly, for this example, it translates to

sims9 = simobs(onecomp, pop, ext_param_array)

where, onecomp is the model, pop represents the Population, pop. The parameters of the model are passed in as a vector of the NamedTuples ext_param_array as defined earlier represents the array of parameter values we want to simulate the model at.

  1. Population level simulation with passing in an array of parameters and random effects

The two usage patterns above can be combined such that one can pass in an array of parameters and corresponding random effects with the usage pattern below

simobs(model::PumasModel, population::AbstractArray{T,1} where T<:Subject, vparam::AbstractArray{T,1} where T, vrandeffs::Union{Nothing, AbstractArray{T,1} where T}, args...; rng, kwargs...)
ext_randeffs_array = [(η = rand(2),) for i in 1:length(ext_param_array)*length(pop)]
sims10 = simobs(onecomp, pop, ext_param_array, ext_randeffs_array)
  1. Simulation using FittedPumasModel

The usage pattern is represented below

simobs(fpm::Pumas.FittedPumasModel)

Any model that has been fitted results in a FittedPumasModel. simobs can simulate from the object where it utilizes the final parameter estimates of the fit and the design of the data that was fitted.