Repeated measures ANOVA: Standard errors of estimated marginal means produced by Stata's margins were different from that of SPSS's EMMEANS

[Keowmani Thamron]

Hi,

Why the SE's of estimated marginal means differ in this case? And which one should be used if we want to estimate the factor level means?
(Note: EMMEANS is equivalent to margins, asbalanced when called after fitting other fixed-factor anova models without any repeated factors.)
Stata/SE 15.1 vs IBM SPSS 22

Stata code:
anova y subject time, repeated(time)
margins time

SPSS code:
GLM y.1 y.2 y.3 y.4
/WSFACTOR=j 4 Polynomial
/METHOD=SSTYPE(3)
/EMMEANS=TABLES(time)
/CRITERIA=ALPHA(.05)
/WSDESIGN=time.

* Example generated by -dataex-. To install: ssc install dataex
clear
input float(y subject time)
20 1 1
24 1 2
28 1 3
28 1 4
15 2 1
18 2 2
23 2 3
24 2 4
18 3 1
19 3 2
24 3 3
23 3 4
26 4 1
26 4 2
30 4 3
30 4 4
22 5 1
24 5 2
28 5 3
26 5 4
19 6 1
21 6 2
27 6 3
25 6 4
end

[Bruce Weaver]

Hi Keowmani. I don't (think I) have a complete answer for you. But your post did entice me to play around with this problem for a while, and here's what I've come up with thus far!

* Example generated by -dataex-. To install: ssc install dataex
clear *
input float(y subject time)
20 1 1
24 1 2
28 1 3
28 1 4
15 2 1
18 2 2
23 2 3
24 2 4
18 3 1
19 3 2
24 3 3
23 3 4
26 4 1
26 4 2
30 4 3
30 4 4
22 5 1
24 5 2
28 5 3
26 5 4
19 6 1
21 6 2
27 6 3
25 6 4
end

anova y subject time, repeated(time)
margins time

* Q1. Can I get different SEs for each time point via -mixed-?
mixed y i.time || subject:, noconstant residuals(un, t(time))
margins time
* A1. Yes, I can.  But they do not match the SEs from SPSS GLM.

* Q2. What if I use REML rather than MLE?
mixed y i.time || subject:, noconstant residuals(un, t(time)) reml
margins time
* A2.  The SEs now match those from SPSS GLM.
* But note that -mixed- is using the critical z-value when
* computing the 95% CI, whereas SPSS used the critical t-value.
* I believe it used t-crit with df = n-1, or 5.
* Let's see if we can duplicate the CIs reported by SPSS GLM.
matrix m = r(table)'
svmat m
rename (m1 m2) (yhat se)
generate double lower = yhat - se*invt(5,(1-.05/2))
generate double upper = yhat + se*invt(5,(1-.05/2))
list yhat se lower upper if !missing(yhat)
* These CI's match those from SPSS GLM.

Output from the last -list- command above:

. list yhat se lower upper if !missing(yhat)

     +---------------------------------------------+
     |     yhat         se       lower       upper |
     |---------------------------------------------|
  1. |       20   1.527525   16.073372   23.926628 |
  2. |       22   1.290994   18.681394   25.318606 |
  3. | 26.66667   1.085255    23.87693   29.456402 |
  4. |       26   1.064581   23.263407   28.736593 |
     +---------------------------------------------+

. * These CI's match those from SPSS GLM.

HTH.

[Marcos Almeida]

That seems to be related to the delta method, explained here.

[Bruce Weaver]

The -mean- command can also duplicate the SEs and 95% CIs that SPSS GLM gives. But you have to avoid use of the over() option. If you specify over(time), -mean- uses a critical t-value for df = the total N minus 1. (There's an older thread on this, but I can't find it right now.) What is needed here is df = group n minus 1. This can be achieved by issuing one -mean- command for each level of time with an appropriate -if-. E.g.,

* Using the same data as in #2...
forvalues t = 1/4 {
 display "time = "`t'
 mean y if time==`t'
}

Here's the output:

time = 1

Mean estimation                   Number of obs   =          6

--------------------------------------------------------------
             |       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
           y |         20   1.527525      16.07337    23.92663
--------------------------------------------------------------
time = 2

Mean estimation                   Number of obs   =          6

--------------------------------------------------------------
             |       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
           y |         22   1.290994      18.68139    25.31861
--------------------------------------------------------------
time = 3

Mean estimation                   Number of obs   =          6

--------------------------------------------------------------
             |       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
           y |   26.66667   1.085255      23.87693     29.4564
--------------------------------------------------------------
time = 4

Mean estimation                   Number of obs   =          6

--------------------------------------------------------------
             |       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
           y |         26   1.064581      23.26341    28.73659
--------------------------------------------------------------

Compare these means, SEs and CIs to the SPSS output in #1. They're the same.

[Keowmani Thamron]

Hi Bruce, thank you for your illustration. So, would it be correct for me to say that SPSS assumes unstructured within-subject residuals variance-covariance matrix to estimate the factor level means, whereas Stata assumes exchangeable structure? If it is correct, then, in your opinion, which approach is more appropriate given the repeated-measures study design? Thanks in advance!

[Joseph Coveney]

Why the SE's of estimated marginal means differ in this case? And which one should be used if we want to estimate the factor level means?

Stata's margins is correct.

SPSS is giving you the unconditional standard errors of the means, as if they were independent standalone means. Stata's standard errors account for the between-subjects factor in the repeated-measures ANOVA. Under compound symmetry, you get a common residual variance and if you want to estimate the factor level means (i.e., model-based), then you should use the standard error from Stata's margins and not the unconditional ones that SPSS reports.

.ÿversionÿ15.1

.ÿ
.ÿclearÿ*

.ÿ
.ÿprogramÿdefineÿshowem
ÿÿ1.ÿÿÿÿÿÿÿÿÿversionÿ15.1
ÿÿ2.ÿÿÿÿÿÿÿÿÿsyntax
ÿÿ3.ÿ
.ÿÿÿÿÿÿÿÿÿtempnameÿVarianceÿsem
ÿÿ4.ÿÿÿÿÿÿÿÿÿscalarÿdefineÿ`Variance'ÿ=ÿe(rss)ÿ/ÿr(df_r)
ÿÿ5.ÿÿÿÿÿÿÿÿÿdisplayÿinÿsmclÿasÿtextÿ"varianceÿofÿresidualsÿ=ÿMSerrorÿ=ÿ"ÿ%09.7fÿ`Variance'
ÿÿ6.ÿÿÿÿÿÿÿÿÿdisplayÿinÿsmclÿasÿtextÿ"between-subjectÿcountÿ=ÿ"ÿe(N_bse)
ÿÿ7.ÿ
.ÿÿÿÿÿÿÿÿÿscalarÿdefineÿ`sem'ÿ=ÿsqrt(`Variance'ÿ/ÿe(N_bse))
ÿÿ8.ÿÿÿÿÿÿÿÿÿdisplayÿinÿsmclÿasÿtextÿ_newline(1)ÿ"semÿ(modelÿbased)ÿ=ÿ"ÿ%08.6fÿ`sem'
ÿÿ9.ÿend

.ÿ
.ÿquietlyÿinputÿbyte(yÿsubjectÿtime)

.ÿ
.ÿanovaÿyÿsubjectÿtime,ÿrepeated(time)

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿ=ÿÿÿÿÿÿÿÿÿ24ÿÿÿÿR-squaredÿÿÿÿÿ=ÿÿ0.9571
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿRootÿMSEÿÿÿÿÿÿ=ÿÿÿÿÿ1.0328ÿÿÿÿAdjÿR-squaredÿ=ÿÿ0.9343

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿSourceÿ|ÿPartialÿSSÿÿÿÿÿÿÿÿÿdfÿÿÿÿÿÿÿÿÿMSÿÿÿÿÿÿÿÿFÿÿÿÿProb>F
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-----------+----------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿModelÿ|ÿÿ357.33333ÿÿÿÿÿÿÿÿÿÿ8ÿÿÿ44.666667ÿÿÿÿÿ41.88ÿÿ0.0000
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿsubjectÿ|ÿÿ173.33333ÿÿÿÿÿÿÿÿÿÿ5ÿÿÿ34.666667ÿÿÿÿÿ32.50ÿÿ0.0000
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿtimeÿ|ÿÿÿÿÿÿÿÿ184ÿÿÿÿÿÿÿÿÿÿ3ÿÿÿ61.333333ÿÿÿÿÿ57.50ÿÿ0.0000
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿResidualÿ|ÿÿÿÿÿÿÿÿÿ16ÿÿÿÿÿÿÿÿÿ15ÿÿÿ1.0666667ÿÿ
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-----------+----------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿTotalÿ|ÿÿ373.33333ÿÿÿÿÿÿÿÿÿ23ÿÿÿ16.231884ÿÿ


Between-subjectsÿerrorÿterm:ÿÿsubject
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿLevels:ÿÿ6ÿÿÿÿÿÿÿÿÿ(5ÿdf)
ÿÿÿÿÿLowestÿb.s.e.ÿvariable:ÿÿsubject

Repeatedÿvariable:ÿtime
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿHuynh-Feldtÿepsilonÿÿÿÿÿÿÿÿ=ÿÿ0.9270
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿGreenhouse-Geisserÿepsilonÿ=ÿÿ0.6038
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿBox'sÿconservativeÿepsilonÿ=ÿÿ0.3333

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ------------ÿProbÿ>ÿFÿ------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿSourceÿ|ÿÿÿÿÿdfÿÿÿÿÿÿFÿÿÿÿRegularÿÿÿÿH-FÿÿÿÿÿÿG-GÿÿÿÿÿÿBox
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-----------+----------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿtimeÿ|ÿÿÿÿÿÿ3ÿÿÿÿ57.50ÿÿÿ0.0000ÿÿÿ0.0000ÿÿÿ0.0000ÿÿÿ0.0006
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿResidualÿ|ÿÿÿÿÿ15
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿ----------------------------------------------------------------

.ÿmarginsÿtime

PredictiveÿmarginsÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿÿÿÿ24

Expressionÿÿÿ:ÿLinearÿprediction,ÿpredict()

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿDelta-method
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿMarginÿÿÿStd.ÿErr.ÿÿÿÿÿÿtÿÿÿÿP>|t|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿtimeÿ|
ÿÿÿÿÿÿÿÿÿÿ1ÿÿ|ÿÿÿÿÿÿÿÿÿ20ÿÿÿÿ.421637ÿÿÿÿ47.43ÿÿÿ0.000ÿÿÿÿÿÿ19.1013ÿÿÿÿÿ20.8987
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿÿÿÿÿÿÿÿ22ÿÿÿÿ.421637ÿÿÿÿ52.18ÿÿÿ0.000ÿÿÿÿÿÿ21.1013ÿÿÿÿÿ22.8987
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿÿ26.66667ÿÿÿÿ.421637ÿÿÿÿ63.25ÿÿÿ0.000ÿÿÿÿÿ25.76797ÿÿÿÿ27.56536
ÿÿÿÿÿÿÿÿÿÿ4ÿÿ|ÿÿÿÿÿÿÿÿÿ26ÿÿÿÿ.421637ÿÿÿÿ61.66ÿÿÿ0.000ÿÿÿÿÿÿ25.1013ÿÿÿÿÿ26.8987
------------------------------------------------------------------------------

.ÿ
.ÿshowem
varianceÿofÿresidualsÿ=ÿMSerrorÿ=ÿ1.0666667
between-subjectÿcountÿ=ÿ6

semÿ(modelÿbased)ÿ=ÿ0.421637

.ÿ
.ÿexit

endÿofÿdo-file


.

SPSS's EMMEANS ignores the model that it just fit. You should ignore SPSS.

[Joseph Coveney]

SPSS's EMMEANS ignores the model that it just fit.

You don't show the output, but are you fitting a MANOVA with SPSS and not a repeated-measures ANOVA? If so, then I take that back.

[Keowmani Thamron]

You don't show the output, but are you fitting a MANOVA with SPSS and not a repeated-measures ANOVA? If so, then I take that back.

In SPSS, when repeated measures ANOVA is requested the output include both MANOVA and univariate ANOVA results. My interest is in fitting repeated measures ANOVA, not MANOVA. I will ignore the MANOVA results.

[Bruce Weaver]

In #6, Joseph Coveney wrote:

Stata's margins is correct.

SPSS is giving you the unconditional standard errors of the means, as if they were independent standalone means. Stata's standard errors account for the between-subjects factor in the repeated-measures ANOVA. Under compound symmetry, you get a common residual variance and if you want to estimate the factor level means (i.e., model-based), then you should use the standard error from Stata's margins and not the unconditional ones that SPSS reports.

Joseph, this suggests (to me) that one ought to be able to generate the same marginal means and SEs after using -mixed- to estimate the model. But my attempts (thus far) to do so have failed. For example, this code....

* Joseph's post (#6) suggests that -mixed- with compound symmetry
* should yield the same results (more or less) as RM ANOVA.
* Stata has no cs option, but I believe exchangeable is the same thing.
mixed y i.time || subject:, noconstant residuals(exch, t(time)) reml
margins time

gives this table of marginal means:

Adjusted predictions                            Number of obs     =         24

Expression   : Linear prediction, fixed portion, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        time |
          1  |         20   1.256096    15.92   0.000      17.5381     22.4619
          2  |         22   1.256096    17.51   0.000      19.5381     24.4619
          3  |   26.66667   1.256096    21.23   0.000     24.20476    29.12857
          4  |         26   1.256096    20.70   0.000      23.5381     28.4619
------------------------------------------------------------------------------

Do you have any thoughts on how to tweak the -mixed- command to make -margins- yield SEs of .421637, matching what Keowmani got using -anova-?

Thanks,
Bruce

[Keowmani Thamron]

Hi Bruce,

Stata's

margins

is correct.

SPSS is giving you the unconditional standard errors of the means, as if they were independent standalone means. Stata's standard errors account for the between-subjects factor in the repeated-measures ANOVA. Under compound symmetry, you get a common residual variance and if you want to estimate the factor level means (i.e., model-based), then you should use the standard error from Stata's

margins

and not the unconditional ones that SPSS reports.

To get the SE similar to that produced by anova, you can use the following code below. However, please note that the CI is calculated based on z-distribution not t-distribution.

. xtset subject
       panel variable:  subject (balanced)

. quietly xtreg y i.time, fe

. margins time

Adjusted predictions                            Number of obs     =         24
Model VCE    : Conventional

Expression   : Linear prediction, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        time |
          1  |         20    .421637    47.43   0.000     19.17361    20.82639
          2  |         22    .421637    52.18   0.000     21.17361    22.82639
          3  |   26.66667    .421637    63.25   0.000     25.84027    27.49306
          4  |         26    .421637    61.66   0.000     25.17361    26.82639
------------------------------------------------------------------------------

[Joseph Coveney]

Do you have any thoughts on how to tweak the -mixed- command to make -margins- yield SEs of .421637, matching what Keowmani got using -anova-?

Think of repeated-measures ANOVA as like an extended paired t-test. It has the subject effect subtracted out, leaving the conditionally independent residuals. If you want to think of the situation in terms of a variance-covariance matrix, you'd have the residual variance σ²_e on the diagonal and zeroes off-diagonal. It’s basically a scalar, the one (MS_error) that you see at the bottom of the ANOVA table. (You compute the subject variance σ²_u via expected mean squares.) The standard errors for margins after repeated-measures ANOVA uses this σ²_e, as illustrated in #6.

When you use mixed to fit the same model, the compound-symmetric variance-covariance matrix has on the diagonals not the residual variance σ²_e, but rather σ²_e plus the subject variance σ²_u, and has σ²_u off-diagonal as the covariance. (That "var(e)" that mixed prints out after you use the noconstant option to force the random effects equation to be empty isn't what you'd otherwise think it is—it's not the variance of the residuals.) The standard errors for margins uses this combined variance as shown below. To recover the residual variance, you have to subtract it out, also shown below.

To answer your question, I don’t know how to coax margins after mixed to use this fixed effects residual variance, other than to cynically fit a fixed effects model (shown last).

version 15.1

clear *

quietly input byte(y subject time)
20 1 1
24 1 2
28 1 3
28 1 4
15 2 1
18 2 2
23 2 3
24 2 4
18 3 1
19 3 2
24 3 3
23 3 4
26 4 1
26 4 2
30 4 3
30 4 4
22 5 1
24 5 2
28 5 3
26 5 4
19 6 1
21 6 2
27 6 3
25 6 4
end

// Repeated-measures ANOVA
anova y subject time, repeated(time)
margins time , df(`=e(df_r)')

// -mixed- with exchangeable residual correlation (compound symmetric variance-covariance matrix)
mixed y i.time || subject: , noconstant ///
    reml dfmethod(satterthwaite) residuals(exchangeable) ///
    nolrtest nolog
margins time , df(`=e(ddf_m)')

/* Recovering residual variance for time (subject variance removed) */
estat wcorrelation, covariance
tempname R n res_var
matrix define `R' = r(Cov)
matrix define `n' = e(N_g)
scalar define `res_var' = `R'[1, 1] - `R'[1, 2] // <= here is the crux

/* Compare below with Post #6 above in the thread */
display in smcl as text "Residual variance = " %09.7f `res_var'
display in smcl as text "Standard error based on residual variance alone = " %08.6f sqrt(`res_var' / `n'[1, 1])

/* Compare below with -margins- after -mixed- just above */
display in smcl as text "Standard error based on sum residual and subject variances, " ///
    "i.e., -margins- after -mixed- = " %08.6f sqrt(`R'[1, 1] / `n'[1, 1])

// How to use -mixed- to do repeated-measures ANOVA
mixed y i.time i.subject, reml dfmethod(residual) nolog
testparm i.subject, df(`=e(ddf_m)')
testparm i.time, df(`=e(ddf_m)')

margins time , df(`=e(ddf_m)')

exit

[Joseph Coveney]

It has the subject effect subtracted out

The residual actually has both terms factored out.

[Bruce Weaver]

Thanks Joseph. I can duplicate the output from your final -margins- command in SPSS using either of the following commands:

UNIANOVA Y BY time id
  /METHOD=SSTYPE(3)
  /EMMEANS=TABLES(time)
  /DESIGN=time id.

MIXED Y BY time id
  /FIXED=time id | SSTYPE(3)
  /METHOD=REMML
  /EMMEANS=TABLES(time).

Here are the EMMEANS results from those models (UNIANOVA first, then MIXED):

time                
Dependent Variable:   Y
time    Mean    SE    95% Confidence Interval    
            Lower    Upper
1    20.000    .421637    19.1013    20.8987
2    22.000    .421637    21.1013    22.8987
3    26.667    .421637    25.7680    27.5654
4    26.000    .421637    25.1013    26.8987

time(a)                    
time    Mean    SE    df    95% Confidence Interval    
                Lower     Upper
1    20.000    .421637    15    19.1013    20.8987
2    22.000    .421637    15    21.1013    22.8987
3    26.667    .421637    15    25.7680    27.5654
4    26.000    .421637    15    25.1013    26.8987
a Dependent Variable: Y.

And here is the output from Joseph's final -margins- command in #11:

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        time |
          1  |         20    .421637    47.43   0.000      19.1013     20.8987
          2  |         22    .421637    52.18   0.000      21.1013     22.8987
          3  |   26.66667    .421637    63.25   0.000     25.76797    27.56536
          4  |         26    .421637    61.66   0.000      25.1013     26.8987
------------------------------------------------------------------------------

[Bruce Weaver]

Hi Keowmani. I missed your -xtreg- example in #10 until now. Thanks.

[Bruce Weaver]

One consequence of using -margins- following repeated measures -anova- (as in #1) is that any pairwise contrasts of means will also use a pooled error term. Some authors (e.g., Maxwell, 1980) recommend using a separate error term for each pairwise contrast (i.e., an ordinary paired t-test), and applying a Bonferroni correction. (Nowadays, one might prefer to use Bonferroni-Holm, I suppose.) One of the reasons, I gather, is that the nature of the Treatment x Subjects interaction can vary quite a lot across all possible pairs of treatments.

Here are a couple links for Maxwell (1980).

• https://journals.sagepub.com/doi/10....69986005003269
• https://www.jstor.org/stable/1164969...n_tab_contents

See also this excerpt from Thom Baguley's book, Serious Stats: A Guide to Advanced to Advanced Statistics for the Behavioral Sciences (p. 642).

HTH.