Emax Model for Drug Combination

Modeling
1

Please share your Stan program and accompanying data if possible.

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use

data {
int<lower=1> N; // total number of subjects // 
int<lower=0> D;  // number of all unique combinations
vector[N] dose1; // list of dose1 for each subject
vector[N] dose2;
real<lower=0> Dose_set[D,2]; // dose amount by combination all unique
vector[N] y; //response EFF 0.12,.13....

  real<lower=0> a1;
  real<lower=0> b1;
  real<lower=0> a2;
  real<lower=0> b2;
  real<lower=0> a3;
  real<lower=0> b3;
  real<lower=0> a4;
  real<lower=0> b4;
  real          a5;
  real<lower=0> b5;
  real<lower=0> a6;
  real<lower=0> b6;
   real          a7;
   real<lower=0> b7;
  int<lower=0, upper=1> Alpha_Dist; 
}

parameters {
 
  //real <lower=0> alpha; 
  real<lower=0> Emax;		// positive
  real<lower=0> theta3;	// IC50 for drug A
  real<lower=0> theta4;	// IC50 for drug B
  real<lower=0> tau;		// positive (gamma, slope of the sigmoid)
  real<lower=0> phi;		// positive precision
 // real<lower=0> sigma;
 real alpha0; // interaction parameter //normal/lognormal
 real B;
}


transformed parameters{
real sigma;
real alpha;
vector[N] X;
vector[D] X2;
for (i in 1:N){ 
 
  X[i]  = (((dose1[i]/theta3)+(dose2[i]/theta4)+alpha*(dose1[i]/theta3)*(dose2[i]/theta4))^tau)/(1+((dose1[i]/theta3)+(dose2[i]/theta4)+alpha*(dose1[i]/theta3)*(dose2[i]/theta4))^tau);
  
  }
  
for (i in 1:D) {
     X2[i] =((((Dose_set[i,1]/theta3)+(Dose_set[i,2]/theta4)+alpha*(Dose_set[i,1]/theta3)*(Dose_set[i,2]/theta4))^tau)/(1+((Dose_set[i,1]/theta3)+(Dose_set[i,2]/theta4)+alpha*(Dose_set[i,1]/theta3)*(Dose_set[i,2]/theta4))^tau));

}
  sigma = 1 / sqrt(phi);
  
  if (Alpha_Dist == 0){ 
    // 0:Normal, 1:LogNormal
    alpha = alpha0;
  } 
  else {
    alpha = exp(alpha0);
  }
}

model {
  
 
real mu[N];
   
for( i in 1:N){
  

mu[i] = B+ (Emax*X[i]);
//y[i] ~ normal(mu[i],sigma);

}

y ~ normal(mu,sigma);  ////// phi
    // prior
  Emax ~ gamma(a1, b1);
  theta3 ~ gamma(a2, b2);  //ic50A
  theta4 ~ gamma(a3, b3);  // IC50B
  tau ~ gamma(a4, b4);    // GAMMA SHAPE OF SIGMOID
alpha0 ~ normal(a5,b5);  // INTERACTION
 phi ~ gamma(a6, b6);   //  gamma
 B~normal(a7,b7);

}


generated quantities {

 real Y_mean[D]; 
 real Y_pred[D]; 
  
for (i in 1:D) {

	
  Y_mean[i] = (Emax*X2[i])+B;
    // Posterior predictive distribution
  Y_pred[i] = normal_rng(Y_mean[i], sigma);   
}
}

I am trying to get a bayesian version for the Greco model (1990) however getting this error:

Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Location parameter[1] is nan, but must be finite! (in ‘model15703f8f6e1f_comboemaxdata’ at line 82)

Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Location parameter[1] is nan, but must be finite! (in ‘model15703f8f6e1f_comboemaxdata’ at line 82)

Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”

Can you please help?

2

Put

vector<lower=0>[N] X;

in the transformed parameters block and it will show you when X[i] is NaN (presumably due to raising a negative number to a real power).

3

This is what I get.
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: validate transformed params: X[1] is nan, but must be greater than or equal to 0 (in ‘model16103575424e_comboemaxdata’ at line 46)

Did you mean that the initial values of tau is making the difference?

4

No, I mean you are taking a negative value and raising it to a positive real power and hence get a NaN for X[1] (and probably a bunch of other observations).