Drug concentration as a function of sex

Questions
1

I’m using Pymc3 to fit a hierarchical model for concentration of a drug. The model should include sex as a covariate.

Our model posits that concentration of drug, conditioned on parameters, should look like…

C(t) = \dfrac{2.5}{V} \dfrac{k}{k-k_a}\left( e^{-kt} - e^{k_at}\right).

Here, V,k, k_a are all parameters which are modeled heirarchically. I would like to model the mean of each of these parameters as a function of sex.

The data I have consist of 36 patients sampled at 8 time points. I am currently ignoring the fact that measurements come from different patients (that will come at a later time) and instead imagining I have 2 patients (one male, one female) who have provided these measurements. Here is my data.

Here is how I construct my model

df = pd.read_csv('Data/ApixibanExperimentData.csv').query('Time>0')

#Rescaled to mg/L
Time= df.Time.values
Concentration_hat = df.Concentration.values*10e-3

sample_times = np.sort(np.unique(Time))

Sex = df.Sex.apply(lambda x: 1 if x=='Male' else 0).values
Sex = np.c_[np.ones_like(Sex),Sex]


with pm.Model() as model1:
    
    D = 2.5 #Dose in mg
    
    #Sex Effects
    beta_k_a_sex = pm.Normal('ß_k_a_sex', 0,1, shape = 2)
    beta_k_sex = pm.Normal('ß_k_sex', 0,1, shape = 2)
    beta_V_sex = pm.Normal('V_sex', 0,1, shape = 2)
    
    #PKPD parameters
    k_a = pm.TruncatedNormal('k_a',pm.math.exp(pm.math.dot(Sex,beta_k_a_sex)),1, lower = 0)
    k = pm.TruncatedNormal('k',pm.math.exp(pm.math.dot(Sex,beta_k_sex)),1, lower = 0)
    
    #Volume of the absorptive compartment
    V = pm.TruncatedNormal('V',pm.math.exp(pm.math.dot(Sex,beta_V_sex)),1,lower = 0)
    sigma = pm.HalfCauchy('sigma',1)
    
    #Analytic solution to the ODE
    solution = (D/V)*(k_a/(k-k_a))*(tt.exp(-k_a*Time) - tt.exp(-k*Time))
    #Convienience
    solution_samp = (D/V)*(k_a/(k-k_a))*(tt.exp(-k_a*sample_times) - tt.exp(-k*sample_times))
    
    C = pm.Deterministic('C',solution)
    Cs = pm.Deterministic('Cs', solution_samp)
        
    Y = pm.Lognormal('Y', mu = tt.log(C), sd = sigma, observed = Concentration_hat )
    
    trace = pm.sample(2000, tune = 500, chains = 4)
    ppc = pm.sample_ppc(trace)

Problem

After sampling, the male and female curves look exactly the same, which leads me to believe I have improperly specified my model.

Have I set up my model correctly? Will solution give me two curves, depending on the sex of the patient? If not, how can I correct my model?

2

Thanks @dpananos for the detailed description and code!

I think the problem with the two identical sex curves lies in the way you combine Sex and the ODE parameters. I would do as follows, and this could be done equally for patients:

(1) instead of apply/lambda, which is correct, you could use pandas.Categorical to get a vector of sex indices:

df['Sex'] = pd.Categorical(df['Sex'].values)
sex_index = df['Sex'].cat.codes.values
print (sex_index)

(2) Then, you should index the hyperpriors (instead of using dot product):

k_a = pm.TruncatedNormal('k_a',pm.math.exp(beta_k_a_sex[sex_index]),1, lower = 0)
k = pm.TruncatedNormal('k',pm.math.exp(beta_k_sex[sex_index]),1, lower = 0)
V = pm.TruncatedNormal('V',pm.math.exp(beta_V_sex[sex_index]),1,lower = 0)

… because each observation is either female or male. Through the Sex matrix structure and then the dot product, you got two columns, and each was inferred via the same observations. Now, you index, which means a vector is assembled in which either the female version of the hyperparameter or the male one are used, depending on the data.

This still retrieves eight distributions for C, which I think is meaningless. This is the eight time points. There is no need to look at C by creating a deterministic.

Cheers,

Falk