Can -margins- after -mixed- with small sample inference show t-tests rather than z-tests?

[Bruce Weaver]

Hello everyone. After using one of the dfmethod() options for -mixed- to obtain small-sample inference, I see t-tests on the coefficients and an F-test for the overall model, as expected. But when I follow up with a -margins- command, I see z-tests on the fitted values and CIs that have been computed using the critical z-value rather than the critical t-value. A bit of digging led me to one of Joseph Coveney's old posts from which I learned about the small option for -contrast-. But as far as I can tell, there is no similar option for -margins- following -mixed-. I hope I am wrong about that, and that someone can show me how to get t-tests in the -margins- output. I've pasted below code for a variation on Example 13 in the documentation. And I am using Stata 16 (for Windows), by the way.

Thanks,
Bruce

// Example 13 in documentation for -mixed-.
clear
use https://www.stata-press.com/data/r16/t43
// Use dfemethod(repeated) option to show t-and F-tests
// rather than z- and Chi-square tests
mixed score i.drug || person:, reml dfmethod(repeated)
// Replaying the model will show the default large-sample
// z and Chi-square tests
mixed
// The -contrast- command has a small option:
contrast drug
contrast drug, small
// But apparently, -margins- does not have a small option:
margins drug
*margins drug, small
// The previous line is commented out because it causes an error.
// ----------------------------------------------------
// Q. Is there an option to make -margins- show t-tests
// and CIs computed with critical t-values?  
// ----------------------------------------------------

[Andrew Musau]

I think the -df()- option achieves what you want.

df(#) use t distribution with # degrees of freedom for computing p-values and confidence intervals

clear
use https://www.stata-press.com/data/r16/t43
// Use dfemethod(repeated) option to show t-and F-tests
// rather than z- and Chi-square tests
mixed score i.drug || person:, reml dfmethod(repeated)
local df_m= e(df_m)
margins drug, df(`df_m')

Res.:

. mixed score i.drug || person:, reml dfmethod(repeated)

Performing EM optimization: 

Performing gradient-based optimization: 

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

Computing standard errors:

Computing degrees of freedom:

Mixed-effects REML regression                   Number of obs     =         20
Group variable: person                          Number of groups  =          5

                                                Obs per group:
                                                              min =          4
                                                              avg =        4.0
                                                              max =          4
DF method: Repeated                             DF:           min =       4.00
                                                              avg =      10.00
                                                              max =      12.00

                                                F(3,    12.00)    =      24.76
Log restricted-likelihood = -49.640099          Prob > F          =     0.0000

------------------------------------------------------------------------------
       score |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        drug |
          2  |        -.8   1.939072    -0.41   0.687    -5.024874    3.424874
          3  |      -10.8   1.939072    -5.57   0.000    -15.02487   -6.575126
          4  |        5.6   1.939072     2.89   0.014     1.375126    9.824874
             |
       _cons |       26.4   3.149604     8.38   0.001      17.6553     35.1447
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
person: Identity             |
                  var(_cons) |   40.20004   30.10272      9.264606    174.4319
-----------------------------+------------------------------------------------
               var(Residual) |   9.399997   3.837532       4.22305    20.92325
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 15.03         Prob >= chibar2 = 0.0001

. 
. local df_m= e(df_m)

. 
. margins drug, df(`df_m')

Adjusted predictions                            Number of obs     =         20

Expression   : Linear prediction, fixed portion, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        drug |
          1  |       26.4   3.149604     8.38   0.004     16.37655    36.42345
          2  |       25.6   3.149604     8.13   0.004     15.57655    35.62345
          3  |       15.6   3.149604     4.95   0.016     5.576554    25.62345
          4  |         32   3.149604    10.16   0.002     21.97655    42.02345
------------------------------------------------------------------------------

.

[Bruce Weaver]

Thanks Andrew Musau. I had not noticed the df() option, and that is certainly a step in the right direction. But notice that the p-value and CI for the constant in the model output do not match the p-value and CI for the drug 1 margins.

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        drug |
          1  |     26.400      3.150    8.382   0.0036       16.377      36.423
       _cons |     26.400      3.150    8.382   0.0011       17.655      35.145
------------------------------------------------------------------------------

But with a bit of digging, I discovered that if I use e(df_min) rather than e(df_m), the results match.

. clear

. use https://www.stata-press.com/data/r16/t43
(T4.3 -- Winer, Brown, Michels)

. mixed score i.drug || person:, reml dfmethod(repeated)

Performing EM optimization:

Performing gradient-based optimization:

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

Computing standard errors:

Computing degrees of freedom:

Mixed-effects REML regression                   Number of obs     =         20
Group variable: person                          Number of groups  =          5

                                                Obs per group:
                                                              min =          4
                                                              avg =        4.0
                                                              max =          4
DF method: Repeated                             DF:           min =       4.00
                                                              avg =      10.00
                                                              max =      12.00

                                                F(3,    12.00)    =      24.76
Log restricted-likelihood = -49.640099          Prob > F          =     0.0000

------------------------------------------------------------------------------
       score |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        drug |
          1  |      0.000  (base)
          2  |     -0.800      1.939   -0.413   0.6872       -5.025       3.425
          3  |    -10.800      1.939   -5.570   0.0001      -15.025      -6.575
          4  |      5.600      1.939    2.888   0.0136        1.375       9.825
             |
       _cons |     26.400      3.150    8.382   0.0011       17.655      35.145
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
person: Identity             |
                  var(_cons) |     40.200     30.103         9.265     174.432
-----------------------------+------------------------------------------------
               var(Residual) |      9.400      3.838         4.223      20.923
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 15.03         Prob >= chibar2 = 0.0001

. local df_m = e(df_m)

. local df_min = e(df_min)

. margins drug, df(`df_m')

Adjusted predictions                            Number of obs     =         20

Expression   : Linear prediction, fixed portion, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        drug |
          1  |     26.400      3.150    8.382   0.0036       16.377      36.423
          2  |     25.600      3.150    8.128   0.0039       15.577      35.623
          3  |     15.600      3.150    4.953   0.0158        5.577      25.623
          4  |     32.000      3.150   10.160   0.0020       21.977      42.023
------------------------------------------------------------------------------

. // The p-value & CI from -margins- do not match those for the constant
. // Try with df_min rather than df_m
. margins drug, df(`df_min')

Adjusted predictions                            Number of obs     =         20

Expression   : Linear prediction, fixed portion, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        drug |
          1  |     26.400      3.150    8.382   0.0011       17.655      35.145
          2  |     25.600      3.150    8.128   0.0012       16.855      34.345
          3  |     15.600      3.150    4.953   0.0077        6.855      24.345
          4  |     32.000      3.150   10.160   0.0005       23.255      40.745
------------------------------------------------------------------------------

Here's the code.

clear
use https://www.stata-press.com/data/r16/t43
mixed score i.drug || person:, reml dfmethod(repeated)
local df_m = e(df_m)
local df_min = e(df_min)
margins drug, df(`df_m')
// The p-value & CI from -margins- do not match those for the constant
// Try with df_min rather than df_m
margins drug, df(`df_min')

Thanks so much for pointing me in the right direction, Andrew.

[Bruce Weaver]

Here's a little more code that might help someone perusing the archives some day.

// Confirm that e(df_min) works with other dfmethod settings
local method "repeated satterthwaite kroger"
foreach m of local method {
  display "`m'"
  mixed score i.drug || person:, reml dfmethod("`m'")
  local df_min = e(df_min)
  margins drug, df(`df_min')
}

[Tiago Pereira]

Thank you, Bruce. Very useful.