webreg <- "http://www.sagepub.co.uk//wrightandlondon//"
maleptsd <- read.table(paste(webreg,"maleptsd.dat",sep=""), header=T)

attach(maleptsd)
hist(PTSD)
hist(log(PTSD)+1)

ptsd <- log(PTSD + 1)
preds <- cbind(OVER2,OVER3,OVER5,BOND,POSIT,NEG,CONTR,SUP,CONS,AFF)
ptsdmat <- cor(preds)
print(round(ptsdmat,2))
print(apply(preds,2,sd),digits=2)

library(car)

allvars <- lm(ptsd~preds)
summary(allvars)

#####################################################################
#Linear Regression: Selection using leaps
install.packages("leaps")
library(leaps)
x1 <- regsubsets(preds,ptsd)
summary(x1)

par(mfrow=c(2,2))
plot(x1,scale=c("r2"))
title("Default for leaps")
plot(x1,scale=c("adjr2"))
title("Default for leaps")
x2 <- leaps(preds,ptsd,nbest=1, method="r2")
x3 <- leaps(preds,ptsd,nbest=1, method="adjr2")
m2 <- max(x3$adjr2)
m1 <- which.max(x3$adjr2)

plot(x2$size-1,x2$r2,xlab="Number of predictors",ylab=expression(R^2))
lines(spline(x2$size-1,x2$r2))
plot(x3$size-1,x3$adjr2, xlab="Number of predictors",ylab=expression(adj.R^2))
lines(spline(x3$size-1,x3$adjr2))

par(mfrow=c(1,1))
plot(x3$size-1,x3$adjr2,xlab="# of predictors not including intercept", ylab="Adj. R-squared")
lines(spline(x3$size-1,x3$adjr2))

arrows(m1,m2-.02,m1,m2)
text(m1,m2-.023,expression(paste("max ",r[adj]^2 == 0.2)))
text(m1,m2-.031, paste("with",format(m1,digits=1),"variables"),pos=3)

par(mfrow=c(1,1))

ivars<-c(1:10)
select.vars<-ivars[x3$which[which.max(x3$adjr2),]]



labels(preds[1,])

summary(lm(ptsd~preds[,ivars]))
summary(lm(ptsd~preds[,select.vars]))

#####################################################################

library(MASS)

lam.seq<-seq(0,100,by=0.1)

ptsd.ridge<-lm.ridge(ptsd~preds,lambda=lam.seq)

plot(lam.seq,ptsd.ridge$GCV,type="l",xlab=expression(lambda),ylab="GCV")
abline(v=lam.seq[which.min(ptsd.ridge$GCV)],col="red",lty=2)

plot(ptsd.ridge)
title("Ridge Regression: Estimates",xlab=expression(lambda))
abline(h=0)
abline(v=50,lty=3)

ptsd.ridge$coef[,402]

#####################################################################

#Lasso
install.packages("lars")
library("lars")
ptsd.lasso1 <- lars(preds,ptsd)
plot(ptsd.lasso1)

predict(ptsd.lasso1,s=7, mode="step", type="coefficient")$coefficients

ptsd.lasso1.cv <- cv.lars(preds,ptsd,plot.it=T)
plot(ptsd.lasso1.cv$index,ptsd.lasso1.cv$cv)

ptsd.frac <- ptsd.lasso1.cv$index[which.min(ptsd.lasso1.cv$cv)]

predict(ptsd.lasso1,fraction=ptsd.frac,mode="fraction", type="coefficient")$coefficients

detach(maleptsd)