How to write a program to calculate bootstrap SEs for a multilevel logistic regression

[Luis Mijares Castaneda]

My question is in the title of this post. I am interested in writing a program to calculate the bootstrap SE, to get a smaller CI, for a multilevel logistic regression model. I wrote some code, but it failed, and I'm not certain if it's correct. Attached below is my model, any help would be appreciated.

local contrls i.long_guns_law c.ppage_c i.biosex i.race i.income i.education ///
                i.employment i.maritalstatus_alt i.children i.political_party 
                

melogit firearm_locked `contrls' [pw=weight] || states:, or

[Andrew Musau]

You have a two-level random-intercept model, so rather than reinventing the wheel, you can estimate it using xtlogit, which has a -vce(bootstrap)- option.

webuse bangladesh, clear
melogit c_use i.urban age i.children || district:, nolog
*EQUIVALENT
xtset district
xtlogit c_use i.urban age i.children, re nolog
*BOOTSTRAP SEs
xtlogit c_use i.urban age i.children, vce(bootstrap, seed(12122025) reps(50)) nolog

Res.:

. melogit c_use i.urban age i.children || district:, nolog

Mixed-effects logistic regression               Number of obs     =      1,934
Group variable: district                        Number of groups  =         60

                                                Obs per group:
                                                              min =          2
                                                              avg =       32.2
                                                              max =        118

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(5)      =     109.60
Log likelihood = -1206.8322                     Prob > chi2       =     0.0000
-------------------------------------------------------------------------------------
              c_use | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
--------------------+----------------------------------------------------------------
              urban |
             Urban  |   .7322765   .1194857     6.13   0.000     .4980888    .9664641
                age |  -.0264981   .0078916    -3.36   0.001    -.0419654   -.0110309
                    |
           children |
           1 child  |   1.116001   .1580921     7.06   0.000     .8061465    1.425856
        2 children  |   1.365895   .1746691     7.82   0.000      1.02355     1.70824
3 or more children  |   1.344031   .1796549     7.48   0.000     .9919139    1.696148
                    |
              _cons |   -1.68929   .1477591   -11.43   0.000    -1.978892   -1.399687
--------------------+----------------------------------------------------------------
district            |
          var(_cons)|    .215618   .0733222                      .1107208    .4198954
-------------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 43.39       Prob >= chibar2 = 0.0000

. 
. *EQUIVALENT

. 
. xtset district

Panel variable: district (unbalanced)

. 
. xtlogit c_use i.urban age i.children, re nolog

Random-effects logistic regression                   Number of obs    =  1,934
Group variable: district                             Number of groups =     60

Random effects u_i ~ Gaussian                        Obs per group:
                                                                  min =      2
                                                                  avg =   32.2
                                                                  max =    118

Integration method: mvaghermite                      Integration pts. =     12

                                                     Wald chi2(5)     = 109.60
Log likelihood = -1206.8322                          Prob > chi2      = 0.0000

-------------------------------------------------------------------------------------
              c_use | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
--------------------+----------------------------------------------------------------
              urban |
             Urban  |   .7322764   .1194857     6.13   0.000     .4980887    .9664641
                age |  -.0264981   .0078916    -3.36   0.001    -.0419654   -.0110309
                    |
           children |
           1 child  |   1.116001   .1580922     7.06   0.000     .8061466    1.425856
        2 children  |   1.365895   .1746691     7.82   0.000      1.02355     1.70824
3 or more children  |   1.344031   .1796549     7.48   0.000     .9919139    1.696148
                    |
              _cons |   -1.68929   .1477592   -11.43   0.000    -1.978892   -1.399687
--------------------+----------------------------------------------------------------
           /lnsig2u |  -1.534244   .3400608                     -2.200751   -.8677373
--------------------+----------------------------------------------------------------
            sigma_u |   .4643475   .0789532                      .3327461    .6479973
                rho |   .0615089   .0196302                      .0325591    .1131878
-------------------------------------------------------------------------------------
LR test of rho=0: chibar2(01) = 43.39                  Prob >= chibar2 = 0.000

. 
. 
. 
. *BOOTSTRAP SEs

. 
. xtlogit c_use i.urban age i.children, vce(bootstrap, seed(12122025) reps(50)) nolog
(running xtlogit on estimation sample)

Bootstrap replications (50): .........10.........20.........30.........40.........50 done

Random-effects logistic regression                   Number of obs    =  1,934
                                                     Replications     =     50
Group variable: district                             Number of groups =     60

Random effects u_i ~ Gaussian                        Obs per group:
                                                                  min =      2
                                                                  avg =   32.2
                                                                  max =    118

Integration method: mvaghermite                      Integration pts. =     12

                                                     Wald chi2(5)     =  61.31
Log likelihood = -1206.8322                          Prob > chi2      = 0.0000

                                      (Replications based on 60 clusters in district)
-------------------------------------------------------------------------------------
                    |   Observed   Bootstrap                         Normal-based
              c_use | coefficient  std. err.      z    P>|z|     [95% conf. interval]
--------------------+----------------------------------------------------------------
              urban |
             Urban  |   .7322764   .1650622     4.44   0.000     .4087604    1.055792
                age |  -.0264981   .0068093    -3.89   0.000    -.0398441   -.0131522
                    |
           children |
           1 child  |   1.116001     .18062     6.18   0.000     .7619928     1.47001
        2 children  |   1.365895   .1973953     6.92   0.000     .9790075    1.752783
3 or more children  |   1.344031   .2219872     6.05   0.000     .9089443    1.779118
                    |
              _cons |   -1.68929   .1927767    -8.76   0.000    -2.067125   -1.311454
--------------------+----------------------------------------------------------------
           /lnsig2u |  -1.534244    .292176                     -2.106899   -.9615898
--------------------+----------------------------------------------------------------
            sigma_u |   .4643475   .0678356                      .3487328    .6182917
                rho |   .0615089    .016866                      .0356486    .1041037
-------------------------------------------------------------------------------------
LR test of rho=0: chibar2(01) = 43.39                  Prob >= chibar2 = 0.000