¿Qué desviación está usando glmnet para comparar valores de ?

Un criterio para seleccionar el valor óptimo de con una regresión penalizada neta elástica o similar es examinar una gráfica de la desviación contra el rango de y seleccionar cuando la desviación se minimiza (o dentro de un error estándar del mínimo).$\lambda$$\lambda$$\lambda$$\lambda$

Sin embargo, estoy teniendo dificultades para entender con qué se glmnetmuestra , precisamente, plot.cv.glmnetporque la trama que se muestra no se parece en absoluto a los resultados de trazar la desviación contra .$\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))

Parece que el segundo gráfico no incorpora la penalización neta elástica, y también se escala verticalmente incorrectamente. Baso la afirmación sobre la base de que la forma de la curva para valores más grandes de parece a la de la salida. Sin embargo, cuando intenté calcular la penalización por mi cuenta, mi intento también parece ser muy inexacto.$\lambda$glmnet

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)

Mi pregunta es: ¿cómo se calcula manualmente la desviación informada en el plot.cv.glmnetdiagrama predeterminado ? ¿Cuál es su fórmula y qué he hecho mal en mi intento de calcularla?

Solo quería agregar a la entrada, pero por el momento no tengo una respuesta concisa y es demasiado largo para un comentario. Esperemos que esto dé más información.

Parece que la función de interés está en la biblioteca glmnet desempaquetada, y se llama cv.lognet.R Es difícil rastrear explícitamente todo, ya que mucho está en el código S3 / S4, pero la función anterior aparece como una 'función glmnet interna , 'utilizado por los autores y parece coincidir con la forma en que cv.glmnet calcula la desviación binomial.

Si bien no lo vi en ninguna parte del documento, desde el rastreo del código glmnet hasta cv.lognet, lo que deduzco es que está usando algo llamado la desviación binomial limitada descrita aquí .

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

predmat es una matriz de la salida de valores de probabilidad limitados (E, 1-E) para cada lambda, que se comparan con los valores de complemento de y e y que resultan en lp. Luego se colocan en la forma de desviación 2 * (ly-lp) y se promedian sobre pliegues de retención con validación cruzada para obtener cvm - El error medio de validación cruzada - y rangos de cv, que usted ha trazado en la primera imagen.

Creo que la función de desviación manual (2º gráfico) no se calcula de la misma manera que esta función interna (1º gráfico).

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

Entonces visité el sitio de CRAN y descargué lo que creo que es la fuente del paquete glmnet . En ./glmnet/R/plot.cv.glmnet.R parece que encontrarás el código fuente que estás buscando. Es bastante breve, así que lo pegaré aquí, pero probablemente sea mejor si lo revisa usted mismo para asegurarse de que realmente se esté ejecutando el código.

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()
}