### R code from vignette source 'lasso.Rnw'

###################################################
### code chunk number 1: lasso.Rnw:61-72
###################################################
library(alr3)
set.seed(1)
test <- sample(1:202,150,replace=FALSE)
ais0 <- with(ais,data.frame(LBM=LBM,logSSF=log(SSF),logWt=log(Wt),Sex=Sex,
           logHg=log(Hg),logHt=log(Ht),logRCC=log(RCC),logHc=log(Hc),
           logWCC=log(WCC),logFerr=log(Ferr)))
ais0 <- as.data.frame(scale(ais0))
m0 <- lm(LBM~1,ais0, subset=test)
m1 <- update(m0,LBM~logWt+logSSF+logRCC+Sex+logHg+logHt+logWCC+logHc+logFerr)
(f <- as.formula(paste("~",paste(names(coef(m1))[-1],collapse="+"))))
step(m0, scope=list(lower=~1, upper=f), direction="forward")


###################################################
### code chunk number 2: lasso.Rnw:79-92
###################################################
library(MASS)
FScoefs <- function(m0, m1, data, trace=FALSE) {
   keepCoef <- function(m, aic) {
    all <- names(coef(m1))
    new <- names(coef(m))
    ans <- rep(0,length(all))
    ans[match(new,all)] <- coef(m)
    ans
   }
   out <- with(data,stepAIC(m0, scope=list(lower=~1, upper=f), k=0,
      trace=trace, keep=keepCoef, direction="forward"))
   rownames(out$keep) <- names(coef(m1))
   out$keep}


###################################################
### code chunk number 3: lasso.Rnw:94-97
###################################################
print(coefs <- FScoefs(m0,m1,ais0), digits=2)
n <- length(coef(m1))-1
steps <- 0:(dim(coefs)[2]-1)


###################################################
### code chunk number 4: fs
###################################################
matplot(steps,t(coefs[-1,]),lty=1,type="l",xlim=c(0,n+2),
  xlab="Step number",ylab="Coef est")
xpos = rep(rev(steps)[1],4)
ypos = coefs[2:5,xpos[1]]
text(xpos,ypos,rownames(coefs)[2:5],cex=0.6,pos=4)


###################################################
### code chunk number 5: lasso1
###################################################
library(glmnet)
X <- as.matrix(ais0[ , -1])
Y <- ais0[, 1]
ans <- glmnet(X[test,], Y[test], standardize=FALSE)
plot(ans,"norm",label=TRUE)


###################################################
### code chunk number 6: lasso.Rnw:188-200
###################################################
cv.glmnet <- function(m,X,Y,blocksize=10,nblocks=10,show=TRUE,...) {
 lam <- m$lambda
 error <- NULL
 for (i in 1:nblocks) {
  j <- sample(1:dim(X)[1], blocksize, replace=FALSE)
  m1 <- glmnet(X[-j,],Y[-j], lambda=lam, standardize=TRUE)
  error <- rbind(error,Y[j]-predict(m1,X[j,]))}
  sds <- apply(error,2,sd)
  if(show==TRUE) plot(lam,sds,...)
  which <- order(sds)[1]
  list(lam=lam, sds=sds, winner=which)
  }


###################################################
### code chunk number 7: lasso1a
###################################################
out <- cv.glmnet(ans,X[test,], Y[test], log="x")
ans$beta[,winner<-out$winner,drop=FALSE]


###################################################
### code chunk number 8: lasso.Rnw:210-214
###################################################
rss.lasso <- sum( (Y[-test] - predict(ans,X[-test,])[,winner])^2)
coef(befit <- stepAIC(m0,scope=list(lower=~1,upper=f),trace=F,direction="forward"))
rss.be <- sum( (Y[-test] - predict(befit,ais0[-test,]))^2)
data.frame(rss.lasso=rss.lasso, rss.be=rss.be)


###################################################
### code chunk number 9: swiss1
###################################################
pairs(banknote[,-7],col=rainbow(2)[banknote$Y+1])
mtext("Red are genuine",3,1.5,col="red",adj=0)


###################################################
### code chunk number 10: lasso4
###################################################
set.seed(44)
sel <- sample(1:200, 150, replace=FALSE)
fit <- glmnet(as.matrix(banknote[sel, -7]),banknote[sel, 7],"binomial")
plot(fit,label=TRUE)