# Slide 7
adj <- 0.07
par(mar=c(2,0,0,2))
plot(0, 0, type="l", xaxt="n", yaxt="n", xlim=c(-0.1,1), ylim=c(-0.1,0.9), bty="n", xlab="", ylab="")
arrows(0,0,1,0)
text(1+adj, -adj, expression(bold(x)[1]), xpd=TRUE)
arrows(0,0,1/2,sqrt(3)/2)
text(1/2+adj, sqrt(3)/2-adj, expression(bold(x)[2]), xpd=TRUE)
points(0, 0, pch=19)
text(-adj, -adj, expression(bold(0)))
points(0.8,0,pch=19)
text(0.8-adj, -adj, expression(bar(bold(y))[1]))
points(0.4+sqrt(3)/3, 1/3, pch=19)
text(0.4+sqrt(3)/3+adj, 1/3-adj, expression(bar(bold(y))[2]), xpd=TRUE)
col <- "#008DFFFF"
lines(c(0, 0.4), c(0,0), col=col, lwd=3)
lines(c(0.4, 0.4+1/2), c(0, sqrt(3)/2), col="gray", lty=2, lwd=3)
lines(c(0.4, 0.4+sqrt(3)/3), c(0, 1/3), col=col, lwd=3)

# Pollution data
load(url('http://myweb.uiowa.edu/pbreheny/data/pollution.RData'))
source('http://myweb.uiowa.edu/pbreheny/7600/s16/notes/fun.R')
require(glmnet)
X <- std(X)

# Slide 20
fit <- glmnet(X, y)
plot(fit, label=TRUE)
lam <- 1.84            # From cross-validation (next lecture)

## Slide 21
par(mar=c(5, 5, 5, 5))
fit <- glmnet(X, y)
l <- fit$lambda
b <- coef(fit)[-1,]
nv <- sapply(predict(fit, type="nonzero"), length)
matplot(log(l), t(b), col=pal(ncol(X)), lwd=3, bty="n", xlim=rev(range(log(l))), type="l", lty=1, xaxt="n", xlab=expression(log(lambda)), ylab="", las=1)
mtext(expression(hat(beta)), 2, 3, las=1)
at <- c(40, 4, 0.4, 0.04)
axis(1, at=log(at), labels=at)
ind <- abs(b[,ncol(b)]) > 15
text(x=log(min(l)/2.5), y=b[which(ind),ncol(b)], colnames(X)[ind], xpd=TRUE)
text(x=log(min(l)/2.5), y=10, "SO2", xpd=TRUE)
ind <- seq(1, length(l), length=5)
axis(3, at=log(l)[ind], labels=round(nv[ind], 1))
mtext("Number of nonzero coefficients", 3, 3)
abline(v=log(lam), col="gray")

## Slide 21, original (i.e., not log) axis
par(mar=c(5, 5, 5, 5))
b <- coef(fit)[-1,]
nv <- sapply(predict(fit, type="nonzero"), length)
matplot(l, t(b), col=pal(ncol(X)), lwd=3, bty="n", xlim=rev(range(l)), type="l", lty=1, xlab=expression(lambda), ylab="", las=1)
mtext(expression(hat(beta)), 2, 3, las=1)
ind <- abs(b[,ncol(b)]) > 15
text(x=-5, y=b[which(ind),ncol(b)], colnames(X)[ind], xpd=TRUE)
text(x=-5, y=10, "SO2", xpd=TRUE)
ind <- sapply(seq(40, 0, -10), function(x) which.min(abs(l-x)))
axis(3, at=l[ind], labels=round(nv[ind], 1))
mtext("Number of nonzero coefficients", 3, 3)
abline(v=log(lam), col="gray")

# What are the nonzero coefficients?
predict(fit, s=lam, type="nonzero")
colnames(X)[unlist(predict(fit, s=lam, type="nonzero"))]

# Slide 23
fit.ridge <- ridge(X,y)
b.lasso <- as.numeric(coef(fit, s=lam))[-1]
b.ridge <- coef(fit.ridge, which=which.min(fit.ridge$GCV))[-1]
boxplot(b.ridge, b.lasso, col='gray', las=1, names=c('Ridge', 'Lasso'),
        frame.plot=FALSE)
mtext(expression(hat(beta)), 2, 3, las=1)