Verwendung von glm () als Ersatz für einen einfachen Chi-Quadrat-Test

Ich bin daran interessiert, die Nullhypothesen glm()in R zu ändern .

Beispielsweise:

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

testet die Hypothese, dass . Was ist, wenn ich die Null in = einen beliebigen Wert innerhalb von ändern möchte ? $p = 0.5$$p$glm()

Ich weiß, dass dies auch mit prop.test()und möglich ist chisq.test(), aber ich möchte die Idee untersuchen, glm()mit allen Hypothesen, die sich auf kategoriale Daten beziehen, einen Test durchzuführen.

Sie können einen Versatz verwenden : glmmit family="binomial"Schätzparametern auf der Log-Odds- oder Logit-Skala, sodass einer Log-Odds von 0 oder einer Wahrscheinlichkeit von 0,5 entspricht. Wenn Sie mit einer Wahrscheinlichkeit von vergleichen möchten , möchten Sie, dass der Grundlinienwert . Das statistische Modell ist jetztp q = logit ( p ) = log ( p / ( 1 - p ) $\beta_0=0$$p$$q = \textrm{logit}(p)=\log(p/(1-p))$

wobei sich nur die letzte Zeile gegenüber der Standardeinstellung geändert hat. Im R-Code:

• Verwenden Sie offset(q)in der Formel
• Die Logit / Log-Odds-Funktion ist qlogis(p)
• Etwas ärgerlich ist, dass Sie für jedes Element in der Antwortvariablen einen Versatzwert angeben müssen - R repliziert für Sie nicht automatisch einen konstanten Wert. Dies wird im Folgenden durch das Einrichten eines Datenrahmens durchgeführt, aber Sie könnten es einfach verwenden 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"))

Sehen Sie sich das Konfidenzintervall für die Parameter Ihres GLM an:

> 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 

Dies ist ein Konfidenzintervall für Log-Odds.

$p=0.5$$\log(odds) = \log \frac{p}{1-p} = \log 1 = 0$$p=0.5$

$p$

Es ist nicht (vollständig) korrekt / genau, die p-Werte, die auf den z- / t-Werten in der Funktion glm.summary basieren, als Hypothesentest zu verwenden.

1. Das ist verwirrende Sprache. Die angegebenen Werte werden als z-Werte bezeichnet. In diesem Fall wird jedoch der geschätzte Standardfehler anstelle der tatsächlichen Abweichung verwendet. In Wirklichkeit sind sie daher näher an t-Werten . Vergleichen Sie die folgenden drei Ausgaben:
   1) summary.glm
   2) t-test
   3) z-test

   > 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. Sie sind keine exakten p-Werte. Eine genaue Berechnung des p-Wertes mit der Binomialverteilung würde besser funktionieren (mit der Rechenleistung ist dies heutzutage kein Problem mehr). Die t-Verteilung unter der Annahme einer Gaußschen Verteilung des Fehlers ist nicht genau (sie überschätzt p, ein Überschreiten des Alpha-Niveaus tritt in der "Realität" seltener auf). Siehe folgenden Vergleich:

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

   

   Die schwarze Kurve steht für Gleichheit. Die rote Kurve ist darunter. Dies bedeutet, dass wir für einen gegebenen berechneten p-Wert durch die glm-Zusammenfassungsfunktion diese Situation (oder einen größeren Unterschied) in der Realität seltener finden, als der p-Wert angibt.