%let location=C:\Documents and Settings\woodph\My Documents\glm2\kutnersolutions\chapter18;
PROC IMPORT OUT= WORK.RUST 
            DATAFILE= "&location\chapter18.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="RustInhibitorData$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "Listing of Rust data (taken from Chapter 17). Refer to Table 18.1, page 777";run;
proc glm;class brand;
model inhibit=brand;
output out=rustr rstudent=rstudent p=predict r=residual;
symbol1 i=none v=dot;
axis1 offset=(5) minor=none;
proc gplot;
plot residual*predict/haxis=axis1 noframe;
title "Residual plot against Y predict. Refer to Table 18.1a, page 777";run;
proc univariate;var residual;class brand;
histogram residual/normal noframe;
title "Aligned residual plots. Refer to Figure 18.1b, page 777";run;
proc capability data=rustr noprint;
  qqplot rstudent/noframe;
title 'normal probability plot. Refer to Figure 18.1c, page 777';run;
title "studentized deleted residuals. Refer to Table 18.1, page 777";run;
PROC IMPORT OUT= WORK.solder 
            DATAFILE= "&location\chapter18.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH18TA04$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;var fluxtype joint strength;
title "solder joint pull strengths- ABT Electronics example. Refer to Table 18.2, page 783";run;
proc means mean std var ;var strength;class fluxtype;
title "summary statistics. Refer to Table 18.2, page 783";run;
proc univariate;var strength;class fluxtype;
histogram strength/normal noframe;
title "Plots of pull strengths. Refer to Figure 18.6, page 784";run;
proc glm;class fluxtype;
model strength=fluxtype;
means fluxtype/hovtest=bf;
title "Brown-Forsythe test, Refer to page 785";run;

proc means data=solder var ;class fluxtype;var strength;output out=wgt var=rwgt;
title 'computing and saving the variances of each class';run;
data wgt;set wgt;if fluxtype>.;
wgt=1/rwgt;keep wgt fluxtype;*Refer to Equation 18.14, page 786;
proc sort data=wgt out=wgt;by fluxtype;
proc sort data=solder out=solder;by fluxtype;
data solder;merge solder wgt;by fluxtype;
x=1;
proc boxplot data=solder;
plot strength*fluxtype/boxconnect=mean;
title 'boxplots of the data revealing heterogeniety';run;
proc glm data=solder;class fluxtype;
model strength=fluxtype/clparm;
means fluxtype/hovtest=bf;
contrast '123 vs 45' fluxtype 2 2 2 -3 -3;
title 'regular glm with contrast';
title2 'note the brown-forsythe results';
run;
proc model data=solder;
parms b0 b1 b2 b3 b4;
      strength = b0 + b1*x1+b2*x2+b3*x3+b4*x4;
      fit strength / white pagan=(1 x1 x2 x3 x4);
title 'Easy Way: SAS modified Breusch-Pagan Test';run;

proc glm data=solder;class fluxtype;
weight wgt;
model strength=fluxtype/clparm;
means fluxtype/tukey;
contrast '123 vs 45' fluxtype 2 2 2 -3 -3;
title 'weighted least squares regression weighted by reciprocal of variance. Refer to page 789';
run;
proc reg;weight wgt;
model strength=x1-x5/noint;
model strength=x/noint;
title 'analysis run using two regression models as shown in text. Refer to Figure 18.7, page 788';run;
PROC IMPORT OUT= WORK.Servo 
            DATAFILE= "&location\chapter18.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH18TA05$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "Time between Computer failures. Refer to Table 18.5, page 791";run;
%include "C:\Documents and Settings\woodph\My Documents\glm2\kutnersolutions\chapter18\boxglm.sas";
title "Box-Cox transformation using Friendly's Boxglm macro. Refer to Table 18.6 & discussion, page 792";
title2 "Note that this macro will tell you about the one significantly different residual via the leverage plot";
title3 "Which you see in Figure 18.8 later";
%boxglm(data=servo, resp=failure,model=location,class=location,gplot=effect rmse infl);run;

/*Alternatively, using transreg*/
proc transreg data = servo ss2;
	model boxcox(failure/lambda=-1 to 1 by .05)=class(location);
title "Box-Cox transformation using Transreg. Refer to Table 18.6 & discussion, page 792";
run;
data servo;set servo;
lfailure=log(failure);
proc glm;class location;
model failure lfailure=location/clparm;
lsmeans location/adjust=bon alpha=.1 pdiff cl;
output out=lservo r=rawresid logresid;
title "Bonferroni comparisons of log failures. Refer to Page 793";run;quit;
proc capability data=lservo noprint;
  qqplot rawresid logresid/noframe;
title "normal probability plots for original and transformed data. Refer to Figure 18.8a and 18.8b, page 792";run;
proc rank data=servo out=servor;var failure;
proc glm;class location;
model failure=location;
title "rank transformed data. Refer to Page 796";run;
proc npar1way data=servo;
class location;
exact wilcoxon;
var failure;
ods output KruskalWallisTest=Kruskalwallis anova=anova;
title "rank transform via npar1way. Refer to Page 796";run;quit;
data _null_;set anova;
if source='Among' then call symput('dfbetween',df);
if source='Within' then call symput('dfwithin',df);
data kruskalwallis;set kruskalwallis;
file print;
if label1='Chi-Square';
alpha=.1;
f=(&dfwithin*cvalue1)/(&dfbetween*(&dfwithin+&dfbetween-cvalue1));
critf=finv(1-alpha,&dfbetween,&dfwithin);
pval=1-cdf('F',f,&dfbetween,&dfwithin);
put 'Test statistic is ' f best12. ' with p value ' pval best12./
 'critical f for alpha of ' alpha 5.2 ' is ' critf 8.2;run;
proc glm data=servor;class location;
model failure=location/clparm;
lsmeans location/adjust=bon cl alpha=.1 pdiff;
title "Bonferroni adjusted pairwise comparisons. Refer to Page 797";run;
PROC IMPORT OUT= WORK.heart 
            DATAFILE= "&location\chapter18.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH18TA07$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;

proc means;run;
proc format;
value mismf
1='1=low'
2='2=medium'
3='3=high';
data heart;set heart;
format tissue mismf.;
proc sort data=heart out=heart;by tissue;
proc boxplot;
plot survival*tissue/boxconnect=mean;
title 'original boxplot of the data';run;
proc univariate;var survival;class tissue;
histogram survival/normal;inset mean std max min;
title 'histograms by tissue mismatch. Refer to Figure 18.9';run;
proc glm data=heart;class tissue;
model survival=tissue;
output out=resids r=residh rstudent=studentd;
title 'raw glm';run;
proc boxplot;
plot survival*tissue/boxconnect=mean;
title "Box plot of studentized deleted residuals. Refer to Figure 18.9b";
title2 "figure in book is of studentized residuals only- to get those, you would have specified student= instead of rstudent= in output stmt above";run;
proc capability data=resids noprint;
  qqplot residh;
title 'raw probability plot';run;
proc capability data=resids noprint;
  qqplot studentd;
title 'Studentized deleted residuals normal probability plot. refer to Figure 18.9c, page 799';run;

proc gplot;
plot survival*tissue;
title 'observed scores by tissue';run;
axis1 offset=(5);
symbol1 i=none v=dot;
proc gplot;
plot studentd*tissue/haxis=axis1 noframe;
title 'studentized deleted residuals by tissue. Refer to Figure 18.10, page 800';run;
proc transreg data = heart;
	model boxcox(survival)=identity(tissue);
title "Box-Cox transformation of heart data using SAS";
run;
%include 'C:\Documents and Settings\woodph\My Documents\glm2\kutnersolutions\chapter18\boxglm.sas';
title "Box-Cox transformation of heart data using Friendly's macro";
%boxglm(data=heart, resp=survival,model=tissue,class=tissue,gplot=effect rmse infl);run;
data heart;set heart;
survival=log(survival);

proc glm data=heart;class tissue;
model survival=tissue;
output out=residh r=residf rstudent=studentd;
title 'Anova of log transformed heart data Refer to Figure 18.10c, Page 800';run;
proc capability data=residh noprint;
  qqplot residf;
title 'log transformed probability plot. Refer to Figure 18.10b, page 800';run;
proc gplot;
plot survival*tissue;
title 'observed scores by tissue';run;
proc gplot;
plot studentd*tissue;
title 'studentized deleted residuals by tissue. Refer to Figure 10.10a, page 800';run;

proc capability data=residh noprint;
  qqplot residf studentd;
title 'q-q plot of residuals and studentized deleted residuals, Refer to Figure 18.10a, page 800';
run;

proc rank data=heart out=heartr;var survival;
proc glm;class tissue;
model survival=tissue;
lsmeans tissue/t adjust=bon cl pdiff tdiff alpha=.1;run;
title 'rank transformed data';run;