# R/tuneLasso.R In LINselect: Selection of Linear Estimators

# ---------------------------------------------------------------
#  LINselect R package
#  INRA, UR1404, Research Unit MaIAGE
#  F78352 Jouy-en-Josas, France.
# ---------------------------------------------------------------

tuneLasso <- function
### tune the  lasso parameter in the
### regression model : \eqn{Y= X \beta + \sigma N(0,1)}
### using the lasso or the gauss-lasso method
##references<<  See Baraud et al. 2010
## \url{http://hal.archives-ouvertes.fr/hal-00502156/fr/} \cr
## Giraud et al., 2013,
## \url{http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.ss/1356098553}
(Y, ##<<vector with n components : response variable.
X, ##<<matrix with n rows and p columns : covariates.
normalize=TRUE, ##<<logical : corresponds to the input \code{normalize}
## of the functions \code{\link[elasticnet]{enet}} and \code{\link[elasticnet]{cv.enet}}. \cr If TRUE the variates \code{X} are
## normalized.
method=c("lasso","Glasso"), ##<<vector of characters whose components are subset of
## (\dQuote{lasso}, \dQuote{Glasso})
dmax=NULL,##<<integer : maximum number of variables in the lasso
##estimator. \code{dmax} \eqn{\le} D where \cr
##  D = min (3*p/4 , n-5) if  p\eqn{ \ge }n
## \cr D= min(p,n-5) if
## p < n. \cr Default : \code{dmax} = D.
Vfold=TRUE,##<<logical : if TRUE the tuning is done by Vfold-CV
V=10,##<<integer. Gives the value of V in the Vfold-CV procedure
LINselect=TRUE,##<<logical : if TRUE the tuning is done by LINselect
a=0.5,##<<scalar : value of the parameter \eqn{\alpha} in the LINselect criteria
K=1.1,##<<scalar : value of the parameter \eqn{K} in the LINselect criteria
verbose=TRUE,##<<logical : if TRUE a trace of the current process is displayed in real time.
max.steps=NULL##<<integer : maximum number of steps in the lasso
##procedure. \cr Corresponds to the input \code{max.steps} of the function
##Default :
##\code{max.steps} = 2*min(p,n)
)
{
#
#
# packages needed : elasticnet
calc.proj <- FALSE
calc.pen <- FALSE
res <- NULL
n=nrow(X)
p=ncol(X)
Xmean <- apply(X,2,mean)
mu=mean(Y)
if (is.null(max.steps)) max.steps=2*min(p,n)
if (p>=n) dmax <- floor(min(c(3*p/4,n-5,dmax)))
if (p<n) dmax <- floor(min(c(p,n-5,dmax)))
#   library(elasticnet)

res.lasso <- try(enet(X,Y,lambda=0, intercept=TRUE, normalize=normalize,
max.steps=max.steps))
if (inherits(res.lasso, "try-error")) {
print("error  when calling enet")
res <- res.lasso[1]
}
if (!inherits(res.lasso, "try-error")) {
lesBeta <- array(0,c(dim(res.lasso$beta.pure)[1],p)) lesBeta[,res.lasso$allset] <- res.lasso$beta.pure lesDim <- apply(lesBeta!=0,1,sum) if (max(lesDim) <= dmax) Nmod.lasso <- length(lesDim) if (max(lesDim)>dmax) Nmod.lasso<-(1:dim(lesBeta)[1])[(lesDim>=dmax)][1] f.lasso <- matrix(0,nrow=n,ncol=Nmod.lasso) Intercept <- mu-lesBeta[1:Nmod.lasso,]%*%Xmean for (l in 1:Nmod.lasso) { f.lasso[,l] <- rep(Intercept[l],n)+X%*%lesBeta[l,] } if (is.element("lasso",method)) { res<- list(lasso=NULL) if (Vfold) { # lasso + V fold CV if (verbose) { print(paste("Tuning Lasso with Vfold CV: V=",V,", dmax=",dmax)) } s <- 1:(Nmod.lasso+1) res.cv.enet=0 while (!is.list(res.cv.enet)) { s <- s[-length(s)] res.cv.enet <- try(cv.enet(X, Y, K = V, lambda=0, s=s, mode="step", normalize=normalize, intercept=TRUE, plot.it=FALSE,se=FALSE),silent=TRUE) } l.lasso <- which.min(res.cv.enet$cv)
f.pred <- f.lasso[,l.lasso]
coeff <- c(Intercept[l.lasso],lesBeta[l.lasso,][lesBeta[l.lasso,]!=0])
names(coeff) <- c("Intercept",which(lesBeta[l.lasso,]!=0))
res$lasso$CV <- list(support=which(lesBeta[l.lasso,]!=0),
coeff=coeff, fitted=f.pred, crit=res.cv.enet$cv, crit.error=res.cv.enet$cv.error,
lambda=res.lasso$penalty[s]) } if (LINselect) { # lasso + LinSelect calc.proj <- TRUE if (verbose) print(paste("Tuning Lasso with LINselect: K=",K,", a=",a,", dmax=",dmax)) # pen <- penalite(p, n, K=K, dmax+1) D <- 0:dmax dm <- pmin(D,rep(p/2,dmax+1)) Delta <- lgamma(p+1)-lgamma(dm+1)-lgamma(p-dm+1)+2*log(D+1) pen <- penalty(Delta, n, p, K) calc.pen <- TRUE I.lasso <- list(NULL) ProjMod <- array(0,c(n,n,Nmod.lasso)) A <- array(Inf,c(Nmod.lasso,Nmod.lasso)) B=A SCR <- rep(0,Nmod.lasso) penSCR <- rep(0,Nmod.lasso) sumY2 <- sum((Y-mu)**2) un <- rep(1,n) for (l in 1:Nmod.lasso) { I.lasso[[l]] <- (1:p)[lesBeta[l,]!=0] if (length(I.lasso[[l]])>0) { ProjM <- Proj(cbind(un,X[,I.lasso[[l]]]),length(I.lasso[[l]])+1) ProjMod[,,l] <- ProjM$Proj
SCR[l] <- sum((Y-ProjMod[,,l]%*%Y)**2)
penSCR[l] <- SCR[l]*pen[ProjM$rg+1]/(n-ProjM$rg)
}
if (length(I.lasso[[l]])==0) {
ProjMod[,,l] <-  un%*%t(un)/n
ProjM <- list(Proj=ProjMod[,,l],rg=1)
SCR[l] <- sumY2
penSCR[l] <- sumY2*pen[ProjM$rg+1]/(n-ProjM$rg)
}
}
Ind <- rep(1:Nmod.lasso,rep(Nmod.lasso,Nmod.lasso))
YY<-Y%*%t(rep(1,Nmod.lasso))
for (m in 1:Nmod.lasso) {
Proj.f<-ProjMod[,,m]%*%f.lasso
B[m,] <- apply((YY-Proj.f)**2,2,sum)+a*apply((f.lasso-Proj.f)**2,2,sum)
A[m,] <- B[m,]+penSCR[m]
}
l.lasso <- Ind[which.min(A)]
f.pred <- f.lasso[, l.lasso]
coeff <- c(Intercept[l.lasso],lesBeta[l.lasso,][lesBeta[l.lasso,]!=0])
names(coeff) <- c("Intercept",which(lesBeta[l.lasso,]!=0))
crit <- apply(A,2,min)
res$lasso$Ls <- list(support=which(lesBeta[l.lasso,]!=0),
coeff=coeff, fitted=f.pred, crit=crit,
lambda=res.lasso$penalty[1:Nmod.lasso]) } } if (is.element("Glasso",method)) { if (!is.list(res)) res <- list(Glasso=NULL) if (!calc.pen) { D <- 0:dmax dm <- pmin(D,rep(p/2,dmax+1)) Delta <- lgamma(p+1)-lgamma(dm+1)-lgamma(p-dm+1)+2*log(D+1) pen <- penalty(Delta, n, p, K) # pen <- penalite(p, n, K=K, dmax+1) } if (!calc.proj) { I.lasso <- list(NULL) SCR <- rep(0,Nmod.lasso) penSCR <- rep(0,Nmod.lasso) sumY2 <- sum((Y-mu)**2) un <- rep(1,n) for (l in 1:Nmod.lasso) { I.lasso[[l]] <- (1:p)[lesBeta[l,]!=0] if (length(I.lasso[[l]])>0) { ProjM <- Proj(cbind(un,X[,I.lasso[[l]]]),length(I.lasso[[l]])+1) SCR[l] <- sum((Y-ProjM$Proj%*%Y)**2)
penSCR[l] <- SCR[l]*pen[ProjM$rg+1]/(n-ProjM$rg)
}
if (length(I.lasso[[l]])==0) {
SCR[l] <- sumY2
penSCR[l] <- SCR[l]*pen[2]/(n-1)
}
}
}
if (Vfold) {
# Gauss lasso + V fold CV
if (verbose) {
print(paste("Tuning Gauss Lasso with Vfold CV: V=",V,", dmax=",dmax))
}
I.CV <- list(NULL)
IM.CV <- list(NULL)
for (iv in 1:V) {
I.CV[[iv]] <- (floor((iv-1)*n/V)+1):(iv*n/V)
IM.CV[[iv]] <- (1:n)[-I.CV[[iv]]]
}
active_set<-list(NULL)
all_lambda<-c(-1)
lambda <- list(NULL)
for (i in 1:V) {
in_s<-enleve_var_0(X,IM.CV[[i]])
res.i<-enet(X[IM.CV[[i]],in_s],Y[IM.CV[[i]]],lambda=0,
intercept=TRUE,normalize=normalize,max.steps=max.steps)
lesBeta.i <- array(0,c(dim(res.i$beta.pure)[1],length(in_s))) lesBeta.i[,res.i$allset] <- res.i$beta.pure lesDim.i <- apply(lesBeta.i!=0,1,sum) if (max(lesDim.i) <= dmax) Nmod.lasso.i <- length(lesDim.i) if (max(lesDim.i)>dmax) Nmod.lasso.i<-(1:dim(lesBeta.i)[1])[(lesDim.i>=dmax)][1] lambda[[i]]<-res.i$penalty[1:Nmod.lasso.i]
active_set[[i]]<-active(res.i,in_s)[1:Nmod.lasso.i]
all_lambda<-c(all_lambda,lambda[[i]])
}
#         all_lambda<-all_lambda[2:(length(all_lambda))]
#         all_lambda<-sort(all_lambda,decreasing=TRUE)
all_lambda<-all_lambda[2:(length(all_lambda))]
L1 <- length(all_lambda)
all_lambda<-c(all_lambda,res.lasso$penalty[1:Nmod.lasso]) L<-length(all_lambda) rank.all_lambda<-L-rank(all_lambda)+1 sort.all_lambda<-sort(all_lambda,decreasing=TRUE) erreur<-matrix(0,V,L) for (i in 1:V) { est<-estime2(Y,X,I.CV[[i]],IM.CV[[i]],active_set[[i]]) ind<-1 for (l in 1:L) { if (lambda[[i]][ind]>sort.all_lambda[l]) { if (ind < length(lambda[[i]])) { ind<-ind+1 } } erreur[i,l]<-est[[ind]] } } resultat<-apply(erreur,2,mean) err.resultat<-sqrt(apply(erreur,2,var))/sqrt(V) resultat[resultat<10**(-10)] <- 0 lambda_cv<-sort.all_lambda[which.min(resultat)] ind<-indice(lambda_cv,res.lasso$penalty[1:Nmod.lasso])
rlm <- lm(Y~X[,I.lasso[[ind]]])
beta <- rlm$coef names(beta) <- c("Intercept",I.lasso[[ind]]) resultatF <- c(min(resultat), resultat[rank.all_lambda][(L1+1):L]) err.resultatF <- c(err.resultat[which.min(resultat)], err.resultat[rank.all_lambda][(L1+1):L]) les.lambda <- c(lambda_cv,res.lasso$penalty[1:Nmod.lasso])
crit <- resultatF[order(les.lambda,decreasing=TRUE)]
crit.error <- err.resultatF[order(les.lambda,decreasing=TRUE)]
#         ii <- length(resultat)
#         if (min(resultat)==0) ii <- which(resultat==0)[1]
res$Glasso$CV <- list(support=I.lasso[[ind]],
coeff=beta, fitted=rlm$fitted, crit=crit, crit.error=crit.error, lambda=sort(les.lambda,decreasing=TRUE)) } if (LINselect) { # Gauss lasso + LinSelect if (verbose) print(paste("Tuning Gauss Lasso with LinSelect: K=",K)) if (!calc.pen) { D <- 0:dmax dm <- pmin(D,rep(p/2,dmax+1)) Delta <- lgamma(p+1)-lgamma(dm+1)-lgamma(p-dm+1)+2*log(D+1) pen <- penalty(Delta, n, p, K) # pen <- penalite(p, n, K=K, dmax) } crit <- SCR + penSCR ind <- which.min(crit) rlm <- lm(Y~X[,I.lasso[[ind]]]) beta <- rlm$coef
names(beta) <- c("Intercept",I.lasso[[ind]])
res$Glasso$Ls <- list(support=I.lasso[[ind]],
coeff=beta, fitted=rlm$fitted+mu, crit=crit, lambda=res.lasso$penalty[1:Nmod.lasso] )
}
}
}
result <- res
return(result)
### A list with one or two components according to
### \code{method}. \cr
### \code{lasso} if \code{method} contains "lasso" is a list with  one or two components
### according to \code{Vfold} and \code{LINselect}.
###
### \itemize{
###   \item{\code{Ls} if \code{LINselect}=TRUE. A list with components
###     \itemize{
###       \item{\code{support}: vector of integers. Estimated support of the
###              parameter vector \eqn{\beta}.}
###       \item{\code{coef}: vector whose first component is the estimated
###             intercept. \cr The other components are the estimated non zero
###             coefficients.}
###       \item{\code{fitted}: vector with length n. Fitted value of
###       the response.}
###       \item{\code{crit}: vector containing the values of the criteria for
###       each value of \code{lambda}.}
###       \item{\code{lambda}: vector containing the values of the tuning
###       parameter of the lasso algorithm.}
###     }
###   }
###   \item{\code{CV} if \code{Vfold}=TRUE. A list with components
###     \itemize{
###       \item{\code{support}: vector of integers. Estimated support of the
###              parameter vector \eqn{\beta}.}
###       \item{\code{coef}: vector whose first component is the estimated
###             intercept. \cr The other components are the estimated non zero
###             coefficients.}
###       \item{\code{fitted}: vector with length n. Fitted value of
###       the response.}
###       \item{\code{crit}: vector containing the values of the criteria for
###       each value of \code{lambda}.}
###       \item{\code{crit.err}: vector containing the estimated
###       standard-error of the criteria.}
###       \item{\code{lambda}: vector containing the values of the tuning
###       parameter of the lasso algorithm.}
###     }
###   }
### }
### \code{Glasso} if \code{method} contains "Glasso". The same as \code{lasso}.
}


## Try the LINselect package in your browser

Any scripts or data that you put into this service are public.

LINselect documentation built on May 2, 2019, 4:09 a.m.