library(MASS)
library(Design)
library(psych)
library(Rcmdr)
library(lmtest)
library(robustbase)
library(simpleboot)
library(MBESS)
library(relaimpo)

mireault=read.table("http://www.unt.edu:8080/rss/class/mike/data/Mireault.txt", header=T)


#Model will only use those who have experienced loss
mireault2=na.omit(mireault[,c("DEPRESST","AGEATLOS","PVLOSS","SUPPTOTL")])
attach(mireault2)
describe(mireault2)  
scatterplot.matrix(mireault2)  
cor(mireault2)  
cov.rob(mireault2, cor=T)$cor 

mymodel <- lm(DEPRESST~., data=mireault2)  
summary(mymodel)  
plot(mymodel)  
influencePlot(mymodel) #make sure you 'stop' the graph when you're done identifying outliers (right click on the graph)

#Statistical tests of assumptions
shapiro.test(mymodel$residuals)
vif(mymodel)
bptest(mymodel)
reset(mymodel)  


#Robust comparison.
myrobmodel=rlm(DEPRESST~., data=mireault2) 
summary(myrobmodel)  
myrobmodel$w 

AIC(mymodel)
AIC(myrobmodel)

#Robust bootstrapped intervals 
mybootmod= lm.boot(myrobmodel, R=500)  
mybootmod   
perc.lm(mybootmod$rsquare, p = c(0.025, 0.975))
confint(mymodel)  

#Validate ols model
mymodelval <- ols(DEPRESST~AGEATLOS+ PVLOSS + SUPPTOTL, x=T, y=T,data=mireault2)  
mymodelval
validate(mymodelval) 

#Interval on R2
ci.R2(yourbiasadjR2,p=3,N=135)

#Variable Importance
varimportance= calc.relimp(mymodel,rela=T)
varimportance
varimportance2=boot.relimp(DEPRESST~AGEATLOS+ PVLOSS + SUPPTOTL, rela=TRUE, b=500)
booteval.relimp(varimportance2,typesel="lmg",bty="bca")