0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ If your hypothesis requires that all 5 variables are used, you may not need to reconcile the significance of individual variables. It may be obvious, but most of your variables are social constructs and as such, are unlikely to represent true biological or non-biological determinants of disease. $\endgroup$
    – Todd D
    Apr 5, 2022 at 2:08

1 Answer 1

0
+50
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Thank you for your answer. I have edited the question with a comparison of the coefficient estimates from the two models. As you can see, there is the expected correlation between race/ethnicity predictor and the language predictor. $\endgroup$ Mar 29, 2022 at 22:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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