Bagaimana R menghitung nilai p untuk regresi binomial ini?

Pertimbangkan regresi binomial berikut:

# Create some data
set.seed(10)
n <- 500
x <- runif(n,0,100)
y <- x + rnorm(n,sd=100) < 0
 
# Fit a binomial regression model
model <- glm(y ~ x, family="binomial")

summary(model)

The summarymengembalikan fungsi p-nilai 1.03e-05. Saat menggunakan anova.glm, satu mendapat nilai p sedikit lebih ekstrim terlepas dari metode apa yang digunakan untuk menghitung nilai p.

anova(model, test="Rao")   # p.value = 7.5e-6
anova(model, test="LRT")   # p.value = 6.3e-6
anova(model, test="Chisq") # p.value = 6.3e-6

Apakah nilai p dari summaryfungsi berlaku untuk hipotesis yang sama dengan yang dikembalikan oleh anovafungsi? Jika ya, bagaimana summarymenghitung nilai-p ini dan apakah mungkin untuk melakukan perhitungan yang sama secara langsung anova?

Ini dapat membantu Anda membaca jawaban saya di sini: Mengapa nilai-p saya berbeda antara output regresi logistik, uji chi-square, dan interval kepercayaan untuk OR? Pertanyaan Anda di sini hampir merupakan duplikat dari itu, tetapi ada beberapa elemen tambahan dalam pertanyaan Anda yang dapat diatasi.

Seperti yang dicatat oleh @CliffAB, nilai-p dalam summary.glm()output berasal dari tes Wald. Ini analog dengan$t$-menguji koefisien untuk model linier dalam hal mereka adalah perbedaan antara nilai yang pas dari koefisien dan nilai referensi (dianggap sebagai $0$), dibagi dengan kesalahan standar. Perbedaannya adalah bahwa ini diambil untuk didistribusikan sebagai standar normal, bukan$t$. Di sisi lain, ini berlaku untuk sampel besar dan kami belum tentu tahu apa yang dimaksud dengan 'sampel besar' dalam kasus apa pun.

Menggunakan anova.glm()memberi Anda akses ke berbagai tes. Ketika Anda mengatur test="Rao", itu memberi Anda nilai p dari tes skor. Dan ketika Anda mengatur salah satu test="Chisq"atau test="LRT"(mereka sama), itu memberi Anda nilai p dari tes rasio kemungkinan.

The anova.glm()Fungsi tidak menguji hipotesis nol sama dengan uji Wald di summary()keluaran dalam kasus ini . Itu hanya karena model Anda hanya memiliki satu variabel. The anova.glm()fungsi akan melakukan tes berurutan, yang analog dengan 'tipe I SS' dalam pengaturan linear, sedangkan tes Wald dari summary()analog ke 'ketik III SS' dalam pengaturan linear (lihat jawaban saya di sini: Bagaimana menafsirkan tipe I, tipe II, dan tipe III ANOVA dan MANOVA? ). Mempertimbangkan:

x2 = rnorm(n)
m2 = glm(y~x+x2, family="binomial")
summary(m2)$coefficients
#                Estimate  Std. Error    z value     Pr(>|z|)
# (Intercept) -0.05906436 0.186876339 -0.3160612 7.519561e-01
# x           -0.01567551 0.003537183 -4.4316372 9.352029e-06
# x2          -0.05967796 0.099093504 -0.6022388 5.470152e-01
anova(m2, test="LRT")
# Terms added sequentially (first to last)
# 
#      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
# NULL                   499     619.10              
# x     1  20.3841       498     598.72 6.335e-06 ***
# x2    1   0.3627       497     598.35     0.547    
m3 = glm(y~x2+x, family="binomial")  # I just switched the order of x & x2 here
summary(m3)$coefficients
#                Estimate  Std. Error    z value     Pr(>|z|)
# (Intercept) -0.05906436 0.186876339 -0.3160612 7.519561e-01
# x2          -0.05967796 0.099093504 -0.6022388 5.470152e-01  # these are the same
# x           -0.01567551 0.003537183 -4.4316372 9.352029e-06  #  as above
anova(m3, test="LRT")
# Terms added sequentially (first to last)
# 
#      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
# NULL                   499     619.10              
# x2    1   0.1585       498     618.94    0.6906      # these differ from the
# x     1  20.5883       497     598.35 5.694e-06 ***  #  anova output above

Anda dapat memilih anova.glm()fungsi untuk memberi Anda skor dan kemungkinan rasio tes variabel individu dalam model regresi logistik ganda yang analog dengan 'tipe III SS', tetapi itu membosankan. Anda harus tetap memperbaiki model Anda sehingga setiap variabel pada gilirannya terdaftar terakhir dalam formula yang disediakan untuk glm()panggilan. Nilai p terakhir yang tercantum dalam anova.glm()output adalah yang akan dianalogikan dengan 'tipe III SS'.

Untuk mendapatkan tes rasio skor atau kemungkinan variabel individual lebih nyaman, gunakan drop1()saja. Mempertimbangkan:

drop1(m3, test="LRT")
# Single term deletions
# 
# Model:
# y ~ x2 + x
#        Df Deviance    AIC     LRT  Pr(>Chi)    
# <none>      598.35 604.35                      
# x2      1   598.72 602.72  0.3627     0.547      # the same as when x2 is last above
# x       1   618.94 622.94 20.5883 5.694e-06 ***  # the same as when x  is last above

Dalam R, summaryfungsi untuk glmmenghitung nilai p menggunakan statistik Wald sederhana, yaitu

$2 \times \Phi \left(\frac{-|\hat \beta|} { SE(\hat \beta) } \right)$

dimana $\hat \beta$ adalah parameter regresi yang menarik, $SE(\hat \beta)$ adalah kesalahan standar estimasi parameter regresi ini dan $\Phi$ adalah CDF dari distribusi normal standar.

Untuk membuat ulang ini dari output Anda, coba

 beta = coef(model)[2]
 # getting estimate
 B_SE = sqrt(vcov(model)[2,2])
 # extracting standard error
 pvalue =  pnorm(-abs(beta) / B_SE)  * 2
 # pvalue = 1.027859e-05