# How to handle several random effects with pairwise correlation between some

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[1])
Vc = tvvc*exp(η_vc_cl[2])
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 `DataFrame` constructor

``````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!

2 Likes