/*  Fit gee with linear estimating equation for beta and
    quadratic equation for covariance parameters with gaussian
    working assumption using macro NLINMIX  */

/*  This program uses the newest version of the macro availble on the
    sas web site */

options ps=59 ls=80 nodate;
run;

/*  Include the SAS macro code -- the user should of course 
    substitute the appropriate path  */

%inc "nlinmix_macro.sas";

/*  Read in the data -- Seizure data from Thall and Vail (1990)

  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;

title "Seizure data of Thall and Vail (1990)";

/*proc print; run;*/

/*  Invoke the NLINMIX macro 

    -  use the variable name "predv" to define the mean function

    -  if you do not wish to include analytical derivatives (d_b1 through
       d_b5 here) SAS will calculate them automatically

    -  the parms statement gives starting values
 
    -  the stmts portion basically contains elements of a  call to 
       proc mixed with the linearized version of the mean model.  
       We use the repeated statement to specify the working correlation
       structure

    -  the procopt statement allows us to specify "empirical"
       which will produce the robust sandwich standard errors

*/

title2 "Unstructured working correlation";

%nlinmix(data=seizure,
  model=%str(
     predv = exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
  ),
  derivs=%str(
     d_b1 = exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b2 = logbase*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b3 = logage*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b4 = trt*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b5 = visit4*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b6 = basetrt*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     wt=1/predv;
  ),
  parms=%str(b1=-2.0 b2=1 b3=1 b4=-1 b5=-0.1 b6=0.5),
  stmts=%str(
     class subject;
     model pseudo_seize = d_b1 d_b2 d_b3 d_b4 d_b5 d_b6/ noint notest solution;
     repeated / subject=subject type=un rcorr; 
     weight wt;
  ), 
  expand=zero,
  procopt=%str(empirical method=ml)
)

title2 "Exchangeable working correlation";

%nlinmix(data=seizure,
    model=%str(
     predv = exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
  ),
  derivs=%str(
     d_b1 = exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b2 = logbase*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b3 = logage*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b4 = trt*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b5 = visit4*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b6 = basetrt*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     wt=1/predv;
  ),
  parms=%str(b1=-2.0 b2=1 b3=1 b4=-1 b5=-0.1 b6=0.5),
  stmts=%str(
     class subject;
     model pseudo_seize = d_b1 d_b2 d_b3 d_b4 d_b5 d_b6/ noint notest solution;
     repeated / subject=subject type=cs rcorr; 
     weight wt;
  ), 
  expand=zero,
  procopt=%str(empirical method=ml)
)


title2 "AR(1) working correlation matrix";


%nlinmix(data=seizure,
  model=%str(
     predv = exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
  ),
  derivs=%str(
     d_b1 = exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b2 = logbase*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b3 = logage*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b4 = trt*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b5 = visit4*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     d_b6 = basetrt*exp(b1+b2*logbase+b3*logage+b4*trt+b5*visit4+b6*basetrt);
     wt=1/predv;
  ),
  parms=%str(b1=-2.0 b2=1 b3=1 b4=-1 b5=-0.1 b6=0.5),
  stmts=%str(
     class subject;
     model pseudo_seize = d_b1 d_b2 d_b3 d_b4 d_b5 d_b6/ noint notest solution;
     repeated / subject=subject type=ar(1) rcorr; 
     weight wt;
  ), 
  expand=zero,
  procopt=%str(empirical method=ml)
)