lincom with random-effects parameters

[Stefan Kreisel]

I have no idea if this is statistically viable, but:

Post mixed of the following sort

[...]

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
_all: Identity               |
                var(R._time) |    .000703   .0006863      .0001038    .0047634
-----------------------------+------------------------------------------------
idpatient: Identity          |
                  var(_cons) |   1.38e-25   3.76e-25      6.70e-28    2.86e-23
-----------------------------+------------------------------------------------
               var(Residual) |   .0233767   .0012138      .0211148     .025881
------------------------------------------------------------------------------

one could be tempted to add variances together for further calculation - such as when one is interested in intraclass correlation coefficients as in this Statalist post here.

Adding 0.000703 + 0.0233767 = 0.0240797 is fun and easy; but adding the confidence intervals!?

lincom is the obvious option. However, Stata doesn't make it easy for you to get the correct naming of the random-effect parameters. Doing a

mixed outcome  || _all:R.time || patient: , reml var coeflegend

doesn't do the trick - there's no specific legending output for the random effects part.

But

matrix list e(b)

e(b)[1,4]
         Cwalk:   lns1_1_1:   lns2_1_1:    lnsig_e:
         _cons       _cons       _cons       _cons
y1   .00131434  -3.6300429  -28.619886  -1.8780069

does - but that's unintelligible.

1. What does lns2_1_1 stand for?
2. And: one can't related to the output'ed parameters - they are nowhere to be found in the orginal mixed result.

https://stats.oarc.ucla.edu/stata/fa...using-_diparm/ tells you why:

You can access these values using the undocumented command _diparm (which stands for display parameter). [...] To use _diparm you have to understand how Stata computes the random effects. Stata computes the variances as the log of the standard deviation (ln_sigma) and computes covariances as the arc hyperbolic tangent of the correlation.

So that means one gets to do the following, that only the most tenured Stata'ist will understand:

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

That actually gives you the original output - Eureka!

_diparm lns1_1_1, f(exp(@)^2) d(2*exp(@)^2)
/lns1_1_1 |    .000703   .0006863                      .0001038    .0047634

And now the final act:

lincom [lns1_1_1]_cons + [lnsig_e]_cons

 ( 1)  [lns1_1_1]_cons + [lnsig_e]_cons = 0

------------------------------------------------------------------------------
       Cwalk |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   -5.50805   .4884959   -11.28   0.000    -6.465484   -4.550615
------------------------------------------------------------------------------

But that get's me into a dead end and I quit...
I have no idea how to get from here back to the "real" sum of variances (ie 0.0240797) and the actual goal: what are my confidence intervals?

Cheers

[Clyde Schechter]

Because the variance components are not linear in the parameters that Stata has actually used for the estimation, you cannot get a sum of variances with -lincom-. Use -nlcom- instead, with the appropriate exponential transformations:

. webuse pig, clear
(Longitudinal analysis of pig weights)

.
. mixed weight week || id:

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:  Log likelihood = -1014.9268  
Iteration 1:  Log likelihood = -1014.9268  

Computing standard errors ...

Mixed-effects ML regression                        Number of obs    =      432
Group variable: id                                 Number of groups =       48
                                                   Obs per group:
                                                                min =        9
                                                                avg =      9.0
                                                                max =        9
                                                   Wald chi2(1)     = 25337.49
Log likelihood = -1014.9268                        Prob > chi2      =   0.0000

------------------------------------------------------------------------------
      weight | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        week |   6.209896   .0390124   159.18   0.000     6.133433    6.286359
       _cons |   19.35561   .5974059    32.40   0.000     18.18472    20.52651
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   14.81751   3.124225      9.801716    22.40002
-----------------------------+------------------------------------------------
               var(Residual) |   4.383264   .3163348      3.805112     5.04926
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 472.65        Prob >= chibar2 = 0.0000

.
. matrix list e(b)

e(b)[1,4]
       weight:    weight:  lns1_1_1:   lnsig_e:
         week      _cons      _cons      _cons
y1  6.2098958  19.355613  1.3479048  .73889679

.
. nlcom exp(2*_b[lns1_1_1:_cons]) + exp(2*_b[lnsig_e:_cons])

       _nl_1: exp(2*_b[lns1_1_1:_cons]) + exp(2*_b[lnsig_e:_cons])

------------------------------------------------------------------------------
      weight | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |   19.20077   3.136657     6.12   0.000     13.05304    25.34851
------------------------------------------------------------------------------

[Stefan Kreisel]

Crass! Thanks Clyde Schechter, I thought it could be something like that, but wouldn't dare imagine how ...

So, to complete

nlcom exp(2*[lns1_1_1]_cons) + exp(2*[lnsig_e]_cons)

       _nl_1:  exp(2*[lns1_1_1]_cons) + exp(2*[lnsig_e]_cons)

------------------------------------------------------------------------------
       Cwalk |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _nl_1 |   .0240798   .0013877    17.35   0.000       .02136    .0267996
------------------------------------------------------------------------------

is exactly 0.000703 + 0.0233767 = 0.0240797 as in post #1!

Just out of curiosity and astonishment, a couple of follow-ups:

• Wonder if I'll ever need _diparm lns1_1_1, f(exp(@)^2) d(2*exp(@)^2)again!? It was kind of fun discovering "that" it exists ...
• Clyde's way of writing the parameters as in _b[lns1_1_1:_cons] works just like [lns1_1_1]_cons : wasn't aware that there are differenct syntaxes for that; albeit: the whole behind-the-scenes way how Stata labels things is a bit of a mystery, or is it?
• Clyde Schechter: While I get why one has to exp(parameter), I'm a bit lost on why it is exp(2*parameter)? Just out of curiosity from a non-statistician ...

Regards!

[Clyde Schechter]

While I get why one has to exp(parameter), I'm a bit lost on why it is exp(2*parameter)? Just out of curiosity from a non-statistician ...

So lns (and lnsig) stand for ln(standard deviation). The variance is the square of the standard deviation. And by the standard laws of logarithms and exponents, to get from the log of the standard deviation to the (untransformed) variance, you do exp(2*lns).

_b[lns1_1_1:_cons] works just like [lns1_1_1]_cons : wasn't aware that there are differenct syntaxes for that;

Some commands, but not all, allow simplified syntax for coefficients. But the _b[...] notation is the standard one that works in all commands that use coefficients. Rather than trying to remember which commands allow the shorter version(s), I just routinely use _b[...].