Menggunakan glm () sebagai pengganti uji chi square sederhana

Saya tertarik untuk mengubah hipotesis nol menggunakan glm()R.

Sebagai contoh:

x = rbinom(100, 1, .7)  
summary(glm(x ~ 1, family = "binomial"))

menguji hipotesis bahwa . Bagaimana jika saya ingin mengubah null ke p = beberapa nilai arbitrer, di dalam ? $p = 0.5$$p$glm()

Saya tahu ini dapat dilakukan juga dengan prop.test()dan chisq.test(), tetapi saya ingin mengeksplorasi ide menggunakan glm()untuk menguji semua hipotesis yang berkaitan dengan data kategorikal.

Anda dapat menggunakan offset : glmdengan family="binomial"estimasi parameter pada log-odds atau skala logit, jadi sesuai dengan log-odds 0 atau probabilitas 0,5. Jika Anda ingin membandingkan dengan probabilitas p , Anda ingin nilai baseline menjadi q = logit ( p ) = log ( p / ( 1 - p ) ) . Model statistik sekarang$\beta_0=0$$p$$q = \textrm{logit}(p)=\log(p/(1-p))$

di mana hanya baris terakhir telah berubah dari pengaturan standar. Dalam kode R:

• gunakan offset(q)dalam formula
• fungsi logit / peluang-log adalah qlogis(p)
• sedikit mengganggu, Anda harus memberikan nilai offset untuk setiap elemen dalam variabel respons - R tidak akan secara otomatis mereplikasi nilai konstan untuk Anda. Ini dilakukan di bawah ini dengan mengatur kerangka data, tetapi Anda bisa menggunakannya rep(q,100).

x = rbinom(100, 1, .7)
dd <- data.frame(x, q = qlogis(0.7)) 
summary(glm(x ~ 1 + offset(q), data=dd, family = "binomial"))

Lihatlah interval kepercayaan untuk parameter GLM Anda:

> set.seed(1)
> x = rbinom(100, 1, .7)
> model<-glm(x ~ 1, family = "binomial")
> confint(model)
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.3426412 1.1862042 

Ini adalah interval kepercayaan untuk log-odds.

$p=0.5$$\log(odds) = \log \frac{p}{1-p} = \log 1 = 0$. So testing hypothesis that $p=0.5$ is equivalent to checking if confidence interval contains 0. This one does not, so hypothesis is rejected.

Now, for any arbitrary $p$, you can compute log-odds and check if it is inside confidence interval.

It is not (entirely) correct/accurate to use the p-values based on the z-/t-values in the glm.summary function as a hypothesis test.

1. Ini adalah bahasa yang membingungkan. Nilai yang dilaporkan bernama nilai-z. Tetapi dalam kasus ini mereka menggunakan estimasi kesalahan standar sebagai pengganti penyimpangan yang sebenarnya. Karena itu dalam kenyataannya mereka lebih dekat dengan nilai-t . Bandingkan tiga output berikut:
   1) ringkasan.glm
   2) uji-
   3 3) uji-z

   > set.seed(1)
   > x = rbinom(100, 1, .7)
   
   > coef1 <- summary(glm(x ~ 1, offset=rep(qlogis(0.7),length(x)), family = "binomial"))$coefficients
   > coef2 <- summary(glm(x ~ 1, family = "binomial"))$coefficients
   
   > coef1[4]  # output from summary.glm
   [1] 0.6626359
   > 2*pt(-abs((qlogis(0.7)-coef2[1])/coef2[2]),99,ncp=0) # manual t-test
   [1] 0.6635858
   > 2*pnorm(-abs((qlogis(0.7)-coef2[1])/coef2[2]),0,1) # manual z-test
   [1] 0.6626359

2. Mereka bukan nilai p yang tepat. Perhitungan tepat dari nilai-p menggunakan distribusi binomial akan bekerja lebih baik (dengan kekuatan komputasi saat ini, ini bukan masalah). Distribusi t, dengan asumsi distribusi kesalahan Gaussian, tidak tepat (terlalu tinggi p, melebihi tingkat alpha terjadi lebih jarang dalam "kenyataan"). Lihat perbandingan berikut:

   # trying all 100 possible outcomes if the true value is p=0.7
   px <- dbinom(0:100,100,0.7)
   p_model = rep(0,101)
   for (i in 0:100) {
     xi = c(rep(1,i),rep(0,100-i))
     model = glm(xi ~ 1, offset=rep(qlogis(0.7),100), family="binomial")
     p_model[i+1] = 1-summary(model)$coefficients[4]
   }
   
   
   # plotting cumulative distribution of outcomes
   outcomes <- p_model[order(p_model)]
   cdf <- cumsum(px[order(p_model)])
   plot(1-outcomes,1-cdf, 
        ylab="cumulative probability", 
        xlab= "calculated glm p-value",
        xlim=c(10^-4,1),ylim=c(10^-4,1),col=2,cex=0.5,log="xy")
   lines(c(0.00001,1),c(0.00001,1))
   for (i in 1:100) {
     lines(1-c(outcomes[i],outcomes[i+1]),1-c(cdf[i+1],cdf[i+1]),col=2)
   #  lines(1-c(outcomes[i],outcomes[i]),1-c(cdf[i],cdf[i+1]),col=2)
   }
   
   title("probability for rejection as function of set alpha level")

   

   Kurva hitam mewakili kesetaraan. Kurva merah di bawahnya. Itu berarti bahwa untuk nilai-p hitung yang diberikan oleh fungsi ringkasan glm, kami menemukan situasi ini (atau perbedaan yang lebih besar) lebih jarang dalam kenyataan daripada yang ditunjukkan oleh nilai-p.