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

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

  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 nlmixed documentation for information on syntax,
  required statements, and options.

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

proc nlmixed data=seizure;
  parms b0=-2.0 b1=1.0 b2=1.0 b3=-1 b4=-0.1 b5=0.5 sb2=0.2;
  eta = b0+b1*logbase + b2*logage + b3*trt + b4*visit4 + b5*basetrt + bi;
  f = exp(eta);
  model seize ~ poisson(f);
  random bi ~ normal(0,sb2) subject=subject;
run;


proc nlmixed data=seizure;
  parms b0=-2.0 b1=1.0 b2=1.0 b3=-1 b4=-0.1 b5=0.5 
        sb12=0.3 cb12=-0.02 sb22=0.009;
  eta = b0+b1*logbase + b2*logage + b3*trt + b4*visit4 + b5*basetrt 
        + b1i + b2i*visit4;
  f = exp(eta);
  model seize ~ poisson(f);
  random b1i b2i ~ normal([0,0],[sb12,cb12,sb22]) subject=subject;
run;