Multilevel model of drug efficacy with low sized clusters does not run. What if I aggregate clusters?

[Gianfranco Di Gennaro]

Dear all,
I am conducting an analysis where my outcome is effectiveness (effective/not effective) of a drug.
My patients come from different hospitals, and each of them has been observed repeatedly over time. My problem is that I have clusters (hospitals) often of very low sample size (1, 2... maximum 5 patients), and when I include the hospital clustering in my model, the model either does not converge or cannot estimate appropriate initial values .

The model is:
xtmelogit efficacy....... || Hospital: || patientid:

Anyway, when I remove the hospital clustering, the model runs smoothly.

xtmelogit efficacy....... || patientid:

In your opinion, does it make more sense to keep the latter model without hospital-level clustering or to group the low-sized clusters into a big cluster named "other" ?

[Erik Ruzek]

First things first, do not create a "big cluster" with all the low-sized clusters put into it. The structure of the data is meaningful. When you start assigning patients to clusters based on some arbitrary distinction, all meaning is lost.

With that out of the way, this requires a bit of investigation. It's not clear that the problem is that you have too small of sample sizes within hospitals or if there is little variation in the outcome between hospitals. Can you post the results of the following? I would like to know if the model will run and the ICC that comes from a model with just hospitals as the clustering variable

*Does the model estimate without predictors?
melogit efficacy || Hospital: || patientid: // drop the xt before melogit to use the most up to date version
estat icc
*If not, what do you get from this?
melogit efficacy || Hospital: /
estat icc

[Gianfranco Di Gennaro]

Thanks for your time.
I have a very high icc by considering the hospital clustering only.

Please be aware: "EventiYN" is my outcome (efficacy); "Prefisso_Centro" is my hospital.

. melogit EventiYN  || Prefisso_Centro:

Fitting fixed-effects model:

Iteration 0:   log likelihood = -136.36114  
Iteration 1:   log likelihood = -136.03671  
Iteration 2:   log likelihood = -136.03608  
Iteration 3:   log likelihood = -136.03608  

Refining starting values:

Grid node 0:   log likelihood =  -116.1713

Fitting full model:

Iteration 0:   log likelihood =  -116.1713  
Iteration 1:   log likelihood = -112.22082  
Iteration 2:   log likelihood = -111.53104  
Iteration 3:   log likelihood = -111.49579  
Iteration 4:   log likelihood = -111.49539  
Iteration 5:   log likelihood = -111.49543  
Iteration 6:   log likelihood = -111.49544  

Mixed-effects logistic regression               Number of obs     =        363
Group variable: Prefisso_Cen~o                  Number of groups  =         19

                                                Obs per group:
                                                              min =          4
                                                              avg =       19.1
                                                              max =         57

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(0)      =          .
Log likelihood = -111.49544                     Prob > chi2       =          .
---------------------------------------------------------------------------------
       EventiYN | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
----------------+----------------------------------------------------------------
          _cons |  -2.723959   .6719897    -4.05   0.000    -4.041034   -1.406883
----------------+----------------------------------------------------------------
Prefisso_Centro |
      var(_cons)|   5.781188   3.755498                      1.618344    20.65206
---------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 49.08       Prob >= chibar2 = 0.0000

. estat icc

Intraclass correlation

------------------------------------------------------------------------------
                       Level |        ICC   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
             Prefisso_Centro |   .6373225   .1501517      .3297217    .8625896
------------------------------------------------------------------------------

Here's when I go with the id (here I have to use xtmelogit because with melogit the model does not converge):

. xtmelogit EventiYN  || Prefisso_Centro: || id:

Refining starting values:

Iteration 0:   log likelihood = -109.45546  
Iteration 1:   log likelihood = -99.436292  
Iteration 2:   log likelihood = -98.981539  

Performing gradient-based optimization:

Iteration 0:   log likelihood = -98.981539  
Iteration 1:   log likelihood = -97.397885  
Iteration 2:   log likelihood = -97.381487  
Iteration 3:   log likelihood = -97.381482  

Mixed-effects logistic regression               Number of obs     =        363

----------------------------------------------------------------------------
                |     No. of       Observations per group       Integration
 Group variable |     groups    Minimum    Average    Maximum      points
----------------+-----------------------------------------------------------
   Prefisso_C~o |         19          4       19.1         57           7
             id |        115          1        3.2          4           7
----------------------------------------------------------------------------

                                                Wald chi2(0)      =          .
Log likelihood = -97.381482                     Prob > chi2       =          .

------------------------------------------------------------------------------
    EventiYN | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |  -5.020491    1.31831    -3.81   0.000    -7.604332   -2.436651
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
Prefisso_C~o: Identity       |
                   sd(_cons) |   2.526879   1.258755      .9518347    6.708221
-----------------------------+------------------------------------------------
id: Identity                 |
                   sd(_cons) |   3.576373   .9672312      2.104917    6.076458
------------------------------------------------------------------------------
LR test vs. logistic model: chi2(2) = 77.31               Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

. estat icc

Residual intraclass correlation

------------------------------------------------------------------------------
                       Level |        ICC   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
             Prefisso_Centro |   .2842197   .2053643      .0520479    .7417124
          id|Prefisso_Centro |   .8535586   .0666755      .6720135    .9431212
------------------------------------------------------------------------------

Thank you very much Erik!

[Erik Ruzek]

You do have a fair bit of variance between hospitals, however the bulk, about 67%, of the variance at the higher levels of the data lies between patients (3.6²/(3.6²+2.53²)). At what level of the hierarchy is your key predictor? If it is at the measurement occasion or at the patient level, then you might consider including indicators for each of the hospitals. This treats hospitals as more of a nuisance than something of substantive interest. As a result, all coefficients give you log/log odds in which the contrast represents mean log/log odds differences when comparing patients to other patients in the same hospital.

If your key predictor is at the hospital level, then you are in more of a difficult position. I would carefully look at how the covariates at the measurement occasion and at the patient level influence the between hospital variance. It could be that one or more of your lower-level predictors have an adverse impact on between hospital variance. To help with this, you could include group level means for each of the lower level predictors. This helps to appropriately partition variance in those predictors and could help w/ estimation, although no guarantee:

*For a predictor(s) at the measurement occasions
egen id_mn_var = mean(patient_variable), by(id)
egen hosp_mn_var = mean(hospital_variable), by(Hospital)

[Gianfranco Di Gennaro]

Thanks again, Erik.
I'm sorry for taking up your time. But... I don't understand what you mean by "you might consider including indicators for each of the hospitals." "Hospital" in my dataset is a string variable. How, in practice, should I create indicators?

All my predictors of interest are at the patient level (Age, gender, drug taken, number of previous drugs, months since the onset of the disease). How would you specify the model?

melogit efficacy var1 var2 ... || ...?

Feel free not to respond if I'm taking advantage of your availability.
Gianfranco

[Clyde Schechter]

To create indicators from a string variable, you first have to make a corresponding numeric variable. If there are fewer than 65,536 hospitals you can do this with -encode-:
[code
encode hospital, gen(n_hospital)
[/code]
Then in your model you include it as

melogit efficacy var1 var2 ... i.n_hospital || patient_id:

In the unlikely event you have more than 65,536 hospitals, use -egen n_hospital = group(hospital)- instead of -encode-. (But if you have that many hospitals, I think that many indicators will cause -melogit- to blow through matrix size limits, so the approach isn't really viable in this situation.)

[Gianfranco Di Gennaro]

Thank you Clyde Schechter.
Unfortunately I don't have enough data to properly estimate a coefficient for each hospital.

PS: I succeded in running the two-levels model by using a different integration method (gh, nonadaptive Gauss-Hermite quadrature) but I obtained very large (and hardly trustable) Odds Ratios.

As far as you know, is there a way to determine if the results obtained from these quadrature methods are reliable?
Thank you for your time, as always

[Clyde Schechter]

The details of the different integration methods is beyond my scope of knowledge and I can't advise you about that. But let me just challenge you on your concern that the odds ratios are too large to trust. What is the scale of the variables whose odds ratios seem too large? If the variables are dichotomous 0/1 variables, then, yes you can form a judgment that an odds ratio beyond a threshold (I usually use 4) is suspiciously large and warrants further investigation. But if the variable is continuous, remember that the coefficient scales inversely with the scale of the variable itself. So if the variable has a very small range of values but is highly influential on the outcome, then the coefficient will necessarily be very large, and that will in turn translate into a large odds ratio. Also, if the outcome probabilities are all very close to 0 or all very close to 1, bear in mind that even a huge odds ratio will correspond go a minuscule change in the outcome probability. So you really need to judge these things based on the associated impact on outcome probabilities; you cannot judge the odds ratios in isolation when the regressor is continuous.

Remember that the odds ratio you are shown is the odds ratio associated with a unit change in the variable--and if the variable's scale is such that unit changes are too large to occur, or can only occur very rarely, then a very large odds ratio is perfectly fine. I think a reasonable diagnostic is to get the 10th and 90th percentiles of the variable itself and then calculate the corresponding expected difference in the outcome probabilities (use -margins-) and see if that seems reasonable or not.

[Gianfranco Di Gennaro]

Dear Clyde Schechter , thank you.
Unfortunately in my variable of interest is categorical (my drug compared to a reference, "AED_aggiu" in my output) and an OR of 35K really tells me that there is something wrong.
I don't have empty cells or other similar problems and I'm having a lot of trouble understanding where the problem comes from.

Mixed-effects logistic regression               Number of obs     =        359

        Grouping information
        -------------------------------------------------------------
                        |     No. of       Observations per group
         Group variable |     groups    Minimum    Average    Maximum
        ----------------+--------------------------------------------
           Prefisso_C~o |         18          4       19.9         57
                     id |        114          1        3.1          4
        -------------------------------------------------------------

Integration method: ghermite                    Integration pts.  =          7

                                                Wald chi2(6)      =      16.97
Log likelihood = -86.625242                     Prob > chi2       =     0.0094
--------------------------------------------------------------------------------------------
                  EventiYN | Odds ratio   Std. err.      z    P>|z|     [95% conf. interval]
---------------------------+----------------------------------------------------------------
Mesi_da_diagnosi_epilessia |   .9956514   .0023775    -1.83   0.068     .9910025    1.000322
                       Età |   1.407398   .1718527     2.80   0.005     1.107846    1.787946
                 _ISesso_1 |   538.1197   1061.272     3.19   0.001     11.27589     25680.7
             _IAED_aggiu_2 |   79.53746   166.1264     2.10   0.036     1.326493    4769.125
             _IAED_aggiu_3 |   875.1254   1860.568     3.19   0.001     13.56333    56464.35
             _IAED_aggiu_4 |   35270.55   108073.7     3.42   0.001      86.9391    1.43e+07
                     _cons |   3.06e-20   4.23e-19    -3.26   0.001     5.60e-32    1.68e-08
---------------------------+----------------------------------------------------------------
Prefisso_Centro            |
                 var(_cons)|   47.84136   23.85362                      18.00511    127.1192
---------------------------+----------------------------------------------------------------
Prefisso_Centro>id         |
                 var(_cons)|   51.05243    26.5831                      18.39899    141.6572
--------------------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation to odds ratios.
Note: _cons estimates baseline odds (conditional on zero random effects).
LR test vs. logistic model: chi2(2) = 68.95               Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

As regards the size of the ORs in the case of continuous variables, I absolutely agree with you and I also use margins to evaluate the probabilities at different percentiles (or at values of the variable of my interest).

[Clyde Schechter]

I'm going to infer that your variables _IAED_aggiu* are indicators for levels 2 through 4 of some variable _IAED_aggiu, and that level 1 is the omitted reference value. I look at those odds ratios and they are indeed eye-popping. As is the one for -ISesso_1, which I guess is the 1 level of a 0/1 dichotomous variable. I'm also guessing that Mesi_da_diagnosi_epilessia and Età are some continuous variables (though what I have to say next doesn't depend crucially on this--their odds ratios being close to 1, my reasoning would go through regardless).

Now I look at another eye-popping result: a constant term of 3.06e-20. The implication of all of these, unless the estimates are somehow grossly wrong, is that for an observation in which ISesso == 0 and _IAED_aggiu = 1, the probability of EventiYN being 1 is excruciatingly close to zero. Thus it seems like you are modeling an event that is extremely rare, at least in these base-value conditions. Maximum likelihood estimation is not good at estimating very rare events (nor the complementary situation where the event almost always occurs). This may also explain why initially you encountered convergence problems. If this were not a multi-level model, I would tell you to use -firthlogit-, which uses penalized maximum likelihood, instead. But I am not aware of any implementation of penalized maximum likelihood estimation for multi-level logistic models. (And the variance components at the higher levels are clearly not negligible, so we can't just go back to a single level.)

You can confirm that this is your problem with:

tab EventiYN if ISesso == 0 & I_AED_aggiu == 1 & !missing(Mesi_da_diagnosi_epilessia, Età, id,Prefisso_Centro)

I expect you will find that EventiYN is nearly always zero in this table.

If I am right about this, the first thing you should do is check to see if this represents an error in your data set. Perhaps somehow during data management, a whole bunch of EventiYN = 1 observations with these values of ISesso and I_AED_aggiu got wiped out, or overwritten.

Assuming, however, that the this rare outcome behavior is not due to data errors, I think if I were in your shoes, I would remove the observations where ISesso = 0 and _IAED_aggiu = 1 and try again. You might have to go farther and remove the observations where ISesso = 0 or _IAED_aggiu = 1 (which would entail removing ISesso from the model altogether, and would restructure _IAED_aggiu so that 2 is the base value.) And you would have to qualify any presentation of your results by saying that they are not generalizable to these excluded situations. But I think you will probably get more reasonable results if you do this. Also, while it will not affect the results you get, I recommend you redo these variables with factor-variable notation: it's a good habit to get into, and it's actually easier than creating your own indicator variables. Crucially, it makes it possible to use the -margins- command after estimation, which you may or may not need to do here, but which is just so often really helpful. (See -help fvvarlist- if you are not familiar with factor-variable notation.)

[Joseph Coveney]

I don't have empty cells or other similar problems and I'm having a lot of trouble understanding where the problem comes from.

According to the intercept-only hierarchical logistic regression models that you successfully fitted with melogit (two-level) and xtmelogit (three-level) above in #3, you have a total of 363 observations distributed among a total of 115 patients who in turn are seen at or admitted to a total of 19 hospitals.

And according to their intercept regression coefficients, among these 363 observations you have somewhere between about 2 and about 22 total positive outcomes marginally for your EventiYN efficacy variable.

. display invlogit(-5.020491) * 363
2.3805516

. display invlogit(-2.723959) * 363
22.351473

Are those in the ballpark as to what you've got?

If the Statistical Considerations section of your clinical study protocol will sanction it, then you might want to consider a far less ambitious regression model—or even an alternative statistical method that's not so subject to the risk of nonconvergence—to address your research question.

[Gianfranco Di Gennaro]

Dear Erik Ruzek, Clyde Schechter and Joseph Coveney, thanks for your support.
I hope you had a happy Christmas time and I wish you a great 2024.

I'm running more parsimonious models, as Joseph suggested.
I will report all the limits of the definitive model in the discussion of my study.

Thanks again for your time and patience.