data {

  int N;      // number of observations
  int J;      // number of  patients

  int y[N];   // This should be the "outcome" variable
              // in the "toenail" dataframe, i.e., 0's an 1's

  real month[N];

  int ID[N];  // convert ID in the dataframe to successive integers
              // with no gaps

  real treatment[N];

}
  
parameters{

  real b0;

  real b1;

  real b2;

  real eta0[J];

  real eta1[J];

  real<lower=0> tau0;

  real<lower=0> tau1;

}
  
transformed parameters {

  real logit_p[N];

  real a0[J];

  real a1[J];

  for (j in 1:J) {

    a0[j] <- b0 + eta0[j];

    a1[j] <- b1 + eta1[j];

  }

  for (i in 1:N) {

    logit_p[i] <- a0[ID[i]] + a1[ID[i]]*month[i] + b2*treatment[ID[i]];
  
  }

}
  
model {

  for (i in 1:N) {

    y[i] ~ bernoulli_logit(logit_p[i]) ;

  }

  for (j in 1:J) {

    eta0[j] ~ normal(0,tau0);

    eta1[j] ~ normal(0,tau1);

  }

  b0 ~ normal(0,1e6);

  b1 ~ normal(0,1e6);

  b2 ~ normal(0,1e6);

  tau0 ~ uniform(0,50);

  tau1 ~ uniform(0,50);
  
}
 
generated quantities {

  real yrep[N];

  for (i in 1:N) {
  
    yrep[i] <- bernoulli_logit_rng(logit_p[i]);

  }

}

CXXFLAGS = -O2 -Wall $(DEBUGFLAG) -mfpmath=sse -msse2 -mstackrealign $(LTO)