Interpreting the output of a multilevel (xtmelogit) model

[Vian Rajabzadeh]

Hi there

I just wanted to understand more about interpreting the random effects section of a multilevel model below:

------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
usualgpsor~e: Identity |
sd(_cons) | .1897804 .0254148 .1459692 .2467411
------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 36.02 Prob >= chibar2 = 0.0000

I am looking at statin prescribing (the probability of), including several a priori predictors and have clustered by practice level (usualgpsor~e). How can I tell that including practice level in the model has helped to explain some of the variation.

And in this scenario how do I present this part of the output in a table. Happy to give more context to this if required.

Look forward to hearing from any STATA gurus

Thanks
Vian

[Vian Rajabzadeh]

I have used the postestimation estat icc command to calculate how strongly statin prescribing in the same practice is. Is there any way of looking at this as a pre-estimation.

I calculated that patient-level characteristics explained 76% of the variation, by generating a risk score from a logisitic regression model, to then include the risk score into a multilevel (by practice) model. Does this make it worthwhile to conduct a multilevel model given that it could explain the remaining variation?

[Weiwen Ng]

Vian, it's best to present both code and output using the code delimiters - use the # button on the formatting toolbar (between the " and <> buttons).

That said, the output in post #1 says that the standard deviation of the random intercepts is 0.18878 in logit units, with a 95% CI of 0.146 to 0.247. That excludes 0. The test statistic below is a likelihood ratio test for the model you fit minus random effects (i.e. as if you did a pure logit model without a fixed or random effect for practice) versus the random effects model you fit. The null hypothesis is that the models explain the data equally well. The null hypothesis is rejected, therefore there's evidence your model (with a random intercept) fit the data better.

When I've presented this stuff in tables, I've usually just added a line saying "standard deviation of random intercepts".

If there is clustering, a lot of people in health services research lean toward using a random effects model even if the intra-class correlation is low (e.g. 0.05 or less). You are on solid grounds here, even given that 76% of variation in outcomes is explained by a pure logit model - that said, I need to check how you got that statistic. I am not sure that logistic pseudo-R^2s can be interpreted as proportion of variation explained. The AUC definitely cannot be interpreted as such.

[Vian Rajabzadeh]

xtmelogit, or

Mixed-effects logistic regression               Number of obs     =     21,426
Group variable: usualgpsorga~e                  Number of groups  =        167

                                                Obs per group:
                                                              min =          6
                                                              avg =      128.3
                                                              max =        545

Integration points =   7                        Wald chi2(14)     =    3012.87
Log likelihood = -10354.662                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
     statins | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         ccg |
          1  |   1.946596    .140914     9.20   0.000     1.689108    2.243336
          2  |   1.164944   .0752334     2.36   0.018      1.02644    1.322138
          3  |    1.34448   .0974707     4.08   0.000     1.166392    1.549757
             |
   ethnicity |
          2  |   1.333998   .0651296     5.90   0.000     1.212264    1.467957
          3  |   .6851677   .0372073    -6.96   0.000     .6159894     .762115
          4  |   .8237944   .0737609    -2.16   0.030     .6912004    .9818242
          5  |   .6508426   .0534885    -5.23   0.000     .5540148    .7645935
             |
    age_cats |
          1  |   2.465576   .1379175    16.13   0.000     2.209553    2.751264
          2  |   3.423772   .1810174    23.28   0.000     3.086748    3.797594
          3  |   2.287052   .1138616    16.62   0.000      2.07443    2.521468
             |
         sex |   1.327581   .0472624     7.96   0.000     1.238106    1.423522
          ht |   2.266378   .0880347    21.06   0.000     2.100238    2.445661
        diab |   4.667528   .1990217    36.13   0.000     4.293307    5.074366
    ckdcoded |   1.317028    .055296     6.56   0.000     1.212989     1.42999
       _cons |   .2469519   .0174913   -19.75   0.000     .2149428    .2837278
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
usualgpsor~e: Identity       |
                   sd(_cons) |   .1897804   .0254148      .1459692    .2467411
------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 36.02       Prob >= chibar2 = 0.0000

[Vian Rajabzadeh]

Thank you for your response Weiwen - this has helped me a lot.

I calculated the 76%, by running a logistic regression including all of the above predictors and then used the 'predict' function to calculate the log(odds of being prescribed statins) - essentially a risk score. I then included this as a predictor in my multilevel model as the only predictor(model 1). I then also run the multilevel model without any predictors (model 2). I subtract the practice level variance of model 1 from model 2, then divide by model 2 variance (mod2var-mod1var/mod2var).