# Simultaneously estimating baseline and treatment parameters: hierarchical ODE modelling

I have a hierarchical modelling puzzle that I would like your collective wisdom on.

For clarity, my model is an ODE model of PK/PD dynamics tracking the both the pathogen dynamics (both in the absence and presence of treatment) and drug dynamics. For this example, let N be a set of subjects: N_{untreated} and N_{treated} are sets of untreated and treated subjects, respectively. Similarly, \theta is a set of parameters, \theta_{baseline} is a set of baseline parameters that is shared among untreated and treated subjects, and \theta_{treated} contains PK/PD parameters that are only informed by N_{treated}. I’m then interested in estimating the mean (the entire parameter set, \theta) and subject-level effects (including correlations among subject-level effects).

Based on the above premises, I present my proposed model structure below. But I am afraid there is some issues with subject-level effects in this formulation.

    u = diag_pre_multiply(sd_u, L_u) * z_u;

for(j in 1:size_N) {
theta[j][1] = theta_baseline_1 + u[1,j];
theta[j][2] = theta_baseline_2 + u[2,j];
theta[j][3] = theta_treated_1  + u[3,j];
theta[j][4] = theta_treated_2  + u[4,j];
}


where L_u ~ lkj_corr_cholesky(5), sd_u ~ normal(0,1) and to_vector(z_u) ~ normal(0,1). Specifically the issue is that the model would estimate  u[3,j] and u[4,j]  for N_{untreated}, which seems nonsensical. But in order to keep the dimension of  u as \theta \times \theta, I felt forced into this formulation.

Do you think I can just ignore  u[3,j] and u[4,j]  from the output for N_{untreated}? Or, is this formulation going to cause a problem for estimation of subject-level parameter correlations, L_u.

I am wondering an alternative might be a two-step fitting: first to fit the baseline model describing only the pathogen dynamics to N_{untreated} and plug in the posteriors as priors in the second stage fitting of N_{treated}. However, my concern is that there may be interactions (or subject-level correlations) between M_{baseline} and M_{treated} so that all effects should be estimated simultaneously. Also it would be a problem for interpretation if M_{baseline} estimated from N_{treated} differs from that estimated from N_{untreated} due to potential interactions.

Sorry, don’t have to delve into this now, but @wds15 is quite knowledgeable about PK/PD models and hopefully isn’t too busy.

Thank you, Martin, for the reply. Hopefully, @wds15 has some time to look at this post.

In principle there is no issue with your formulation. I presume that the theta 3 & 4 don’t endup in some likelihood term, right? So that’s fine - you will for these subjects just sample the hierarchical prior and you can ignore these estimates (they will anyway be the same for all sibjects).

What does bother me though, is that you would usually express the treated in relation to untreated ones (and essentially only fit a delta to the untreated ones for the treated) - in case this is a randomized treatment assignment then this would be the way to go if you want to learn about the treatment effect.

I hope this helps a little bit to move on.

(this is not a PK/PD question per se so far)

1 Like

I know this post is a little old now, but I have noticed that subject-level effects for baseline parameters (i.e., u[1,j] and u[2,j]) get biased with respect to treatment, such that “mixed effects” (i.e., theta[j][1] and theta[j][2]) are not independent of treatment. What I think is happening is that samples, u[1,j] and u[2,j], are compensating for the lack of flexibility in theta_baseline_1 and theta_baseline_2. But, I have no reason to believe that theta_baseline_1 and theta_baseline_2 should be different between treatments, and I would like theta[j][1] and theta[j][2] to be, on average, independent of treatments.
In the attached figure, I am showing the aggregate theta[j][1] for control (no outline) versus treated subjects (outlined). You see that theta[j][1] for treated subjects are higher on average than control subjects.