I'm trying to use machine learning to model the risk of healthcare-associated infections (HAI) for patients in a number of hospitals. I have variables both at the patient, ward and hospital level.

In particular, some variables at the hospital level are also present at the ward level. For example, milliliters of alcoholic hand gel used per patient day are collected both in each ward and for the whole hospital. Of course, those two variables may be correlated, but in a complex way.

Should I include both in the model or just the more proximal one (i.e. at the ward level) since the effect of the hospital variable is actually mediated by the ward one?


By modeling a predictor's effect at each level you get a more precise estimate of its association with the outcome at that level. Remember that a multilevel model splits the variance in the outcome across the N levels you specify (sounds like 3 in your case). So you have separate intercepts for patients, wards, and hospitals. Predictors at each of those levels explain variation at their respective level. If you have a patient-level predictor and you do not account for its mean (or other measure of it) at the higher levels, then your coefficient for that predictor will be a so-called conflated estimate. It is called such because it represents a blend of the predictor's association with the outcome at each of the levels of the hierarchy.

TL:DR, add the predictors at each level; this is especially useful if you are interested in predicting the outcome at any of the levels (i.e., the empirical Bayes prediction of the random intercept).

  • $\begingroup$ Thank you for the explanation. I have to add a level of complexity here: I'm using machine learning (xgboost) to model the outcome, to account for the non-linearity and interaction of the variables at various levels. My hope is that I have enough information (variables) at each level to account for the increased variability in effects and intercepts and provide a more generalizable model. Indeed I get a test set AUC of almost 90%. But my main goal is not prediction, but interpretation, possibly causal. I know that adding more variables increases fit, but what about interpretation? $\endgroup$ – Bakaburg Dec 7 '19 at 13:05
  • $\begingroup$ PS: for main effect extraction from the xgboost black box I use ALE plots, that try to extract the effect of a simulated increase in a small predictor space defined on the predictor range itself christophm.github.io/interpretable-ml-book/ale.html orhttps://arxiv.org/pdf/1612.08468.pdf $\endgroup$ – Bakaburg Dec 7 '19 at 13:08
  • $\begingroup$ You can model non-linearity and interactions in a traditional multilevel model and use model comparison (likelihood ratio tests, AIC, and BIC) to determine whether those provide a better fit to the data. These model should be quite interpretable. You can use the R package ggeffects() to plot non-linear effects and interaction terms. But the choice of model type is up to you. $\endgroup$ – Erik Ruzek Dec 9 '19 at 15:15
  • $\begingroup$ Yep, If I were modeling the data with a linear model I would definitely use random effects, but in this case I'm using a machine learning model (too many variables and possible complex interaction and non-linear effects to use a glm) $\endgroup$ – Bakaburg Dec 18 '19 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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