%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter2;
PROC IMPORT OUT= WORK.toluca 
            DATAFILE= "&location\chapter2.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH01TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
/* Let's start by using Proc Reg to compute and save residuals - This will
   help us understand what the numbers in Table 2.1 refer to in terms of these
   calculates variables.*/
Proc reg data=toluca;
model workhrs=lotsize/clb alpha=.05;
output out=tolucavars p=pworkhrs r=rworkhrs;
title "using Proc Reg to calculate predicted and residual scores";
title2 "Refer to Figure 2.2, page 46";run;
proc means mean std sum uss data=tolucavars;
title "Basic Statistics. Refer to Table 2.1 page 46 & Proc Reg output above";run;

/*The book refers at several points to Tables B.2 and B.4 in the back of the book
and asks you to look up the critical value. This is cumbersome. We actually do
this now just using software when we need these values:*/
data test;
a=tinv(1-.025,23);*note the division of alpha by 2;
proc print data=test;format a 20.16;
title "Toluca critical t value: two-sided test alpha=.05 (p. 47 ""we require...."")";run;
data test;
a=tinv(1-.05,23);
proc print data=test;format a 20.16;
title "critical t value: one-sided test alpha=.05 (p. 47-48 ""we require..."")";run;

/*Some times you might be interested in calculating the p level
  associated with a given observed statistic*/
data test;
a=1-probt(10.29,23);
proc print;format a 20.16;run;
title "p level associated with test statistic (p. 47, ""Many statistical calculators...."")";run;
/*Note that with the t-statistic, you'll have to keep track of which side of zero
  you are on- You don't have to do this with F, for example*/

proc reg data=toluca;
model workhrs=lotsize/clb alpha=.05;
title "regression with confidence interval around the regression weight p. 47";run;
proc reg data=toluca;
model workhrs=lotsize/clb alpha=.10;
title "regression with confidence interval around the regression weight p. 49";run;

/*When you plan a regression, it's a good idea to assess your odds of finding
  a statistically significant effect under various scenarios. This can get
  complicated for other studies, but for the case of a single predictor, here is the
  code for determining power. Use this presentation in lieu of the material
  presented on pages 50-51. We will present other software for testing specific 
  values of the unstandardized regression weight, but the following is sufficient
  to get us started*/
proc power;
      multreg
         model = fixed
         nfullpredictors = 1
         nreducedpredictors = 0
		 rsquarediff=.3 .5 .7 .8215
         ntotal = 23
         power = .;
      plot x=n min=10 max=50;
title "statistical power for statistical significance test when rsquare is .3, .5, .7, .8215";
run;

symbol1 i=rlcli90 v=dot;
proc gplot data=toluca;
plot workhrs*lotsize;
title "individual prediction 90% confidence intervals page 59";
title2 "refer to Figure 2.6, page 62";run;
symbol1 i=rlclm90 v=dot;
proc gplot data=toluca;
plot workhrs*lotsize;
title "mean level prediction 90% confidence intervals page 55";run;
proc reg data=toluca alpha=.10;
model workhrs=lotsize;
output out=outvals lcl=loweri lclm=lowerm uclm=upperi ucl=upperm;
title "example analysis in which alpha is changed";run;
proc print data=outvals;
title "upper and lower confidence intervals for observed data";
title2 "Refer to pages 55 and 59";run;
/*How to calculate confidence intervals for data points not in your original
  data: Basically, you add an observation which has predictor, but no criterion
  values*/
data newobs;
lotsize=100;output;
lotsize=65;output;
/*Concatenate the newobs data set onto the old toluca data set*/
data toluca;set newobs toluca;
proc reg data=toluca alpha=.10;
model workhrs=lotsize;
output out=outvals lcl=loweri lclm=lowerm uclm=upperm ucl=upperi;
proc print data=outvals;
title "upper and lower confidence intervals for observed data x=100 Refer to page 55";
title2 "a disadvantage is that this shows you all the data";
title3 "upper and lower confidence intervals for mean of 100, Refer to page 59";
data outval1;set outvals;
if _n_=1 or _n_=2;
proc print;
title "confidence intervals for lotsizes of 65 and 100 Refer to pages 54-55 & 59";run;

data regadmit;
file print;
hsgpa=3.5;
sigma=.12;
predy=.10+(0.95*hsgpa);
lowery=predy-(3*.12);
uppery=predy+(3*.12);
put "For predicted HS GPA of " hsgpa 5.3 " predicted grade is " predy 5.3/
 "with Lower CI of " lowery 5.3 " and upper CI of " uppery 5.3;
title "Calculations based on summary statistics. Refer to page 56";run;

/*You can also use SAS to write the output in your printout to a data set and then
  work with that- This is useful in keeping you from having to type things
  when you need ot make a table or chart for your research*/
ods listing close;

proc reg data=toluca outest=estims tableout;
model workhrs=lotsize/all alpha=.10;

output out=stats stdi=stdi stdp=stdp pred=yhat stdr=stdr;
title "regression with supplementary statistics";
ods output anova=out;
ods listing;
data stats;set stats;if workhrs=.;*keep only the value of interest;
proc print data=stats;
title "statistics in a data set";run;
data out;set out;
if source="Error";*keep only the ss & ms error;
proc print data=out;
title "ss and mse in a data set";run;
/*If you wanted, you coulod calculate the sigma of prediction, and the confidence interval
  by hand, for example. Why bother? Sometimes your research question (such as page 61)
  isn't directly available from your printout*/
data both;merge out stats;
sdpred=sqrt((ms/3)+stdp*stdp);
lower=yhat-tinv(1-.05,23)*sdpred;
upper=yhat+tinv(1-.05,23)*sdpred;
proc print data=both;
title "confidence interval for average of size 3 (page 61)";run;
data both;set both;
lower=lower*3;
upper=upper*3;
proc print data=both;
title "confidence interval for the total time required (page 61)";
proc glm data=toluca; 
model workhrs=lotsize/cli  clparm p inverse xpx solution;
title "analysis done in glm individual confidence interval";quit;
run;
symbol1 v=none i=rlcli90 c=red;
axis1 order=50 to 500 by 50 minor=none;
axis2 order=20 to 110 by 10 minor=none;
proc gplot data=toluca;
plot workhrs*lotsize/noframe vaxis=axis1 haxis=axis2;
title "Confidence Band for Regression Line. Refer to Figure 2.6, page 62";run;
proc glm data=toluca;
model workhrs=lotsize;
title "Overall model F refer to page 71";run;
proc reg data=toluca;
model workhrs=lotsize;
test lotsize=0;
title "same analysis comparing reduced and full models";
title2 "also refer page 75";run;
data oil;
file print;
r=.52;
n=23;
df=n-1-1;
criticalt=tinv(1-.05,df);
obst=(r*sqrt(df))/sqrt(1-(r*r));
probt=1-probt(obst,df);
fisherz=.5*log((1+r)/(1-r));
sigma=sqrt(1/(n-3));
alpha=.05;
fzlow=fisherz-(probit(alpha/2)*sigma);
fzhigh=fisherz+(probit(alpha/2)*sigma);
rlow=(exp(2*fzlow)-1)/(1+exp(2*fzlow));
rhigh=(exp(2*fzhigh)-1)/(1+exp(2*fzhigh));
put "For r " r 5.3 " Fisher z is " fisherz 6.3 " lower 90% is " fzlow 6.3 " upper 90% is " fzhigh 6.3
/ " with probability " probt 5.3
/ " lower bound on correlation is therefore " rlow 6.3 " and upper bound is " rhigh 6.3;
title "Fisher z refer to page 85";run;

data food;
file print;
r=-.61;
n=200;
df=n-1-1;
criticalt=tinv(1-.05,df);
obst=(r*sqrt(df))/sqrt(1-(r*r));
probt=1-probt(obst,df);
fisherz=.5*log((1+r)/(1-r));
sigma=sqrt(1/(n-3));
alpha=.05;
fzlow=fisherz-(probit(alpha/2)*sigma);
fzhigh=fisherz+(probit(alpha/2)*sigma);
rlow=(exp(2*fzlow)-1)/(1+exp(2*fzlow));
rhigh=(exp(2*fzhigh)-1)/(1+exp(2*fzhigh));
put "For r " r 5.3 " Fisher z is " fisherz 6.4 " lower 90% is " fzlow 6.4 " upper 90% is " fzhigh 6.4
/ " with probability " probt 6.4 " and sigma of " sigma 7.5
/ " lower bound on correlation is therefore " rlow 6.4 " and upper bound is " rhigh 6.4;
title "Fisher z refer to page 86";run;
PROC IMPORT OUT= WORK.food 
            DATAFILE= "&location\chapter2.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH02TA04$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title 'Refer to Table 2.4, page 88. Population and expenditures';run;
proc corr data=food fisher(alpha=.05);
title "Fisher z based on data.";run;
proc corr spearman fisher(alpha=.05);
title 'spearman correlation (page 88)';run;
/*Or, for fun, you can also rank the data and run a proc corr on that*/
proc rank data=food out=foodr;var pop expend;
proc corr;var pop expend;
title 'Alternatively, calculated as a correlation of ranks (also page 88)';run;
/*Rank normalizing is different than the rank transform*/
proc rank normal=blom data=food out=foodr;var pop expend;
proc corr;var pop expend;
title 'rank normalized correlation (not discussed in text)';run;
proc univariate data=food;
histogram pop expend/normal;
inset mean median skewness ;
title 'raw data';run;
proc univariate data=foodr;
histogram pop expend/normal;
inset mean median skewness ;
title 'rank normalized data';run;