%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter4;
PROC IMPORT OUT= WORK.toluca 
            DATAFILE= "&location\Chapter4.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH01TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "toluca data";run;
/*Bonferroni correction applied to these data (alpha=.1) means we use alpha=.05
  in these data, which is the default (i.e., alpha/2)).*/
proc reg;
model WorkHrs=Lotsize/clb covb;
title "90% confidence intervals on regression weights. Refer Eq. 4.3, p. 156";
footnote "note that square root of diagonal of cov are the standard errors of b on page 156";run;
/*If you"re too lazy to do that, you can also have SAS output the standard errors to a data set and print it*/
proc reg data=toluca outest=seb outseb;
model WorkHrs=Lotsize/clb ;
title "p. 156";
footnote ; run;
proc print data=seb;
title "standard errors of b. Refer to Page 156";run;

/*Notice the formula for the covariance of intercept and slope in eq. 4.5. We can get
  SAS to confirm this calculation by outputting the variance/covariance matrix of b's*/
proc reg data=toluca outest=covb covout ;
model WorkHrs=Lotsize/clb corrb covb;
title "Covariance of b0 and b1 Refer to p. 156";run;
proc print data=covb;
title "variance/covariance matrix of b's p. 156";run;
proc means data=toluca;
title "Summary statistics for toluca data";run;
data eq4_5;set covb;
file print;
if _name_="LotSize";
covb0b1=-70*LotSize;
put "observed covariance between beta''s is: " Intercept;
put "and calculated by equation 4.5 it is: " covb0b1;
put "Isn't Science wonderful?";run;
/*To calculate the confidence intervals for particular values, create a data set in which the
X values are present and for which no criterion values are present as follows:*/
data newobs;
lotsize=30;output;
lotsize=65;output;
lotsize=100;output;
data toluca2;set toluca newobs;
/* Now, for the Working-Hotelling procedure, there is no "canned" way of doing this in SAS that
I know of. We'll need the mean square error and degrees of freedom error from the regression.
One straightforward way to do this is to capture this in sas:
Example on page 158:
*/
proc glm data=toluca2 outstat=stats;
model workhrs=lotsize;
output out=working p=pworkhrs stdi=stdi;
title "GLM model for securing error SS Refer to page 158";run;
proc print data=stats;
title "output statistics data set";run;
data stats;set stats;
drop f prob;
if _type_="ERROR";
mse=ss/df;
grp=1;
data working;set working;grp=1;
data workingm;merge working stats;by grp;
drop grp _source_ _Type_ workhrs _name_;
if workhrs=.;
syhat = sqrt((stdi)**2-mse);
  alpha=.1;*operant alpha;
  conf=1-alpha;*Family confidence coefficient;
  f = finv(1-alpha, 2, df); 
  w = sqrt( 2*f );
  lower = pworkhrs - w*syhat;
  upper = pworkhrs + w*syhat;
  label alpha="alpha level"
        conf="family confidence coefficient"
        lower="Lower Working-Hotelling bound"
        upper="Upper Working-Hotelling bound";
proc print;run;
data _null_;set workingm;
file print;
put "For a Predicted Low Size of " Lotsize ":"/
@5 " The Predicted Work Hours are: " pworkhrs 7.2 " with an alpha interval of " alpha 3.2 /
@10 " ranging from " lower 7.2 " to " upper 7.2;
title "Working-Hotelling procedure. Refer to page 158";
run;
/*Bonferroni adjustments to the predicted value are quite easy and require only
  calculating the adjusted alpha*/
proc reg data=toluca2 ;
model WorkHrs=Lotsize/clm corrb covb alpha=.0333333;
title "Bonferroni interval refer to p. 159";run;
data newobs;
lotsize=80;output;
lotsize=100;output;
data toluca2;set newobs toluca;

proc glm data=toluca2 outstat=stats;
model workhrs=lotsize;
output out=scheffeb p=pworkhrs stdi=stdi;
run;
data stats;set stats;
if _source_="ERROR";
drop _name_ _type_ f prob;
grp=1;
data scheffeb;set scheffeb;grp=1;
data scheffeb;merge scheffeb stats;by grp;
drop _source_ grp SS;
if workhrs=.;
*What are our conditions?;
alpha=.05;
predicts=2;*number of predictions;
Label
smult= "Scheffe Multiple"
bmult= "Bonferroni Multiple"
slow= "Scheffe Lower Bound"
shigh= "Scheffe Upper Bound"
blow= "Bonferroni Lower Bound"
bhigh= "Bonferroni Upper Bound";

smult=sqrt(2*finv(1-alpha,2,df));
bmult=tinv(1-alpha/(2*predicts),df);
slow=pworkhrs-smult*stdi;
shigh=pworkhrs+smult*stdi;
blow=pworkhrs-bmult*stdi;
bhigh=pworkhrs+bmult*stdi;
proc print data=scheffeb;run;
data _null_;set scheffeb;
file print;
put "Assuming a family confidence coefficient of 1-" alpha 4.2/
    @5 "The Lower and Upper Scheffe intervals are: " slow 8.2 shigh 8.2 " respectively, and"/
	@5 "The Lower and Upper Bonferroni intervals are: " blow 8.2 bhigh 8.2 " respectively"/
	@10 "(Assuming a Scheffe multiple of " smult 12.8 " and a Bonferroni multiple of " bmult 12.8 ")";
title "Scheffe and Bonferroni Bounds on multiple new observations. Refer to Page 160";run;
PROC IMPORT OUT= WORK.Warehous 
            DATAFILE= "&location\Chapter4.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH04TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "Table 4.2 page 163, Warehouse Example";run;
proc reg sscp;
model labor=workunit/noint alpha=.05 clb covb rmse;
title "Refer to page 164. Note SS present in Table 4.2 on page 163";run;
data sugar;
n=12;
b0=-.1;
b1=1.017;
mse=.0272;
sb1=.0142;
meanX=5.5;
meanY=5.492;
SSx=135;
t=b1/sb1;
critt=tinv(.975,10);
probt=1-probt(t,n-1);

patient=6.52;
xhat=(patient-b0)/b1;
varxhat=(mse/(b1**2.))*(1+(1/n)+(xhat-meanX)**2./SSx);
sxhat=sqrt(varxhat);
low=xhat-sxhat*critt;
high=xhat+sxhat*critt;
file print;
put "test statistic is " t= "critical value is " critt= "probability of t is " probt 
"predicted X is " xhat= "variance of predicted value is " varxhat= 
"standard deviation of predicted value is " sxhat= 
"lower bound is " low= " upper bound is " high=;
title "refer to pages 169-170";
run;