Modelling PKPD parameters as a function of covariates

Modeling
1

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.
apixaban.data = read_csv('Data/ApixibanExperimentData.csv') %>% 
                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?

2
  • 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.