source('http://myweb.uiowa.edu/pbreheny/7600/s16/notes/fun.R')
require(ncvreg)
require(glmnet)

# TCGA: Adaptive lasso (naive), slides 2-4-----------------------------------

load(url("http://myweb.uiowa.edu/pbreheny/data/bcTCGA.RData"))

fit <- ncvreg(X, y, penalty='lasso')
b <- coef(fit, which=which.min(BIC(fit)))[-1]
w <- abs(b)^(-1)
w <- pmin(w, 1e10)  # cv.glmnet does not allow infinite penalty weights
cvfit <- cv.glmnet(X, y, penalty.factor=w, lambda.min=5e-5)
plot(cvfit, bty='n', ylab=expression(CV(lambda)), xlab=expression(log(lambda)),
     las=1, xlim=rev(range(log(cvfit$lambda))))

# TCGA: Adaptive lasso (cross-validating the entire procedure) --------------

# Fit
fit.init <- ncvreg(X, y, penalty='lasso')
b <- coef(fit, which=which.min(BIC(fit)))[-1]
w <- abs(b)^(-1)
w <- pmin(w, 1e10)  # cv.glmnet does not allow infinite penalty weights
fit.full <- ncvreg(X, y, penalty.factor=w, lambda.min=5e-5, penalty='lasso')
cvfit <- cv.ncvreg(X, y, penalty.factor=w, lambda.min=5e-5, penalty='lasso')

set.seed(1)
nfolds <- 10
n <- length(y)
cv.ind <- ceiling(sample(1:n)/n * nfolds)

E <- matrix(NA, n, length(fit.full$lambda))
pb <- txtProgressBar(1, nfolds, style=3)
for (i in 1:nfolds) {
  fit0 <- ncvreg(X[cv.ind!=i,], y[cv.ind!=i], penalty='lasso')
  b <- coef(fit0, which=which.min(BIC(fit0)))[-1]
  w <- abs(b)^(-1)
  w <- pmin(w, 1e10)  # cv.glmnet does not allow infinite penalty weights
  fit1 <- ncvreg(X[cv.ind!=i,], y[cv.ind!=i], penalty.factor=w, lambda=fit.full$lambda, penalty='lasso')
  yhat <- predict(fit1, X[cv.ind==i,])
  E[cv.ind==i,] <- ncvreg:::loss.ncvreg(y[cv.ind==i], yhat, 'gaussian')
  setTxtProgressBar(pb, i)
}

min(apply(E, 2, mean))
ll <- log(fit.full$lambda)
cve <- apply(E, 2, mean)
cvse <- apply(E, 2, sd)/sqrt(nrow(E))
lwr <- cve - cvse
upr <- cve + cvse
ylim <- c(0.16, 0.54)

## Plot (slide 6)
# Biased
plot(cvfit, bty='n', ylab=expression(CV(lambda)), xlab=expression(log(lambda)),
     las=1, ylim=ylim, selected=FALSE, vertical.line=FALSE)

# Unbiased
plot(ll, cve, type='n', xlim=rev(range(ll)), ylim=ylim, bty='n', las=1,
     ylab=expression(CV(lambda)), xlab=expression(log(lambda)))
arrows(x0 = ll, x1 = ll, y0 = lwr, y1 = upr, code = 3, angle = 90, col = "gray80", length = 0.05)
points(ll, cve, pch=19, col='red', cex=0.5)

# Original lasso
cvfit0 <- cv.ncvreg(X, y, penalty='lasso')
plot(cvfit0, bty='n', ylab=expression(CV(lambda)), xlab=expression(log(lambda)),
     las=1, ylim=ylim, selected=FALSE, vertical.line=FALSE)

# TCGA: MCP/SCAD ----------------------------------------------------------

# MCP (gam=3, gam=7)
set.seed(1)
cvfit3 <- cv.ncvreg(X, y)
fit3 <- cvfit3$fit
set.seed(1)
cvfit7 <- cv.ncvreg(X, y, gam=7)
fit7 <- cvfit7$fit

# Slide 9
ylim=c(-0.21, 0.51)
par(mfrow=c(1,2))
xlim <- log(c(fit3$lambda[1], cvfit3$lambda.min))
plot(fit3, xlim=xlim, log.l=TRUE, bty='n', ylim=ylim)
mtext(expression(gamma==3))
xlim <- log(c(fit7$lambda[1], cvfit7$lambda.min))
plot(fit7, xlim=xlim, log.l=TRUE, bty='n', ylim=ylim)
mtext(expression(gamma==7))

# Slide 10
ylim=c(0.18, 0.551)
par(mfrow=c(1,2))
plot(cvfit3, bty='n')
mtext(expression(gamma==3), line=3)
plot(cvfit7, bty='n')
mtext(expression(gamma==7), line=3)

# Slides 11-12
summary(cvfit3)
summary(cvfit7)

# SCAD (slides 14-15)
set.seed(1)
cvfit <- cv.ncvreg(X, y, penalty='SCAD', gam=8)
fit <- cvfit$fit
xlim <- log(c(fit$lambda[1], cvfit$lambda.min))
par(mfrow=c(1,2))
plot(fit, xlim=xlim, log.l=TRUE, bty='n')
plot(cvfit, bty='n')
summary(cvfit)

# Lasso (for comparison)
set.seed(1)
cvfit <- cv.ncvreg(X, y, penalty='lasso')
fit <- cvfit$fit
xlim <- log(c(fit$lambda[1], cvfit$lambda.min))
par(mfrow=c(1,2))
plot(fit, xlim=xlim, log.l=TRUE, bty='n')
plot(cvfit, bty='n')
summary(cvfit)

# WHO-ARI -----------------------------------------------------------------

load(url('http://myweb.uiowa.edu/pbreheny/data/whoari.RData'))
XX <- std(X)

# Slide 17
set.seed(2)
cvfit.mcp <- cv.ncvreg(XX, y, gam=6)
fit.mcp <- cvfit.mcp$fit
set.seed(2)
cvfit.las <- cv.ncvreg(XX, y, penalty="lasso")
fit.las <- cvfit.las$fit

# Slide 18
ylim <- c(0, 0.4)
par(mfrow=c(1,2))
plot(cvfit.mcp, type='rsq', bty='n', ylim=ylim)
mtext('MCP', line=3)
plot(cvfit.las, type='rsq', bty='n', ylim=ylim)
mtext('Lasso', line=3)

# Slide 19
ylim <- c(-0.5, 0.5)
par(mfrow=c(1,2))
xlim <- log(c(fit.mcp$lambda[1], cvfit.mcp$lambda.min))
plot(fit.mcp, bty='n', xlim=xlim, log.l=TRUE, ylim=ylim)
mtext('MCP')
xlim <- log(c(fit.las$lambda[1], cvfit.las$lambda.min))
plot(fit.las, bty='n', xlim=xlim, log.l=TRUE, ylim=ylim)
mtext('Lasso')

# Slides 20-21
summary(cvfit.mcp)
summary(cvfit.las)