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

  Fit a loglinear mixed model to the epileptic seizure data
  of Thall and Vail (1990) using the SAS macro glimmix.

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

  First access the glimmix macro.

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

options ls=80 ps=59 nodate; run;

%inc 'glmm800.sas'  / nosource;

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

  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 (visit) to allow the dependence with time to be 
    subject-specfic.  

  See the macro for information on syntax, required statements, 
  and options.  See the proc mixed documentation for information
  on specifying the model using the "stmts=" statement.

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

title "RANDOM INTERCEPT ONLY";
%glimmix(data=seizure,
  procopt=method=ml,
  stmts=%str(
     class subject;
     model seize = logbase logage trt visit4 basetrt / solution;
     random intercept / subject=subject type=un;
  ),
  error=poisson,
  link=log
);
run;


title "RANDOM INTERCEPT AND TIME";
%glimmix(data=seizure,
  procopt=method=ml,
  stmts=%str(
     class subject;
     model seize = logbase logage trt visit4 basetrt / solution;
     random intercept visit / subject=subject type=un g gcorr;
  ),
  error=poisson,
  link=log
);
run;