Reconciling LRT of two mixed-effects models with the Wald tests of individual predictors

Question

I hypothesize that a clinical outcome can be predicted by 5 categorical predictors and one numeric predictor. One of the predictors is the physician that cares for the patient. I have the data in the following structure in R:

str(ModelData)
#tibble [775 × 6] (S3: tbl_df/tbl/data.frame)
# $ Outcome       : Factor w/ 2 levels "No","Yes"
# $ Doctor        : Factor w/ 11 levels 
# $ Race          : Factor w/ 4 levels "Non-Hispanic White"
# $ Language      : Factor w/ 2 levels "English","Other"
# $ GestationalAge: num [1:694] 34.7 38 40 39 36.3 ...
# $ GeneticTesting: Factor w/ 2 levels "No","Yes"

I have created a mixed effects logistic regression model using the GLMMadaptive package with the physician as the random effect.

library(GLMMadaptive)
MM1 <- mixed_model(fixed = Outcome ~ Race + Language + GestationalAge + GeneticTesting,
                   random = ~ 1 | Doctor, data = ModelData, family = binomial(link = "logit"))

I hypothesize the language the patient's family prefers to speak is an important predictor of the outcome.

I created a second mixed effects model without language and find that the first model is significantly better at predicting the outcome:

MM2 <- mixed_model(fixed = Outcome ~ Race + GestationalAge + GeneticTesting,
                   random = ~ 1 | Doctor, data = ModelData, family = binomial(link = "logit"))
anova(MM1,MM2)
#       AIC    BIC log.Lik  LRT df p.value
#MM1 871.46 874.64 -427.73                
#MM2 874.81 877.59 -430.40 5.34  1  0.0208

However, the summary of the fixed effects does not indicate that Language = "Other" is significant.

summary(MM1)
#Call:
#mixed_model(fixed = Outcome ~ Race + Language + GestationalAge + 
#    GeneticTesting, random = ~1 | Doctor, data = ModelData, family = binomial(link = #"logit"))
#
#Data Descriptives:
#Number of Observations: 694
#Number of Groups: 11 
#
#Model:
# family: binomial
# link: logit 
#
#Fit statistics:
#   log.Lik      AIC      BIC
# -427.7306 871.4611 874.6443
#
#Random effects covariance matrix:
#               StdDev
#(Intercept) 0.7172494
#
#Fixed effects:
#                       Estimate Std.Err z-value p-value
#(Intercept)             -0.3230  0.9904 -0.3261 0.74436
#RaceHispanic or Latino   0.1602  0.2533  0.6322 0.52723
#RaceNon-Hispanic Black  -0.3452  0.2492 -1.3851 0.16603
#RaceOther                0.0227  0.2461  0.0921 0.92658
#LanguageOther           -0.4167  0.3344 -1.2461 0.21272
#GestationalAge          -0.0021  0.0254 -0.0843 0.93282
#GeneticTestingYes        1.3040  0.2099  6.2140 < 1e-04
#
#Integration:
#method: adaptive Gauss-Hermite quadrature rule
#quadrature points: 11
#
#Optimization:
#method: hybrid EM and quasi-Newton
#converged: TRUE 

Likewise, the marginal interpretation of the fixed effects are also not significant:

marginal_coefs(MM1, std_errors = TRUE)
#                       Estimate Std.Err z-value p-value
#(Intercept)             -0.2842  0.8980 -0.3164 0.75166
#RaceHispanic or Latino   0.1440  0.2204  0.6535 0.51343
#RaceNon-Hispanic Black  -0.3090  0.2224 -1.3896 0.16465
#RaceOther                0.0192  0.2273  0.0846 0.93254
#LanguageOther           -0.3760  0.2988 -1.2587 0.20814
#GestationalAge          -0.0020  0.0224 -0.0905 0.92792
#GeneticTestingYes        1.1727  0.1801  6.5100 < 1e-04

As noted in the answer by @EdM, it might be useful to explore the estimates between the two models:

map(list(MM1=MM1,MM2=MM2),fixef) %>%
   bind_rows(.id = "Model") %>% 
   pivot_longer(-"Model",names_to = "Predictor") %>% 
   pivot_wider(names_from="Model")
#  Predictor                   MM1      MM2
#  <chr>                     <dbl>    <dbl>
#1 (Intercept)            -0.323   -0.371  
#2 RaceHispanic or Latino  0.160    0.0535 
#3 RaceNon-Hispanic Black -0.345   -0.341  
#4 RaceOther               0.0227  -0.0495 
#5 LanguageOther          -0.417   NA      
#6 GestationalAge         -0.00214 -0.00103
#7 GeneticTestingYes       1.30     1.34   

My question is how do I reconcile the differences among the LRT test between the two models and the Wald and Monte Carlo integration? More importantly, how do I effectively present these nuances in a publication?

I have reviewed this Q&A, this Q&A, and this answer, but I am still unclear on the best practice.

Answer 1

I suspect that this has to do in part with the inherent omitted-variable bias in logistic regression. Leaving out any predictor associated with outcome, even if it isn't correlated with the included predictors, tends to bias the coefficient estimates of included predictors toward 0.

Although the coefficient for Language itself doesn't pass the arbitrary p < 0.05 criterion for "statistical significance," its magnitude (-0.42) is sufficient to support some association with outcome. Including it might improve the estimates for the other predictors and thus leads to a model that is "statistically significantly" better overall. It would be interesting to examine the coefficient estimate differences for the other predictors between the two models.

The difference between Wald and likelihood-ratio tests, with the first used for the coefficient estimate and the second for the difference between the models, might also come into play here. If the likelihood profile doesn't have the ideal quadratic shape, the likelihood-ratio test is more reliable.