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
: thePumasModel
, either defined via the@model
DSL or the function-based interface.data
: either aSubject
or aPopulation
.param
: aNamedTuple
of parameters which conform to theParamSet
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:
alg
: the type for which differential equation solver method to use. For example,alg=Rodas5()
specifies the usage of the 5th order Rosenbrock method for ODEs described in the DifferentialEquations.jl solver documentation page. Defaults to an automatic stiffness detection algorithm for ODEs.ensemblealg
: the parallel algorithm to use internally for simulating aPopulation
. The options are derived from DifferentialEquations.jl. The default isEnsembleThreads()
.- Any keyword argument in the DifferentialEquations.jl common solver arguments. These are documented on the DifferentialEquations.jl common solver options page.
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
: theSubject
used to generate the observationtime
: the times associated with the observationsobservations
: 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 SimulatedObservation
s, and when indexed like sim[i]
it returns the SimulatedObservation
of the i
th 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
- 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.
- 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.
- 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 NamedTuple
s ext_param_array
.
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
.
- 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
- 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.
- 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.
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.
- 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 NamedTuple
s ext_param_array
as defined earlier represents the array of parameter values we want to simulate the model at.
- 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)
- 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.