%let location=C:\Documents and Settings\woodph\My Documents\glm2\kutnersolutions\chapter27;
PROC IMPORT OUT= WORK.wine 
            DATAFILE= "&location\chapter27.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH27TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means mean;class wine judge;
title "Wine Judging- Refer to Table 27.2, page 1132";run;
proc means mean;class wine;var rating;
title "Marginal means by Wine. Refer to Table 27.2, page 1132";run;

proc glm;class judge wine;
model rating=judge wine;
random judge/test;
means judge wine;
output out=resids r=rrating p=prating;
title "simple ss- basic GLM model Refer to Figure 27.3, page 1133";run;
proc capability data=resids noprint;var rrating;
  qqplot rrating;
title 'qq plot for overall model';run;
data resids;set resids;
rnormal=rrating;
proc rank normal=blom data=resids out=resids;var rnormal;
proc corr;var rnormal rrating;
title 'correlation between actual and normal refence page 1136';run;
proc gplot data=resids;
plot wine*rating=judge;
title 'dot plots by judge reference page 1136';run;


symbol1 i=join v=dot color=black l=1;
symbol2 i=join v=dot color=black l=2;
symbol3 i=join v=dot color=black l=3;
symbol4 i=join v=dot color=black l=26;
symbol5 i=join v=dot color=black l=33;
symbol6 i=join v=dot color=black l=41;
proc gplot;
plot rating*wine=judge/vaxis=15 to 30 by 5 haxis=1 to 4 by 1 vminor=0 hminor=0;
title "Plot of Wine Scores for Each Judge- Wine Tasting - Refer to Figure 27.2, page 1133";run;proc glm data=wine;class judge wine;
model rating=judge wine;
lsmeans wine/adjust=tukey cl pdiff ;
title 'least squares confidence intervals- reference p. 1137';run;
proc sort data=wine;by judge wine;
proc transpose data=wine out=winet prefix=wine;by judge;var rating;
proc corr data=winet cov;var wine1-wine4;
title 'within-subjects variance/covariance matrix reference page 1134, table 27.3';run;
proc glm;class judge;
model wine1 wine2 wine3 wine4=judge;
repeated wine 4;
title 'wine tasting as a repeated measures model note error term is judge*wine';run;
proc glm data=wine;class judge wine;
model rating=judge wine;
random judge/test;
lsmeans wine/adjust=tukey alpha=.05 pdiff tdiff cl;
title "Confidence intervals on Wies- Refer to Page 1137";run;
PROC IMPORT OUT= WORK.sweet
            DATAFILE= "&location\chapter27.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH27TA04$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "data means Refer to Table 27.4, page 1138";run;
/*I know that the data are ranked already, but maybe you'd be interested in getting the ranks
  computed for your data- here's now*/
proc sort data=sweet out=sweet;by subject;
proc rank data=sweet out=sweet;by subject;var prefer;
proc means;
title 'data means';run;
proc glm;class sweeten subject;
model prefer=sweeten subject;
lsmeans sweeten/alpha=.2 adjust=bon pdiff cl;
title 'Bonferroni corrected lsmeans reference page 1139';run;quit;
/*You can, if you want, get Kendall's concordance for such data
using the F* from the model. Here's how */
ods listing close;
ods trace on;
ods output OverallANOVA=anova FitStatistics=fit modelanova=model;
proc glm data=sweet;class sweeten subject;
model prefer=sweeten subject;
title 'concordance coefficient reference page 1139';
run;quit;
ods trace off;

ods listing;
proc print data=model;run;
data model;set model;
retain n_1 w;
if hypothesistype=1 and source='subject' then n_1=df;
*n_1 is n-1 which is df for subject;
w=fvalue/(fvalue+N_1);
proc print;run;
data model;set model;
if hypothesistype=3 and source='sweeten';
proc print;
title "kendall's concordance. Refer to Page 1139";run;
PROC IMPORT OUT= WORK.Shoe 
            DATAFILE= "&location\chapter27.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH27TA07$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;class campaign time;
title "Shoes Sales Example. Refer to Table 27.7, page 1146";run;
symbol1 i=join v=dot color=black l=1;
symbol2 i=join v=dot color=black l=2;
symbol3 i=join v=dot color=black l=3;
symbol4 i=join v=dot color=black l=26;
symbol5 i=join v=dot color=black l=33;
symbol6 i=join v=dot color=black l=41;
proc sort data=shoe out=shoe;by campaign;
proc gplot data=shoe;plot sales*time=market/vaxis=200 to 1200 by 100 haxis=1 to 3 hminor=0 vminor=0;by campaign;
title "Plots of sales Data by Test Market & Campaign. Refer to Figure 27.6, page 1146";run;
proc sort data=shoe out=shoe;by campaign market;
proc glm;class campaign market time;
model sales=campaign|time market(campaign);
random market(campaign);
test h=campaign e=market(campaign);
lsmeans time/adjust=tukey alpha=.01 pdiff tdiff cl;
title "Analysis of Variance Athletic Shoes Example. Refer to Figure 27.7, and LSMeans, Page 1148";
run;
proc transpose data=shoe out=shoet prefix=Period;by  campaign market;var sales;
run;
proc print data=shoet;run;
/* you're going to be looking mostly at the statistical tests for the repeated measures part, 
The univariate tests can be useful, but they require that you meet the covariance matrix
to meet the sphericity assumption relative to the hypothesis being tested. People frequently use
Mauchley's test for this, which is what you get via the printe option in the repeated measures
statement*/
/*For these data, sphericity is rejected, so you have the choice of using the Huynh-Feldt adjustment
or the Greenhouse-Geisser adjustment. The latter is conservative for smallsamples,
while the Huynh-Feldt is a bit more liberal. Generally, if sphericity is modeslty large (.15-p-.05)
one would use Huynh-Feldt. If small, use G-G. If it's very small, multivariate tests would be appropriate
and that's generally covered in a course on Manova.*/
/*Multivariate tests make no assumptions about within-subject covariance, and are generally valid
for all covariance structures. Univariate tests are, however, more powerful provided sphericity is met.*/
/*Proc GLM with repeated option will give you all tests except the test of market within campaign*/
/*The following code demonstrates various alternative contrasts within repeated*/
proc glm data=shoet;class campaign;
model period1 period2 period3=campaign;
repeated time 3 polynomial/mean printe;
title 'repeated measures anova with polynomial constrasts';run;
proc glm data=shoet;class campaign;
model period1 period2 period3=campaign;
repeated time 3 helmert/mean printe;
title 'repeated measures anova with helmert constrasts';run;
proc glm data=shoet;class campaign;
model period1 period2 period3=campaign;
repeated time 3 profile /mean printe;
title 'repeated measures anova with profile';run;
proc glm data=shoet;class campaign;
model period1 period2 period3=campaign;
repeated time 3 contrast/mean printe;
title 'repeated measures anova with contrast @ time 3 constrasts';run;
proc glm data=shoet;class campaign;
model period1 period2 period3=campaign;
repeated time 3 contrast(1)/mean printe;
title 'repeated measures anova with contrast @ time 1 constrasts';run;
/* Blood Flow Example*/
PROC IMPORT OUT= WORK.exercise 
            DATAFILE= "&location\chapter27.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH27TA08$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value bdyfmt
1="1=Bone  "
2="2=Brain"
3="3=Skin"
4="4=Muscle"
5="5=Heart";
data exercise;set exercise;
format bodypart bdyfmt.;
proc means mean;class exercise bodypart subject;
title 'Blood Flow Data Refer to Table 27.8, p. 1150';run;
/*Data can be transposed to take advantage of repeated option in GLM-
  We do it in this program both ways*/
proc sort data=exercise out=exercise;by exercise subject bodypart;
proc transpose data=exercise out=exerciset prefix=bodypart;by exercise subject;
data exerciset;set exerciset;
format bodypart bdyfmt.;

proc print data=exerciset;
title 'transposed data';run;
symbol1 i=join v=dot;
proc sort data=exercise out=exercise;by exercise bodypart subject;
proc gplot data=exercise;by exercise;
plot flow*bodypart=subject/vaxis=0 to 25 by 5 vminor=0 hminor=0;
title "line plots of data- reference page 1150, Figure 27.8";run;
proc glm data=exercise;class exercise subject bodypart;
model flow=exercise subject(exercise) bodypart bodypart*exercise bodypart*subject(exercise) ;
random subject(exercise);
test h=bodypart*exercise e=bodypart*subject(exercise);
lsmeans bodypart*exercise/slice=exercise adjust=tukey alpha=.1 e=bodypart*subject(exercise) pdiff cl adjust=tukey;
title 'repeated measures SS decomposition- reference p. 1150 Table 27.9';
title2 "Spells out All Error Terms to Match Book";run;


proc glm data=exercise;class exercise subject bodypart;
model flow= exercise subject(exercise) bodypart bodypart*exercise ;
random subject(exercise);
test h=exercise e=subject(exercise);
lsmeans bodypart*exercise/slice=exercise adjust=tukey alpha=.1 pdiff tdiff cl;
title "analysis of least square example using one error term";
title2 "refer to Table 27.9 and page 1151 for the LS comparisons";
title3 "Notice you're getting all the pairwise comparisons here";
run;
/*The following is probably "The Way to Do It" if you're lazy.
  Oddly, you have to run the model twice to get the within comparisons. If you run mixed
  with both LSmeans statements, you get all ten groups compared with all ten groups*/
proc mixed data=exercise method=reml asycov cl covtest;class subject bodypart;
model flow=exercise bodypart bodypart*exercise/cl ddfm=satterth;
random subject;
lsmeans bodypart/at exercise=1 cl diff pdiff adjust=tukey;
title 'alternative analysis using Proc Mixed ref 1151 & 1152';
run;

proc mixed data=exercise method=reml asycov cl covtest;class subject bodypart;
model flow=exercise bodypart bodypart*exercise/cl ddfm=satterth;
random subject;
lsmeans bodypart/at exercise=2 cl diff pdiff adjust=tukey;
title 'alternative analysis using Proc Mixed ref. 1151 & 1152';
run;
/*Blood Flow example. I did this analysis two ways: Once "straight" as the book shows it
  and then using the doubly repeated option in GLM*/
PROC IMPORT OUT= WORK.drugflow 
            DATAFILE= "&location\chapter27.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH27TA12$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means mean;class a b subject;
title "Drug Flow Example. Refer to Table 27.12, page 1158";run;
proc glm;class a b subject;
model flow=subject a b a*b subject*a subject*b;
random subject;
test h=a e=subject*a;
test h=b e=subject*b;
lsmeans a*b/pdiff cl adjust=bon cl;
title 'analysis of variance summary table- reference p. 1159, table 27.11';run;
/*If the correct error term is what's specified in the model, you don't have to go
  to all that work:*/
proc glm;class a b subject;
model flow=subject a|b;
random subject;
lsmeans a*b/pdiff cl adjust=bon cl;
title 'analysis of variance output. Refer Figure 27.11, page 1159';run;
/*Same analysis, but using repeated statement*/
/*Make transpose of data so repeated statement can be used*/
proc sort data=drugflow out=drugflow;by subject a b;

proc transpose data=drugflow out=drugflowt;by subject;var flow;run;

data drugflowt;set drugflowt;
rename col1=a1b1 col2=a1b2 col3=a2b1 col4=a2b2;
proc print;run;
proc glm;
model a1b1 a1b2 a2b1 a2b2=;
repeated a 2, b 2;
title "Doubly repeated analysis of variance.same as p. 1159";
title2 "Probably the way they did it from the looks of the printout in the book";
run;
proc means data=drugflow;class a b;var flow;
output out=outmeans mean=mflow;
data outmeans;set outmeans drugflow;if a>. and b>.;
symbol1 i=join v=dot;
symbol2 i=join v=square h=3;
symbol3 i=none v=circle h=3;
symbol4 i=none v=square h=3;
symbol5 i=none v=square h=3;
symbol6 i=none v=square h=3;
symbol7 i=none v=square h=3;
symbol8 i=none v=square h=3;
symbol9 i=none v=square h=3;
symbol10 i=none v=square h=3;
symbol11 i=none v=square h=3;
symbol12 i=none v=square h=3;
proc print data=outmeans;run;
proc gplot;
plot mflow*a=b/vaxis=-5 to 30 by 5 vminor=0;
plot2 flow*a=b/vaxis=-5 to 30 by 5 vminor=0;

title 'simple plots for the design. Refer to Figure 27.12, page 1160';
run;