Lo que la desviación es glmnet utiliza para comparar los valores de $\lambda$?

Question

Uno de los criterios para seleccionar el valor óptimo de $\lambda$ con una red elástica o similar penalizado de regresión es examinar una parcela de la desviación en contra de la variedad de $\lambda$ y seleccione $\lambda$ cuando la desviación se minimiza (o $\lambda$ dentro de un error estándar de la mínima).

Sin embargo, estoy teniendo dificultades para entender lo que, precisamente, glmnet muestra con plot.cv.glmnet, debido a que el gráfico muestra no se parece en nada a los resultados de trazado de la desviación en contra de $\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 la segunda trama no se incorpora a la red elástica de la pena, y es también la escala incorrecta verticalmente. Me la base de la demanda sobre la base de que la forma de la curva para valores más grandes de $\lambda$ asemeja a la de la glmnet de salida. Sin embargo, cuando he intentado calcular la multa en el mío propio, mi intento asimismo parece ser imprecisos.

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 hace manualmente calcular las desviaciones reportadas en el valor predeterminado plot.cv.glmnet diagrama? ¿Cuál es su fórmula, y lo que he hecho mal en mi intento de calcular?

Answer 1

Yo sólo quería añadir a la entrada, pero no en el momento de tener una respuesta sucinta y es demasiado largo para un comentario. Esperemos que esto le da una visión más clara.

Parece que la función de interés está en el desempaquetado glmnet de la biblioteca, y se llama cv.lognet.R Es difícil explícitamente un seguimiento de todo, como en S3/S4 código, pero la función anterior está en la lista interna de glmnet función,' utilizado por los autores y parece coincidir cómo el cv.glmnet es calcular el binomio de la desviación.

Mientras no me vea en cualquier lugar en el documento de seguimiento de la glmnet código de cv.lognet, lo que deduzco es que es usar algo que se llama el máximo de binomio de la desviación se describe aquí.

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

predmat es una matriz de una de las cubiertas de los valores de la probabilidad (E, 1-E) de salida para cada lambda, que son comparados con los de la y y de y complementar los valores resultantes de la lp. Luego se ponen en la 2*(ly-lp) la desviación forma y promediados de la cruz-validado mantener los pliegues para conseguir cvm - La media de la cruz-validado error - y el cv de los rangos, que se han trazado en la primera imagen.

Creo que el manual de la desviación de la función (2ª parcela) no se calcula de la misma manera interno de uno (1 de parcela).

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

Answer 2

Así que me metí en el CRAN sitio y descargar lo que yo creo que es la fuente de la glmnet paquete. En ./glmnet/R/parcela.cv.glmnet.R parece que te gustaría encontrar en el código fuente que usted está después. Es muy breve, así que me voy a pegar aquí, pero es probablemente mejor si usted echa un vistazo a ti mismo para asegurarse de que es de hecho el código que se está ejecutando.

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