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

  CHAPTER 11, EXAMPLE 1
  
  Fit a loglinear regression model to the epileptic seizure data.
  These are count data, thus we use the Poisson mean/variance
  assumptions.  This model is fitted with different working 
  correlation matrics.

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

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;

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

  Fit the loglinear regression model using PROC GENMOD and
  three different working correlation matrix assumptions:

  -  unstructured
  -  compound symmetry (exchangeable)
  -  AR(1)

  Subject 207 has what appear to be very unusual data -- for
  this subject, both baseline and study-period numbers of seizures
  are huge, much larger than any other subject.  We leave this subject
  in for our analysis; in some published analyses, this subject
  is deleted.  See Diggle, Liang, and Zeger (1995) and 
  Thall and Vail (1990) for more on this subject.

  We fit a mean model similar to that in Thall and Vail (1990), with
  covariates
  
  -  logbase = log (base/4)
  -  logage = log(age) 
  -  trt
  
  The DIST=POISSON option in the model statement specifies
  that the Poisson requirement that mean = variance, be used.  
  The LINK=LOG option asks for the loglinear model.  Other 
  LINK= choices are available.

  The REPEATED statement specifies the "working" correlation
  structure to be assumed.    The CORRW option in the REPEATED 
  statement prints out the estimated working correlation matrix
  under the assumption given in the TYPE= option.  The COVB
  option prints out the estimated covariance matrix of the estimate
  of beta -- both the usual estimate and the "robust" version
  are printed.  The MODELSE option specifies that the standard
  error estimates printed for the elements of betahat are based
  on the usual theory.  By default, the ones based on the "robust"
  version of the sampling covariance matrix are printed, too.  

  The dispersion parameter phi is estimated rather then being held 
  fixed at 1 -- this allows for the possibility of "overdispersion"

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

data seizure; set seizure;
*  if subject=207 then delete;
  logbase=log(base/4);
  logage=log(age);
run;

title "UNSTRUCTURED CORRELATION";
proc genmod data=seizure;
  class subject;
  model seize = logbase logage trt visit /  dist = poisson link = log;
  repeated subject=subject / type=un corrw covb modelse;
run;

title "EXCHANGEABLE (COMPOUND SYMMETRY) CORRELATION";
proc genmod data=seizure;
  class subject;
  model seize = logbase logage trt visit /  dist = poisson link = log;
  repeated subject=subject / type=cs corrw covb modelse;
run;

title "AR(1) CORRELATION";
proc genmod data=seizure;
  class subject;
  model seize = logbase logage trt visit /  dist = poisson link = log;
  repeated subject=subject / type=ar(1) corrw covb modelse;
run;

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

  We fit a modified model exactly like the one in Thall and Vail (1990).
  Here, visit4 is an indicator of the last visit and basetrt is
  the "interaction" term for logbase and trt, allowing the possibility
  that how baseline seizures affect the response may be different
  for treated and placebo subjects.  Visit4 allows that the mean
  response changed over the study period, but not gradually; rather,
  the change only showed up toward the end of the study;

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

data seizure; set seizure;
  basetrt=logbase*trt;
  if visit<4 then visit4=0;
  if visit=4 then visit4=1;
run;

title "MODIFIED MODEL WITH UNSTRUCTURED CORRELATION";
proc genmod data=seizure;
  class subject;
  model seize = logbase logage trt visit4 basetrt /  
                 dist = poisson link = log;
  repeated subject=subject / type=un corrw covb modelse;
run;