# Modelling PKPD parameters as a function of covariates

I have given a drug to 36 patients and sampled their blood at 8 time points to measure the concentration of that drug. We posit that the pharmacokinetic curve conditioned on parameters should look like

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

I have recorded each patient’s sex and would like to model parameters (V, k, and k_a) as a function of sex.

I would like to ignore the fact that the data come from different people from now as I am slowly building model complexity. Thusly, I would like to pretend that the data come from one male and one female patient.

I’m having some troubles with the modelling. The current Stan program I have written yields chains with Rhat much larger than 1. Here is my Stan model:

//adapted from Michael Betancourt.  Thank you
functions {
real C_anal(real t, real D, real V, real k_a, real k) {
return (D / V) * (k_a / (k_a - k))
* ( exp(- k * t) - exp(-k_a * t) );
}
}

data {
real t0;    // Initial time in days;
real C0[1]; // Initial concentration at t0 in mg/L

real D;   // Total dosage in mg

int<lower=1> N_subjects;
int<lower=0> sex[N_subjects];

int<lower=1> N_t;
real times[N_t];   // Measurement times in days

// Measured concentrations in effect compartment in mg/L
real C_hat[N_subjects,N_t];
}

parameters {

real beta_k_a[2];
real<lower=0> sigma_k_a;

real beta_k[2];
real<lower=0> sigma_k;

real beta_V[2];
real<lower=0> sigma_V;

real<lower=0> sigma;

}

transformed parameters{

real<lower=0> C[N_subjects, N_t];
real<lower=0> k_a;
real<lower=0> k;
real<lower=0> V;

for( n in 1:N_subjects){

k_a = exp(beta_k_a[1] + beta_k_a[2]*sex[n]);
k = exp(beta_k[1] + beta_k[2]*sex[n]);
V = exp(beta_V[1] + beta_V[2]*sex[n]);

for(t in 1:N_t)
C[n,t] = C_anal(times[t], D, V, k_a, k);
}
}

model {
// Priors
beta_k_a ~ normal(0,1);
beta_k ~ normal(0,1);
beta_V ~ normal(0,1);
sigma ~ cauchy(0,1);

// Likelihood
for (n in 1:N_subjects)
for (t in 1:N_t)
C_hat[n, t] ~ lognormal(log(C[n, t]), sigma);
}



Here is a link to the data. Shown below is some code to do the sampling.

library(tidyverse)
library(magrittr)
library(rstan)
options(mc.cores = parallel::detectCores())

# Need to rescale the concentration from ng/ML to mg/L since dose is in mg.
arrange(Subject,Time) %>%
filter(Time>0) %>%
mutate(Concentration_orig = Concentration,
Concentration = 10e-3*Concentration_orig) #resalced from ng/ml to mg/L

#Prepare the data to be passed into stan
C_hat = apixaban.data$Concentration C0 = array(c(0), dim = 1) t0 = 0 times = apixaban.data$Time
times = sort(unique(times))
N_t = length(times)
N_ut = length(utimes)

D = 2.5 #mg/L

subject.covariates = apixaban.data %>% select(Subject, Sex, Group) %>% distinct()

sex = subject.covariates$Sex %>% factor %>% as.numeric group = subject.covariates$Group %>% factor %>% as.numeric()

subject = subject.covariates\$Subject %>% factor %>% as.numeric

N_subjects= length(unique(subject))

dim(C_hat)<-c(N_subjects,N_t)


What is wrong with the model I have written? How can I model these parameters as if they had come from two individuals?

• Does it work for a single patient?
• Does it work for multiple patients who all have sex==0?
• Are all the units correct and do priors makes sense (I am too tired to check myself right now)?
• Running the stan program with fixed params and initials which makes sense - does this yield sound entries in C?

… but now I am seeing that you do not put a prior on all those sigma_*… this is probably a problem here since in addition you do not even associate with any data. So please put a prior here (otherwise Stan’s adaptation will go out of control, I think).

my 2 cents.