*Import SplitCov.xls as WORK.splitcov;
proc import out=splitcov
datafile="/home/grego1/STAT 706/SplitCov.xlsx"
dbms=xlsx
replace;
run;
/* Create a Whole Plot Covariate */
proc sort data=splitcov; by wp;
proc means data=splitcov noprint; by wp;
var spcov;
output out=b (where=(_stat_='MEAN') rename=(spcov=wpcov) drop=_TYPE_ _FREQ_); proc print; run;
/* Append the Whole Plot Covariate to the data */
data c; merge splitcov b; by wp; drop _STAT_ run; proc print; run;
/* Create a Whole Plot Factor Covariate; this */
/* term isn't really used                     */
proc sort data=splitcov; by wpf;
proc means data=splitcov noprint; by wpf;
var spcov;
output out=d mean=wpfcov; run; 
/* Append the WP Factor Cov to the data as well */
proc sort data=c; by wpf; proc print; run;
/* Inspect data set e, though it is not used for analysis */
data e; merge d c; by wpf;
/* Create a Split Plot Factor Covariate */
data f; set c;
spfcov=spcov-wpcov; run; 
/* Data set F has a split plot covariate and whole plot covariate, but no WPF covariate */
proc print data=f; run;
*Plot the covariates against the response;
proc sgplot data=f;
scatter x=spcov y=response/group=wpf;
loess x=spcov y=response/group=wpf nomarkers;
xaxis label="Covariate";
yaxis label="Response";
run;

proc sgplot data=f;
scatter x=wpcov y=response/group=wpf;
loess x=wpcov y=response/group=wpf nomarkers;
xaxis label="Whole plot covariate";
yaxis label="Response";
run;

proc sgplot data=f;
scatter x=spfcov y=response/group=wpf;
loess x=spfcov y=response/group=wpf nomarkers;
xaxis label="Split plot covariate";
yaxis label="Response";
run;


 proc print; run;
/* Traditional ANOVA */
proc glm data=f;
class Blck SPF WPF;
model response=Block WPF WPF*Block SPF WPF*SPF SPF*Block SPF*WPF*Block;
random Block/test; 
test h=WPF e=wpf*block;
test h=SPF e=spf*block;
test h=spf*wpf e=spf*wpf*block;
run;
/* Yandell ANOVA */
proc glm data=f;
class WP SPF WPF;
model response=WPF WP(WPF) SPF WPF*SPF;
random WP(WPF)/test;
run;
/* Yandell ANOVA with Whole Plot covariates */
proc glm data=f;
class SPF WPF WP;
model response=wpcov wpf wp(wpf) spf wpf*spf;
/* This next statement won't work */
*random wp(wpf)/test;
test h=wpf e=wp(wpf);
test h=wpcov e=wp(wpf)/htype=1 etype=1;
run;
/* Yandell ANOVA with Split Plot covariates */
proc glm data=f;
class SPF WPF WP;
model response=wpcov wpf wp(wpf) spfcov spf wpf*spf;
test h=wpf e=wp(wpf);
test h=wpcov e=wp(wpf)/htype=1 etype=1;
run;