Two-parameter Bayesian logistic regression model (BLRM) model

Modeling
1

I am working on a Two-parameter Bayesian logistic regression model (BLRM) model, published by Neuenschwander et al. 2008, i.e. a Bayesian approach to finding the maximum-tolerated dose in phase I cancer trials. The model represents dose-toxicity curve with P(toxicity) specified for each dose level similar to
dose_toxicity_curve

The orignal winBUGS code reads

model{
# prior covariance matrix
cova[1,1] <- Prior[3]*Prior[3]
cova[2,2] <- Prior[4]*Prior[4]
cova[1,2] <- Prior[3]*Prior[4]*Prior[5]
cova[2,1] <- cova[1,2]
prec[1:2,1:2] <- inverse(cova[,])
log.alphabeta[1:2] Ëś dmnorm(Prior[1:2],prec[1:2,1:2])
alphabeta[1] <- exp(log.alphabeta[1])
alphabeta[2] <- exp(log.alphabeta[2])
# sampling model
for (j in 1:Ncohorts){
    logit(Pr.Tox1[j])<- log.alphabeta[1]+alphabeta[2]*log(DosesAdm[j]/DoseRef)
    Ntox[j] Ëś dbin(Pr.Tox1[j],Npat[j])
}

# for each dose: probabilities of toxicity, category probabilities, risks
for (i in 1:Ndoses) {
   lin[i]<-log.alphabeta[1]+alphabeta[2]*log(Doses[i]/DoseRef)
   logit(Pr.Tox[i]) <- lin[i]
   Pr.Cat[i,1] <- step(Pint[1]-Pr.Tox[i])
   Pr.Cat[i,2] <- step(Pint[2]-Pr.Tox[i])-step(Pint[1]-Pr.Tox[i])
   Pr.Cat[i,3] <- step(Pint[3]-Pr.Tox[i])-step(Pint[2]-Pr.Tox[i])
   Pr.Cat[i,4] <- step(1-Pr.Tox[i])-step(Pint[3]-Pr.Tox[i])
   Risk1[i] <- inprod(Pr.Cat[i,],LossFunction1[1:4]) 
   Risk2[i] <- inprod(Pr.Cat[i,],LossFunction2[1:4])
   Risk3[i] <- inprod(Pr.Cat[i,],LossFunction3[1:4])
}
} # data
list( Ncohorts=5, DoseRef=250,
DosesAdm=c(1,2.5,5,10,25),Npat=c(3,4,5,4,2),
Ntox=c(0,0,0,0,2), Ndoses=10,Doses=c(1,2.5,5,10,15,20,25,30,40,50),
Prior=c(2.15,0.52,0.84,0.78,0.20),Pint=c(0.2,0.35,0.6),
LossFunction1=c(1,0,1,1),LossFunction2=c(1,0,1,2),LossFunction3=c(1,0,2,4))

The RStan ChatGPT-translated code reads:

data {
  int<lower=1> Ncohorts;  // Number of cohorts
  real<lower=0> DoseRef;  // Reference dose level
  int<lower=1> Ndoses;  // Number of dose levels
  int<lower=1> Npat[Ncohorts];  // Number of patients in each cohort
  int<lower=0> Ntox[Ncohorts];  // Number of toxicities observed in each cohort
  real<lower=0> DosesAdm[Ncohorts];  // Doses administered to each cohort
  real<lower=0> Doses[Ndoses];  // List of all possible dose levels
  vector[5] Prior;  // Prior mean and covariance parameters
  vector[3] Pint;  // Probability intervals for categorization
  vector[4] LossFunction1;  // Loss function for risk calculation 1
  vector[4] LossFunction2;  // Loss function for risk calculation 2
  vector[4] LossFunction3;  // Loss function for risk calculation 3
}

parameters {
  vector[2] log_alphabeta;  // Log-transformed alpha and beta parameters
}

transformed parameters {
  real alphabeta[2];  // Exponentiated alpha and beta parameters
  real<lower=0, upper=1> Pr_Tox1[Ncohorts];  // Toxicity probability for administered doses
  real<lower=0, upper=1> Pr_Tox[Ndoses];  // Toxicity probability for all dose levels
  real lin[Ndoses];  // Linear transformation for dose-response modeling
  real Pr_Cat[Ndoses,4];  // Probability categories for toxicity risk levels
  real Risk1[Ndoses];  // Risk calculation using LossFunction1
  real Risk2[Ndoses];  // Risk calculation using LossFunction2
  real Risk3[Ndoses];  // Risk calculation using LossFunction3
  
  // Convert log-transformed parameters to actual values
  alphabeta[1] = exp(log_alphabeta[1]);
  alphabeta[2] = exp(log_alphabeta[2]);
  
  // Compute probability of toxicity for administered doses
  for (j in 1:Ncohorts) {
    Pr_Tox1[j] = inv_logit(log_alphabeta[1] + alphabeta[2] * log(DosesAdm[j] / DoseRef));
  }
  
  // Compute toxicity probabilities and risk categories for all dose levels
  for (i in 1:Ndoses) {
    lin[i] = log_alphabeta[1] + alphabeta[2] * log(Doses[i] / DoseRef);
    Pr_Tox[i] = inv_logit(lin[i]);
    
    // Compute probability category for each toxicity threshold
    Pr_Cat[i,1] = step(Pint[1] - Pr_Tox[i]);  // Low risk
    Pr_Cat[i,2] = step(Pint[2] - Pr_Tox[i]) - step(Pint[1] - Pr_Tox[i]);  // Moderate risk
    Pr_Cat[i,3] = step(Pint[3] - Pr_Tox[i]) - step(Pint[2] - Pr_Tox[i]);  // High risk
    Pr_Cat[i,4] = step(1 - Pr_Tox[i]) - step(Pint[3] - Pr_Tox[i]);  // Very high risk
    
    // Compute risk scores using loss functions
    Risk1[i] = dot_product(Pr_Cat[i], LossFunction1);
    Risk2[i] = dot_product(Pr_Cat[i], LossFunction2);
    Risk3[i] = dot_product(Pr_Cat[i], LossFunction3);
  }
}

model {
  vector[2] mu = Prior[1:2];  // Mean vector for prior distribution
  matrix[2,2] cova;  // Covariance matrix for prior distribution
  
  // Define covariance matrix using prior parameters
  cova[1,1] = square(Prior[3]);
  cova[2,2] = square(Prior[4]);
  cova[1,2] = Prior[3] * Prior[4] * Prior[5];
  cova[2,1] = cova[1,2];
  
  // Prior distribution for log(alpha) and log(beta)
  log_alphabeta ~ multi_normal(mu, cova);
  
  // Likelihood: Binomial model for toxicity observations
  for (j in 1:Ncohorts) {
    Ntox[j] ~ binomial(Npat[j], Pr_Tox1[j]);
  }
}

// library(rstan)
// data_list <- list(
//   Ncohorts = 5, DoseRef = 250, Ndoses = 10,
//   Npat = c(3, 4, 5, 4, 2), Ntox = c(0, 0, 0, 0, 2),
//   DosesAdm = c(1, 2.5, 5, 10, 25),
//   Doses = c(1, 2.5, 5, 10, 15, 20, 25, 30, 40, 50),
//   Prior = c(2.15, 0.52, 0.84, 0.78, 0.20),
//   Pint = c(0.2, 0.35, 0.6),
//   LossFunction1 = c(1, 0, 1, 1),
//   LossFunction2 = c(1, 0, 1, 2),
//   LossFunction3 = c(1, 0, 2, 4)
// )
// fit <- stan("model.stan", data = data_list, iter = 2000, chains = 4)
// print(fit)

My question is how to plot

  • Toxicity probability for all dose levels, Pr_Tox
  • Probability category for each toxicity threshold, Pr_Cat[i,1], …, Pr_Cat[i,4]
2

Have you validated that you’re getting the same answers as BUGS? Or have you done other tests that the model works like posterior predictive checks? If not, I’d strongly suggest doing that as a first step.

I’m not sure what you want to plot. What you’ll get out is a num_draws x num_parameters matrix of posterior draws, one per row. You can extract the posterior draws for any of the individual variables such as Pr_Tox[1] or Pr_Tox[2]. But how to plot these usually depends on your goals. For example, a common way to plot more than two variables but less than 20 would be a pairs plot that scatterplots each pair of plots with a histogram on the diagonal. That works for any small set of parameters. We don’t otherwise have generic plotting. You can plot the histogram for each parameter. Or you can make 3D scatterplots you can rotate with 3 of the parameters.

I’d suggest focusing on posterior predictive checks for quantities of interest.

Warning: GPT is not very good at Stan. Having said that, it’s getting better—o1 is a lot better than the older models were. The problem is that there’s so much inefficient and/or outdated Stan code out there that it’s prone to generate highly inefficient code.

Note: There’s not “RStan code” per se—Stan code is portable across all the interfaces. You might find cmdstanr easier to use (and certainly easier to install) than RStan. It also keeps up with Stan better because it’s not bound by CRAN policy.

Some notes on coding efficiency

For example, dealing with log_alphabeta, the way to write this most efficiently is as follows:

transformed data {
  cov_matrix[2] cova = [[square(Prior[3]), Prior[3] * Prior[4] * Prior[5]],
                        [Prior[3] * Prior[4] * Prior[5], scale(Prior[4])]];
  cholesky_factor_cov[2] cova_chol = cholesky_factor(cova);
  ...
parameters {
  vector[2] log_alphabeta_std;  // Log-transformed alpha and beta parameters
  ...
transformed parameters {
  vector[2] alphabeta = exp(cova_chol * log_alphabeta_std);
  ...

The binomial should be vectorized,

Ntox ~ binomial(Neat, Pr_Tox1);

I would advise against giving parameters generic names like Prior and would prefer to see something more meaningful here. I also found this odd:

alphabeta[1] = exp(log_alphabeta[1]);
alphabeta[2] = exp(log_alphabeta[2]);
...
inv_logit(log_alphabeta[1] + alphabeta[2] * log(DosesAdm[j] / DoseRef));

You can just use alphabeta = exp(log_alphabeta);, but why do you have a logistic regression with the slope constrained to be positive here?

Pretty much everything else can be vectorized.

I wouldn’t name a vector LossFunction as that implies it’s a function and not a vector.

1 Like
3

The logistic curve is a dose-toxicity curve - higher dose means higher toxicity, so we are looking for positive slope here. Hope that makes sens.