R
似然比檢驗和 Wald 檢驗為 R 中的 glm 提供了不同的結論
我正在從Generalized, Linear, and Mixed Models複製一個示例。我的 MWE 如下:
Dilution <- c(1/128, 1/64, 1/32, 1/16, 1/8, 1/4, 1/2, 1, 2, 4) NoofPlates <- rep(x=5, times=10) NoPositive <- c(0, 0, 2, 2, 3, 4, 5, 5, 5, 5) Data <- data.frame(Dilution, NoofPlates, NoPositive) fm1 <- glm(formula=NoPositive/NoofPlates~log(Dilution), family=binomial("logit"), data=Data) summary(object=fm1)
輸出
Call: glm(formula = NoPositive/NoofPlates ~ log(Dilution), family = binomial("logit"), data = Data) Deviance Residuals: Min 1Q Median 3Q Max -0.38326 -0.20019 0.00871 0.15607 0.48505 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 4.174 2.800 1.491 0.136 log(Dilution) 1.623 1.022 1.587 0.112 (Dispersion parameter for binomial family taken to be 1) Null deviance: 8.24241 on 9 degrees of freedom Residual deviance: 0.64658 on 8 degrees of freedom AIC: 6.8563 Number of Fisher Scoring iterations: 6
代碼
anova(object=fm1, test="Chisq")
輸出
Analysis of Deviance Table Model: binomial, link: logit Response: NoPositive/NoofPlates Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 9 8.2424 log(Dilution) 1 7.5958 8 0.6466 0.00585 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
代碼
library(aod) wald.test(b=coef(object=fm1), Sigma=vcov(object=fm1), Terms=2)
輸出
Wald test: ---------- Chi-squared test: X2 = 2.5, df = 1, P(> X2) = 0.11
估計的係數與書中給出的結果完全匹配,但 SE 相距甚遠。基於 LRT 檢驗,斜率顯著,但基於 Wald 和 Z 檢驗斜率係數不顯著。我想知道我是否錯過了一些基本的東西。在此先感謝您的幫助。
主要問題是,如果您要使用比率作為響應變量,則應該使用
weights
參數。您一定忽略了有關“二項式 glm 中的非整數 #successes”的警告…Dilution <- c(1/128, 1/64, 1/32, 1/16, 1/8, 1/4, 1/2, 1, 2, 4) NoofPlates <- rep(x=5, times=10) NoPositive <- c(0, 0, 2, 2, 3, 4, 5, 5, 5, 5) Data <- data.frame(Dilution, NoofPlates, NoPositive) fm1 <- glm(formula=NoPositive/NoofPlates~log(Dilution), family=binomial("logit"), data=Data, weights=NoofPlates) coef(summary(fm1)) ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 4.173698 1.2522190 3.333042 0.0008590205 ## log(Dilution) 1.622552 0.4571016 3.549653 0.0003857398 anova(fm1,test="Chisq") ## Df Deviance Resid. Df Resid. Dev Pr(>Chi) ## NULL 9 41.212 ## log(Dilution) 1 37.979 8 3.233 7.151e-10 ***
LRT和Wald測試結果還是有很大區別的(-值對比),但出於實際目的,我們可以繼續說它們都非常重要……(在這種情況下(使用單個參數),
aod::wald.test()
給出與 完全相同的 p 值summary()
。)Wald vs profile 置信區間也有適度的不同,但 (0.7,2.5) (Wald) 和 (0.9, 2.75) (LRT) 的 CI [如下所示] 是否實際上不同取決於具體情況。
沃爾德:
confint.default(fm1) ## 2.5 % 97.5 % ## (Intercept) 1.7193940 6.628002 ## log(Dilution) 0.7266493 2.518455
輪廓:
confint(fm1) ## 2.5 % 97.5 % ## (Intercept) 2.2009398 7.267565 ## log(Dilution) 0.9014053 2.757092