OPTIONS LS=70;
proc import out=fly datafile='\\cocasrwebsrv2.ds.sc.edu\people\stat\grego\courses\stat706\bibd_fly.txt'
replace;
run;

/* MEANS VS LSMEANS  */
PROC SORT data=fly; BY CHEMICAL;
PROC MEANS data=fly; BY CHEMICAL;
VAR COUNT; RUN;
/* GRAND MEAN FOR COMPUTING EFFECTS */
PROC MEANS; VAR COUNT; 
output out=gmean;
RUN;

data gmean; set gmean; if _stat_='MEAN';
/* Clean up the data set a little                   */
xmean=count;
drop _stat_ _type_ _freq_ count;
run;
/* GLM  */
ods graphics on;
PROC GLM data=fly plots=all;
CLASS DAY CHEMICAL;
MODEL COUNT=DAY CHEMICAL/SOLUTION;
LSMEANS CHEMICAL/STDERR PDIFF;
run;
proc glimmix data=fly plots=all;
CLASS DAY CHEMICAL;
MODEL COUNT=DAY CHEMICAL/SOLUTION;
ods output LSMeans=LSM;
LSMEANS CHEMICAL/ PDIFF adjust=tukey cl;
run;
ods graphics;
*Combine the grand mean with LSMeans and compute alpha_hat;
data lsmalpha; 
if _N_=1 then set gmean;
set lsm; 
alpha_hat=estimate-xmean;
proc print data=lsmalpha;
run;