Uji rasio kemungkinan dalam R

25

Misalkan saya akan melakukan regresi logistik univariat pada beberapa variabel independen, seperti ini:

mod.a <- glm(x ~ a, data=z, family=binominal("logistic"))
mod.b <- glm(x ~ b, data=z, family=binominal("logistic"))

Saya melakukan perbandingan model (uji rasio kemungkinan) untuk melihat apakah model lebih baik daripada model nol dengan perintah ini

1-pchisq(mod.a$null.deviance-mod.a$deviance, mod.a$df.null-mod.a$df.residual)

Lalu saya membangun model lain dengan semua variabel di dalamnya

mod.c <- glm(x ~ a+b, data=z, family=binomial("logistic"))

Untuk melihat apakah variabel signifikan secara statistik dalam model multivariat, saya menggunakan lrtestperintah dariepicalc

lrtest(mod.c,mod.a) ### see if variable b is statistically significant after adjustment of a
lrtest(mod.c,mod.b) ### see if variable a is statistically significant after adjustment of b

Saya bertanya-tanya apakah pchisqmetode dan lrtestmetode ini setara untuk melakukan tes loglikelihood? Seperti yang saya tidak tahu cara menggunakan lrtestuntuk model logistik univate.

lokheart
sumber
@ Gavin terima kasih telah mengingatkan saya, karena membandingkan dengan stackoverflow, saya perlu meluangkan lebih banyak waktu untuk "mencerna" jawabannya sebelum memutuskan apakah jawabannya sesuai atau tidak, bagaimanapun, terima kasih lagi.
lokheart
Saya tidak akan merekomendasikan menggunakan waldtest dari lmtest. Gunakan paket aod untuk pengujian model. Jauh lebih mudah. cran.r-project.org/web/packages/aod/aod.pdf
Mr. Nobody
epicalctelah dihapus ( sumber ). Sebuah alternatif bisa jadi lmtest.
Martin Thoma

Jawaban:

21

Pada dasarnya, ya, asalkan Anda menggunakan perbedaan yang benar dalam kemungkinan log:

> library(epicalc)
> model0 <- glm(case ~ induced + spontaneous, family=binomial, data=infert)
> model1 <- glm(case ~ induced, family=binomial, data=infert)
> lrtest (model0, model1)
Likelihood ratio test for MLE method 
Chi-squared 1 d.f. =  36.48675 , P value =  0 
> model1$deviance-model0$deviance
[1] 36.48675

dan bukan penyimpangan untuk model nol yang sama dalam kedua kasus. Jumlah df adalah jumlah parameter yang berbeda antara dua model bersarang, di sini df = 1. BTW, Anda dapat melihat kode sumber lrtest()hanya dengan mengetik

> lrtest

di prompt R.

chl
sumber
terima kasih, dan saya baru saja menemukan bahwa saya dapat menggunakan glm (output ~ NULL, data = z, keluarga = binomial ("logistik")) untuk membuat model NULL, dan jadi saya dapat menggunakan lrtest sesudahnya. FYI, terima kasih lagi
lokheart
2
@lokheart anova(model1, model0)akan bekerja juga.
chl
5
@lokheart glm(output ~ 1, data=z, family=binomial("logistic"))akan menjadi model nol yang lebih alami, yang mengatakan bahwa outputdijelaskan oleh suku konstan (intersep) / intersep tersirat dalam semua model Anda, jadi Anda menguji efek asetelah memperhitungkan intersep.
Pasang kembali Monica - G. Simpson
Atau Anda dapat melakukannya "secara manual": nilai p dari uji LR = 1-pchisq (penyimpangan, dof)
Umka
22

Alternatif adalah lmtestpaket, yang memiliki lrtest()fungsi yang menerima model tunggal. Berikut adalah contoh dari ?lrtestdalam lmtestpaket, yang untuk LM tetapi ada metode yang bekerja dengan GLM:

> require(lmtest)
Loading required package: lmtest
Loading required package: zoo
> ## with data from Greene (1993):
> ## load data and compute lags
> data("USDistLag")
> usdl <- na.contiguous(cbind(USDistLag, lag(USDistLag, k = -1)))
> colnames(usdl) <- c("con", "gnp", "con1", "gnp1")
> fm1 <- lm(con ~ gnp + gnp1, data = usdl)
> fm2 <- lm(con ~ gnp + con1 + gnp1, data = usdl)
> ## various equivalent specifications of the LR test
>
> ## Compare two nested models
> lrtest(fm2, fm1)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ gnp + gnp1
  #Df  LogLik Df  Chisq Pr(>Chisq)    
1   5 -56.069                         
2   4 -65.871 -1 19.605  9.524e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 
>
> ## with just one model provided, compare this model to a null one
> lrtest(fm2)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ 1
  #Df   LogLik Df  Chisq Pr(>Chisq)    
1   5  -56.069                         
2   2 -119.091 -3 126.04  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 
Pasang kembali Monica - G. Simpson
sumber
+1 Sangat menyenangkan untuk mengetahui (dan sepertinya saya lupa tentang paket itu).
chl
2
@ GavinSimpson Ini mungkin tampak konyol, tetapi bagaimana Anda menafsirkan hasil 'lrtest (fm2, fm1)'? Model 2 secara signifikan berbeda dari model 1 dan oleh karena itu penambahan variabel con1 berguna? Atau lrtest (fm2) mengatakan bahwa model 2 berbeda secara signifikan dari model 1? Tapi model mana yang lebih baik?
Kerry
5
@Kerry fm1memiliki kemungkinan log yang lebih rendah dan karenanya lebih buruk daripada fm2. LRT memberi tahu kami bahwa sejauh mana kami membuat fm1model yang lebih buruk daripada yang fm2tidak terduga, besar jika istilah yang berbeda di antara model itu berguna (dijelaskan tanggapannya). lrtest(fm2)tidak dibandingkan dengan fm1sama sekali, model fm2dibandingkan dengan dalam kasus itu jika, seperti yang dinyatakan dalam output, ini: con ~ 1. Model itu, model nol, mengatakan bahwa prediktor terbaik conadalah mean sampel dari con(istilah intersep / konstanta).
Pasang kembali Monica - G. Simpson