options ls=72 ps=72;

/*------------------------------------------------------*/
/*                                                      */
/* Proc Glimmix to fit random effcts model to the       */
/* epileptic seizure count data                         */
/*                                                      */
/*------------------------------------------------------*/


data seizure; 
  infile "seize.dat";
  input id seize visit trt age;
  nobs=_n_;
  interval = 2;
  if visit=0 then interval=8;
  logtime = log(interval);
  assign = (visit>0);
  assigntrt = assign*trt;
run;


title "Random intercept model for seizure data using proc GLIMMIX";
title2 "assuming conditional Poisson model";
proc glimmix data=seizure;
  class id;
  model seize = assign assigntrt / link=log dist=poisson offset=logtime s;
  random int / subject=id;
run;


title "Random intercept model for seizure data using proc GLIMMIX";
title2 "assuming conditional overdispersed-Poisson model";
proc glimmix data=seizure;
  class id;
  model seize = assign assigntrt / link=log dist=poisson offset=logtime s;
  random int / subject=id;
  random _residual_;
run;


/* Equivalent program using Proc Nlmixed if Proc Glimmix is not available */
title "Random intercept model for seizure data using proc NLMIXED";
title2 "assuming conditional Poisson model";
proc nlmixed; * qpoints=15;
  parms beta0=-1.4 beta1=0.12 beta2=-0.12 theta=0.1;
  bound theta>0;
  eta = beta0 + beta1*assign + beta2*assigntrt + b;
  mu = interval*exp(eta);
  model seize ~ poisson(mu);
  random b ~ normal(0, theta) subject=id;
run;