8
$\begingroup$

One criterion for selecting the optimal value of $\lambda$ with an elastic net or similar penalized regression is to examine a plot of the deviance against the range of $\lambda$ and select $\lambda$ when deviance is minimized (or $\lambda$ within one standard error of the minimum).

However, I'm having difficulty understanding what, precisely, glmnet displays with plot.cv.glmnet, because the plot displayed does not at all resemble the results of plotting the deviance against $\lambda$.

set.seed(4567)
N       <- 500
P       <- 100
coefs   <- NULL
for(p in 1:P){
    coefs[p]    <- (-1)^p*100*2^(-p)
}
inv.logit <- function(x) exp(x)/(1+exp(x))
X   <- matrix(rnorm(N*P), ncol=P, nrow=N)
Y   <- rbinom(N, size=1, p=inv.logit(cbind(1, X)%*%c(-4, coefs)))
plot(test   <- cv.glmnet(x=X, y=Y, family="binomial", nfolds=10, alpha=0.8))
plot(log(test$lambda), deviance(test$glmnet.fit))

enter image description here enter image description here

It appears that the second plot does not incorporate the elastic net penalty, and is also incorrectly scaled vertically. I base the claim on the basis that the shape of the curve for larger values of $\lambda$ resembles that of the glmnet output. However, when I've attempted to compute the penalty on my own, my attempt likewise appears to be wildly inaccurate.

penalized.dev.fn    <- function(lambda, alpha=0.2, data, cv.model.obj){
    dev <- deviance(cv.model.obj$glmnet.fit)[seq_along(cv.model.obj$lambda)[cv.model.obj$lambda==lambda]]
    beta <- coef(cv.model.obj, s=lambda)[rownames(coef(cv.model.obj))!="(Intercept)"]
    penalty <- lambda * ( (1-alpha)/2*(beta%*%beta) + alpha*sum(abs(beta)) )
    penalized.dev <- penalty+dev
    return(penalized.dev)
}

out <- sapply(test$lambda, alpha=0.2, cv.model.obj=test, FUN=penalized.dev.fn)
    plot(log(test$lambda), out)

My question is: how does one manually compute the deviance reported in the default plot.cv.glmnet diagram? What is its formula, and what have I done wrong in my attempt to compute it?

$\endgroup$
2
  • $\begingroup$ You are aware that cv.glmnet is performing 10-fold cross validation, right? So it's plotting the mean +/- 1 standard error of the deviance on the 10% hold-out data? $\endgroup$ – Andrew M Jan 26 '15 at 3:43
  • $\begingroup$ I am aware of this, yes. $\endgroup$ – Sycorax Jan 26 '15 at 12:51
6
+50
$\begingroup$

I just wanted to add to the input, but don't at the moment have a concise answer and it's too long for a comment. Hopefully this gives more insight.

It seems that the function of interest is in the unpacked glmnet library, and is called cv.lognet.R It's hard to explicitly trace everything, as much is in S3/S4 code, but the above function is listed as an 'internal glmnet function,' used by the authors and seems to match how the cv.glmnet is calculating the binomial deviance.

While I didn't see it anywhere in the paper, from tracing the glmnet code to cv.lognet, what I gather is that it is using something called the capped binomial deviance described here.

$-[Y\log_{10}(E) + (1-Y)\log_{10}(1-E)]$

predmat is a matrix of the capped probability values (E, 1-E) output for each lambda, that are compared to the y and y's complement values resulting in lp. They are then put in the 2*(ly-lp) deviance form and averaged over cross-validated hold out folds to get cvm - The mean cross-validated error - and cv ranges, that you have plotted in the first image.

I think the manual deviance function (2nd plot) is not calculated the same way this internal one (1st plot) is.

    # from cv.lognet.R

    cvraw=switch(type.measure,
    "mse"=(y[,1]-(1-predmat))^2 +(y[,2]-predmat)^2,
    "mae"=abs(y[,1]-(1-predmat)) +abs(y[,2]-predmat),
    "deviance"= {
      predmat=pmin(pmax(predmat,prob_min),prob_max)
      lp=y[,1]*log(1-predmat)+y[,2]*log(predmat)
      ly=log(y)
      ly[y==0]=0
      ly=drop((y*ly)%*%c(1,1))
      2*(ly-lp)

   # cvm output
   cvm=apply(cvraw,2,weighted.mean,w=weights,na.rm=TRUE)
$\endgroup$
1
  • $\begingroup$ Thank you for the answer, pat. This addresses all of the questions that I had about how the procedure works, and the underlying statistical concepts, not just the software. $\endgroup$ – Sycorax Jan 31 '15 at 17:08
3
$\begingroup$

So I visited the CRAN site and downloaded what I think is the source of the glmnet package. In ./glmnet/R/plot.cv.glmnet.R it seems that you'd find the source code you're after. It's pretty brief so I'll paste here but it's probably best if you check it out yourself to be sure that it is indeed the code that is running.

plot.cv.glmnet=function(x,sign.lambda=1,...){
  cvobj=x
  xlab="log(Lambda)"
  if(sign.lambda<0)xlab=paste("-",xlab,sep="")
  plot.args=list(x=sign.lambda*log(cvobj$lambda),y=cvobj$cvm,ylim=range(cvobj$cvup,cvobj$cvlo),xlab=xlab,ylab=cvobj$name,type="n")
      new.args=list(...)
      if(length(new.args))plot.args[names(new.args)]=new.args
    do.call("plot",plot.args)
    error.bars(sign.lambda*log(cvobj$lambda),cvobj$cvup,cvobj$cvlo,width=0.01,col="darkgrey")
  points(sign.lambda*log(cvobj$lambda),cvobj$cvm,pch=20,col="red")
axis(side=3,at=sign.lambda*log(cvobj$lambda),labels=paste(cvobj$nz),tick=FALSE,line=0)
abline(v=sign.lambda*log(cvobj$lambda.min),lty=3)
    abline(v=sign.lambda*log(cvobj$lambda.1se),lty=3)
  invisible()
}
$\endgroup$
2
  • 1
    $\begingroup$ S3 methods are slightly hidden in R, but to see exactly what's being run you can type getS3method('plot', 'cv.glmnet') without going to the trouble of downloading the source package. (Internally, glmnet has just defined a function called plot.cv.glmnet but hasn't exported it. You can still see it by peeking inside the namepace with the ::: operator: glmnet:::plot.cv.glmnet). $\endgroup$ – Andrew M Jan 26 '15 at 3:47
  • $\begingroup$ (+1) Thank you for the answer, Diego. This points me in the right direction, and implicitly points out where I went wrong. However, I'm going to hold-off on accepting for the time being because this doesn't answer the specific statistical (vice programming) question, which is stated at the bottom of my post. $\endgroup$ – Sycorax Jan 26 '15 at 12:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.