Two way repeated measures ANOVA (experimental data)

[Fabian Hernandez]

Hello everyone,

I hope you are doing fine. I am currently assessing an experiment. I had 90 subjects performing a task over 5 rounds. They were divided into three groups (between subjects), and next to rounds and scores I also had their gender (here coded Female). I want to look at the interaction between female and round, to see whether there are differences in round effects depending on gender. E.g. maybe there is a round (time) effect for male subjects, but not for female.
I found this code and adapted it to my needs, however, I am unsure about its interpretation.

Below you can find the code I used (I hope it is in the right format). You can also see the comments on the effects I made.
I have the following understanding of what I am doing:
- First, the anova is performed, which indicates significant interaction. Contrast shows this for the specific round.
- Through the contrast command afterwards I see the effect of time for each gender and if it is significant
- the margins in the next one is checking separatededly for both genders and all pairs
- the next part checks for the singular effect of the gender at each time

My questions are:
- Are my comments correct?
- Is my understanding correct?
- Am I missing anything?

anova Score Female / ParticipantID|Female Round Female#Round, repeated(Round)
** -> interaction effect is  significant (gender and round interact with each other)
** -> Round effect is significant (some tests) and indicates that there is a difference between the rounds

contrast a.Round#Female
** -> Jointly, only effect from second to third round is significant

****** additional aspects
contrast Round@Female, effect
** -> effects only for male: 3 to base, 4 to base, and 5 to base
margins Round, at(Female=1) pwcompare(effects) noestimcheck
** -> no effects
margins Round, at(Female=0) pwcompare(effects) noestimcheck
** -> effects: 3 to 1, 4 to 1, 5 to 1, 4 to 2, and 5 to 2,
anova Score Female##Round
contrast Round@Female, effect
** -> only significance in comparing 1 to 5
margins Female#Round
marginsplot, x(Round)

* Example generated by -dataex-. To install: ssc install dataex
clear
input byte(ParticipantID Group Round) double Score byte Female
 1 1 1   -.7910944306858633 1
 1 1 2  -1.0269490744398917 1
 1 1 3   -.8738091030329711 1
 1 1 4   -.5729675484954635 1
 1 1 5    -.400363384534748 1
 2 1 1    .5927726206979186 0
 2 1 2    .8297738353289045 0
 2 1 3    .5762526177174433 0
 2 1 4    .6744897501960817 0
 2 1 5    .2967378382598981 0
 3 1 1   -.8495855373182969 1
 3 1 2   -.8376590430971385 1
 3 1 3   -.7910944306858633 1
 3 1 4   -.6195477803879315 1
 3 1 5   -.8180417759124541 1
 4 1 1   -.4614936942181583 1
 4 1 2  -1.2627844162947808 1
 4 1 3  -1.6448536269514726 1
 4 1 4   -1.763728040633686 1
 4 1 5   -.8738091030329711 1
 5 1 1    .7242880874234008 0
 5 1 2    .2967378382598981 0
 5 1 3   .48329133890772974 0
 5 1 4    .8985565917048498 0
 5 1 5    .9944578832097528 0
 6 1 1  -1.8807936081512509 1
 6 1 2  -1.8191477690437625 1
 6 1 3  -1.8807936081512509 1
 6 1 4  -2.0098747721953845 1
 6 1 5  -2.0098747721953845 1
 7 1 1  -.18486702761944254 1
 7 1 2   -.4246235657803619 1
 7 1 3  -.32891197566572217 1
 7 1 4  -.08924287887565444 1
 7 1 5 .0027851461281113033 1
 8 1 1   1.4757910281791709 0
 8 1 2   1.3271942183604362 0
 8 1 3   1.2324587197131789 0
 8 1 4  -.49900228007558534 0
 8 1 5   -.6433454053929168 0
 9 1 1  -1.1476547672158675 0
 9 1 2  -1.3407550336902165 0
 9 1 3   -.9990349151802846 0
 9 1 4  -1.4757910281791702 0
 9 1 5  -1.3829941271006387 0
10 1 1   1.0364333894937894 1
10 1 2   .33479759110257196 1
10 1 3   .08365173390712907 1
10 1 4  -.06968492031845566 1
10 1 5   .22183228666917584 1
11 1 1   -.5927726206979184 1
11 1 2   -.9853657550352516 1
11 1 3    .5084880591093566 1
11 1 4     .735192571938582 1
11 1 5    .1368983904280161 1
12 1 1  -.22754497664114945 1
12 1 2  -.06968492031845566 1
12 1 3  .061311526679252146 1
12 1 4  -.32891197566572217 1
12 1 5    .5084880591093566 1
13 1 1    .7098812513288475 0
13 1 2     .612812991016627 0
13 1 3   .11444014881194711 0
13 1 4    .8025718880593095 0
13 1 5  -.18486702761944254 0
14 1 1  -1.4925479547163985 1
14 1 2   1.5184980572589322 1
14 1 3   1.2627844162947803 1
14 1 4   1.1263911290388005 1
14 1 5   .35845879325119373 1
15 1 1    .5469043537825481 0
15 1 2    .5469043537825481 0
15 1 3    .3853204664075676 0
15 1 4    .8297738353289045 0
15 1 5    .6263107889358445 0
16 1 1    .6920771366140133 1
16 1 2    .2418596857990917 1
16 1 3    -.400363384534748 1
16 1 4    .1368983904280161 1
16 1 5  -.08924287887565444 1
17 1 1   -1.287908610584834 0
17 1 2  -1.0460118625082266 0
17 1 3  -.18486702761944254 0
17 1 4  -.06968492031845566 0
17 1 5  -.08924287887565444 0
18 1 1  -1.0954184989297098 0
18 1 2   -1.238432987203541 0
18 1 3  -1.1694503792916737 0
18 1 4   -.9718782705519045 0
18 1 5   -.9853657550352516 0
19 1 1  -.34956302147213447 1
19 1 2   -.4246235657803619 1
19 1 3   -.5308026326187639 1
19 1 4  -.31717461426306087 1
19 1 5  -.22754497664114945 1
20 1 1  -.49900228007558534 0
20 1 2  -.31717461426306087 0
20 1 3  -.37932781842312313 0
20 1 4   -.5308026326187639 0
20 1 5   .11444014881194711 0
end

[Joseph Coveney]

Wow, so much precision. Anyway, I cannot answer your first two questions, but here's my stab at your third: below is how I'd approach an exploratory analysis, assuming that the experiment's written protocol didn't spell out a formal confirmatory-style analysis in its Statistical Considerations section. (I've attached it as a do-file, as well, if you care to execute it against all 90 participants' worth of data.)

version 17.0

clear *

quietly input byte(ParticipantID Group Round) double Score byte Female
/* Elided here, but they're in the do-file attached */
end

// Cleanup
rename *, lower
drop group
rename participantid pid
rename round tim
rename score sco
rename female sex

*
* Primary exploration
*
mixed sco i.sex##i.tim || pid: , ///
    noconstant residuals(unstructured, t(tim)) ///
        reml dfmethod(kroger) ///
            nolrtest nolog
estimates store Manova

quietly mixed sco i.sex##i.tim || pid: , ///
    noconstant residuals(exchangeable) ///
        reml dfmethod(kroger) ///
            nolrtest nolog
local line_size `c(linesize)'
set linesize 80
lrtest Manova
set linesize `line_size'

/* You can -ssc describe mauchly-, too, if you want.
   Anyway, MANOVA seems to fit better, and so we'll go with that. */

quietly estimates restore Manova

// Diagnostics
set more on
predict double res, residuals
pnorm res , ///
    mcolor(black) msize(vsmall) ///
        rlopts(lcolor(black) lpattern(dash)) ///
            ylabel( , angle(horizontal) nogrid)
more
qnorm res , ///
    mcolor(black) msize(vsmall) ///
        rlopts(lcolor(black) lpattern(dash)) ///
            ylabel( , angle(horizontal) nogrid)
more

predict double xb, xb
graph twoway scatter res xb, ///
    mcolor(black) msize(vsmall) ///
        yline(0, lcolor(black) lpattern(dash)) ///
            ylabel( , angle(horizontal) nogrid)
more

*
* Primary analysis
*
quietly margins sex#tim
marginsplot , ///
    xdimension(tim) plotdimension(sex) ///
        level(50) ///
            plot1opts(mcolor(blue) msymbol(O) mfcolor(white)) ///
            plot2opts(mcolor(red) msymbol(O) mfcolor(white)) ///
                title("") ylabel( , angle(horizontal) nogrid) legend(off)


exit

[Fabian Hernandez]

Thanks a lot Joseph!

Could you elaborate on what you are doing? I have never seen manova with margins. What exactly is the outcome we are getting?

And how would you describe the outcome? (see picture below)
Is it ok to say that while there does not seem to be a significant effect of time for females (red), there seems to be a significant effect of time for males?
(as the confidence intervals are still incorporating previos rounds for red, but not for blue?
-> As it is the effect of time/rounds, I would say there is learning (or similar) for men but not for women?

Thank you very much in advance, I appreciate your help a lot!

[Joseph Coveney]

Is it ok to say that while there does not seem to be a significant effect of time for females (red), there seems to be a significant effect of time for males?
. . . I would say there is learning (or similar) for men but not for women?

In the absence of any context, if those are 50% CIs and all 90 participants, then I'd say that the profile plot doesn't suggest much that's earth-shattering going on for either sex. But you'd also want to look at the pattern of temporal (i.e., within-participant) changes and their variability, pattern of dropout, the measurement qualities of the score that you're using, and what makes biological sense in the context of the phenomenon that you're trying to measure.

Also helpful in interpretation: (1) What did the diagnostic plots look like? (2) What was the relative magnitude of the estimates of variances and covariances, and how different were they between time points?

[Fabian Hernandez]

Thank you Joseph!

To give a little bit of context. As I said earlier, participants perform a task over five rounds. They compete in a tournament in each round against two other players (anonymously). I also checked for a simple round effect, and I found nothing.
There was no dropout, the measurement is a simple task (so not of concern), and the biological sense is not given because it is an economic experiment. Theory suggests differences between gender, but there is almost nothing on gender effects over time. Does that help in interpreting? I would agree the effects are not strong, however, there is something. (and yes it is all 90 participants and 50% CIs)

What is a diagnostic plot? And how can I look at the relative magnitude?

[Joseph Coveney]

To give a little bit of context. As I said earlier, participants perform a task over five rounds. They compete in a tournament in each round against two other players (anonymously).

Again, I don't have the subject matter knowledge to be able to address your bigger questions of interpretation, but I can try below to mention a few observations that struck me, if that will help.

There was no dropout

Good, but consider that there may be more than one way to drop out of a tournament (behavioral study) & still make mere presence appear as participation.

the measurement is a simple task (so not of concern)

That a task is simple would not obviate my concern over its measurement quality—the usual psychometric considerations, signal-to-noise ratio, whether I'm measuring what I intend to (what I think I am) etc.

Theory suggests differences between gender, but there is almost nothing on gender effects over time. Does that help in interpreting?

Yes, it does, to me. I would try to let the absence of theoretical justification (prediction) be cautionary, to temper my interpretative enthusiasm. (By the way, do you observe a statistically significant main effect of sex that theory predicts? If not, then I would take that in itself as cautionary.)

What is a diagnostic plot?

The diagnostic plots that I referred to are those identified by a comment in the code. (See #2 above.)

And how can I look at the relative magnitude?

You can look at the relative magnitude of the variances and temporal pattern of error covariances by inspection of their parameter estimates that are printed out beneath the table of those for fixed effects in the fitted the regression model's output.

[Fabian Hernandez]

Hello Joseph!
Thanks a lot, I got bugged down with some other stuff, so it took me some time to get back.

[Fatimah Gh]

I wanted to do this for my data but there seems to be a lot of errors.
my data has 4 variables:
record id
pain score
day
primary anesthesia

first I have two group in the primary anesthesia
and for the pain score it is repeated each day (Day) variable

* Example generated by -dataex-. For more info, type help dataex
clear
input long primaryanesthesia byte(painscore Day recordid)
1 8 0 28
1 5 0 24
1 5 0  5
1 9 0 10
1 4 0 45
1 5 0 30
1 5 0 29
1 4 0 21
1 6 0 17
1 6 0 27
1 5 0 42
1 6 0 35
1 7 0 34
1 9 0 13
1 0 0 19
1 0 0 26
1 7 0 36
1 6 0 39
1 5 0 16
1 7 0 23
1 8 0 12
1 2 0  6
2 0 0 25
2 4 0 33
2 4 0 31
2 0 0  9
2 0 0 40
2 4 0 37
2 0 0  7
2 0 0 18
2 2 0 44
2 0 0 43
2 0 0 22
2 0 0  4
2 0 0  8
2 0 0 38
2 0 0 32
2 . 0 11
2 0 0 41
. . 0  1
2 0 0 15
2 0 0 14
. 0 0 20
2 . 0  2
. . 0  3
1 5 1 28
1 4 1 24
1 3 1  5
1 5 1 10
1 2 1 45
1 4 1 30
1 3 1 29
1 3 1 21
1 4 1 17
1 5 1 27
1 5 1 42
1 5 1 35
1 6 1 34
1 9 1 13
1 6 1 19
1 4 1 26
1 5 1 36
1 4 1 39
1 3 1 16
1 5 1 23
1 6 1 12
1 2 1  6
2 4 1 25
2 3 1 33
2 0 1 31
2 2 1  9
2 3 1 40
2 3 1 37
2 1 1  7
2 4 1 18
2 0 1 44
2 3 1 43
2 1 1 22
2 1 1  4
2 1 1  8
2 3 1 38
2 3 1 32
2 . 1 11
2 2 1 41
. . 1  1
2 3 1 15
2 2 1 14
. 6 1 20
2 . 1  2
. . 1  3
1 3 2 28
1 5 2 24
1 3 2  5
1 3 2 10
1 2 2 45
1 3 2 30
1 1 2 29
1 2 2 21
1 2 2 17
1 3 2 27
end
label values primaryanesthesia poat
label def poat 1 "General Anesthsia", modify
label def poat 2 "Spinal Anesthsia", modify

[Joseph Coveney]

I wanted to do this for my data but there seems to be a lot of errors.

You don't show any code for your attempt; what kinds of error are you encountering?

Anyway, the following lines illustrate one approach and maybe you can try something along them. They would begin right after the dataset listing that you give.

I'm not sure what the purpose of your study is, but the code below includes a couple of optional model-fitting exercises (they're equivalent in this case) that yield test statistics. The main theme, though, is exploration and display. (Complete do-file and log file are attached.)

// Preliminary data management
_strip_labels primaryanesthesia
rename (primaryanesthesia painscore Day recordid) (trt sco tim rid)

* Optional for this purpose:
label variable trt "Primary Anesthsia (1 General 2 Spinal)"
label variable sco "Pain Score (Probably  0–10)"
label variable tim "Postoperative Day"
label variable rid "Study Participant (Record ID)"

order rid trt tim sco
sort rid tim

// Preliminary examination
list, noobs sepby(rid) // Note Partiicpant 20

tabulate sco

version 16.1: table trt tim, contents(mean sco sd sco min sco max sco ///
    n sco) format(%03.1f)

version 16.1: table trt tim, contents(median sco iqr sco) format(%03.1f)


// Because data are balanced through first postoperative day, this is doable
anova sco trt/rid|trt tim trt#tim if tim < 2

// Alternatively and with the same results
mixed sco i.trt##i.tim if tim < 2 || rid: , reml dfmethod(kroger) nolrtest nolog
contrast trt tim trt#tim, small

// For display
quietly mixed sco i.trt##i.tim || rid: , reml dfmethod(kroger) nolrtest
quietly predict double xb, xb
quietly predict double res, res

quietly margins trt#tim, df(`e(df_min)')
quietly marginsplot , ///
    xdimension(tim) ///
    plotopts(lcolor(black) mcolor(black) mfcolor(white)) ///
    plot1opts(msymbol(O)) plot2opts(msymbol(T) msize(medlarge)) ///
    ciopts(lcolor(black)) level(67) ///
    ytitle(Pain Score) ylabel(0(2)6, angle(horizontal) nogrid) ///
    xtitle(Day) ///
    title("") scheme(s2color) legend(off)
quietly graph export anesthesia.png

// seedem
graph twoway ///
    scatter res xb, ///
    mcolor(black) msize(vsmall) jitter(5) jitterseed(374789367) ///
    yline(0, lcolor(black) lpattern(dash)) ///
    scheme(s2color) ylabel( , angle(horizontal) nogrid)
quietly graph export rvf.png

qnorm res, ///
    mcolor(black) msize(vsmall) ///
    rlopts(lcolor(black)) ///
    scheme(s2color) ylabel( , angle(horizontal) nogrid)
quietly graph export qq.png

The missing treatment group assignment for Participant 20 is a little curious. You might want to go over the raw-data records one more time.

[Fatimah Gh]

i didnt notice that dataex only listed 100obs out of 180obs
can i run the do file for my data or i need to do some changes?

[Joseph Coveney]

i didnt notice that dataex only listed 100obs out of 180obs
can i run the do file for my data or i need to do some changes?

Just attach the dataset to your post and we can take a look.

[Fatimah Gh]

here is the file

[Joseph Coveney]

here is the file

Yes, you can use run the do-file, substituting the new data. The if tim < 2 condition in the ANOVA and mixed-model specifications is no longer necessary and you can remove that if condition from them.

You might want to use the fitted prediction (predict double fit, fitted) for the residuals-versus-fitted plot instead of the linear prediction for the fixed portion of the model only (predict double xb, xb) that I show above. The coarseness of the pain measurement makes that graph a little ugly, and the floor effect that's present in the data gives rise to a couple of apparent outliers and to a little sway in the distribution of the residuals that's evident in their Q-Q plot.