/******************************************************************

  Fit a loglinear mixed model to the epileptic seizure data
  of Thall and Vail (1990) using the SAS procedure glimmix as
  in Version 9.2 of SAS.

  These are count data, so we use the Poisson mean/variance
  assumptions.  

 ******************************************************************/

options ls=80 ps=59 nodate; run;

/******************************************************************

  The data look like (first 8 records on first 2 subjects)

       104 5 1  0 11 31
       104 3 2  0 11 31
       104 3 3  0 11 31
       104 3 4  0 11 31
       106 3 1  0 11 30
       106 5 2  0 11 30
       106 3 3  0 11 30
       106 3 4  0 11 30
          .
          .
          .

  column 1      subject
  column 2      number of seizures
  column 3      visit (1--4 biweekly visits)
  column 4      =0 if placebo, = 1 if progabide
  column 5      baseline number of seizures in 8 weeks prior to study
  column 6      age

******************************************************************/

data seizure; infile 'seize.dat';
  input subject seize visit trt base age;
run;

data seizure; set seizure;
*  if subject=207 then delete;
  logbase=log(base/4);
  logage=log(age);
  basetrt=logbase*trt;
  if visit<4 then visit4=0;
  if visit=4 then visit4=1;
run;

/*****************************************************************

  Fit two models: 

  - a model that contains only a random "intercept" in the linear
    predictor to capture inter-subject variation and induce
    correlation

  - a fancier model with random effects both for "intercept" and
    time (visit4) to allow the dependence with time to be 
    subject-specfic.  

  Use the results from running the glimmix macro to deduce starting
  values.

  See the proc glimmix documentation for information on syntax,
  required statements, and options.

******************************************************************/

proc glimmix data=seizure method=quad;
  class subject;
  model seize = logbase logage trt visit4 basetrt / solution link=log 
        dist=poisson;
  random int / subject=subject;
run;

proc glimmix data=seizure method=quad;
  class subject;
  model seize = logbase logage trt visit4 basetrt / solution link=log 
        dist=poisson;
  random int visit4  / type=un subject=subject;
run;