options ls=80 ps=1000 nodate;

/*------------------------------------------------------*/
/*                                                      */
/* Proc Glimmix to fit random intercept 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);
  agn_trt = assign*trt;
  int = 1;
run;


title "Random intercept model for seizure data with conditional overdispersion";
proc glimmix data=seizure;
  class id;
  model seize = assign agn_trt / dist=poisson link=log offset=logtime s;
  random int / subject=id type=vc;
  random _residual_;  *for conditional overdispersion; 
run; 

title "Random intercept model for seizure data without cond. overdispersion";
proc glimmix data=seizure;
  class id;
  model seize = assign agn_trt / dist=poisson link=log offset=logtime s;
  random int / subject=id type=vc;
run; 


title "Random intercept model for seizure data";
title2 "Using over-dispersed Poisson quasi-likelihood";
proc nlmixed qpoints=15;
  parms beta0=-1.4 beta1=0.12 beta2=-0.12 theta=0.1 phi=1;
  eta = beta0 + beta1*assign + beta2*agn_trt + b;
  mu = interval*exp(eta);
  if seize=0 then
    l = -(mu - seize)/phi - log(phi)/2;
  else 
    l = (seize*(log(mu) - log(seize)) - (mu - seize))/phi - log(phi)/2;
  model seize ~ general(l);
  random b ~ normal(0, theta) subject=id;
run;

title "Random intercept model for seizure data using proc NLMIXED";
title2 "assuming conditinal generalized Poisson model";
proc nlmixed; * qpoints=15;
  parms beta0=-1.4 beta1=0.12 beta2=-0.12 theta=0.1 xi=0.5;
  bound theta>0, xi>-1, xi<1;

  eta = beta0 + beta1*assign + beta2*agn_trt + b;
  mu = interval*exp(eta);
  mu1 = (1-xi)*mu;
  mu2 = mu1 + xi*seize;

  l = log(mu1) + (seize-1)*log(mu2) - mu2 - lgamma(seize+1);

  model seize ~ general(l);
  random b ~ normal(0, theta) subject=id;
run;
/* E(y|b) = mu; var(y|b) = mu/(1-xi)^2 */