random effect estimates for all multilevel level regression model

[Melih Bozkurt]

Hello,

My question is on estimating coefficients for random effects in mixed effect models. I am new to mixed effect models. One reviewer insisted on using mixed effect models with AIC and BIC results for my paper. I managed to do that so far. I am building a 3 level mixed effect model. The model includes 2367 students clustered within 93 classes which are clustered within 24 schools. I made 6 different models and the most complex one as follows;

Fitting fixed-effects model:

Iteration 0:   log likelihood = -1363.5138  
Iteration 1:   log likelihood =  -1085.804  
Iteration 2:   log likelihood = -1059.1874  
Iteration 3:   log likelihood = -1058.9712  
Iteration 4:   log likelihood = -1058.9711  

Refining starting values:

Grid node 0:   log likelihood = -1056.6141

Fitting full model:

Iteration 0:   log likelihood = -1056.6141  (not concave)
Iteration 1:   log likelihood =   -1046.24  (not concave)
Iteration 2:   log likelihood = -1040.4289  
Iteration 3:   log likelihood = -1038.3414  
Iteration 4:   log likelihood = -1038.1282  
Iteration 5:   log likelihood = -1038.1257  
Iteration 6:   log likelihood = -1038.1257  

Mixed-effects ologit regression                 Number of obs     =      2,637

-------------------------------------------------------------
                |     No. of       Observations per Group
 Group Variable |     Groups    Minimum    Average    Maximum
----------------+--------------------------------------------
   school_num~r |         24         24      109.9        160
     class_code |         93         10       28.4         41
-------------------------------------------------------------

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(15)     =     332.82
Log likelihood = -1038.1257                     Prob > chi2       =     0.0000
----------------------------------------------------------------------------------------------
        __14_Going_to_Parks1 | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------------------+----------------------------------------------------------------
              Age_Smaller_10 |   3.386472   1.080845     3.82   0.000     1.811646     6.33026
           Age_between_10_14 |   3.232854   .8032719     4.72   0.000     1.986492    5.261207
           Age_between_15_18 |          1  (omitted)
                GenderFemale |   1.157021   .1446025     1.17   0.243     .9056479    1.478165
  Involving Physical activity|   1.006707   .1186667     0.06   0.955      .799037    1.268351
            Overweight_Obeze |   .7651281   .1032076    -1.98   0.047      .587376    .9966716
       Mother_Father_Uni_Gra |   .8212992   .1248524    -1.30   0.195      .609682    1.106368
  Living_in_Gated_community1 |   1.769015   .3537469     2.85   0.004     1.195413     2.61785
             No parks around |    .192004   .0236017   -13.42   0.000     .1508961    .2443106
             Overall quality |   .2330365   .0336746   -10.08   0.000     .1755587    .3093326
 
-----------------------------+----------------------------------------------------------------
                       /cut1 |  -1.568161   .2280446    -6.88   0.000     -2.01512   -1.121202
-----------------------------+----------------------------------------------------------------
school_number                |
                   var(_cons)|   .1506769   .1114597                      .0353501     .642248
-----------------------------+----------------------------------------------------------------
school_number>class_code     |
                   var(_cons)|   .2778238   .1027409                      .1345829    .5735204
----------------------------------------------------------------------------------------------
LR test vs. ologit model: chi2(2) = 41.69                 Prob > chi2 = 0.0000

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

. estat ic

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |      2,637         .  -1038.126      18    2112.251   2218.045
-----------------------------------------------------------------------------
               Note: N=Obs used in calculating BIC; see [R] BIC note.

I have a minor correction for the paper and she wants to see random effect estimates for all multilevel level regression models. However, as I am new to multi-effect models. I could not figure out how to get random effect estimates. I would appreciate help from experienced stata users.

Best Regards,

[Clyde Schechter]

I would ask her to clarify what exactly she wants. Typically when we speak of the random effect estimates from a model we mean:

-----------------------------+----------------------------------------------------------------
school_number |
var(_cons)| .1506769 .1114597 .0353501 .642248
-----------------------------+----------------------------------------------------------------
school_number>class_code |
var(_cons)| .2778238 .1027409 .1345829 .5735204
----------------------------------------------------------------------------------------------

Speaking carefully, one would call those the estimates of the variance components at each level. But she may mean something else, like some display of the actual values of the random intercepts themselves for some particular observations. Anyway, you need to find out for sure what she means, if it isn't what I just showed above.

[Melih Bozkurt]

Hi Clyde,

Thanks for reply. My understanding is exactly like you have described. School number and school_number>class_code shows variance components. All literature I refer to shows me variance components as random effect estimates. I have reported them at the table. I shared the table showing all model summaries below.

Reviewers' whole comments as follows "I recommend reporting random effect estimates for all multilevel level regression models and maintaining consistency in reporting decimal place values for all estimates. It is not clear if the mixed-effects models use random intercept or random slope—authors need to include the information and relevant estimates. These are minor revisions that can be easily addressed using the model results."

Should I understand something different and look at some other results? I am a bit confused. If I misunderstood the reviewer's comment then please can you tell me what to do and how to do it?

Regards,

[Clyde Schechter]

I think perhaps the reviewer did not understand your table row stubs School and School > Class to refer to the variance components of the random effects. I suggest that you simply rename those rows to "Variance of Random Intercepts at School Level" and "Variance of Random Intercepts at Class Level" so it will be unmistakable. As for the decimal points, I would suggest showing only one, or at most two, decimal places. Looking at the standard errors in your Stata output, that is all the precision you can really claim to have. In addition, for the age variable, I recommend showing the Age 15 to 18 group as 1 [ref.], so that it is clear that this is the omitted reference category. And when you describe your analytic methods, explicitly state that you incorporated random intercepts, but no random slopes, at the class and school levels, with school nested in class.

[Melih Bozkurt]

Dear Clyde,

Thank you very much for your clarification. That sorted a lot of things for me. I will definitely state this is random intercepts, but no random slopes. However, I am still struggling to understand whether I could include "relevant estimates" of random intercept. The referee asked for this too as you can see in the comments. Is that a different estimate than what I already showed in the table? If so, how can I calculate those estimates?

[Clyde Schechter]

I don't think the referee is asking for anything more.

It is possible to calculate the actual estimated values of the random intercepts themselves for each class and school.

predict random_intercept*, reffects

will create two new variables, random_intercept1 and random_intercept2 that contain the empirical bayes estimates of the school and class level intercepts, respectively. But there will be 24 such estimates for schools and 93 such for classes. I suppose you could make a table of these results, but it would be rather long and it is hard to imagine anybody would really want to see it. Also, if you were to name the schools and classes in that table by anything other than an encoded study ID, you might encounter angry pushback from the schools or class teachers themselves, or possibly even litigation.

I don't think the reviewer really wants that. If you are unsure, I suggest you contact the journal editor and ask him or her to clarify what the reviewer is asking for.

[Melih Bozkurt]

Dear Clyde,

Thank you very much for your detailed explanation. That really clarified everything. This is really the best help I have ever received about statistics because every time I dig down to understand those random effects, random intercepts I have come across something really confusing. So at some point, I thought, ok I cannot do this until I got a reply from you. It is not very far from what I understood in the first place but many tutorials on the net made me so confused. Thanks again.

[Allan Angeli]

Dear all,

I am working with panel data.

I did the tests and they suggest fixed effect adjustment. However, analyzing the model and context, I cannot safely discard the variable effect.

1) With these tests, is it possible to safely say that I only have a fixed effect?

2) How is the code to see the mixed effect (fixed and random effect at the same time)? So how do I do xtreg fe and re, what is the "xtreg" for the mixed effect? xtmixed?

I appreciate any help. Thanks a lot!

Banco = Bank = id
Data = date/time (semesters)
$dependente = dependent
$explicativas = independent variables

xtreg $dependente $explicativas, fe

Fixed-effects (within) regression               Number of obs      =      3534
Group variable: banco                           Number of groups   =       186

R-sq:  within  = 0.5596                         Obs per group: min =         1
       between = 0.5303                                        avg =      19.0
       overall = 0.5155                                        max =        38

                                                F(22,3326)         =    192.07
corr(u_i, Xb)  = -0.1482                        Prob > F           =    0.0000

         rho |  .67080754   (fraction of variance due to u_i)
------------------------------------------------------------------------------
F test that all u_i=0:     F(185, 3326) =    15.88           Prob > F = 0.0000

xtreg $dependente $explicativas, re theta

Random-effects GLS regression                   Number of obs      =      3534
Group variable: banco                           Number of groups   =       186

R-sq:  within  = 0.5572                         Obs per group: min =         1
       between = 0.7092                                        avg =      19.0
       overall = 0.6631                                        max =        38

                                                Wald chi2(22)      =   4680.29
corr(u_i, X)   = 0 (assumed)                    Prob > chi2        =    0.0000

------------------- theta --------------------
  min      5%       median        95%      max
0.2381   0.3606     0.7094     0.8126   0.8126

 rho |  .41958416   (fraction of variance due to u_i)

quietly xtreg $dependente $explicativas, re

xttest0

Breusch and Pagan Lagrangian multiplier test for random effects

        inter[banco,t] = Xb + u[banco] + e[banco,t]

        Estimated results:
                         |       Var     sd = sqrt(Var)
                ---------+-----------------------------
                   inter |   480.5217        21.9208
                       e |   84.10534       9.170897
                       u |   60.79998       7.797434

        Test:   Var(u) = 0
                             chibar2(01) =  5715.62
                          Prob > chibar2 =   0.0000

 quietly xtreg $dependente $explicativas, fe
estimates store fixed
quietly xtreg $dependente $explicativas, re
estimates store random

Note: the rank of the differenced variance matrix (19) does not equal the number of coefficients being tested (22); be sure this is what you
        expect, or there may be problems computing the test.  Examine the output of your estimators for anything unexpected and possibly
        consider scaling your variables so that the coefficients are on a similar scale.

                           b = consistent under Ho and Ha; obtained from xtreg
            B = inconsistent under Ha, efficient under Ho; obtained from xtreg

    Test:  Ho:  difference in coefficients not systematic

                 chi2(19) = (b-B)'[(V_b-V_B)^(-1)](b-B)
                          =  -992.36    chi2<0 ==> model fitted on these
                                        data fails to meet the asymptotic
                                        assumptions of the Hausman test;
                                        see suest for a generalized test

hausman fixed random, sigmamore

 b = consistent under Ho and Ha; obtained from xtreg
            B = inconsistent under Ha, efficient under Ho; obtained from xtreg

    Test:  Ho:  difference in coefficients not systematic

                 chi2(19) = (b-B)'[(V_b-V_B)^(-1)](b-B)
                          =      123.82
                Prob>chi2 =      0.0000

hausman fixed random, sigmaless

 b = consistent under Ho and Ha; obtained from xtreg
            B = inconsistent under Ha, efficient under Ho; obtained from xtreg

    Test:  Ho:  difference in coefficients not systematic

                 chi2(19) = (b-B)'[(V_b-V_B)^(-1)](b-B)
                          =      127.54
                Prob>chi2 =      0.0000

xi: xtreg $dependente $explicativas, re

xtoverid

Test of overidentifying restrictions: fixed vs random effects
Cross-section time-series model: xtreg re   
Sargan-Hansen statistic 127.820  Chi-sq(22)   P-value = 0.0000

Testing for time-fixed effects

xtreg $dependente $explicativas i.data, fe

Fixed-effects (within) regression               Number of obs      =      3534
Group variable: banco                           Number of groups   =       186

R-sq:  within  = 0.5679                         Obs per group: min =         1
       between = 0.5260                                        avg =      19.0
       overall = 0.5098                                        max =        38

                                                F(57,3291)         =     75.87
corr(u_i, Xb)  = -0.1572                        Prob > F           =    0.0000

 rho |  .67554064   (fraction of variance due to u_i)
------------------------------------------------------------------------------
F test that all u_i=0:     F(185, 3291) =    15.95           Prob > F = 0.0000


testparm i.data

F( 36,  3291) =    2.46
            Prob > F =    0.0000

Wald test

xttest3

Modified Wald test for groupwise heteroskedasticity
in fixed effect regression model

H0: sigma(i)^2 = sigma^2 for all i

chi2 (186)  =   1.5e+32
Prob>chi2 =      0.0000

Wooldridge test

xtserial $dependente $explicativas, output

Linear regression                                      Number of obs =    3294
                                                       F( 21,   177) =   44.75
                                                       Prob > F      =  0.0000
                                                       R-squared     =  0.5291
                                                       Root MSE      =  7.7268

                                (Std. Err. adjusted for 178 clusters in banco)

Wooldridge test for autocorrelation in panel data
H0: no first-order autocorrelation
    F(  1,     170) =     14.070
           Prob > F =      0.0002

. 
. 
. 
. xtserial $dependente $explicativas

Wooldridge test for autocorrelation in panel data
H0: no first-order autocorrelation
    F(  1,     170) =     14.070
           Prob > F =      0.0002

Trying MIXED EFFECT MODEL

 xtmixed $dependente $explicativas|| banco:

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -13101.186  
Iteration 1:   log likelihood = -13101.186  

Computing standard errors:

Mixed-effects ML regression                     Number of obs      =      3534
Group variable: banco                           Number of groups   =       186

                                                Obs per group: min =         1
                                                               avg =      19.0
                                                               max =        38


                                                Wald chi2(22)      =   4628.55
Log likelihood = -13101.186                     Prob > chi2        =    0.0000


------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
banco: Identity              |
                   sd(_cons) |   9.789867   .5879408       8.70276    11.01277
-----------------------------+------------------------------------------------
                sd(Residual) |   9.159925   .1122329      8.942573     9.38256
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) =  1512.54 Prob >= chibar2 = 0.0000

. 
. 
. 
. estat ic

-----------------------------------------------------------------------------
       Model |    Obs    ll(null)   ll(model)     df          AIC         BIC
-------------+---------------------------------------------------------------
           . |   3534           .   -13101.19     25     26252.37    26406.63
-----------------------------------------------------------------------------
               Note:  N=Obs used in calculating BIC; see [R] BIC note