estimates of variances in esttab for all the levels

[Charles Ngiendo]

I'm running a multilevel mixed model with:
eststo: xtmixed y x1 x2 || country: || household:, reml
eststo: xtmixed y x1 x2 x3 x4|| country: || household:, reml

I want to obtain the estimates of variances in all the levels recorded in the esttab
This is the code I'm using now but it does not give me all the variances (from the random effects parameters), both at level 3, level2, and level 1
esttab using file.rtf,se bic aic eqlabels("" "sd(country)" "sd(_cons)" "sd(Residual)", none) title(Regression Output) nonumbers mtitles("Model 1" "Model 2") addnote("Source: Here")

[Andrew Musau]

estout is from SSC (FAQ Advice #12). Note that xtmixed is renamed mixed in recent versions of Stata. The mixed coefficients are the logged standard deviations of the random effects. You need to take the exponent and then square these coefficients to get the variances. You can then save them as scalars in the estimates using estadd.

webuse productivity, clear
mixed gsp private emp hwy water other unemp || region: || state:, reml
foreach v in s1_1_1 s2_1_1 sig_e{
    estadd scalar  v`v'= (exp(`=e(b)[1, "ln`v':_cons"]'))^2
}
esttab ., keep(`e(depvar)':) scalars(vs1_1_1 vs2_1_1 vsig_e)

Res.:

. mixed gsp private emp hwy water other unemp || region: || state:, reml

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood =  1404.7101  
Iteration 1:   log restricted-likelihood =  1404.7101  

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =        816

-------------------------------------------------------------
                |     No. of       Observations per Group
 Group Variable |     Groups    Minimum    Average    Maximum
----------------+--------------------------------------------
         region |          9         51       90.7        136
          state |         48         17       17.0         17
-------------------------------------------------------------

                                                Wald chi2(6)      =   18382.38
Log restricted-likelihood =  1404.7101          Prob > chi2       =     0.0000

------------------------------------------------------------------------------
         gsp |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     private |   .2660308   .0215471    12.35   0.000     .2237993    .3082623
         emp |   .7555059   .0264556    28.56   0.000     .7036539     .807358
         hwy |   .0718857   .0233478     3.08   0.002     .0261249    .1176465
       water |   .0761552   .0139952     5.44   0.000     .0487251    .1035853
       other |  -.1005396   .0170173    -5.91   0.000    -.1338929   -.0671863
       unemp |  -.0058815   .0009093    -6.47   0.000    -.0076636   -.0040994
       _cons |   2.126995   .1574865    13.51   0.000     1.818327    2.435663
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
region: Identity             |
                  var(_cons) |   .0018964   .0016225      .0003545    .0101438
-----------------------------+------------------------------------------------
state: Identity              |
                  var(_cons) |   .0064439   .0015334      .0040419    .0102732
-----------------------------+------------------------------------------------
               var(Residual) |   .0013543   .0000695      .0012247    .0014976
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 1162.40               Prob > chi2 = 0.0000

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

. esttab ., keep(`e(depvar)':) scalars(vs1_1_1 vs2_1_1 vsig_e)

----------------------------
                      (1)  
                      gsp  
----------------------------
gsp                        
private             0.266***
                  (12.35)  

emp                 0.756***
                  (28.56)  

hwy                0.0719**
                   (3.08)  

water              0.0762***
                   (5.44)  

other              -0.101***
                  (-5.91)  

unemp            -0.00588***
                  (-6.47)  

_cons               2.127***
                  (13.51)  
----------------------------
N                     816  
vs1_1_1           0.00190  
vs2_1_1           0.00644  
vsig_e            0.00135  
----------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

.

Also see the -transform- option of estout which may allow you to directly transform the coefficients to variances.

[Charles Ngiendo]

Thank you so much Andrew. Attached is the output. Would you kindly help me interpret the variance part: random effects parameters / part just to double check.

[Charles Ngiendo]

Also, the independent variable values ranges from zero (€0) to €700000. I wanted to convert them into log for ease of interpretation but the zero values are rendered missing, how do I handle that? Is there any way out?

[Andrew Musau]

Would you kindly help me interpret the variance part: random effects parameters / part just to double check.

As I said, the coefficients that mixed outputs are the logged standard deviations of the estimated random effects parameters. You need to take the exponent and square them to obtain the variances. In a three-level random intercept model, e.g. #2 above, 1_1_1 is the variance of level 3 errors (region) \(\approx\) 0.00190, 2_1_1 is the variance of level 2 errors (state) \(\approx\) 0.00644, and sig_e is the variance of the residuals (or level 1 errors) \(\approx\) 0.00135.

Also, the independent variable values ranges from zero (€0) to €700000. I wanted to convert them into log for ease of interpretation but the zero values are rendered missing, how do I handle that? Is there any way out?

Just rescale the variable by dividing by 1000. Then the new range will be 0-700 and the amounts will be in thousands of euros.

[Charles Ngiendo]

I get you, thanks

[Charles Ngiendo]

Andrew, let me ask again. I did the scaling on the variable but how then do I make the translation be in percentages. For instance a change in x contributes to change in y by x per cent, if not by converting log form?

[Andrew Musau]

The log transformation is only valid for a variable which takes on strictly positive values. After estimation, see

help margins

and the -eydx()- option.

[Charles Ngiendo]

Ok. Thanks. Once again, If I also need to add the likelihood-ration test: LR chi2(2) and prob>chi2, what code do I have to add on the esttab?

[Andrew Musau]

webuse productivity, clear
mixed gsp private emp hwy water other unemp || region: || state:, reml
estadd local Chi2c = cond(e(p_c)<0.01,"`:di %5.3f `=e(chi2_c)''***", ///
cond(e(p_c)<0.05,"`:di %5.3f `=e(chi2_c)''**", ///
cond(e(p_c)<0.1,"`:di %5.3f `=e(chi2_c)''*",  "`:di %5.3f `=e(chi2_c)'''")))
esttab ., keep(`e(depvar)':) stats(N Chi2c, labels(Observations "Chi2-LR stat.")) modelwidth(20) noabbrev

Res.:

. mixed gsp private emp hwy water other unemp || region: || state:, reml

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log restricted-likelihood =  1404.7101  
Iteration 1:   log restricted-likelihood =  1404.7101  

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =        816

-------------------------------------------------------------
                |     No. of       Observations per Group
 Group Variable |     Groups    Minimum    Average    Maximum
----------------+--------------------------------------------
         region |          9         51       90.7        136
          state |         48         17       17.0         17
-------------------------------------------------------------

                                                Wald chi2(6)      =   18382.38
Log restricted-likelihood =  1404.7101          Prob > chi2       =     0.0000

------------------------------------------------------------------------------
         gsp |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     private |   .2660308   .0215471    12.35   0.000     .2237993    .3082623
         emp |   .7555059   .0264556    28.56   0.000     .7036539     .807358
         hwy |   .0718857   .0233478     3.08   0.002     .0261249    .1176465
       water |   .0761552   .0139952     5.44   0.000     .0487251    .1035853
       other |  -.1005396   .0170173    -5.91   0.000    -.1338929   -.0671863
       unemp |  -.0058815   .0009093    -6.47   0.000    -.0076636   -.0040994
       _cons |   2.126995   .1574865    13.51   0.000     1.818327    2.435663
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
region: Identity             |
                  var(_cons) |   .0018964   .0016225      .0003545    .0101438
-----------------------------+------------------------------------------------
state: Identity              |
                  var(_cons) |   .0064439   .0015334      .0040419    .0102732
-----------------------------+------------------------------------------------
               var(Residual) |   .0013543   .0000695      .0012247    .0014976
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 1162.40               Prob > chi2 = 0.0000

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

. estadd local Chi2c = cond(e(p_c)<0.01,"`:di %5.3f `=e(chi2_c)''***", ///
> cond(e(p_c)<0.05,"`:di %5.3f `=e(chi2_c)''**", ///
> cond(e(p_c)<0.1,"`:di %5.3f `=e(chi2_c)''*",  "`:di %5.3f `=e(chi2_c)'''")))

added macro:
              e(Chi2c) : "1162.403***"

. esttab ., keep(`e(depvar)':) stats(N Chi2c, labels(Observations "Chi2-LR stat.")) modelwidth(20) noabbrev

------------------------------------
                              (1)   
                              gsp   
------------------------------------
gsp                                 
private                     0.266***
                          (12.35)   

emp                         0.756***
                          (28.56)   

hwy                        0.0719** 
                           (3.08)   

water                      0.0762***
                           (5.44)   

other                      -0.101***
                          (-5.91)   

unemp                    -0.00588***
                          (-6.47)   

_cons                       2.127***
                          (13.51)   
------------------------------------
Observations                  816   
Chi2-LR stat.          1162.403***   
------------------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

[Charles Ngiendo]

awesome Musa!

[Charles Ngiendo]

but wait now the random effects components are missing

[Andrew Musau]

Append the codes. In #10, I am just providing an answer to what you asked in #9.

[Charles Ngiendo]

Thanks man

[Charles Ngiendo]

I'm getting incomplete output, where is the error?

foreach v in s1_1_1 s2_1_1 sig_e{
estadd scalar v`v'= (exp(`=e(b)[1, "ln`v':_cons"]'))^2
} ,
local Chi2c = cond(e(p_c)<0.01,"`:di %5.3f `=e(chi2_c)''***", ///
cond(e(p_c)<0.05,"`:di %5.3f `=e(chi2_c)''**", ///
cond(e(p_c)<0.1,"`:di %5.3f `=e(chi2_c)''*", "`:di %5.3f `=e(chi2_c)'''")))

esttab using table_1ji.rtf, keep(`e(depvar)' scalars(vs1_1_1 vs2_1_1 vsig_e) stats(N Chi2c, labels(Observations "Chi2-LR stat.")) modelwidth(20) noabbrev