# Compare estimates between non-nested mixed-effects models

I have two mixed-effect linear regression models, both fit on the same same data and same outcome, but with slightly different predictors. Is there a way to compare (i.e., generate a p-value for the comparison of) either (A) the model fits or (B) the regression coefficients of the predictor that differs between the models?

An example in R:

# load in data and packages
library(palmerpenguins)
library(lme4)
library(lmerTest)
penguins = na.omit(palmerpenguins::penguins)

# scale and center
penguins = penguins %>% mutate_at(.vars = c("bill_length_mm","bill_depth_mm","flipper_length_mm","body_mass_g"),.funs = scale,center=T,scale=T)

# fit mixed effect models
mod1 = lmer(flipper_length_mm  ~ bill_length_mm + body_mass_g + (1|species) + (1|island),data = penguins )
mod2 = lmer(flipper_length_mm  ~ bill_depth_mm  + body_mass_g + (1|species) + (1|island),data = penguins )


Model summaries

    # model summary
################################

summary(mod1)

Linear mixed model fit by REML. t-tests use Satterthwaites method ['lmerModLmerTest']
Formula: flipper_length_mm ~ bill_length_mm + body_mass_g + (1 | species) +      (1 | island)
Data: penguins

REML criterion at convergence: 316

Scaled residuals:
Min      1Q  Median      3Q     Max
-3.6437 -0.5624 -0.0317  0.6390  2.7633

Random effects:
Groups   Name        Variance Std.Dev.
species  (Intercept) 0.336373 0.57998
island   (Intercept) 0.005493 0.07412
Residual             0.140220 0.37446
Number of obs: 333, groups:  species, 3; island, 3

Fixed effects:
Estimate Std. Error        df t value Pr(>|t|)
(Intercept)    2.671e-03  3.386e-01 2.009e+00   0.008    0.994
bill_length_mm 1.914e-01  4.646e-02 3.273e+02   4.120  4.8e-05 ***
body_mass_g    3.877e-01  4.410e-02 3.274e+02   8.793  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) bll_l_
bll_lngth_m -0.031
body_mass_g  0.021 -0.581

##########################
summary(mod2)

Linear mixed model fit by REML. t-tests use Satterthwaites method ['lmerModLmerTest']
Formula: flipper_length_mm ~ bill_depth_mm + body_mass_g + (1 | species) +      (1 | island)
Data: penguins

REML criterion at convergence: 321.2

Scaled residuals:
Min      1Q  Median      3Q     Max
-2.7062 -0.5982 -0.0221  0.6285  3.3089

Random effects:
Groups   Name        Variance Std.Dev.
species  (Intercept) 0.695135 0.83375
island   (Intercept) 0.005245 0.07242
Residual             0.141861 0.37664
Number of obs: 333, groups:  species, 3; island, 3

Fixed effects:
Estimate Std. Error        df t value Pr(>|t|)
(Intercept)     0.03265    0.48383   1.94725   0.067 0.952480
bill_depth_mm   0.16046    0.04626 323.43083   3.469 0.000594 ***
body_mass_g     0.39210    0.04609 323.20298   8.507 6.77e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) bll_d_
bll_dpth_mm -0.008
body_mass_g  0.007 -0.620



Edit:

A couple of half-baked ideas I've had:

1. Compute a Z score the way we might if these were independent samples. The problem with this is that we can't necessarily assume that the covariance of the effects are independent, as these aren't independent samples. Would it for example be valid to put variables in a single model and pull the covariance from there?

  mod3 = lmer(flipper_length_mm  ~ bill_length_mm + bill_depth_mm+ body_mass_g + (1|species) + (1|island),data = penguins )

compare_z = function(b1,se1,b2,se2,cv){
z = (b1 - b2)/sqrt(se1^2 + se2^2 - 2*cv )
p =  2*pnorm(q = abs(z),mean = 0,sd = 1,lower.tail = FALSE)
return(p)
}

# compare effects of the two
compare_z(b1 = summary(mod1)$$coefficients[2,1], se1 = summary(mod1)$$coefficients[2,2],
b2 = summary(mod2)$$coefficients[2,1], se2 =summary(mod2)$$coefficients[2,2] ,
cv = mod3@vcov_beta[2,3] )
[1] 0.07955337

2. Another idea would be to compute the $$R^2$$ of the fixed effects, and associated confidence intervals, and then try to use the confidence intervals to infer a p-value. But, I'm not sure how to move from a CI to a p-value (besides p<0.05).

library(r2glmm)

r2glmm::r2beta(model = mod1,partial = FALSE,method = 'kr')
Effect   Rsq upper.CL lower.CL
1  Model 0.376     0.45    0.304
r2glmm::r2beta(model = mod2,partial = FALSE,method = 'kr')
Effect   Rsq upper.CL lower.CL
1  Model 0.371    0.446    0.299

• Is your goal interpretation or prediction? You may be able to use something as simple as percent of variance explained or compare model accuracy. Commented Jul 26 at 19:00
• The focus is more on interpretation, and ideally getting p-values for the comparison of the two. Commented Jul 26 at 19:34
• Everything in this dataset is quite connected. I could interpret the coefficients in mod3, but the ones in mod1/2, seem to me, to only make sense as prediction tools. What do you think, they mean? Commented Jul 29 at 14:30
• Please clarify the null hypothesis for which you want to calculate a p-value. The null hypothesis of your compare_z() function seems to be that the two coefficients have 0 difference in their values. I don't think that's what you really want to test. If your null hypothesis is that one or both of bill_length_mm or bill_depth_mm has no association with outcome, you can test that via summary(mod3). Also, why you would risk omitted-variable bias by using a model that omits a potentially outcome-associated predictor?
– EdM
Commented Jul 29 at 15:07
• Totally agree regarding the example. Again, no, in my actual data there are hundreds of groups (9,000 + individuals, longitudinal data, with random intercepts for participant ID, family ID, and data collection site ID). Commented Aug 1 at 17:16

You could use the model comparison function in anova for nested models and ask if each model is significantly better than the base model. In this case, both are significantly better, so you would pick the one with the smaller AIC or smaller BIC. They all point to model 1 with the bill length.

# load in data and packages
library(palmerpenguins, quietly = TRUE, warn.conflicts = FALSE)
library(lme4, quietly = TRUE, warn.conflicts = FALSE)
library(lmerTest, quietly = TRUE, warn.conflicts = FALSE)
library(dplyr, quietly = TRUE, warn.conflicts = FALSE)

penguins = na.omit(palmerpenguins::penguins)

# scale and center
penguins = penguins |> dplyr::mutate_at(.vars = c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g"),
.funs = scale, center=TRUE, scale=TRUE)

# fit mixed effect models
mod0 = lmer(flipper_length_mm  ~ body_mass_g + (1|species) + (1|island), data = penguins )
mod1 = lmer(flipper_length_mm  ~ bill_length_mm + body_mass_g + (1|species) + (1|island), data = penguins )
mod2 = lmer(flipper_length_mm  ~ bill_depth_mm  + body_mass_g + (1|species) + (1|island), data = penguins )

anova(mod1, mod0)
#> refitting model(s) with ML (instead of REML)
#> Data: penguins
#> Models:
#> mod0: flipper_length_mm ~ body_mass_g + (1 | species) + (1 | island)
#> mod1: flipper_length_mm ~ bill_length_mm + body_mass_g + (1 | species) + (1 | island)
#>      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
#> mod0    5 333.02 352.06 -161.51   323.02
#> mod1    6 318.36 341.21 -153.18   306.36 16.665  1   4.46e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(mod2, mod0)
#> refitting model(s) with ML (instead of REML)
#> Data: penguins
#> Models:
#> mod0: flipper_length_mm ~ body_mass_g + (1 | species) + (1 | island)
#> mod2: flipper_length_mm ~ bill_depth_mm + body_mass_g + (1 | species) + (1 | island)
#>      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
#> mod0    5 333.02 352.06 -161.51   323.02
#> mod2    6 324.27 347.12 -156.14   312.27 10.749  1   0.001043 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


## EDIT:

Although still not a p-value, you could ask in what percentage of bootstrap samples would you make the same decision, or make a histogram or confidence interval on the AIC differences. I do not have any academic references for an approach like this.

# load in data and packages
library(palmerpenguins, quietly = TRUE, warn.conflicts = FALSE)
library(lme4, quietly = TRUE, warn.conflicts = FALSE)
library(lmerTest, quietly = TRUE, warn.conflicts = FALSE)
library(dplyr, quietly = TRUE, warn.conflicts = FALSE)

penguins = na.omit(palmerpenguins::penguins)

# scale and center
penguins = penguins |> dplyr::mutate_at(.vars = c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g"),
.funs = scale, center=TRUE, scale=TRUE)

# fit mixed effect models
boot_fun <- function(X, i) {
mod0 = lmer(flipper_length_mm  ~ body_mass_g + (1|species) + (1|island), data = X[i,], REML = FALSE)
mod1 = lmer(flipper_length_mm  ~ bill_length_mm + body_mass_g + (1|species) + (1|island), data = X[i,], REML = FALSE )
mod2 = lmer(flipper_length_mm  ~ bill_depth_mm  + body_mass_g + (1|species) + (1|island), data = X[i,], REML = FALSE )

c(AIC(mod0), AIC(mod1), AIC(mod2))
}
# multiple singluar fits
suppressMessages({
b <- boot::boot(penguins, statistic = boot_fun, R = 200, strata = as.factor(paste0(penguins$$species, penguins$$island)))
})

length(which(b$$t[,2] < b$$t[,3])) / 200
#> [1] 0.69


In 69% of the bootstrap samples, AIC for model 1 is lower, and you would choose the bill_length_mm variable. You would fail to reject the Null hypothesis that the variables are equivalent based on this decision criterion.

• Yes I agree! If the models are fit with ML, not REML, then we can compare the AIC. Ideally I would like a p-value for the comparison of the two models. Commented Jul 29 at 13:42
• @David B As you can see, the anova method does switch the fit from REML to ML to produce the statistics and allow for the AIC comparison. If you run AIC(mod1) you get a different AIC for the REML estimates than what is shown in the anova for the ML estimates. Commented Jul 29 at 15:16
• Sorry I wasn't clear, I was agreeing with you that that is the correct way to get the AICs. Commented Jul 29 at 15:33
• Sorry @R Carnell, just saw your edit! What about bootstrapping the difference in the AIC between the two models? e.g., diff = b$t[,2] - b$t[,3]; quantile(diff,probs = c(0.025,0.975))? Commented Aug 2 at 13:49
• Agreed. Bootstrapping the difference works the same way. Commented Aug 2 at 14:12

@R Carnell gave a good answer, but I think I can still add a lot. The difference in AICs is a decently common and understood method for comparing non nested models. See the following links:

Now I'm going to mostly copy R Carnells setup, except I focus on the difference in AICs with $$d_{\text{observed}} = AIC_1 - AIC_2$$. Also I think R = 200 is too small for comfort.


library(palmerpenguins)
library(tidyverse)
library(lme4)

penguins = na.omit(palmerpenguins::penguins)

# scale and center
penguins = penguins |> dplyr::mutate_at(.vars = c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g"),
.funs = scale, center=TRUE, scale=TRUE)

mod0 = lmer(flipper_length_mm  ~ body_mass_g + (1|species) + (1|island), data = penguins )
mod1 = lmer(flipper_length_mm  ~ bill_length_mm + body_mass_g + (1|species) + (1|island), data = penguins )
mod2 = lmer(flipper_length_mm  ~ bill_depth_mm  + body_mass_g + (1|species) + (1|island), data = penguins )
anova(mod0, mod1)
anova(mod0, mod2)
d_obs <- anova(mod0, mod1)[2, "AIC"] - anova(mod0, mod2)[2, "AIC"]

# fit mixed effect models
boot_fun <- function(X, i) {
#mod0 = lmer(flipper_length_mm  ~ body_mass_g + (1|species) + (1|island), data = X[i,], REML = FALSE)
mod1 = lmer(flipper_length_mm  ~ bill_length_mm + body_mass_g + (1|species) + (1|island), data = X[i,], REML = FALSE )
mod2 = lmer(flipper_length_mm  ~ bill_depth_mm  + body_mass_g + (1|species) + (1|island), data = X[i,], REML = FALSE )

c(AIC(mod1), AIC(mod2))
}
# multiple singluar fits
suppressMessages({
# parallel might not work on windows
b <- boot::boot(penguins, statistic = boot_fun, R = 2000, strata = as.factor(paste0(penguins$$species, penguins$$island)),
parallel = "multicore", ncpus = 8)
})
d_boot <- b$$t[,1] - b$$t[,2]

boot_ci <- quantile(d_boot, c(0.025, 0.975))
c(d_obs, boot_ci)


We obtain $$\hat d_{\text{observed}} = -5.9$$ 95%-CI $$(-28.9, 12.5)$$. Now reading the linked answers, we see that AIC-differences bigger than 4 are not exactly small and 10 is considered large. The CI covering 10 in the opposite direction implies we are working with very low power.

Since you are working with same number of parameters in both models you probably don't have to bother with $$AIC_c$$ or BIC. Everything just breaks down to differences in deviance, which brings me to:

### Efficient approximation

Your model sounds large, so bootstrap might not be the thing for you. Differences in deviance like $$d_1 = Deviance_0 - Deviance_1, d_2 = Deviance_0 - Deviance_2$$ (We are flipping signs)and converge towards (non central) $$\chi^2$$-distributions(page 35) and $$d_{\text{observed}} =d_2 - d_1$$. Simulating $$d_1, d_2$$ and taking the difference is extremely cheap and shows very good results:

# flip for deviance
d1 <- diff(anova(mod0, mod1)[, "deviance"])
d2 <- diff(anova(mod0, mod2)[, "deviance"])

d_approx <- rchisq(10^5, df = 1, ncp = abs(d2)) - rchisq(10^5, df = 1, ncp = abs(d1))
lines(density(d_approx), col = 4)
approx_ci <- quantile(d_approx, c(0.025, 0.975))
approx_ci
hist(d_boot, freq = F, breaks = 40)
abline(v = c(d_obs, boot_ci), col = 2)
lines(density(d_approx), col = 4)


Note that the bootstrap isn't necessarily perfect with the convergence problems on this dataset. The alternative CI is $$(- 27.8, 14.7)$$. Basically the same.

We can even use this approach to construct a sampling distribution under the null assuming $$d_1, d_2$$ i.i.d. with $$ncp = (d_1 + d_2)/2$$, but I personally think p-values are inferior to CIs.

d_null <- rchisq(10^5, df = 1, ncp = abs(mean(d1, d2))) - rchisq(10^5, df = 1, ncp = abs(mean(d1, d2)))
p <- 2*min(c(mean(d_obs > d_null), mean(d_obs < d_null)))


PS: Since I see it a lot on this website. Random effects need a lot of groups: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#should-i-treat-factor-xxx-as-fixed-or-random

Thierry Onkelinx has a blog post with some simulations on the impact of the number of levels and concludes with a few recommendations for the number of levels of the grouping variable ns: > - get ns>1000 levels when an accurate estimate of the random effect variance is crucial. E.g. when a single number will be use for power calculations. > - get ns>100 levels when a reasonable estimate of the random effect variance is sufficient. E.g. power calculations with sensitivity analysis of the random effect variance. > - get ns>20 levels for an experimental study > - in case 10<ns<20 you should validate the model very cautious before using the output > - in case ns<10 it is safer to use the variable as a fixed effect.

• Interesting! So knowing the difference in deviance is enough to simulate a valid distribution that converges on what we would get with an actual bootstrap? That seems kind of surprising. Is there somewhere I can read more about it? Is this an approach that other people have used/published on? And yes, I'm well versed in mixed effect models. The example data is just an example, fit with two crossed random effects to get the idea across. Commented Aug 3 at 0:33
• @DavidB Yes and no, I've given a source for $d_1, d_2$ being noncentral $\chi^2$(also the null in the anova-function is $\chi^2$), but using the difference is my own logic. Also I had a mistake in defining $d_1, d_2$ with AIC, because mod0 has fewer parameters. You need to look directly at the deviance. Commented Aug 3 at 11:01
• One part of your answer that surprises me is that the difference between two bootstraps (or simulations), yields a distribution that's only a bit wider than bootstrapping the difference. Is that just a known property of bootstraps that I wasn't aware of? Or is something else going on here? Commented Aug 4 at 0:52
• @DavidB With the convergence issues and the less than infinite sample size the deviance of the bootstrap models won't be exactly $\chi^2$ distributed. If there are some particular properties of the bootstrap, I don't know about them. Commented Aug 4 at 17:27