Bagaimana cara menghitung pseudo-

46

Tulisan Christopher Manning pada regresi logistik di R menunjukkan regresi logistik di R sebagai berikut:

ced.logr <- glm(ced.del ~ cat + follows + factor(class), 
  family=binomial)

Beberapa output:

> summary(ced.logr)
Call:
glm(formula = ced.del ~ cat + follows + factor(class),
    family = binomial("logit"))
Deviance Residuals:
Min            1Q    Median       3Q      Max
-3.24384 -1.34325   0.04954  1.01488  6.40094

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.31827    0.12221 -10.787 < 2e-16
catd          -0.16931    0.10032  -1.688 0.091459
catm           0.17858    0.08952   1.995 0.046053
catn           0.66672    0.09651   6.908 4.91e-12
catv          -0.76754    0.21844  -3.514 0.000442
followsP       0.95255    0.07400  12.872 < 2e-16
followsV       0.53408    0.05660   9.436 < 2e-16
factor(class)2 1.27045    0.10320  12.310 < 2e-16
factor(class)3 1.04805    0.10355  10.122 < 2e-16
factor(class)4 1.37425    0.10155  13.532 < 2e-16
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 958.66 on 51 degrees of freedom
Residual deviance: 198.63 on 42 degrees of freedom
AIC: 446.10
Number of Fisher Scoring iterations: 4

Dia kemudian masuk ke beberapa detail tentang bagaimana menafsirkan koefisien, membandingkan model yang berbeda, dan sebagainya. Cukup berguna.

Namun, berapa banyak varian yang diperhitungkan oleh model? Sebuah halaman Stata pada regresi logistik mengatakan:

Secara teknis, tidak dapat dihitung dengan cara yang sama dalam regresi logistik seperti di regresi OLS. Pseudo- , dalam regresi logistik, didefinisikan sebagai , di mana mewakili kemungkinan log untuk model "hanya-konstan" dan adalah kemungkinan log untuk model penuh dengan konstan dan prediktor.R2 1 - L 1R2 L0L11L1L0L0L1

Saya mengerti ini di level tinggi. Model hanya-konstan akan tanpa parameter (hanya istilah intersep). Log likelihood adalah ukuran seberapa dekat parameter tersebut sesuai dengan data. Bahkan, Manning semacam petunjuk bahwa penyimpangan yang mungkin . Mungkin penyimpangan nol hanya konstan dan penyimpangan residual adalah dari model? Namun, saya tidak sejelas itu.- 2 log L2logL2logL

Dapatkah seseorang memverifikasi bagaimana seseorang benar-benar menghitung pseudo- dalam R menggunakan contoh ini?R2

makan siang
sumber
5
Halaman komputasi statistik UCLA yang biasanya luar biasa telah membuat kesalahan langka di sini - seharusnya tidak ada tanda kurung dalam ekspresi untuk pseudo- , yaitu seharusnya . (Maaf karena tidak menjawab pertanyaan Anda karena saya akan pergi tidur - saya yakin orang lain akan menjawab ini sebelum saya cukup bangun untuk melakukannya.) 1 - L 1 / L 0R21L1/L0
onestop
3
Halaman ini membahas beberapa pseudo-R ^ 2s.
dfrankow
2
Catatan: pertanyaan terkait tidak menyukai pseudo-R ^ 2s, tetapi lebih memilih cross-validation atau prediksi tes penahan.
dfrankow

Jawaban:

49

Jangan lupa paket rms , oleh Frank Harrell. Anda akan menemukan semua yang Anda butuhkan untuk pemasangan dan validasi GLM.

Berikut adalah contoh mainan (dengan hanya satu prediktor):

set.seed(101)
n <- 200
x <- rnorm(n)
a <- 1
b <- -2
p <- exp(a+b*x)/(1+exp(a+b*x))
y <- factor(ifelse(runif(n)<p, 1, 0), levels=0:1)
mod1 <- glm(y ~ x, family=binomial)
summary(mod1)

Ini menghasilkan:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.8959     0.1969    4.55 5.36e-06 ***
x            -1.8720     0.2807   -6.67 2.56e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 258.98  on 199  degrees of freedom
Residual deviance: 181.02  on 198  degrees of freedom
AIC: 185.02

Sekarang, menggunakan lrmfungsinya,

require(rms)
mod1b <- lrm(y ~ x)

Anda segera mendapatkan banyak indeks kesesuaian model, termasuk Nagelkerke , dengan :R2print(mod1b)

Logistic Regression Model

lrm(formula = y ~ x)

                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       

Obs           200    LR chi2      77.96    R2       0.445    C       0.852    
 0             70    d.f.             1    g        2.054    Dxy     0.705    
 1            130    Pr(> chi2) <0.0001    gr       7.801    gamma   0.705    
max |deriv| 2e-08                          gp       0.319    tau-a   0.322    
                                           Brier    0.150                     


          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept  0.8959 0.1969  4.55  <0.0001 
x         -1.8720 0.2807 -6.67  <0.0001 

Di sini, R2=0.445 dan dihitung sebagai (1exp(LR/n))/(1exp((2L0)/n)) , di mana LR adalah stat χ2 (membandingkan dua model bersarang yang Anda uraikan), sedangkan penyebutnya hanya nilai maks untuk . Untuk model yang sempurna, kita harapkanR2LR=2L0R2=1

Dengan tangan,

> mod0 <- update(mod1, .~.-x)
> lr.stat <- lrtest(mod0, mod1)
> (1-exp(-as.numeric(lr.stat$stats[1])/n))/(1-exp(2*as.numeric(logLik(mod0)/n)))
[1] 0.4445742
> mod1b$stats["R2"]
       R2 
0.4445742 

R2R2c

chl
sumber
Bisakah Anda jelaskan bagaimana Anda mendapatkan 0,445? Saya menggunakan 1-exp (-77,96 / 200) tetapi saya mendapat 0,323. Apa yang saya lakukan salah? Terima kasih.
2
Yang mana Nagelkerke R2?
JetLag
1
@ JetLag Di Bawah Indeks Diskriminasi, Nagelkerke disingkat R2 (yaitu 0,445). Anda dapat memeriksanya menggunakan fungsi NagelkerkeR2 () dari paket fmsb.
Chernoff
11

R2

pengguna48729
sumber
7

R2 :

R2RM2=1lnL^fulllnL^nulllnL^fulllnL^full adalah log-kemungkinan model dengan hanya intercept.

R2

  1. deviance=2ln(Lfull)null.deviance=2ln(Lnull)

    pR2 = 1 - mod$deviance / mod$null.deviance # works for glm

Tetapi pendekatan di atas tidak bekerja untuk Pseudo sampelR2

  1. Gunakan fungsi "logLik" dalam R dan definisi (juga berfungsi untuk sampel)

    mod_null <- glm(y~1, family = binomial, data = insample) 1- logLik(mod)/logLik(mod_null)

R2

Contoh:

out-of-sample pseudo-R

R2

Rp2=1Lest.outLnull.out,
Lest.outLnull.out

Kode:

pred.out.link <- predict(mod, outSample, type = "link") mod.out.null <- gam(Default~1, family = binomial, data = outSample) pR2.out <- 1 - sum(outSample$y * pred.out.link - log(1 + exp(pred.out.link))) / logLik(mod.out.null)

Xiaorui Zhu
sumber
deviance=2ln(Lfull)model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp, data = esoph, family = binomial)model1$deviance-2*logLik(model1)
6

jika penyimpangan proporsional dengan kemungkinan log, dan seseorang menggunakan definisi (lihat misalnya McFadden di sini )

pseudo R^2 = 1 - L(model) / L(intercept)

R21198.63958.66

Pertanyaannya adalah: apakah penyimpangan dilaporkan sebanding dengan kemungkinan log?

dfrankow
sumber
3
Pseudo-R ^ 2 ini sama sekali tidak setuju dengan jawaban Nagelkerke R ^ 2 dari @ chl.
dfrankow
Penyimpangan didefinisikan sebagai -2 * LL ketika saya masih di sekolah.
DWin
@dfrankow tidak setuju, karena Nagelkerke adalah normalisasi Cox dan Snell R2, yang berbeda dari McFaddens R2.
colin
0

R2R2=1llfullllconstantllfullllconstant adalah log-kemungkinan data uji dengan model dengan hanya konstanta yang dipasang pada set pelatihan, dan kemudian menggunakan konstanta yang dipasang untuk memprediksi pada set pengujian menghitung probabilitas dan karenanya mendapatkan kemungkinan log.

R2R2=1i(yiy^i)2i(yiy¯train)2i(yiy¯train)2y¯traini(yiβ0)2β^0=y¯trainR2R2

cthraves
sumber