Sometimes it is useful to have some random effects correlated and others uncorrelated with each other. In Pumas, this is possible by specifying separate random effects for those that are correlated and those that are not correlated. In the below example, we have that the clearance and volume of distribution parameters are correlated through the first and second component of the
η_vc_cl random effect, while the other random effect on the rate of absorption
η_ka is uncorrelated with the other two.
model = @model begin @param begin tvvc ∈ RealDomain(lower = 0) tvcl ∈ RealDomain(lower = 0) tvka ∈ RealDomain(lower = 0) Ωvc_cl ∈ PSDDomain(2) ωka ∈ RealDomain(lower = 1e-5) σ_pk ∈ RealDomain(lower = 1e-5) end @random begin η_vc_cl ~ MvNormal(Ωvc_cl) η_ka ~ Normal(0, ωka) end @pre begin Ka = tvka*exp(η_ka) CL = tvcl*exp(η_vc_cl) Vc = tvvc*exp(η_vc_cl) end @dynamics Depots1Central1 @derived begin conc_pk ~ @. LogNormal(log(Central/Vc), σ_pk) end end doses = DosageRegimen(10, addl = 3, ii = 12) raw_subject = Subject(id="1", events = doses) param = (tvvc = 2.1, tvcl = 2.1, tvka = 0.9, Ωvc_cl = [0.4 0.2; 0.2 0.4], ωka = 0.3, σ_pk = 0.2) subject = Subject(simobs(model, raw_subject, param; obstimes = 0:0.5:72)) using PumasUtilities sim_plot(subject)
The simulation gives a trajectory as follows
To get a confirmation that the two components are indeed correlated, we can simulate many subjects as follows
raw_population = [Subject(id=string(id), events = doses) for id = 1:1000] sim_population = simobs(model, raw_population, param; obstimes = 0:0.5:72) population = Subject.(sim_population) sim_plot(population)
and this gives trajectories as follows
Then, we can look at their tabular representation using the
sim_population_df = DataFrame(sim_population) # First the first row of each subject sim_popuplation_df_first = unique(sim_population_df, [:id])
and from this we can look at correlations between the different random effects using an AlgebraOfGraphics plot
draw(data(sim_population_df) * (visual(AlgebraOfGraphics.Scatter) + AlgebraOfGraphics.linear()*visual(color=:red) * mapping([:η_vc_cl_1, :η_vc_cl_2,:η_ka_1 ], [:η_vc_cl_1 :η_vc_cl_2 :η_ka_1 ], col = dims(1), row = dims(2)))
and we see that the correlations are as expected.
This example had a pair with correlations and a third random effect that was uncorrelated with the other two, but the scalar random effect for rate of absorption that had a
Normal distribution could just as well have been multivariate, or we could have had ever more random effects each with their own correlation structure between the elements such as the following pseudo-example
@model begin @param begin [...] Ω1 ∈ PSDDomain(2) Ω2 ∈ PSDDomain(4) Ω3 ∈ PSDDomain(3) [...] end @random begin η1 ~ MvNormal(Ω1) η2 ~ MvNormal(Ω2) η3 ~ MvNormal(Ω3) end [...] end
Hope this helps!