manipulating stored results: random variance standard error and lower & upper bounds from r(table)? (mixed procedure)

[Wendy Hartogensis]

Hello,
I use Stata version 16 on Windows 10. I'm running models using the mixed procedure, and then trying to work with elements from the r(table) stored results. My outcome is scores, and my predictor is timepoint (I have data from three time points on students nested within classes). The original output table shows estimated variance and standard error of the random variance, but the r(table) results are transformed. I know I need to exponentiate and square the first row of r(table) results to obtain the random variance estimate in the form it appears in the main Stata output. What I can't figure out is how to transform the second row of the r(table) output to obtain the standard error of the random variance estimate in the form it appears in the main Stata output, or the 5th and 6th rows to obtain the confidence bounds for the variance estimate in the form they appears in the main Stata output for the model. Can someone please educate me? Below is some code, with comments showing what works and what does not.

mixed scores i.timepoint || class: || student:
matrix define rtable = r(table)  // r(table) is a 9 x 7 matrix (columns 5,6,7 are random effects, and row 2 is StdErr, row 5&6 are ll&ul of CI)
matlist r(table)
matrix define rtable = r(table)
display exp(2*rtable[1,5])   // this correctly displays the random intercept variance estimate for class
* the above estimate can also be shown as .display exp(2*_b[lns1_1_1:_cons]
display exp(2*rtable[2,5])   // this line does NOT give the estimated standard error corresponding to the random class intercept
display exp(2*rtable[5,5])   // this line does NOT give the estimated lower confidence bound corresponding to the random class intercept
display exp(2*rtable[6,5])   // this line does NOT give the estimated upper confidence bound corresponding to the random class intercept

Can someone please tell me how to use the r(table) results to obtain the values of random variance standard error, and lower and upper confidence bounds for the random variance estimates, from a mixed model (in the form in which they appear in the main model output)? I'm sorry if this should be obvious--and I tried to search for this but was having no luck (if you want to educate me about search terms that would have helped here, I welcome that feedback as well!).

[Scott Merryman]

The standard error used in the confidence intervals can be found in (for the first variance component) _se[lns1_1_1:_cons] . The confidence intervals can be calculated from the results returned _b[] and _se[]. If you want the reported standard error and confidence intervals it is probably easiest to use _diparm or one could construct it from the results of -testnl-

webuse productivity,clear
mixed gsp private || region: || state:, nogroup noheader  nofetable  nolog

display exp(_b[lns1_1_1:_cons] - invnormal(.025)*_se[lns1_1_1:_cons])^2
display exp(_b[lns1_1_1:_cons] + invnormal(.025)*_se[lns1_1_1:_cons])^2

_diparm lns1_1_1, f(exp(@)^2) d(2*exp(@)^2)
return list

qui testnl (exp(_b[lns1_1_1:_cons]))^2 = 0
//G is the derivatives of (exp(_b[lns1_1_1:_cons]))^2) with respect to b
mat a= (r(G)*e(V)*r(G)')
mat b = sqrt(a[1,1])
mat list b

. webuse productivity,clear
(Public Capital Productivity)

. mixed gsp private || region: || state:, nogroup noheader  nofetable  nolog

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
region: Identity             |
                  var(_cons) |   .0441739   .0282398      .0126184    .1546421
-----------------------------+------------------------------------------------
state: Identity              |
                  var(_cons) |   .0688075   .0159875      .0436374    .1084958
-----------------------------+------------------------------------------------
               var(Residual) |   .0038104   .0001947      .0034472    .0042118
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 2196.88               Prob > chi2 = 0.0000

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

. 
. display exp(_b[lns1_1_1:_cons] - invnormal(.025)*_se[lns1_1_1:_cons])^2
.15464207

. display exp(_b[lns1_1_1:_cons] + invnormal(.025)*_se[lns1_1_1:_cons])^2
.01261841

. 
. _diparm lns1_1_1, f(exp(@)^2) d(2*exp(@)^2)
   /lns1_1_1 |   .0441739   .0282398                      .0126184    .1546421

. return list

scalars:
                r(cns) =  0
               r(crit) =  1.959963984540054
                 r(df) =  .
                  r(z) =  .
                  r(p) =  .
                  r(i) =  0
                 r(ub) =  .1546420664868232
                 r(lb) =  .0126184068705908
                r(est) =  .0441739347833051
                 r(se) =  .0282397951523927

. 
. qui testnl (exp(_b[lns1_1_1:_cons]))^2 = 0

. //G is the derivatives of (exp(_b[lns1_1_1:_cons]))^2) with respect to b
. mat a= (r(G)*e(V)*r(G)')

. mat b = sqrt(a[1,1])

. mat list b

symmetric b[1,1]
          c1
r1  .0282398

. 
end of do-file

[Wendy Hartogensis]

Brilliant--thanks so much! This solves my problem several ways and expands my understanding too. I appreciate your help, Scott Merryman!