data {

   int N;
   int y[N];
   real roach1[N];
   real treatment[N];
   real senior[N];
   real exposure2[N];

   int N_y_pred;
   matrix[N_y_pred,5] predictdata;  // read in from an R matrix

}

transformed data {

   real offset[N];

   for (i in 1:N) {
       offset[i] <- log(exposure2[i]);
   }

}

parameters {

   real b0;
   real b_roach1;
   real b_treatment;
   real b_senior;
}

model {

  real loglambda[N];

  for (i in 1:N) {
       loglambda[i] <- b0 + b_roach1*roach1[i] + b_treatment*treatment[i] +
                    b_senior*senior[i] + offset[i];
       y[i] ~ poisson_log(loglambda[i]);
  }

  b0          ~ normal(0,1e6);
  b_roach1    ~ normal(0,1e6);
  b_treatment ~ normal(0,1e6);
  b_senior    ~ normal(0,1e6);

}

generated quantities {

  vector[5] beta;
  
   // there is a subtlety here:
   //
   // "real beta[N]" is a vector for which stan doesn't know
   //    how to do vector & matrix algebra
   // "vector[N] beta" is a vector for which stan knows how to do
   //    vector addition, matrix multiplication, etc...
   //
   // We need the "vector" version in the "generated quantities" block
   // below...

  real y_pred[N_y_pred];

  vector[N_y_pred] pred_loglambda;

  beta[1] <- b0;
  beta[2] <- b_roach1;
  beta[3] <- b_treatment;
  beta[4] <- b_senior;
  beta[5] <- 1;                            // for the logexposure2 offset...  

  pred_loglambda <- predictdata * beta;    // matrix multiplication!

  for (i in 1:N_y_pred) {
      y_pred[i] <- poisson_rng(exp(pred_loglambda[i]));
  }

}

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