I’ve been working on a project trying to model the PK of ertapenem in intermittent hemodialysis patients. The two issues I ran into was: 1) Simulating inter-occasional variability and 2) simulating a population where IOV is expected to occur at certain time points within the data. I’ve included a snapshot of my PUMAS code and the dataset I’ve been working with

Brief background: I have data for ertapenem concentrations during dialysis and off of dialysis. Each participant goes through 4-5 cycles of dialysis (every 48 hours for 6 hours) and we wanted to incorporate IOV on CL for each dialysis cycle.

Please take a look at Defining NLME models in Pumas · Pumas. There are several examples demonstrating between/inter occasion variability. Hopefully that will help you. Please reply here if it requires further explanation.

Please also see this post for more information about IOV in Pumas:

For your data set, you will want to have a single OCC column with 1-4/5 for each cycle of dialysis.

Regarding your 2), if I understand you correctly, you may want to do something like the choose_covariates function here: https://tutorials.pumas.ai/html/simulation/simulating_populations.html . So, you’ll determine each subject’s OCC column based on whatever criteria you choose for that given subject, and then pass OCC as a covariate.

Thanks everyone for your advice. I’ve been looking at all the resources you have provided and tried a couple of methods.

I changed my OCC to a single column to have 1 -6 different occasions

CL = tvcl * exp(tvdialDIAL) * exp(η[1]+DIALη_occ[OCC]) ; this seems to work for me where DIAL is a covariate, and when it is present (i.e. =1), add IOV.
But, I found two different methods to define η_occ in the @random block and both gives different answers.
η_occ ~ MvNormal(Diagonal(fill(Ω_occ,6))) vs.
η_occ ~ MvNormal([Ω_occ, Ω_occ, Ω_occ, Ω_occ, Ω_occ,Ω_occ]).
I was wonder if anyone why this is the case?

The difference is a little subtle. In the first version, you pass a matrix, which the MvNornal constructor interprets as a covariance matrix. In the second version, you pass a vector which the MvNormal constructor interprets as a vector of standard deviations. Compare

i.e. the first argument is the dimension and the second is the standard deviation. Personally, I also prefer to use lower case greek letters, i.e. ω, for standard deviations, squared versions, i.e. ω², for variances and only use the upper case version for the whole covariance matrix. I think that makes it easier to read.

Thanks @donaldlee3 ,
Thanks for sharing the reference. Based on the provided example, you would only be assigning a specific covariate for each simulated subjected (same covariate value at all timepoints for the subject). I was wondering if there was a way to specify OCC to be present at specific time points for each individual (i.e. subject 1: OCC = 0 at t = 0-18h, OCC = 1 at t = 19-24h, OCC = 0 at t = 25-42h, and OCC =2 at t = 43-48h ).
Many thanks!

One way to do this is to convert your Pumas population into a data frame, and then add your OCC column that you made manually. You can then use read_pumas to convert the dataframe back to a Pumas population.

For example:

pop_df = DataFrame(Population(map(j -> Subject(id=j, events=ev), 1:Nsubj)))
pop_df[!, :OCC] .= pre_made_occs
pop = read_pumas(pop_df,
id = :id,
time = :time,
amt = :amt,
cmt = :cmt,
observations = [:dv],
covariates = [:OCC],
covariates_direction = :left,
evid = :evid)