Comparing regression models: Different Tests

[DominiqueMakowski]

We can discuss here potential features regarding tests to compare different models.

Nested models

Two regression models are called nested if one contains all the predictors of the other, and some additional predictors.

• Wald test / F-test

The anova() function by default.

m1 <- lm(Sepal.Length ~ Petal.Width + Species, data=iris)
m2 <- lm(Sepal.Length ~ Species, data=iris)

# Nested
lmtest::waldtest(m1, m2)
#> Wald test
#> 
#> Model 1: Sepal.Length ~ Petal.Width + Species
#> Model 2: Sepal.Length ~ Species
#>   Res.Df Df     F    Pr(>F)    
#> 1    146                       
#> 2    147 -1 22.37 5.251e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2020-12-28 by the reprex package (v0.2.1)}

• LRT

m1 <- lm(Sepal.Length ~ Petal.Width + Species, data=iris)
m2 <- lm(Sepal.Length ~ Species, data=iris)

anova(m1, m2, test="LRT")
#> Analysis of Variance Table
#> 
#> Model 1: Sepal.Length ~ Petal.Width + Species
#> Model 2: Sepal.Length ~ Species
#>   Res.Df    RSS Df Sum of Sq  Pr(>Chi)    
#> 1    146 33.780                           
#> 2    147 38.956 -1   -5.1759 2.248e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2020-12-29 by the reprex package (v0.3.0)}

• Robust LRT

m1 <- lm(Sepal.Length ~ Petal.Width + Species, data=iris)
m2 <- lm(Sepal.Length ~ Species, data=iris)

nonnest2::vuongtest(m1, m2, nested=TRUE)
#> 
#> Model 1 
#>  Class: lm 
#>  Call: lm(formula = Sepal.Length ~ Petal.Width + Species, data = iris)
#> 
#> Model 2 
#>  Class: lm 
#>  Call: lm(formula = Sepal.Length ~ Species, data = iris)
#> 
#> Variance test 
#>   H0: Model 1 and Model 2 are indistinguishable 
#>   H1: Model 1 and Model 2 are distinguishable 
#>     w2 = 0.126,   p = 9.2e-05
#> 
#> Robust likelihood ratio test of distinguishable models 
#>   H0: Model 2 fits as well as Model 1 
#>   H1: Model 1 fits better than Model 2 
#>     LR = 21.384,   p = 1.12e-05

^{Created on 2020-12-28 by the reprex package (v0.2.1)}

• ...

Non-nested models

• BIC-based BF
  • Works on: same outcome variable
• Cox Test
  • Works on: All of the alternative models have the same error structure

The idea of the Cox test is the following: if the first model contains the correct set of regressors, then a fit of the regressors from the second model to the fitted values from first model should have no further explanatory value. But if it has, it can be concluded that model 1 does not contain the correct set of regressors.

m1 <- lm(Sepal.Length ~ Species, data=iris)
m2 <- lm(Sepal.Length ~ Petal.Width, data=iris)
m3 <- lm(Sepal.Length ~ Petal.Width + Species, data=iris)

lmtest::coxtest(m1, m2)
#> Cox test
#> 
#> Model 1: Sepal.Length ~ Species
#> Model 2: Sepal.Length ~ Petal.Width
#>                 Estimate Std. Error z value  Pr(>|z|)    
#> fitted(M1) ~ M2 -18.8530     3.5919 -5.2488 1.531e-07 ***
#> fitted(M2) ~ M1   0.5414     3.9130  0.1384    0.8899    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2020-12-28 by the reprex package (v0.2.1)}

• Davidson-MacKinnon J test

The idea of the J test is the following: if the first model contains the correct set of regressors, then including the fitted values of the second model into the set of regressors should provide no significant improvement. But if it does, it can be concluded that model 1 does not contain the correct set of regressors. Hence, to compare both models the fitted values of model 1 are included into model 2 and vice versa. The J test statistic is simply the marginal test of the fitted values in the augmented model.

m1 <- lm(Sepal.Length ~ Species, data=iris)
m2 <- lm(Sepal.Length ~ Petal.Width, data=iris)

lmtest::jtest(m1, m2)
#> J test
#> 
#> Model 1: Sepal.Length ~ Species
#> Model 2: Sepal.Length ~ Petal.Width
#>                 Estimate Std. Error t value  Pr(>|t|)    
#> M1 + fitted(M2)   1.0319    0.21817  4.7298 5.251e-06 ***
#> M2 + fitted(M1)  -0.0277    0.22551 -0.1228    0.9024    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2020-12-28 by the reprex package (v0.2.1)}

• Encompassing Test of Davidson & MacKinnon

The encompassing test fits an encompassing model which contains all regressors from both models such that the two models are nested within the encompassing model. A Wald test for comparing each of the models with the encompassing model is carried out by wald test.

m2 <- lm(Sepal.Length ~ Species, data=iris)
m3 <- lm(Sepal.Length ~ Petal.Width, data=iris)
lmtest::encomptest(m2, m3)
#> Encompassing test
#> 
#> Model 1: Sepal.Length ~ Species
#> Model 2: Sepal.Length ~ Petal.Width
#> Model E: Sepal.Length ~ Species + Petal.Width
#>           Res.Df Df       F    Pr(>F)    
#> M1 vs. ME    146 -1 22.3705 5.251e-06 ***
#> M2 vs. ME    146 -2  0.0748     0.928    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2020-12-28 by the reprex package (v0.2.1)}

• Non-nested LRT

See also nonnest2::icci() that computes CI around the difference of AIC and BIC, with the caveat that it only works 1) if models are non-nested and 2) if models are distinguishable (as shown by as significant Vuong's test).

m2 <- lm(Sepal.Length ~ Species, data=iris)
m3 <- lm(Sepal.Length ~ Petal.Width, data=iris)

nonnest2::vuongtest(m2, m3)
#> 
#> Model 1 
#>  Class: lm 
#>  Call: lm(formula = Sepal.Length ~ Species, data = iris)
#> 
#> Model 2 
#>  Class: lm 
#>  Call: lm(formula = Sepal.Length ~ Petal.Width, data = iris)
#> 
#> Variance test 
#>   H0: Model 1 and Model 2 are indistinguishable 
#>   H1: Model 1 and Model 2 are distinguishable 
#>     w2 = 0.119,   p = 0.000138
#> 
#> Non-nested likelihood ratio test 
#>   H0: Model fits are equal for the focal population 
#>   H1A: Model 1 fits better than Model 2 
#>     z = -2.509,   p = 0.994
#>   H1B: Model 2 fits better than Model 1 
#>     z = -2.509,   p = 0.006057

^{Created on 2020-12-28 by the reprex package (v0.2.1)}

• ...

Resources:

• Model Specification Tests Against Non-Nested Alternatives

[DominiqueMakowski]

Maybe we should have a function, in insight, that takes an arbitrary number of objects (find_class(...)) and returns the type, with classes like ListNestedModels and ListNonNestedModels. We could have the .default method for each test running this find_class() first, and dispatching on the appropriate method. It would makes it easier to throw warnings (for instance if ListNonNested is passed to a test that requires nested models). @strengejacke thoughts?

[DominiqueMakowski]

with useful attributes about the classes of each models etc.

[strengejacke]

I'm wondering if we should really call this test_performance_xxx()? Why not performance_vuong()? We have other tests that follow this naming pattern. Maybe the "full" test function can just be test_performance(), which internally calls performance_lrt(), performance_vuong() etc.?

[DominiqueMakowski]

yeah I thought about it, and was hesitating between performance_vuong (following performance_lrt) or test_vuong() (which would open a new category of test_*, from which test_performance() would be the wrapper)... I agree test_performance_vuong is not good.

Pros of performance_vuong;

• there's already performance_lrt

Cons:

• It's unclear that it's a test, rather than an index of perf per se

Pros of test_vuong:

• Opens a new test_* group (not especially a good point but juste a note)

Cons:

• Not obvious from the name that it's to test model's performances

[strengejacke]

there's already performance_lrt

Maybe we should indeed think of functions that compute model quality indices, that start with performance_*(), and functions that actually compare models and which could get a different prefix, maybe test_*().

[strengejacke]

Se here: https://easystats.github.io/performance/reference/index.html

We would have the section

Check Model Performance or Quality

and then maybe

Compare Model Fits

or something like that, with performance_lrt() being renamed and moved to that section.

[DominiqueMakowski]

Yes, and under "Comparing Models" there would be two categories like "descriptive comparison" or "relative comparison", which is what compare_performance mostly does, i.e., assemble indices together so what can eye-gauge them, and "comparison testing" which are actual "tests"