how to use test postestimation after contrast with ar. following mixed

[Gunther Meinlschmidt]

Hello,

I am using Stata/SE 15.1.

I have difficulties in running a post-estimation test using "test", comparing two reverse adjacent contrasts estimated using "contrast" following "mixed".

The dependent variable of interest is BMQappre
The variable time has three factors: 0, 1, 2 reflecting t0, t1, t2
The variables studySite, BPTgroup, LFD describe the structure of the data, modelled as random effects.

I want to test, whether the difference in BMQappre from t1 to t2 is equal to the difference from t0 to t1, but I have difficulties in correctly pulling out the reverse adjacent contrast parameters (see output below).

Here is the

*** ==> The mixed command runs the mixed model and works fine
mixed BMQappre i.time if ${all_manuscript01A} || studySite: || BPTgroup: || LFD: , reml dfmethod(kroger, eim)

*** ==> The contrast command contrasts t2 with t1 and t1 with t0 using the reverse adjacent contrast prefix "ar."; this also looks to me as working fine.
contrast ar.time, pveffects post nofvlabel

*** ==> Now, I want to compare these two reverse adjacent contrasts using "test", but apparently I do something wrong, when pulling out the parameters reflecting "2 vs 1" and "1 vs 0"; so there is an error message
test ar2.time == ar1.time

Can anyone help me out, how to correctly pull out the parameters for the "test" comparison? (Alternatively, it may be possible to test them directly within the "contrast" command?)

Thank you very much in advance!
Gunther


HERE ARE SOME SNIPPETS OF SAMPLE DATA:

input double BMQappre float(time studySite) double(BPTgroup LFD)
                29 0 1  1 1
                27 1 1  1 1
                32 2 1  1 1
                32 0 1  1 2
                37 1 1  1 2
                42 2 1  1 2
                36 0 1  1 3
                30 1 1  1 3
                40 2 1  1 3
                36 0 1  1 4
                32 1 1  1 4
                42 2 1  1 4
                34 1 1  2 5
                37 2 1  2 5
                40 0 1  2 5
                41 1 1  2 6
                41 2 1  2 6
                33 0 1  2 6
                17 0 1  2 7
                23 1 1  2 7
                23 2 1  2 7
                35 0 2  5 8
                32 1 2  5 8
                35 2 2  5 8
                18 0 2  5 9
                15 1 2  5 9
                16 2 2  5 9

end
label values time label_time
label def label_time 0 "t0", modify
label def label_time 1 "t1", modify
label def label_time 2 "t2", modify

HERE IS THE OUTPUT FROM THE COMPLETE DATA SET (NOT FROM SAMPLE DATASNIPPETS ABOVE):

. *** ==> The mixed command runs the mixed model and works fine

. mixed BMQappre i.time if ${all_manuscript01A} || studySite: || BPTgroup: || LFD: , reml dfmethod(kroger, eim)



Performing EM optimization:



Performing gradient-based optimization:



Iteration 0: log restricted-likelihood = -350.22414

Iteration 1: log restricted-likelihood = -350.19123

Iteration 2: log restricted-likelihood = -350.19098

Iteration 3: log restricted-likelihood = -350.19098



Computing standard errors:



Computing degrees of freedom:



Mixed-effects REML regression Number of obs = 113



-------------------------------------------------------------

| No. of Observations per Group

Group Variable | Groups Minimum Average Maximum

----------------+--------------------------------------------

studySite | 2 48 56.5 65

BPTgroup | 7 14 16.1 18

LFD | 38 2 3.0 3

-------------------------------------------------------------

DF method: Kenward-Roger DF: min = 1.12

avg = 49.08

max = 73.12



F(2, 73.08) = 18.49

Log restricted-likelihood = -350.19098 Prob > F = 0.0000



------------------------------------------------------------------------------

BMQappre | Coef. Std. Err. t P>|t| [95% Conf. Interval]

-------------+----------------------------------------------------------------

time |

t1 | -1.486842 .8458535 -1.76 0.083 -3.172625 .1989409

t2 | 3.576918 .8541278 4.19 0.000 1.874691 5.279146

|

_cons | 29.90408 1.621449 18.44 0.025 13.97035 45.83781

------------------------------------------------------------------------------



------------------------------------------------------------------------------

Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]

-----------------------------+------------------------------------------------

studySite: Identity |

var(_cons) | 4.32e-19 . . .

-----------------------------+------------------------------------------------

BPTgroup: Identity |

var(_cons) | 6.114891 8.854323 .3579764 104.4535

-----------------------------+------------------------------------------------

LFD: Identity |

var(_cons) | 43.69086 12.24502 25.22485 75.67504

-----------------------------+------------------------------------------------

var(Residual) | 13.5939 2.248831 9.82947 18.8

------------------------------------------------------------------------------

LR test vs. linear model: chi2(3) = 78.57 Prob > chi2 = 0.0000



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



.

. *** ==> The contrast command contrasts t2 with t1 and t1 with t0 using the reverse adjacent contrast prefix "ar."; this also looks to me as working fine.

. contrast ar.time, pveffects post nofvlabel



Contrasts of marginal linear predictions



Margins : asbalanced



------------------------------------------------

| df chi2 P>chi2

-------------+----------------------------------

BMQappre |

time |

(1 vs 0) | 1 3.09 0.0788

(2 vs 1) | 1 35.16 0.0000

Joint | 2 37.00 0.0000

------------------------------------------------



-----------------------------------------------------

| Contrast Std. Err. z P>|z|

-------------+---------------------------------------

BMQappre |

time |

(1 vs 0) | -1.486842 .8458535 -1.76 0.079

(2 vs 1) | 5.063761 .8540279 5.93 0.000

-----------------------------------------------------



.

. *** ==> Now, I want to compare these two reverse adjacent contrasts using "test", but apparently I do something wrong, when pulling out the parameters reflecting "2 vs 1" and "1 vs 0"; so there is an error message

. test ar2.time == ar1.time

variable time not found

r(111);

[Clyde Schechter]

To use the -test- (or -lincom-) commands you need to know exactly what the effects are called. The easy way to do that is to -matrix list e(b)- and then pick it up from there:

. *** ==> The contrast command contrasts t2 with t1 and t1 with t0 using the reverse adjacent contrast prefix "ar."; this also looks to me as working fine.
. contrast ar.time, pveffects post nofvlabel

Contrasts of marginal linear predictions

Margins      : asbalanced

------------------------------------------------
             |         df        chi2     P>chi2
-------------+----------------------------------
BMQappre     |
        time |
   (1 vs 0)  |          1        0.13     0.7175
   (2 vs 1)  |          1        7.17     0.0074
      Joint  |          2        8.44     0.0147
------------------------------------------------

-----------------------------------------------------
             |   Contrast   Std. Err.      z    P>|z|
-------------+---------------------------------------
BMQappre     |
        time |
   (1 vs 0)  |  -.5555556   1.535586    -0.36   0.718
   (2 vs 1)  |   4.111111   1.535586     2.68   0.007
-----------------------------------------------------

.
. matrix list e(b)

e(b)[1,2]
      BMQappre:   BMQappre:
        ar1vs0.     ar2vs1.
          time        time
y1  -.55555556   4.1111111

.
. test _b[BMQappre:ar1vs0.time] = _b[BMQappre:ar2vs1.time]

 ( 1)  [BMQappre]ar1vs0.time - [BMQappre]ar2vs1.time = 0

           chi2(  1) =    3.08

        Prob > chi2 =    0.0793

. lincom _b[BMQappre:ar1vs0.time] - _b[BMQappre:ar2vs1.time] // BETTER THAN TEST!

 ( 1)  [BMQappre]ar1vs0.time - [BMQappre]ar2vs1.time = 0

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -4.666667   2.659713    -1.75   0.079    -9.879608     .546275
------------------------------------------------------------------------------

The American Statistical Association has recommended that the concept of statistical significance be abandoned. See https://www.tandfonline.com/doi/full...5.2019.1583913 for the "executive summary" and
https://www.tandfonline.com/toc/utas20/73/sup1 for all 43 supporting articles. Or https://www.nature.com/articles/d41586-019-00857-9 for the tl;dr. Accordingly, I have shown you the -lincom- command to estimate the difference between these two contrasts, rather than testing for equality. This way you get a since of how different the two contrasts are, and an estimate of the precision with which you have estimated the difference. This is much more useful information than a yes/no significance test in almost any context.

Note: The above results were generated from your example data (and omitting the -if- condition which relies on a global macro not defined in the code you showed), so you need to rerun with your actual data--these are not the correct answers.

[Gunther Meinlschmidt]

Hello Clyde,

Thank you very much for your reply – it solved my problem!
And thank you also for pointing out the lincom command that is indeed producing more useful information – highly appreciated!

Best Gunther