10
$\begingroup$

I'm reproducing an example from Generalized, Linear, and Mixed Models. My MWE is below:

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)

Output


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

Code


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

Output


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

Code


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

Output


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

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

Estimated coefficients are perfectly matching with the results given in the book but SE's are far apart. Based on LRT test the slope is significant but based on Wald and Z-test slope coefficient is insignificant. I wonder if I miss something basic. Thanks in advance for your help.

$\endgroup$
10
$\begingroup$

The main problem is that if you're going to use the ratio as your response variable, you should be using the weights argument. You must have ignored a warning about "non-integer #successes in a binomial glm" ...

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

The LRT and Wald test results are still quite different ($p$-values of $4 \times 10^{-4}$ vs. $7 \times 10^{-10}$), but for practical purposes we can go ahead say they're both strongly significant ... (In this case (with a single parameter), aod::wald.test() gives exactly the same p-value as summary().)

The Wald vs profile confidence intervals are also moderately different, but whether CIs [shown below] of (0.7,2.5) (Wald) and (0.9, 2.75) (LRT) are practically different depends on the particular situation.

Wald:

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

Profile:

confint(fm1)
##                   2.5 %   97.5 %
## (Intercept)   2.2009398 7.267565
## log(Dilution) 0.9014053 2.757092
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.