multilevel, multinomial logistic regression error message with xtmlogit but not with gsem

[Suvi Maatta]

Dear all,

I am interested in conducting a multilevel, multinomial logistic regression model in Stata. The dependent variable is a categorical with 4 distinct groups. The explanatory factors include factors of individual background and characteristics of school and neighborhood settings (as a group average). The data is nested with schools (pupils in schools).

I conducted the analyses with gsem at first, and it worked fine. Then, I found this discussion (https://www.statalist.org/forums/for...random-effects ), and thus, tried to conduct the analysis by using xtmlogit with school random effects. I did get following error message: Cannot compute an improvement -- discontinuous region encountered. I tested it so that I included one predictor at a time without any luck. Now, I wonder what could explain this error message?
Code: xtmlogit expclass i.familyeconomic3 if otos==1, re rrr vce(cluster koulu_koodi)

As I was able to conduct the analyses by using gsem, I wonder if I could still report the results based on gsem?

Thank you for all the help,
Suvi

[Richard Williams]

What was the gsem command? We want to make sure the 2 commands really were equivalent. Also what was your xtset command?

Some generic advice on what to do when you have convergence problems can be found at

https://www3.nd.edu/~rwilliam/xsoc73994/L02.pdf

[Suvi Maatta]

Hi, thank you for the reply.
Gsem command: gsem (i.expclass <- i.familyeconomic3 M[school_code]), mlogit

...and just noticed that I did not translate the measures in English in my previous post, so xtset-command was "xtset koulu_koodi" (in english: "xtset school_code").

Code: xtmlogit expclass i.familyeconomic3 if sample==1, re rrr vce(cluster school_code)

Suvi

[Richard Williams]

Thanks. In the xtmlogit command you are restricting it to cases where sample equals 1, but in the gsem command there is no such restriction.

If I am reading it right, the two commands are not equivalent. Try either limiting the sample the same way in gsem or drop the sample restriction in xtmlogit.

You also have the vce option in xtmlogit but not in gsem, so perhaps that affects things.

[Suvi Maatta]

Hi,
thank you for the reply. I unified the commands (please see below), and I was still able to conduct analyses by using gsem, but not with xtmlogit.
gsem (i.expclass <- i.familyeconomic3 M[school_code]) if sample==1, mlogit vce(cluster school_code)
xtmlogit expclass i.familyeconomic3 if sample==1, re rrr vce(cluster school_code)

I also tested the model so that I did not have vce(cluster school_code) included in xtmlogit, and so that I replaced it with vce(robust), but it did not run through in either of these cases. I also tested the model so that the outcome variable had only three classes (instead of four), but this did not help with xtmlogit either.

Just to clarify one thing, in these analyses, the sample consists of those pupils, who have answered each question that we use in analyses and lives in our city (sample size= 10 680). When I included the pupils, who live outside our city (about 3000 cases), xtmlogit was run through. It is the only time when xtmlogit has worked. The overall data consists of pupils attending schools in our city (i.e. no survey sampling procedure was carried out). As I am interested in the influence of neighborhoods in our city on behavior, I have excluded the pupils living outside the city from the analyses.

I really appreciate, if you have any other suggestions what to do. Thank you.

Suvi

[Andrew Musau]

You will increase your chances of receiving helpful replies by including a data example that reproduces your issue (see FAQ Advice #12 for tips on how to pose questions on the forum). Since you are new, I have included a reproducible example on your behalf. To replicate xtmlogit's output using gsem, you need to specify separate random-effect variances for each non-reference level of the outcome. How you constrain the variances and covariances will determine whether you obtain equivalent results. Below, for example, I request all variances and covariances to be estimated distinctly in xtmlogit using the option -cov(unstructured)-

webuse estatus, clear
xtset id
xtmlogit estatus i.hhchild age hhincome i.hhsigno i.bwinner, cov(unstructured)
tab estatus
gsem (1.estatus <- i.hhchild age hhincome i.hhsigno i.bwinner R1[id]) ///
     (2.estatus <- i.hhchild age hhincome i.hhsigno i.bwinner R2[id]), mlogit

Res.:

. xtset id

Panel variable: id (unbalanced)

. xtmlogit estatus i.hhchild age hhincome i.hhsigno i.bwinner, cov(unstructured) nolog

Random-effects multinomial logistic regression       Number of obs    =  4,761
Group variable: id                                   Number of groups =    800

Random effects u_i ~ Gaussian                        Obs per group:
                                                                  min =      5
                                                                  avg =    6.0
                                                                  max =      7

Integration method: mvaghermite                      Integration pts. =      7

                                                     Wald chi2(10)    = 242.93
Log likelihood = -4438.2887                          Prob > chi2      = 0.0000

------------------------------------------------------------------------------------
           estatus | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------------+----------------------------------------------------------------
Out_of_labor_force |
           hhchild |
              Yes  |   .4924799   .1002988     4.91   0.000      .295898    .6890619
               age |   -.004219   .0070064    -0.60   0.547    -.0179513    .0095133
          hhincome |   -.006046    .001992    -3.04   0.002    -.0099503   -.0021417
                   |
           hhsigno |
              Yes  |   .5036976   .0966982     5.21   0.000     .3141726    .6932225
                   |
           bwinner |
              Yes  |   -.489057   .0745454    -6.56   0.000    -.6351632   -.3429507
             _cons |  -.3930378    .298386    -1.32   0.188    -.9778636     .191788
-------------------+----------------------------------------------------------------
Unemployed         |
           hhchild |
              Yes  |   .0399687   .1238417     0.32   0.747    -.2027565    .2826939
               age |   .0045538   .0085081     0.54   0.592    -.0121219    .0212294
          hhincome |  -.0315377   .0027426   -11.50   0.000    -.0369131   -.0261624
                   |
           hhsigno |
              Yes  |   .1495817   .1214242     1.23   0.218    -.0884053    .3875687
                   |
           bwinner |
              Yes  |  -.2552257   .0968165    -2.64   0.008    -.4449826   -.0654689
             _cons |  -.0417024   .3633406    -0.11   0.909    -.7538368     .670432
-------------------+----------------------------------------------------------------
Employed           |  (base outcome)
-------------------+----------------------------------------------------------------
            var(u1)|   1.132081   .1331468                       .899012    1.425572
            var(u2)|   1.102612   .1698422                      .8152803     1.49121
-------------------+----------------------------------------------------------------
         cov(u1,u2)|   .7871916   .1222148     6.44   0.000      .547655    1.026728
------------------------------------------------------------------------------------
LR test vs. multinomial logit: chi2(3) = 286.41           Prob > chi2 = 0.0000

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

. tab estatus

 Employment status |      Freq.     Percent        Cum.
-------------------+-----------------------------------
Out of labor force |      1,682       35.33       35.33
        Unemployed |        703       14.77       50.09
          Employed |      2,376       49.91      100.00
-------------------+-----------------------------------
             Total |      4,761      100.00

. gsem (1.estatus <- i.hhchild age hhincome i.hhsigno i.bwinner R1[id]) ///
>      (2.estatus <- i.hhchild age hhincome i.hhsigno i.bwinner R2[id]), mlogit nolog

Generalized structural equation model                    Number of obs = 4,761
Response:     estatus    
Base outcome: 3          
Family:       Multinomial
Link:         Logit      
Log likelihood = -4438.2887

 ( 1)  [1.estatus]R1[id] = 1
 ( 2)  [2.estatus]R2[id] = 1
------------------------------------------------------------------------------------
                   | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------------+----------------------------------------------------------------
1.estatus          |
           hhchild |
              Yes  |   .4924799   .1002988     4.91   0.000      .295898    .6890619
               age |   -.004219   .0070064    -0.60   0.547    -.0179513    .0095133
          hhincome |   -.006046    .001992    -3.04   0.002    -.0099503   -.0021417
                   |
           hhsigno |
              Yes  |   .5036976   .0966982     5.21   0.000     .3141726    .6932225
                   |
           bwinner |
              Yes  |   -.489057   .0745454    -6.56   0.000    -.6351632   -.3429507
                   |
            R1[id] |          1  (constrained)
                   |
             _cons |  -.3930378    .298386    -1.32   0.188    -.9778636     .191788
-------------------+----------------------------------------------------------------
2.estatus          |
           hhchild |
              Yes  |   .0399687   .1238417     0.32   0.747    -.2027565    .2826939
               age |   .0045538   .0085081     0.54   0.592    -.0121219    .0212294
          hhincome |  -.0315377   .0027426   -11.50   0.000    -.0369131   -.0261624
                   |
           hhsigno |
              Yes  |   .1495817   .1214242     1.23   0.218    -.0884053    .3875687
                   |
           bwinner |
              Yes  |  -.2552257   .0968165    -2.64   0.008    -.4449826   -.0654689
                   |
            R2[id] |          1  (constrained)
                   |
             _cons |  -.0417024   .3633406    -0.11   0.909    -.7538368     .670432
-------------------+----------------------------------------------------------------
3.estatus          |  (base outcome)
-------------------+----------------------------------------------------------------
        var(R1[id])|   1.132081   .1331468                       .899012    1.425572
        var(R2[id])|   1.102612   .1698422                      .8152803     1.49121
-------------------+----------------------------------------------------------------
 cov(R1[id],R2[id])|   .7871916   .1222148     6.44   0.000      .547655    1.026728
------------------------------------------------------------------------------------