*Lecture I;
PROC IMPORT OUT= WORK.WQ 
            DATAFILE= "W:\courses\stat704\EColi.xls" 
            DBMS=EXCEL REPLACE;
     RANGE="EColi$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;

data lcongwq; set wq; 
where station in ('C-001','C-021','C-068','C-076','C-077');
lEcoli=log(EColi);
lFC=log(FecalColi);
lentero=log(Enterococci);
label lFC='log(Fecal coliform)'
      lEntero='log(Enterococci)'
	  lEcoli='log(Ecoli)';

proc sgscatter data=lcongwq;
matrix lEcoli lFC lEntero/diagonal=(histogram kernel); run;

*To see the data scale;
proc sgscatter data=lcongwq;
matrix lEcoli lFC lEntero; 
run;

proc reg data=lcongwq;
model lEcoli=lFC lEntero/clb; *alpha=0.90 for 90% CI, etc;
run;

proc reg data=lcongwq;
model lEcoli=lFC lEntero/covb corrb; *alpha=0.90 for 90% CI, etc;
run;
       
*Lecture II; 
PROC IMPORT OUT= WORK.senic 
            DATAFILE= "Z:\STAT 704\SENIC.txt" 
            DBMS=TAB REPLACE;
     GETNAMES=YES;
     DATAROW=2; 
RUN;

proc sgscatter data=senic;
matrix stay--facilities /diagonal=(histogram kernel); run;

proc sql;
insert into senic
(stay, facilities, xray)
values(10, 50.0, 100.0)
;
quit;

proc reg data=senic;
model infection=stay facilities xray/clm cli alpha=0.05;
output out=outsenic r=Residuals;
run;

*Diagnostics;
proc sgscatter data=senic;
matrix infection stay facilities xray /diagonal=(histogram kernel); run;

*Linear trend in variance;
proc model data=senic;
parms beta0 betas betaf betax;
infection=beta0+betas*stay+betaf*facilities+betax*xray;
fit infection/breusch=(1 stay facilities xray);
run;

*Normality test;
proc univariate data=outsenic normal;
var Residuals;
run;

*Two-way interactions residuals plot;
data outsenic; set outsenic;
SxF=stay*facilities;
SxX=stay*xray;
FxX=facilities*xray;
run;

proc sgplot data=outsenic;
scatter x=SxF y=Residuals;
reg x=SxF y=Residuals/nomarkers;
loess x=SxF y=Residuals/nomarkers;
refline 0/axis=y;
run;

proc sgplot data=outsenic;
scatter x=SxX y=Residuals;
reg x=SxX y=Residuals/nomarkers;
loess x=SxX y=Residuals/nomarkers;
refline 0/axis=y;
run;

proc sgplot data=outsenic;
scatter x=FxX y=Residuals;
reg x=FxX y=Residuals/nomarkers;
loess x=FxX y=Residuals/nomarkers;
refline 0/axis=y;
run;