I just started playing around with Pumas, so it may very well be that I didn’t find the right convenience function.
I am wondering whether there is a way for NLME models to also get estimates or chain traces for the individual subject parameters. For example for the tutorial model
mymodel = @model begin
@param begin
tvcl ∈ RealDomain(lower=0, init = 1.0)
tvv ∈ RealDomain(lower=0, init = 20)
tvka ∈ RealDomain(lower = 0, init= 1)
Ω ∈ PDiagDomain(init=[0.09,0.09, 0.09])
σ_prop ∈ RealDomain(lower=0,init=0.04)
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
CL = tvcl * (Wt/70)^0.75 * exp(η[1])
Vc = tvv * (Wt/70) * exp(η[2])
Ka = tvka * exp(η[3])
end
@covariates Wt
@dynamics Depots1Central1
#@dynamics begin
# Depot' = -Ka*Depot
# Central' = Ka*Depot - (CL/V)*Central
#end
@derived begin
cp = @. 1000*(Central / Vc)
dv ~ @. Normal(cp, sqrt(cp^2*σ_prop))
end
end
You can sample from the posterior doing something like this
However the Pumas.Chains method only returns the traces for the population parameters. I am interested in the individual η for the Subjects…
Looking into the result.chain it seems that the traces do exist, but it’s not obvious how the traces map to the parameters. Is there any way to clarify that for me or is there maybe even a convenience method to extract also the individual parameters (I was going to try the param_map in the Pumas.Chains, but I am not sure how the Subject ηs would be called)?
Thanks a lot in advance! (Pumas is really great to work with btw!)
Maybe to make it easier to reproduce chains, here is a working code snippet
using Pumas, Plots, CSV, Random
Random.seed!(0)
# Define regimen
single_dose_regimen = DosageRegimen(100)
# Define covariate generator
choose_covariates() = (Wt = rand(55:80),)
# Define Population
pop = Population(
map(i -> Subject(
id = i,
events=single_dose_regimen,
covariates=choose_covariates()),
1:12))
# Define model
mymodel = @model begin
@param begin
tvcl ∈ RealDomain(lower=0, init = 1.0)
tvv ∈ RealDomain(lower=0, init = 20)
tvka ∈ RealDomain(lower = 0, init= 1)
Ω ∈ PDiagDomain(init=[0.09,0.09, 0.09])
σ_prop ∈ RealDomain(lower=0,init=0.04)
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
CL = tvcl * (Wt/70)^0.75 * exp(η[1])
Vc = tvv * (Wt/70) * exp(η[2])
Ka = tvka * exp(η[3])
end
@covariates Wt
@dynamics Depots1Central1
@derived begin
cp := @. 1000*(Central / Vc)
dv ~ @. Normal(cp, sqrt(cp^2*σ_prop))
end
end
# Simulate Population
param = init_param(mymodel)
obs = simobs(mymodel, pop, param, obstimes=0:1:72)
plot(obs)
# Convert to dataset
simdf = DataFrame(obs)
simdf[!, :route] .= "ev"
simdf.cmt = ifelse.(ismissing.(simdf.cmt), 2, simdf.cmt)
est_df = simdf[.!((simdf.dv .== 0.0) .& (simdf.cmt .==2)),:]
data = read_pumas(est_df, covariates=[:Wt], observations=[:dv])
# Use NUTS to infer parameters
result = fit(mymodel, data, param, Pumas.BayesMCMC();
nsamples=2000, nadapts=1000)
# This returns the traces of the population parameters
chains = Pumas.Chains(result)
# Here I can see 43 traces which presumably correspond to the 7 population
# parameters plus the 12 * 3 Subject random effects for CL, Vc and Ka.
result.chain
Thanks for your kind words. It’s true that the values are all in there but we don’t extract them in the Chains constructor by default. This feature is on our list for the next release though.
For now if you poke around the result.chain struct, you can extract the values. Your assumption is correct in that the first 7 are the population parameters and the remaining are the subjects’ random effects. More generally, if nparams is the number of population parameters (7 above) and nrandeffs is the number of random effects per subject (3 above), then the j^{th} sample of the i^{th} subject’s random effects sits in:
