Apa penyimpangan yang digunakan glmnet untuk membandingkan nilai ?

8

Salah satu kriteria untuk memilih nilai optimal dengan jaring elastis atau regresi yang dihukum serupa adalah untuk memeriksa sebidang penyimpangan terhadap rentang dan pilih ketika penyimpangan diminimalkan (atau dalam satu kesalahan standar dari minimum).λλλλ

Namun, saya mengalami kesulitan memahami apa, tepatnya, glmnetmenampilkan dengan plot.cv.glmnet, karena plot yang ditampilkan sama sekali tidak menyerupai hasil dari merencanakan penyimpangan terhadap .λ

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

masukkan deskripsi gambar di sini masukkan deskripsi gambar di sini

Tampaknya plot kedua tidak memasukkan penalti bersih elastis, dan juga diskalakan secara keliru secara vertikal. Saya mendasarkan klaim atas dasar bahwa bentuk kurva untuk nilai yang lebih besar dari menyerupai output. Namun, ketika saya mencoba untuk menghitung penalti sendiri, usaha saya juga tampaknya sangat tidak akurat.λ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)

Pertanyaan saya adalah: bagaimana cara seseorang menghitung penyimpangan yang dilaporkan dalam plot.cv.glmnetdiagram default ? Apa formulanya, dan apa yang saya lakukan salah dalam upaya saya menghitungnya?

Sycorax berkata Reinstate Monica
sumber
Anda sadar cv.glmnetsedang melakukan validasi silang 10 kali lipat, bukan? Jadi itu merencanakan rata-rata +/- 1 kesalahan standar penyimpangan pada data penahan 10%?
Andrew M
Saya sadar akan hal ini, ya.
Sycorax berkata Reinstate Monica

Jawaban:

6

Saya hanya ingin menambahkan input, tetapi tidak pada saat ini memiliki jawaban yang singkat dan terlalu lama untuk berkomentar. Semoga ini memberi lebih banyak wawasan.

Tampaknya fungsi yang menarik adalah di pustaka glmnet yang tidak dibongkar, dan disebut cv.lognet.R Sulit untuk melacak semuanya secara eksplisit, seperti halnya dalam kode S3 / S4, tetapi fungsi di atas terdaftar sebagai 'fungsi glmnet internal , 'digunakan oleh penulis dan tampaknya cocok dengan bagaimana cv.glmnet menghitung penyimpangan binomial.

Sementara saya tidak melihatnya di mana pun di koran, dari menelusuri kode glmnet ke cv.lognet, apa yang saya kumpulkan adalah bahwa ia menggunakan sesuatu yang disebut penyimpangan binomial tertutup yang dijelaskan di sini .

-[Ycatatan10(E)+(1-Y)catatan10(1-E)]

predmat adalah matriks dari nilai probabilitas capped (E, 1-E) output untuk setiap lambda, yang dibandingkan dengan nilai komplemen y dan y yang menghasilkan lp. Mereka kemudian dimasukkan ke dalam bentuk penyimpangan 2 * (ly-lp) dan dirata-rata selama lipatan tahan validasi silang untuk mendapatkan cvm - Kesalahan rata-rata yang divalidasi silang - dan rentang cv, yang telah Anda plot di gambar pertama.

Saya pikir fungsi penyimpangan manual (plot 2) tidak dihitung dengan cara yang sama seperti internal ini (plot 1).

    # 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)
menepuk
sumber
Terima kasih atas jawabannya, pat. Ini menjawab semua pertanyaan yang saya miliki tentang bagaimana prosedur ini bekerja, dan konsep statistik yang mendasarinya, bukan hanya perangkat lunaknya.
Sycorax berkata Reinstate Monica
2

Jadi saya mengunjungi situs CRAN dan mengunduh apa yang saya pikir adalah sumber paket glmnet . Di ./glmnet/R/plot.cv.glmnet.R tampaknya Anda akan menemukan kode sumber yang Anda cari. Ini cukup singkat jadi saya akan menempelkannya di sini, tetapi mungkin lebih baik jika Anda memeriksanya sendiri untuk memastikan bahwa itu memang kode yang sedang berjalan.

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()
}
Diego
sumber
1
Metode S3 sedikit tersembunyi di R, tetapi untuk melihat apa yang sedang berjalan Anda bisa mengetik getS3method('plot', 'cv.glmnet')tanpa harus repot mengunduh paket sumber. (Secara internal, glmnetbaru saja mendefinisikan fungsi yang dipanggil plot.cv.glmnettetapi belum mengekspornya. Anda masih dapat melihatnya dengan mengintip di dalam ruang nama dengan :::operator :) glmnet:::plot.cv.glmnet.
Andrew M
(+1) Terima kasih atas jawabannya, Diego. Ini menunjukkan saya ke arah yang benar, dan secara implisit menunjukkan di mana saya salah. Namun, saya akan menunda menerima untuk saat ini karena ini tidak menjawab pertanyaan statistik khusus (wakil pemrograman), yang dinyatakan di bagian bawah posting saya.
Sycorax berkata Reinstate Monica