Repeated measures ANOVA in Python is useful when you want to analyze data from within-subjects or repeated-measures designs. Although several Python packages now support this analysis, it was not always straightforward to perform using only Python.
At the time this tutorial was written, statsmodels did not include repeated measures ANOVA, and one common solution was to use rpy2 to call R from Python. That approach required a working R installation. An alternative was the pyvttbl package, which made it possible to run repeated measures ANOVA directly in Python using its anova() method.
Note: pyvttbl is no longer maintained, and you should see a new post using Pingouin for carrying out repeated Measures ANOVA: Repeated Measures ANOVA in R and Python using afex & pingouin
Table of Contents
Why Use a Within-Subjects Design?
Within-subjects designs have several advantages over between-subjects designs. One of the main benefits is that each participant is tested under every experimental condition rather than being assigned to only one condition. Because each participant contributes data to all conditions, fewer participants are typically needed to achieve the same statistical power.
Another advantage is that individual differences between participants are accounted for in the analysis. Since each participant serves as their own control, variability due to stable individual characteristics is reduced. This generally decreases the error variance and can increase the sensitivity of the statistical test.
Now, if you are working with data from a study with a between-subjects design and have only two groups you can carry out two-sample t-test with Python or Mann-Whitney U test in Python.
Repeated Measures ANOVA in Python
Throughout this tutorial, we will use the pyvttbl package for performing repeated-measures ANOVA in Python. Although the package is no longer maintained, the examples below remain useful for understanding how the analysis works. First, we need to install this Python package.
Installing pyvttbl
pyvttbl can be installed using pip:
pip install pyvttblCode language: Bash (bash)If you are using Linux, you may need to add ‘sudo’ before the pip command. This method installs pyvttbl and, hopefully, any missing dependencies. Note: if you decide to work with pyvttbl, you also need to use pip to install specific versions of the dependencies.
Creating an Example Dataset
We will begin by simulating a response time dataset for the examples in this tutorial. If you already have your own data, you can skip this step and load your dataset from a CSV file using the read_tbl() method.
from numpy.random import normal
import pyvttbl as pt
from collections import namedtuple
N = 40
P = ["noise","quiet"]
rts = [998,511]
mus = rts*N
Sub = namedtuple('Sub', ['Sub_id', 'rt','condition'])
df = pt.DataFrame()
for subid in xrange(0,N):
for i,condition in enumerate(P):
df.insert(Sub(subid+1,
normal(mus[i], scale=112., size=1)[0],
condition)._asdict())Code language: Python (python)Cunning a repeated measures ANOVA with pyvttbl only requires a few lines of code. Call the anova() method on a pyvttbl DataFrame, specify the dependent variable and the column containing the participant IDs (sub), and then provide the within-subject factor(s) using the wfactors argument. The simulated dataset used in this tutorial contains a single within-subject factor (condition).
Note: If you are using NumPy 1.10 or later, pyvttbl may not work correctly. The easiest solution is to run pyvttbl in a virtual environment with an older NumPy version (see Step-by-step guide for solving the Pyvttbl Float and NoneType error for a detailed solution both for Linux and Windows users).
aov = df.anova('rt', sub='Sub_id', wfactors=['condition'])
print(aov)Code language: Python (python)Tests of Within-Subjects Effects
Measure: rt
| Source | Type III Sum of Squares | ε | df | MS | F | Sig. | η2G | Obs. | SE of x̄ | ±95% CI | λ | Obs. Power | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| condition | Sphericity Assumed | 4209536.428 | – | 1.000 | 4209536.428 | 309.093 | 0.000 | 4.165 | 40.000 | 19.042 | 37.323 | 317.019 | 1.000 |
| Greenhouse-Geisser | 4209536.428 | 1.000 | 1.000 | 4209536.428 | 309.093 | 0.000 | 4.165 | 40.000 | 19.042 | 37.323 | 317.019 | 1.000 | |
| Huynh-Feldt | 4209536.428 | 1.000 | 1.000 | 4209536.428 | 309.093 | 0.000 | 4.165 | 40.000 | 19.042 | 37.323 | 317.019 | 1.000 | |
| Box | 4209536.428 | 1.000 | 1.000 | 4209536.428 | 309.093 | 0.000 | 4.165 | 40.000 | 19.042 | 37.323 | 317.019 | 1.000 | |
| Error(condition) | Sphericity Assumed | 531140.646 | – | 39.000 | 13618.991 | ||||||||
| Greenhouse-Geisser | 531140.646 | 1.000 | 39.000 | 13618.991 | |||||||||
| Huynh-Feldt | 531140.646 | 1.000 | 39.000 | 13618.991 | |||||||||
| Box | 531140.646 | 1.000 | 39.000 | 13618.991 |
The output includes Type III sums of squares, which are also used by statistical software such as SPSS and R packages including afex and ez. The table also reports the Greenhouse-Geisser, Huynh-Feldt, and Box corrections for violations of the sphericity assumption. In this example, the within-subject factor has only two levels, so sphericity is not an issue.
The output also includes generalized eta squared as an effect size measure together with 95% confidence intervals. According to the Anova class documentation, the standard errors and confidence intervals are calculated following Loftus and Masson (1994). Generalized eta squared is useful because it can be compared across both within-subjects and between-subjects designs (Olejnik & Algina, 2003).
If your data require a transformation, pyvttbl provides a transform argument. Supported transformations include log, log10, reciprocal, inverse, sqrt, arcsin, and windsor10. For example, a log transformation can be applied by setting transform="log".
aovlog = df.anova('rt', sub='Sub_id', wfactors=['condition'], transform='log')Code language: Python (python)pyvttbl can also be used for mixed-design (split-plot) ANOVAs. To include between-subject factors, simply pass them to the bfactors argument as a list. If you are looking for a one-way ANOVA for independent groups, see my newer tutorial: Four Ways to Conduct One-Way ANOVAs in Python.
The anova() method returns a dictionary-like object containing the ANOVA results. The helper function below extracts the statistics needed for APA reporting, including the F statistic, mean square error (MSE), generalized eta squared, and the p value. To see the available result keys for your analysis, use aov.keys().
def extract_for_apa(factor, aov, values = ['F', 'mse', 'eta', 'p']):
results = {}
for key,result in aov[(factor,)].iteritems():
if key in values:
results[key] = result
return resultsCode language: Python (python)The ANOVA table shown above was generated using the private _within_html method. Before creating the HTML output, import SimpleHTML:
import SimpleHTML
output = SimpleHTML.SimpleHTML('Title of your HTML-table')
aov._within_html(output)
output.write('results_aov.html')Code language: Python (python)Although pyvttbl makes it possible to run repeated measures ANOVA in Python, it has a few limitations. One is that it uses its own DataFrame implementation rather than working directly with Pandas DataFrames. In practice, this is not a major issue, since pyvttbl provides the functionality needed to run the analysis.
The main limitation today is that pyvttbl is no longer actively maintained. For current projects, I recommend using more recent packages such as Pingouin, which I cover in a newer tutorial on repeated measures ANOVA in Python.
References
Loftus, G.R., & Masson, M.E. (1994). Using confidence intervals in within-subjects designs. The Psychonomic Bulletin & Review, 1(4), 476-490.
Olejnik, S., & Algina, J. (2003). Generalized eta and omega squared statistics: measures of effect size for some common research designs. Psychological Methods, 8(4), 434–47. http://doi.org/10.1037/1082-989X.8.4.434
Dear Erik,
thanks for this excellent blog post! I would like to point out that there is also relatively new, but actively maintained RM-ANOVA support in the mne-python package, see http://martinos.org/mne/stable/generated/mne.stats.f_mway_rm.html#mne.stats.f_mway_rm The data have to be transformed into a certain shape to use this function though; but this should be straightforward using pandas. It would be great if you could cover this approach in future post, as I think it might help many people making the final switch to Python.
All the best,
–Richard
Hey Richard,
Glad you liked my post.
Thank you for mentioning mne-python package. It seems like I’ve missed it.
I will have a look at it later today (travelling atm).
Again, thank you for pointing it out.
Have a nice day,
Erik
Hi Erik,
great you like my suggestion 🙂 I actually discovered mne’s RM-ANOVA through this github issue re statsmodels: https://github.com/statsmodels/statsmodels/issues/749
And by the way, although the docstring of mne-python’s f_mway_rm() states it expects a 3D array as input, the code will handle ordinary 2D (“behavioral”) data just fine: https://github.com/mne-tools/mne-python/blob/master/mne/stats/parametric.py#L295
It would be great to have a side-by-side comparison of the results from this RM-ANOVA implementation and possibly afex or JASP (which, to my knowledge, internally uses afex). Let me know if I can help you with this in any way!
(NB: I’m by no means a stats guru, but simply a neuroscientist/psychologist who happens to be willing to analyze his data 😉 That’s why I’ve been very hesitant to move away from afex so far, because I have the impression that Henrik Singmann, the author of afex — who also happens to be a psychologist! — made very “sane” decision as to how the package should behave, e.g. it allows me to replicate existing SPSS analyses etc., so I can assume that “I’m doing it right” with great confidence.)
Thanks again,
–Richard
Oh and, of course, have a nice day too! 🙂
Great suggestions. I have actually been thinking on comparing R-packages and Python packages. Maybe I’ll reach out to you.
I am no stats guru either… However, I enjoy programming (may sound strange but I find it relaxing to write code) and I’ve found that I understand statistical concepts better when using R (or Python).
I tend to mess things up in GUI interfaces such as SPSS (e.g., forgetting which boxes I checked in last time). Doing stats in R (or Python) make it easier to reproduce the same analysis. In my two first studies of my thesis I use R (i.e., afex and lsmeans). Afex is really great.
Hi Erik,
Do you know where I could find the SimpleHTML module?
Thanks,
John
Hey John,
just import it from pyvttbl:
from pyvttbl import SimpleHTML.Erik
Stellar, thanks for posting this.
Hi,
Thanks for your post.
I’m able to run your example without an issue, but when I try to run it with a second within subject factor, I get the error below;
Exception: (‘data has %d conditions; design only %d’, 2, 4)
I’ve pasted is the script that causes the error below. As you can see, the only change that I have made to your example is to add a second within subject factor (which is simply a duplicate of your single within subject factor), such that now the data contains condition1 and condition2 instead of only condition. Have you been able to run your example with multiple within subject factors?
Here is the script;
from numpy.random import normal
import pyvttbl as pt
from collections import namedtuple
N = 40
P = [‘c1’, ‘c2’]
rts = [998, 511]
mus = rts * N
Sub = namedtuple(‘Sub’, [‘Sub_id’, ‘rt’, ‘condition1’, ‘condition2’])
df = pt.DataFrame()
for subid in xrange(0, N):
for i, condition in enumerate(P):
df.insert(Sub(subid + 1, normal(mus[i], scale=112., size=1)[0], condition, condition)._asdict())
aov = df.anova(‘rt’, sub=’Sub_id’, wfactors=[‘condition1’, ‘condition2’])
print(aov)
Thanks, Phil
Hey Phil,
For simulating data for a 2-way ANOVA for repeated measures look at this example: http://wp.me/p5bFFs-oV
Regards,
Erik
Hello Erik,
thank you for your article!
I am trying to use your guide for running a repeated measures analyisis for my experiment. I have the data in a CSV file which I imported using the method read_table. However, I am getting issues when running the anova method.
Here my script:
df = pt.DataFrame()
df.read_tbl(‘GEQ_DATA_ANALYSIS_CORE_CSV_SUMMARY_p.csv’,delimiter =’;’)
aov = df.anova(‘FLOW’, sub=’P_ID’, wfactors=[‘SETUP’])
print(aov)
The error I get is:
TypeError Traceback (most recent call last)
in ()
1 df = pt.DataFrame()
2 df.read_tbl(‘GEQ_DATA_ANALYSIS_CORE_CSV_SUMMARY_p.csv’,delimiter =’;’)
—-> 3 aov = df.anova(‘FLOW’, sub=’P_ID’, wfactors=[‘SETUP’])
4 print(aov)
C:\Program Files\Anaconda2\lib\site-packages\pyvttbl\base.pyc in anova(self, dv, sub, wfactors, bfactors, measure, transform, alpha)
1973 aov=stats.Anova()
1974 aov.run(self, dv, sub=sub, wfactors=wfactors, bfactors=bfactors,
-> 1975 measure=measure, transform=transform, alpha=alpha)
1976 return aov
1977
C:\Program Files\Anaconda2\lib\site-packages\pyvttbl\stats\_anova.pyc in run(self, dataframe, dv, wfactors, bfactors, sub, measure, transform, alpha)
708
709 if len(wfactors)!=0 and len(bfactors)==0:
–> 710 self._within()
711
712 if len(wfactors)==0 and len(bfactors)!=0:
C:\Program Files\Anaconda2\lib\site-packages\pyvttbl\stats\_anova.pyc in _within(self)
1131 for e in xrange(1,Ne+1):
1132 # code effects so we can build contrasts
-> 1133 cw = self._num2binvec(e,Nf)
1134 efs = asarray(factors)[Nf-1-where(asarray(cw)==2.)[0][::-1]]
1135 r=self[tuple(efs)] # unpack dictionary
C:\Program Files\Anaconda2\lib\site-packages\pyvttbl\stats\_anova.pyc in _num2binvec(self, d, p)
1238 d=floor(d/2.)
1239
-> 1240 return list(array(list(zeros((p-len(b))))+b)+1.)
1241
1242 ## def output2html(self, fname, script=”):
TypeError: ‘float’ object cannot be interpreted as an index
Which are the requirements for the CSV file ? How does it has to be formatted ?
Thank you in advance!
Best,
Luca
Hey Luca,
your data should be stored in long format (See here, for instance.). I am not sure exactly what is going on with your ANOVA but if you could post your the first couple of rows of your data (‘df’, in your code) maybe we can find a solution.
Kind regards,
Erik
Hi – I am getting the same error. It don’t think it is about the data, but instead the packages installed in anaconda. My pyvttbl code will run perfectly until something is upgraded, then I get the identical error:
-> 1240 return list(array(list(zeros((p-len(b))))+b)+1.)
TypeError: ‘float’ object cannot be interpreted as an index
I don’t know how to tell which package is causing the crash. It was previously fixed again *magically* after upgrading some packages. Now I have the same problem again.
If anyone knows which package might be the culprit, or a work around for the error it would be very helpful. This is a great package for ANOVAs and I hope I can keep using it.
Hey Oly,
You are right. After a comment by Damien we know what package is the culprit. It is Numpy. Pyvttbl is not maintained any more and seems only be compatible with Numpy version 1.1.x. I have created a step-by-step guide on how to run Pyvttbl within a virtual environment (https://www.marsja.se/solving-pyvttbl-error-float-nonetype-error/). I may see if I can come up with another solution if I have time. You also have the option to install the beta version of statsmodels. You can run repeated measures ANOVA with this version of statsmodels. However, you will only get an ANOVA table containing the degrees of freedeom, F-values, and p-values.
Great, thank you for your fast response! I will take a look at this information and post my solution in the future for others who may have the same problem.
No worries. It may, of course, not be the same issue.
Hi Erik,
To follow up, I implemented the AnovaRM function in the beta version of statsmodels. It worked great! Results were identical to SPSS output (I’m running two-way, 2×2 levels, repeated measures ANOVAs).
Found your YouTube video very helpful as well: https://www.youtube.com/watch?v=xzET1rpvJ_A
Thank you so much for taking the time to do this.
Best,
Oly
Hey Oly,
Thanks for your comment. That is great! I put together that video in a rush, glad it helped. I did not have time to figure out to get/calculate effect sizes. If, or when, I figure it out, I may put it together in a blog post and/or a YouTube video.
Best,
Erik
That would definitely be great! Should be relatively straight forward if you have the sum of squares, correct? I’m happy to test code on my own data, just let me know.
Looking forward to following your posts. Thanks again.
Best,
Oly
Hi Erik,
I just wanted to let you know that I did find a discrepancy in the AnovaRM function in the beta version of statsmodels when I ran a 1-way RM ANOVA with 2 levels (so akin to a paired t-test). I noticed the degrees of freedom of the error was off, and indeed the function gave different values than in SPSS. So in sum, the two-way RM ANOVAs I tested give identical output as SPSS, but the one-way RM ANOVAs do not. Let me know if you would like me to send the data I used?
Thanks in advance for all your help. The AnovaRM function implementation is great, so I hope this can be sorted out relatively easily.
Best,
Oly
Hey Oly,
Thanks for letting me know. I am not involved in developing the AnovaRM method. Maybe you could open up an issue on their GitHub page (https://github.com/statsmodels/statsmodels). This way it may be fixed.
Erik
Hi Erik,
Ok thanks! Sorry, I thought you were one of the contributors.
Best,
Olympia
Hi Erik! Thank you for your article!
I am a beginner in Python, I’m trying to use your guide for running a split plot anova (my goal is to determine the interaction between two within variables( AGE’, ‘ETHNICGROUP’) and one between variable( ‘SEXP’). and I obtained the output (thank you again).
Now I’m trying to use your function to extract the p-value obtained with the anova method, but I don’t understand how it works.
These are my aov.keys():
[(‘AGE’,),
(‘ETHNICGROUP’,),
(‘AGE’, ‘ETHNICGROUP’),
(‘SEXP’,),
(‘AGE’, ‘SEXP’),
(‘ETHNICGROUP’, ‘SEXP’),
(‘AGE’, ‘ETHNICGROUP’, ‘SEXP’),
(‘SUBJECT’,),
(‘TOTAL’,),
(‘WITHIN’,),
(‘AGE’, ‘SUBJECT’),
(‘ETHNICGROUP’, ‘SUBJECT’),
(‘AGE’, ‘ETHNICGROUP’, ‘SUBJECT’)]
How do I modify your script? ==> def extract_for_apa(factor, aov, values = [‘F’, ‘mse’, ‘eta’, ‘p’]):
results = {}
for key,result in aov[(factor,)].iteritems():
if key in values:
results[key] = result
return results
Hey Leonardo,
I am glad you found it helpful. I realize that I don’t show an example of how to use that function. It can extract the values for each factor:
print(extract_for_apa('condition', aov))Would print F, MSE, eta-squared, and the p-value (based on the example above). In your case you could loop through the 3 variables of interest or just runt the function 3 times with your variables as input. Let me know if it doesn’t work or if you need help.
/Erik
Thank you, it works. These are the results:
print(extract_for_apa(‘AGE’, aov))
{‘p’: 0.38304275607253779, ‘mse’: 0.00071049005495488833, ‘eta’: 0.0047898426784019855, ‘F’: 0.77986179350784657}
print(extract_for_apa(‘ETHNICGROUP’, aov))
{‘p’: 0.019574151602473413, ‘mse’: 0.00069223374419013475, ‘eta’: 0.034659508263642044, ‘F’: 5.9711544610483269}
print(extract_for_apa(‘SEXP’, aov))
‘p’: 0.44921555871908303, ‘mse’: 0.00091254877064392932, ‘eta’: 0.0046182855882132079, ‘F’: 0.58533465964381803}
Due to my lack of knowledge in data analysis, I thought that if you want to explore the interaction between two wfactors and one bfactor you have to obtain just one p value (and not one for every variable.)…but as you said I was wrong!
Sorry if I bother you again Erik, but I was wondering if there’s a way to print the p value of the interaction between these three factors. (Unfortunately you have to obtain just one p value to see if ‘AGE’ * ‘ETHNICGROUP’ * ‘SEXP’ is significant.
I tried to put in your function the “key of the interaction” (print(extract_for_apa((‘AGE’ , ‘ETHNICGROUP’ , ‘SEXP’), aov)) it doesn’t work.