options ls=72 ps=72;

/*------------------------------------------------------*/
/*                                                      */
/* Proc nlmixed 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";
proc nlmixed qpoints=15;
  parms beta0=-1.4 beta1=0.12 beta2=-0.12 theta=0.1;
  eta = beta0 + beta1*assign + beta2*assigntrt + b;
  mu = interval*exp(eta);
  model seize ~ poisson(mu);
  random b ~ normal(0, theta) subject=id;
run;