#SELECTION CRITERIA FOR SERIAL MODELS

#READ STRUCTURED DATA TABLE WITH NUMERIC CODED FACTOR
K=read.table("c:/2008LinearModelsData/SurgicalUnitR.txt")
K
attach(K)
options(digits=6)
Y=log(Y)

#VIEWING DATA
DATA=cbind(lnY,X1,X2,X3,X4)
DATA
#VIEWING CORRELATION MATRIX
cor(DATA)

#FITTING LINEAR MODEL
FM=lm(Y~X1+X2+X3+X4)

#SETTING UP MATRIX ALGEBRA OBJECTS AND HAT MATRIX
n=length(Y)
I=diag(n)
J=matrix(nrow=n,ncol=n,1)
X=model.matrix(FM)
p=ncol(X)
H=X%*%solve(t(X)%*%X)%*%t(X) #HAT MATRIX

#CALCULATING SUMS OF SQUARES
SSR=t(Y)%*%(H-((1/n)*J))%*%Y
SSR
SSE=t(Y)%*%(I-H)%*%Y
SSE
SSTO=t(Y)%*%(I-(1/n)*J)%*%Y
SSTO

#COMPARING WITH anova() OUTPUT
anova(FM)

#CALCULATING CRITERIA FOR MODEL SELECTION

#SUM OF SQUARES ERROR
SSE[1]  

#COEFFICIENT OF MULTIPLE DETERMINATION
Rsq=1-(SSE[1]/SSTO[1])
Rsq     
summary(FM)$r.squared #ALTERNATE CALCULATION

#ADJUSTED COEFFICIEINT OF MULTIPLE DETERMINATION
MSE=SSE[1]/(n-p)
MSTO=SSTO[1]/(n-1)
Rsqa=1-(MSE/MSTO)
Rsqa
summary(FM)$adj.r.squared #ALTERNATE CALCULATION

#MALLOW'S Cp
Cp=SSE[1]/MSE -(n-2*p)
Cp
require(wle) #DOWNLOAD {wle} PACKAGE FROM CRAN WEBSITE
mle.cp(FM)  #DON'T KNOW MUCH ABOUT THIS ONE, BUT INTERESTING!

#AKAIKE'S INFORMATION CRITERION AIC
AIC=n*log(SSE[1])-n*log(n)+2*p
AIC
extractAIC(FM)  #REPORTS: (equivalent df, AIC)  see ?extractAIC

#SCHWARTZ'S BAYESIAN CRITERION SBC
SBC=n*log(SSE[1])-n*log(n)+log(n)*p
SBC
extractAIC(FM,k=log(n))

#PRESS CRITERION
e=residuals(FM)
d=e/(1-diag(H))  #KNNL Eq. 10.21a
diag(H)          #MAIN DIAGONAL OF H AS A VECTOR
hatvalues(FM)    #ALTERNATE CALCULATION USING BUILT-IN FUNCTION & FM
PRESS=sum(d^2)
PRESS

#AUTOMATED STEPWISE REGRESSION
FM=lm(Y~X1+X2+X3+X4)
RM=lm(Y~1)

#STEPWISE REGRESSION USING AIC
step(FM,direction="backward")
step(RM,~X1+X2+X3+X4,direction="forward")
step(RM,~X1+X2+X3+X4,direction="both")

#STEPWISE REGRESSION USING SBC
step(FM,direction="backward",k=log(n))
step(RM,~X1+X2+X3+X4,direction="forward",k=log(n))
step(RM,~X1+X2+X3+X4,direction="both",k=log(n))


influence.measures(FM)