/* This SAS code was used in conjunction with the SAS Macro GLIMMIX to obtain 
standard epidemiologic output from a multilevel analysis of the diet and breast 
cancer data.  We assume here that the SAS dataset "design" already exists, and 
includes observed values for the outcome variable ("disease"), plus fixed and 
random effects.  In our example, fixed effects include the covariates 
("calories", "BMI", "SES", "alcohol", "hormones", "age"), and the product of the 
first- and second-stage design matrices (XZ), given as ("col1" "col2" "col3" 
"col4" "col5" "col6" "col7" "col8").  Random effects are the exposures of 
interest ("caulif" "brocco" "brusse" "cabbrd" "cabbwh" "savoy" "skraut"  
"allium" "tomato" "pizza" "tomsalad" "strawb").  Second-stage variances are 
included in the dataset "tau2".  Results are output in "mmeqsol" as shown in 
figure 1, and then used in SAS proc IML to calculate odds ratios and 95% 
confidence limits.  Note that since we are using a logistic model, a binomial 
error term is specified in GLIMMIX.  Users of this code will need to make minor 
modifications to reflect one's own data structure.  An electronic version of 
this code, and other multilevel modeling software, is available at URL 
http://darwin.cwru.edu/~witte/hm.html */

/* Run the GLIMMIX Macro */
%inc '/glmm612.sas' / nosource;
%glimmix(data=design, procopt=mmeqsol,
  stmts=%str(
    model DISEASE = calories BMI SES alcohol hormones age
                    col1 col2 col3 col4 col5 col6 col7 col8 / solution;
    random          caulif brocco brusse cabbrd cabbwh savoy skraut allium
			  tomato pizza tomsalad strawb 
                    /gdata=tau2 solution;
    make 'mmeqsol' out=mmeqsol;),
    error=binomial)
run;

/* Use SAS proc IML to read estimates from GLIMMIX output into matrix form */ 

proc iml;
use mmeqsol;
read all var {_col1 _col2 _col3 _col4 _col5 _col6 _col7 _col8 _col9 _col10
  _col11 _col12 _col13 _col14 _col15 _col16 _col17 _col18 _col19
  _col20 _col21 _col22 _col23 _col24 _col25 _col26 _col27 _col28}
into mmeqsol;

pi    = mmeqsol[8:15,28];  	/* Fixed coefficient estimates (i.e., of XZ) */ 
covpi = mmeqsol[8:15,8:15];	/* Covariance of fixed coefficient estimates */ 
delta = mmeqsol[16:27,28];	/* Random coefficient estimates		     */	
covd  = mmeqsol[16:27,16:27];   /* Covariance of random coefficient estimates*/ 
cov_pd= mmeqsol[8:15,16:27];	/* Covariance of fixed and random coeff.     */

/* Calculate estimates, odds ratios and CLs. Z needs to be entered here.     */

b     = Z*pi + delta; 		/* Food-effect coefficients, equation (4)    */
orb   = exp(b);  		/* Odds ratios 				     */
varb  = Z*covpi*t(Z) + covd + Z*cov_pd + t(cov_pd)*t(Z); /* Covariance 	     */
 				/* matrix for the betas, equation (5) 	     */
stderr= sqrt(vecdiag(varb));	/* Standard errors of betas 		     */
lower = exp(b-1.96*stderr); 	/* Upper and lower 95% confidence limits     */
upper = exp(b+1.96*stderr);
print orb lower upper; 		/* Print results     			     */