library(hdrm)

# bcTCGA ------------------------------------------------------------------
attachData(bcTCGA)

# Adaptive lasso, biased CV estimate
ylim <- c(0.15, 0.55)
Fig3.7(ylim=ylim)

# CV the entire adaptive lasso procedure
cvfit <- cv.adaptive_lasso(X, y)
plot(cvfit, bty='n', ylab = expression(CV(lambda)), ylim=ylim)

# CV: Regular lasso
cvfit.las <- cv.ncvreg(X, y, penalty='lasso')
plot(cvfit.las, bty='n', ylab = expression(CV(lambda)), ylim=ylim)
plot(cvfit.las$fit, log.l = TRUE, bty = "n")

# CV: MCP (gam=3, gam=7); basically same as Fig3.8()
out <- Fig3.8()
summary(out$cvfit3)
summary(out$cvfit7)

# SCAD
cvfit.scad <- Fig3.9()
summary(cvfit.scad)

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

attachData('whoari')
XX <- std(X)

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

summary(cvfit.mcp)
summary(cvfit.las)

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)

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')

b <- coef(cvfit.mcp)[-1]
sort(b[b!=0])