Hello,

My question is on estimating coefficients for random effects in mixed effect models. I am building a 3 level mixed effect model in order to describe hospital-level rates of testing patients with neuro-imaging. My model includes ~ 45,00 patients clustered within 480 hospitals clustered within 73 geographic regions (health referral regions or HRRs). My goal is to estimate adjusted diagnostic testing rates for each hospital (the second level in the model).

My goal is a graph similar to figure 2 of the following manuscript, which is a dot plot depicting estimated hospital level rates of a binary outcome with 95% CIs:
HTML Code:
 https://www.ncbi.nlm.nih.gov/pubmed/29404575
Below is the mixed effect model I have built to adjust for patient level factors such as age, gender, co-morbidities and to account for clustering within hospitals and regions. The outcome is binary, whether or not a patient received a neuro-imaging test, and thus I used the melogit command:

Code:
 
melogit neuro_imaging i.u_icu i.u_ccu i.obs age i.gender i.pay1 i.aweekend i.anticoag i.neuro_symptom i.head_trauma cm_aids cm_alcohol cm_anemdef cm_arth cm_bldloss cm_chf cm_chrnlung cm_coag cm_depress cm_dm cm_dmcx cm_drug cm_htn_c cm_hypothy cm_liver cm_lymph cm_lytes cm_mets cm_neuro cm_obese cm_para cm_perivasc cm_psych cm_pulmcirc cm_renlfail cm_tumor cm_ulcer cm_valve cm_wghtloss || hrrcode: || hospital:, or intpoints(10)

Fitting fixed-effects model:

Iteration 0:   log likelihood = -27240.211  
Iteration 1:   log likelihood = -27197.166  
Iteration 2:   log likelihood = -27197.154  
Iteration 3:   log likelihood = -27197.154  

Refining starting values:

Grid node 0:   log likelihood = -22450.666

Fitting full model:

Iteration 0:   log likelihood = -22450.666  
Iteration 1:   log likelihood = -22277.979  
Iteration 2:   log likelihood = -22271.953  
Iteration 3:   log likelihood = -22271.818  
Iteration 4:   log likelihood = -22271.818  

Mixed-effects logistic regression               Number of obs     =     45,513

-------------------------------------------------------------
                |     No. of       Observations per Group
 Group Variable |     Groups    Minimum    Average    Maximum
----------------+--------------------------------------------
        hrrcode |         73         11      623.5      7,301
       hospital |        486         10       93.6      1,077
-------------------------------------------------------------

Integration method: mvaghermite                 Integration pts.  =         10

                                                Wald chi2(43)     =    2017.57
Log likelihood = -22271.818                     Prob > chi2       =     0.0000
-------------------------------------------------------------------------------------
      neuro_testing | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
                    |
            1.u_icu |   .9617766   .0441584    -0.85   0.396     .8790077    1.052339
            1.u_ccu |   .9170839   .0490852    -1.62   0.106     .8257529    1.018516
              1.obs |   1.097361   .0564116     1.81   0.071     .9921839    1.213687
                age |   1.006762   .0010272     6.60   0.000      1.00475    1.008777
                    |
             gender |
            Female  |   1.086055   .0273125     3.28   0.001     1.033822    1.140928
                    |
               pay1 |
          Medicaid  |   1.020001   .0497538     0.41   0.685     .9270017     1.12233
 Private Insurance  |    1.07706   .0400239     2.00   0.046     1.001403    1.158433
          Self-Pay  |   1.092733   .0777682     1.25   0.213     .9504636    1.256299
         No-Charge  |   1.134688   .1900918     0.75   0.451     .8171027    1.575711
             Other  |   .9632547   .0924938    -0.39   0.697     .7980074    1.162721
                    |
         1.aweekend |   1.093833   .0296999     3.30   0.001     1.037144    1.153621
         1.anticoag |   .8961464   .0279134    -3.52   0.000     .8430737    .9525601
    1.neuro_symptom |   2.231223   .0697267    25.68   0.000     2.098662    2.372157
      1.head_trauma |   3.775735   .1881069    26.67   0.000      3.42448    4.163018
            cm_aids |   1.497218   .3726232     1.62   0.105     .9192644    2.438538
         cm_alcohol |   1.257131   .0785104     3.66   0.000     1.112298    1.420822
         cm_anemdef |   .9110613   .0304405    -2.79   0.005     .8533106    .9727205
            cm_arth |   1.097065   .0797171     1.27   0.202     .9514388    1.264981
         cm_bldloss |   .7145878   .1033806    -2.32   0.020     .5381593    .9488561
             cm_chf |   .8469047   .0545417    -2.58   0.010     .7464764    .9608443
        cm_chrnlung |   1.062038   .0346219     1.85   0.065     .9963029    1.132111
            cm_coag |   1.048007    .069548     0.71   0.480     .9201882    1.193581
         cm_depress |   1.211311   .0469752     4.94   0.000     1.122653     1.30697
              cm_dm |   1.074146   .0317421     2.42   0.016       1.0137    1.138196
            cm_dmcx |   1.073205   .0643979     1.18   0.239     .9541265    1.207144
            cm_drug |   1.138923   .0853021     1.74   0.082     .9834267    1.319007
           cm_htn_c |   .9917578   .0275891    -0.30   0.766     .9391318    1.047333
         cm_hypothy |   1.009805   .0353911     0.28   0.781     .9427688    1.081609
           cm_liver |    1.03391   .0926915     0.37   0.710     .8673042    1.232521
           cm_lymph |   1.113212   .1431993     0.83   0.404     .8651322    1.432429
           cm_lytes |    .952935   .0257518    -1.78   0.074     .9037757    1.004768
            cm_mets |   1.251564   .1293983     2.17   0.030     1.021993    1.532704
           cm_neuro |    1.42887   .0578913     8.81   0.000     1.319794    1.546962
           cm_obese |   .9985896   .0467823    -0.03   0.976     .9109816    1.094623
            cm_para |    1.33541   .1415068     2.73   0.006     1.084968    1.643661
        cm_perivasc |   .9368681   .0484128    -1.26   0.207     .8466278    1.036727
           cm_psych |   1.402534   .0862985     5.50   0.000     1.243193    1.582297
        cm_pulmcirc |   .9512912    .105854    -0.45   0.654     .7648863    1.183124
        cm_renlfail |   .9259629   .0340846    -2.09   0.037     .8615112    .9952364
           cm_tumor |   .8685913   .0666102    -1.84   0.066     .7473755    1.009467
           cm_ulcer |   1.655414   1.290669     0.65   0.518     .3591321    7.630609
           cm_valve |   .9076863   .0683901    -1.29   0.199     .7830717    1.052132
        cm_wghtloss |   1.148169   .0918705     1.73   0.084      .981515    1.343119
              _cons |   .5648739   .0973378    -3.31   0.001     .4029715     .791824
--------------------+----------------------------------------------------------------
hrrcode             |
          var(_cons)|   1.091473    .302129                      .6344438    1.877728
--------------------+----------------------------------------------------------------
hrrcode>hospital    |
          var(_cons)|   1.713218   .1434292                      1.453954    2.018713
-------------------------------------------------------------------------------------
LR test vs. logistic model: chi2(2) = 9850.67             Prob > chi2 = 0.0000


Next, I use post-estimation commands to obtain empirical bayes means
Code:
 predict hrr_effects hosp_effects, reffects
Finally, I generate an adjusted testing rate accounting for hospital and region random effects
Code:
 gen adj_neurotesting_rate = exp(_cons+hosp_effects+hrr_effects)/(1+exp(_cons+hosp_effects+hrr_effects))
Questions
1. As graphed below, when I compare the adjusted and unadjusted hospital testing rates (unadjusted rates calculated by dividing the # of patients who received a test by the total # of patients seen at the hospital), I notice that all adjusted rates are higher than the crude rates. It seems incorrect that adjustment would only result in higher predicted probabilities for every hospital, but I am unsure where I may have erred in my process, and was hoping the forum might be able to advise me on where I went wrong.
[ATTACH]temp_12910_1545846644215_902[/ATTACH]
2. Do I need to group mean center my variables to obtain accurate estimates of hospital-level rates?
3. I can use the predict, reses command to calculate standard errors of the empirical bayes estimates. How do I transform these estimates into 95% confidence intervals to match the predicted probabilities?
4. If I were to add hospital characteristics in the model as random coefficients, how would I account for this in calculating predicted probabilities of the hospital random effects.
e.g.. a dichotomous marker for profit status as shown in example code below:

melogit neuro_imaging [omitting variables for brevity] || hrrcode: || hospital profit_status:,

Thank you.
Tim