Wie berechnet R den p-Wert für diese binomiale Regression?

Betrachten Sie die folgende binomiale Regression:

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

Die summaryFunktion gibt einen p-Wert von zurück 1.03e-05. Bei Verwendung anova.glmerhält man etwas extremere p-Werte, unabhängig davon, mit welcher Methode der p-Wert berechnet wird.

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

Gilt der p-Wert der summaryFunktion für dieselbe Hypothese wie die von der anovaFunktion zurückgegebenen? Wenn ja, wie wurde summarydieser p-Wert berechnet und ist es möglich, dieselbe Berechnung direkt mit durchzuführen anova?

Es kann Ihnen helfen, meine Antwort hier zu lesen: Warum unterscheiden sich meine p-Werte zwischen der logistischen Regressionsausgabe, dem Chi-Quadrat-Test und dem Konfidenzintervall für den OP? Ihre Frage hier ist fast ein Duplikat davon, aber es gibt ein paar zusätzliche Elemente in Ihrer Frage, die angesprochen werden können.

Wie @CliffAB feststellt, stammen die p-Werte in der summary.glm()Ausgabe aus Wald-Tests. Diese sind insofern analog zu Tests von Koeffizienten für ein lineares Modell, als sie die Differenz zwischen dem angepassten Wert des Koeffizienten und dem Referenzwert (angenommen als ) geteilt durch den Standardfehler sind. Der Unterschied besteht darin, dass diese als Standardnormal anstelle von . Auf der anderen Seite gelten diese für große Stichproben, und wir wissen nicht unbedingt, was in einem bestimmten Fall eine „große Stichprobe“ ausmacht. $t$$0$$t$

Mit können anova.glm()Sie auf verschiedene Tests zugreifen. Wenn Sie einstellen test="Rao", erhalten Sie den p-Wert aus einem Score-Test. Und wenn Sie entweder test="Chisq"oder test="LRT"(sie sind gleich) einstellen , erhalten Sie den p-Wert aus einem Likelihood-Ratio-Test.

Die anova.glm()Funktion summary()testet in diesem Fall dieselbe Nullhypothese wie der Wald-Test in der Ausgabe . Das liegt nur daran, dass Ihr Modell nur eine Variable hat. Die anova.glm()Funktion führt sequentielle Tests durch, die in einer linearen Einstellung analog zu 'Typ I SS' sind, während die Wald-Tests in einer linearen Einstellung summary()analog zu 'Typ III SS' sind (siehe meine Antwort hier: Interpretation von Typ I, Typ II und Typ III ANOVA und MANOVA? ). Erwägen:

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

Sie können die anova.glm()Funktion nutzen, um Score- und Likelihood-Ratio-Tests einzelner Variablen in einem multiplen logistischen Regressionsmodell durchzuführen, das analog zu 'Typ III SS' ist, aber es ist langwierig. Sie müssten Ihr Modell weiter anpassen, damit jede Variable der Reihe nach zuletzt in der Formel aufgeführt wird, die für den glm()Aufruf bereitgestellt wird . Der letzte in der anova.glm()Ausgabe aufgeführte p-Wert ist analog zu 'Typ III SS'.

Verwenden Sie drop1()stattdessen, um die Score- oder Likelihood-Ratio-Tests einzelner Variablen bequemer zu erhalten. Erwägen:

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

In berechnet Rdie summaryFunktion für glmden p-Wert unter Verwendung einer einfachen Wald-Statistik, d. H.

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

$\hat \beta$$SE(\hat \beta)$$\Phi$

Versuchen Sie es, um dies aus Ihrer Ausgabe wiederherzustellen

 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