Wald tests - differences across packages

[strengejacke]

m1 <- glm(vs ~ disp + hp + drat, data = mtcars, family = "binomial")
m2 <- glm(vs ~ disp + hp, data = mtcars, family = "binomial")
m3 <- glm(vs ~ disp, data = mtcars, family = "binomial")

rez <- performance::test_wald(m1, m2, m3)
ref1 <- lmtest::waldtest(m1, m2, m3, test = "F")
ref2 <- anova(m1, m2, m3, test = "F")
#> Warning: using F test with a 'binomial' family is inappropriate

rez
#> Name | Model | df | df_diff |     F |     p
#> -------------------------------------------
#> m1   |   glm | 28 |         |       |      
#> m2   |   glm | 29 |   -1.00 |  8.44 | 0.007
#> m3   |   glm | 30 |   -1.00 | 12.94 | 0.001
#> Models were detected as nested and are compared in sequential order.
ref1
#> Wald test
#> 
#> Model 1: vs ~ disp + hp + drat
#> Model 2: vs ~ disp + hp
#> Model 3: vs ~ disp
#>   Res.Df Df      F  Pr(>F)  
#> 1     28                    
#> 2     29 -1 2.0444 0.16383  
#> 3     30 -1 2.9782 0.09504 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ref2
#> Analysis of Deviance Table
#> 
#> Model 1: vs ~ disp + hp + drat
#> Model 2: vs ~ disp + hp
#> Model 3: vs ~ disp
#>   Resid. Df Resid. Dev Df Deviance      F  Pr(>F)  
#> 1        28     12.869                             
#> 2        29     16.750 -1  -3.8802 3.8802 0.04886 *
#> 3        30     22.696 -1  -5.9462 5.9462 0.01475 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2022-10-30 with reprex v2.0.2}

@DominiqueMakowski @mattansb @bwiernik (or someone else)

Any ideas?

[strengejacke]

Might be due to the binomial family. For linear models, results are the same for test_wald() and anova() (and slightly differ for lmtest()):

m1 <- lm(mpg ~ disp + hp + drat, data = mtcars)
m2 <- lm(mpg ~ disp + hp, data = mtcars)
m3 <- lm(mpg ~ disp, data = mtcars)

rez <- performance::test_wald(m1, m2, m3)
ref1 <- lmtest::waldtest(m1, m2, m3, test = "F")
ref2 <- anova(m1, m2, m3, test = "F")

rez
#> Name | Model | df | df_diff |    F |     p
#> ------------------------------------------
#> m1   |    lm | 28 |         |      |      
#> m2   |    lm | 29 |   -1.00 | 3.33 | 0.079
#> m3   |    lm | 30 |   -1.00 | 3.72 | 0.064
#> Models were detected as nested and are compared in sequential order.
ref1
#> Wald test
#> 
#> Model 1: mpg ~ disp + hp + drat
#> Model 2: mpg ~ disp + hp
#> Model 3: mpg ~ disp
#>   Res.Df Df      F  Pr(>F)  
#> 1     28                    
#> 2     29 -1 3.3319 0.07863 .
#> 3     30 -1 3.4438 0.07368 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ref2
#> Analysis of Variance Table
#> 
#> Model 1: mpg ~ disp + hp + drat
#> Model 2: mpg ~ disp + hp
#> Model 3: mpg ~ disp
#>   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
#> 1     28 253.35                              
#> 2     29 283.49 -1   -30.148 3.3319 0.07863 .
#> 3     30 317.16 -1   -33.665 3.7207 0.06393 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2022-10-30 with reprex v2.0.2}

[mattansb]

Weird indeed... wouldn't a Wald test for a glm use a chi-sq?

[strengejacke]

anova and lmtest still differ for chisq, but the default behaviour of anova is just no test, I think:

m1 <- glm(vs ~ disp + hp + drat, data = mtcars, family = "binomial")
m2 <- glm(vs ~ disp + hp, data = mtcars, family = "binomial")
m3 <- glm(vs ~ disp, data = mtcars, family = "binomial")

lmtest::waldtest(m1, m2, m3, test = "Chisq")
#> Wald test
#> 
#> Model 1: vs ~ disp + hp + drat
#> Model 2: vs ~ disp + hp
#> Model 3: vs ~ disp
#>   Res.Df Df  Chisq Pr(>Chisq)  
#> 1     28                       
#> 2     29 -1 2.0444    0.15276  
#> 3     30 -1 2.9782    0.08439 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1, m2, m3, test = "Chisq")
#> Analysis of Deviance Table
#> 
#> Model 1: vs ~ disp + hp + drat
#> Model 2: vs ~ disp + hp
#> Model 3: vs ~ disp
#>   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
#> 1        28     12.869                       
#> 2        29     16.750 -1  -3.8802  0.04886 *
#> 3        30     22.696 -1  -5.9462  0.01475 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# default behaviour
anova(m1, m2, m3)
#> Analysis of Deviance Table
#> 
#> Model 1: vs ~ disp + hp + drat
#> Model 2: vs ~ disp + hp
#> Model 3: vs ~ disp
#>   Resid. Df Resid. Dev Df Deviance
#> 1        28     12.869            
#> 2        29     16.750 -1  -3.8802
#> 3        30     22.696 -1  -5.9462

^{Created on 2022-10-31 with reprex v2.0.2}