Modelling strategies for analyzing an effect of a predictor through higher hierarchical level What strategies can be considered when a predictor's direct effect can not be measured directly due to unmeasured confounding? However, data has a hierarchical structure (patients within regions) that may solve some of the problem.
We have patients' data including three variables:

*

*complications (outcome),

*rehabilitation (received hours but highly zero-inflated variable),

*region (patient's place of residence).

A model like this gives biased result due to unmeasured confounding (unavailable variables): e.g. very sick patients are excluded from rehabilitation and very fit patients do not need much rehabilitation.
complications ~ rehabilitation + region

However, patients in different regions are relatively similar and I do know that regional disparities in rehabilitation exist. Can I specify the model in a way that I examine rehabilitation's regional variability on complications?

*

*Hierarchical modelling?
complications ~ rehabilitation + (rehabilitation | region)


*Hierarchical modelling analysing correlations between intercepts/rehabilitation?
complications ~ rehabilitation + (rehabilitation |c| region)


*Cross Classification Modelling?
complications ~ 1 + (1 + rehabilitation) + (1 | region)


*Other strategies that I am not aware of.
PS! All ideas, including Bayesian solutions, are also very welcome.
 A: Regarding the first point

A model like this gives biased result due to unmeasured confounding (unavailable variables): e.g. very sick patients are excluded from rehabilitation and very fit patients do not need much rehabilitation.

No modelling solution will be able to "fix" this. You'll want to state this as a limitation of your study and generalize only to the population of patients who "qualify" for rehabilitation.
Edit: Thinking about this more, you might not have "unmeasured confounding", as you state in your post (if the variability is at the participant and regional level, you have that information). These answers really depend a lot on your data and your research question(s), so please consider the implications before running with any of the models, and if you have a confound that you can't model, remember to make that explicit.
You might be able to get at the rehabilitation effect if you allow each participant's intercept and slope to vary and allow a correlation
complications ~ rehabilitation + (rehabilitation|participant) 

Then, model regional differences:
complications ~ rehabilitation*region + (rehabilitation|participant) 

More info:
You could model different (random) base amount of complications for each participant with
complications ~ rehabilitation + (1|participant)

and model different slopes of rehabilitation for each participant with
complications ~ rehabilitation + (rehabilitation|participant)

Is rehabilitation's regional variability important to your research in some way?  If not, and you just want to account for the fact that participants are in different regions, you could use
complications ~ rehabilitation + (1|region)

This accounts for random intercepts (so participants in each region start with different levels of complications).
To look at whether rehabilitation differs by region you could model an interaction, and that does not require random effects:
complications ~ rehabilitation * region

