/*************************************************************************
 ** Program: "Problem 1 xover trial.sas"                                **
 ** Project:  2004 Summer Institute - Multilevel Models Class           **
 ** Author:   Michael Griswold                                          **
 **                                                                     **
 ** Description: SAS code for the 2x2 crossover trial example           **
 **              from DHLZ (2002) pg 148-150 & 180-181.                 **
 **               1) Input data                                         **
 **               2) Marginal Model (Alternating Logistic Regressions)  **
 **               3) Conditional Model (Random Intercepts)              **
 **               4) Comparison                                         **
 ** Notes:                                                              **
 **   Instructors: Francesca Dominici, Scott Zeger, Michael Griswold    **
 *************************************************************************/

***********************************************************************;
******--   1) Input & format data from DHLZ (2002) pg 148.    --*******;
***********************************************************************;
data xover_gps; 
  input trt $ r11 r01 r10 r00 n;
  cards;
  AB 22 0 6 6 34
  BA 18 4 2 9 33
  ;
  run;

**-- expand data to wide format --**;
data xover_wide; set xover_gps;
 do i=1 to n;
   if (i <= r11) then do;                          y1=1; y2=1; end;
   if (r11 < i <= r11 + r01) then do;              y1=0; y2=1; end;
   if (r11 + r01 < i <= r11 + r01 + r10) then do;  y1=1; y2=0; end;
   if (r11 + r01 + r10 < i ) then do;              y1=0; y2=0; end;
   output;
 end;
 run;
 **-- check data --**;
proc freq data=xover_wide;
   tables trt*y1*y2;
   run;

****-- create long-format dataset --****;
**-- pick off first observation (period=0) --**;
data xover1; set xover_wide(keep=trt y1 rename=(y1=y));
  trt_A = (trt="AB");
  period=0;
  id=_n_;
  run;
**-- pick off second observation (period=1) --**;
data xover2; set xover_wide(keep=trt y2 rename=(y2=y));
  trt_A = (trt="BA");
  period=1;
  id=_n_;
  run;
**-- append obs1 and obs2 --**;
data xover; set xover1 xover2; run;
**-- sort by id to get into long-format --**;
proc sort data=xover; by id; run;
 **-- check data --**;
proc freq data=xover;
  tables trt*y*period;
  run;


***********************************************************************;
****-- 2) MARGINAL, Alternating Logistic Regressions (GEE) model --****;
***********************************************************************;
title1 'Marginal Logistic Regression model (GEE-ALR)';
proc genmod data=xover descending;					*descending option models Pr(y=1);
  class id;											*subject identifier in class statement;
  model y = period trt_A / link=logit dist=bin;		*logistic regression MEAN MODEL;
  repeated subject=id / logor=exch;					*logistic regression ASSOC MODEL (exchangeable);
  ods output GEEempPEst=MargEsts;					*save Marginal estimates;
  run;


**********************************************************************;
****-- 3) CONDITIONAL, Random Effects Logistic Regression model --****;
**********************************************************************;

**-- define initial parameter estimates to use in proc nlmixed --**;
data InitEsts; set MargEsts;
  select ; 
    when ( Parm  	 = 'Intercept'	) Parameter = 'alpha0'  	; 
    when ( Parm      = 'period'    	) Parameter = 'alpha1'  	; 
    when ( Parm      = 'trt_A' 		) Parameter = 'alpha2'  	; 
    when ( Parm      = 'Alpha1'     ) Parameter = 'tau'  		;
  end;
  run;
 
title1 'Conditional Logistic Regression model (Random Int)';
   proc nlmixed data=xover;
      parms / data=InitEsts;    *use the marginal estimates as inital values;
      eta = alpha0 + alpha1*period + alpha2*trt_A + a;  *linear predictor specification;
      expeta = exp(eta);								*exponetiating the linear predictor;
      pi_c = expeta/(1+expeta);							*Conditional Pr(Success)=inverse logit;
	  model y ~ binomial(1,pi_c);						*Model Statement;
      random a ~ normal(0,tau*tau) subject=id;			*random intercept specification;
      ods output ParameterEstimates=CondEsts;			*save Conditional estimates;
   run;
title;


**********************************************************************;
******************--      4) Comparison         --********************;
**********************************************************************;
data compare;
  format Parameter $11.;
  merge MargEsts(keep=Parm Estimate Stderr 
				 rename=(Parm=Parameter Estimate=MarginalEstimate Stderr=MarginalStderr))
        CondEsts(keep=Estimate StandardError 
				 rename=(Estimate=ConditionalEstimate StandardError=ConditionalStderr));
  if Parameter="Alpha1" then Parameter="Association";
  run;
title "Marginal & Conditional Model Comparison";
proc print data=Compare; run;
title;