La prueba de razón de probabilidad y la prueba de Wald proporcionan una conclusión diferente para glm en R

11

Estoy reproduciendo un ejemplo de modelos generalizados, lineales y mixtos . Mi MWE está abajo:

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)

Salida


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

Código


anova(object=fm1, test="Chisq")

Salida


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

Código


library(aod)
wald.test(b=coef(object=fm1), Sigma=vcov(object=fm1), Terms=2)

Salida


Wald test:
----------

Chi-squared test:
X2 = 2.5, df = 1, P(> X2) = 0.11

Los coeficientes estimados coinciden perfectamente con los resultados dados en el libro, pero los SE están muy separados. Según la prueba LRT, la pendiente es significativa, pero según el coeficiente de pendiente de Wald y la prueba Z es insignificante. Me pregunto si extraño algo básico. Gracias de antemano por tu ayuda.

MYaseen208
fuente

Respuestas:

12

El principal problema es que si va a usar la razón como su variable de respuesta, debería usar el weightsargumento. Debe haber ignorado una advertencia sobre "#cesos no enteros en una película binomial" ...

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 ***

Los resultados de la prueba LRT y Wald siguen siendo bastante diferentes ( valores de frente a ), pero a efectos prácticos podemos seguir diciendo que ambos fuertemente significativo ... (en este caso (con un solo parámetro), proporciona exactamente el mismo valor p que ).4 × 10 - 4 7 × 10 - 10pag4 4×10-4 47 7×10-10aod::wald.test()summary()

Los intervalos de confianza de Wald vs perfil también son moderadamente diferentes, pero si los IC [mostrados a continuación] de (0.7,2.5) (Wald) y (0.9, 2.75) (LRT) son prácticamente diferentes depende de la situación particular.

Wald:

confint.default(fm1)
##                   2.5 %   97.5 %
## (Intercept)   1.7193940 6.628002
## log(Dilution) 0.7266493 2.518455

Perfil:

confint(fm1)
##                   2.5 %   97.5 %
## (Intercept)   2.2009398 7.267565
## log(Dilution) 0.9014053 2.757092
Ben Bolker
fuente