Dear Statalisters,

I am using -melogit- to compare hospital admission rates between 177 public hospitals while adjusting for patient and hospital factors. The outcome is admission (yes/no), so I'm using -melogit-. I have a set of patient-level and hospital-level variables (eg hospital in a rural or urban area & its capacity/funding - referred to in NSW as 'peer-group') in the fixed effect part of the model, and a random intercept for hospital.

This is the model:
Code:
melogit admitted i.sex age_c i.language_english i.ATS_urgent i.arrival_ambulance ib(last).afterhours i.LHD_metropolitan i.peer_group_simplified || hospital_encoded:,or
This is the output:
Code:
Mixed-effects logistic regression               Number of obs     =    176,610
Group variable:  hospital_enc~d                 Number of groups  =        177

                                                Obs per group:
                                                              min =          3
                                                              avg =      997.8
                                                              max =      5,346

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(13)     =   26383.41
Log likelihood = -75965.437                     Prob > chi2       =     0.0000
---------------------------------------------------------------------------------------
             admitted | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------------+----------------------------------------------------------------
                  sex |
              Female  |   1.162512    .015033    11.64   0.000     1.133418    1.192353
                age_c |   1.031678   .0003643    88.31   0.000     1.030964    1.032392
                      |
     language_english |
             English  |     1.0744   .0251603     3.06   0.002     1.026201    1.124862
                      |
           ATS_urgent |
              Urgent  |   2.263905   .0304137    60.82   0.000     2.205073    2.324307
                      |
    arrival_ambulance |
           Ambulance  |   4.393615     .05869   110.81   0.000     4.280077    4.510164
                      |
           afterhours |
       Working hours  |   1.141963   .0147982    10.24   0.000     1.113324    1.171338
                      |
     LHD_metropolitan |
   Metropolitan LHDs  |   1.210663   .2016394     1.15   0.251     .8734806    1.678005
                      |
peer_group_simplified |
               Major  |   .6380802   .1488773    -1.93   0.054     .4038986    1.008041
            District  |   .2936127   .0686239    -5.24   0.000     .1857076    .4642157
           Community  |   .4790235   .1225161    -2.88   0.004     .2901694    .7907913
       Multi-purpose  |   .4260602   .1066784    -3.41   0.001      .260822    .6959814
           Sub-acute  |   .3520253   .1833226    -2.00   0.045     .1268521    .9769001
           Ungrouped  |   .2191504   .0847144    -3.93   0.000     .1027314    .4674997
                      |
                _cons |   .1309571   .0305555    -8.71   0.000     .0828939    .2068881
----------------------+----------------------------------------------------------------
hospital_encoded      |
            var(_cons)|   .3905498   .0502347                       .303522    .5025307
---------------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.
Note: _cons estimates baseline odds (conditional on zero random effects).
LR test vs. logistic model: chibar2(01) = 6088.61     Prob >= chibar2 = 0.0000

However I am confused about what the next step should be.

Option 1. I understand I run -predict p, mu- to obtain predicted probabilities that will account for both the fixed and the random component of my model. This would allow me to calculate a risk-adjusted admission rate for each hospital following methods described here:
HTML Code:
https://journals.sagepub.com/doi/abs/10.1177/1536867X19854021?journalCode=stja
However there is some literature that suggests that the best way to adjust for case-mix variables is to use shrinkage methods (eg empirical bayes means) eg
HTML Code:
https://pubmed.ncbi.nlm.nih.gov/25742812/
Option 2. I run -predict eb, reffects reses(eb_se) to obtain empirical bayes means and the SE for the random-effects model. However I am having difficulty in interpreting what -eb- and -se- really mean. I'm assuming I cannot interpret them on their own (eg simply plot eb values on a caterpillar plot and claim that 'there is hospital variation') because I do not understand what the metric means. So what would the next recommended step be? Calculate an risk-adjusted admission rate using the eb value?

Thanks in advance,
Giovanni