Sign up ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

The following question is one of those holy grails for me for some time now, I hope someone might be able to offer a good advice.

I wish to perform a non-parametric repeated measures multiway anova using R.

I have been doing some online searching and reading for some time, and so far was able to find solutions for only some of the cases: friedman test for one way nonparametric repeated measures anova, ordinal regression with {car} Anova function for multi way nonparametric anova, and so on. The partial solutions is NOT what I am looking for in this question thread. I have summarized my findings so far in a post I published some time ago (titled: Repeated measures ANOVA with R (functions and tutorials), in case it would help anyone)


If what I read online is true, this task might be achieved using a mixed Ordinal Regression model (a.k.a: Proportional Odds Model).

I found two packages that seems relevant, but couldn't find any vignette on the subject:

So being new to the subject matter, I was hoping for some directions from people here.

Are there any tutorials/suggested-reading on the subject? Even better, can someone suggest a simple example code for how to run and analyse this in R (e.g: "non-parametric repeated measures multiway anova") ?

Thanks for any help, Tal

share|improve this question
Tal, may I ask whether you found a solution for this? I am having the same problem, and the replies below might be helpful in assisting to find an answer, but don't really provide a definitive answer. I have 9 ordinal DVs and 2 time points, and am looking for the same test you wanted to perform. – Torvon Jul 5 '13 at 23:36
Hi Torvon. I have never came by a solution. I guess permutation tests will be the safest route, but I never got the time to sit and make it work. If you do - please come again to post your answer. Cheers, T – Tal Galili Jul 6 '13 at 18:53
Thank you for your quick response. I'll have to work this out and will let you know. – Torvon Jul 6 '13 at 23:46

4 Answers 4

up vote 8 down vote accepted

The ez package, of which I am the author, has a function called ezPerm() which computes a permutation test, but probably doesn't do interactions properly (the documentation admits as much). The latest version has a function called ezBoot(), which lets you do bootstrap resampling that takes into account repeated measures (by resampling subjects, then within subjects), either using traditional cell means as the prediction statistic or using mixed effects modelling to make predictions for each cell in the design. I'm still not sure how "non-parametric" the bootstrap CIs from mixed effects model predictions are; my intuition is that they might reasonably be considered non-parametric, but my confidence in this area is low given that I'm still learning about mixed effects models.

share|improve this answer
Hello Mike. Thank you for the answer, and for your package - it is really great! – Tal Galili Sep 3 '10 at 11:05
@Mike, Your package seems to be the only working for mixed multiple factor designs. The aovp alternative - from the orphaned lmperm package - produces huge variations for the p-values, see this. I have a few questions: Where can I found a bibliographical references for the implementation of ezPerm? How can I interpret that the function probably doesn't do interactions properly? What could be a post-hoc test in this case? Thanks! – toto_tico Nov 11 at 1:10
@Mike, would something like ezPerm( data = DATA, dv = DV, wid = WID, within = interaction(A,B), perms = 1e3) make sense to double check if the interaction is significative? – toto_tico Nov 11 at 1:14

When in doubt, bootstrap! Really, I don't know of a canned procedure to handle such a scenario.

Bootstrapping is a generally applicable way of generating some error parameters from the data at hand. Rather than relying on the typical parametric assumptions, bootstrap procedures capitalize on the characteristics of the sample to generate an empirical distribution against which your sample estimates can be compared.

Google scholar is's been done least once.

Lunneborg, Clifford E.; Tousignant, James P.; 1985 "Efron's Bootstrap with Application to the Repeated Measures Design." Multivariate Behavioral Research; Apr85, Vol. 20 Issue 2, p161, 18p

share|improve this answer
Thank you for the lead Brett! I wonder if someone got to implement it by now in R (I'd guess not). – Tal Galili Aug 5 '10 at 1:43
Right. R has lots of routines for assisting with bootstrap and other randomization methods, but I don't know that you'll find anything specific to this problem. – Brett Magill Aug 5 '10 at 3:15
Very nice first sentence. I hope it is not copyrighted, because I plan to use it :D – gui11aume Jun 28 '12 at 15:55

Consider checking out gmpm, but it is still in development.

share|improve this answer
Thank you for the link. I'll try to contact the author sometime in the future. – Tal Galili Aug 5 '10 at 17:52

There is a "trick" mentioned in some forums and mailing lists - I also found it mentioned in Joop Hox' book "Multilevel Analysis" (second edition, 2010), pp. 189.

The idea is: you reformat your long data into a long long dataset in which you create a new DV that includes all your DV responses, and use an index variable that holds information about the nature of the DVs to predict this outcome.

Let's assume you have 9 depression symptoms (ordinal), 2 measurement points, and 300 subjects. So while you have 300 rows in your normal dataset, and in your long dataset you'll have 600 rows, this new dataset will have 9 (symptoms) x 2 (time) x 300 (subjects) rows.

The new DV variable "symptoms" now contains the symptom severity of participants on 9 symptoms, the variables "index" contains the information about the nature of the symptom (1 to 9), and then there are the two variables "time" and "UserID".

You can now use the ordinal package to run this.

data<-read.csv("data_long_long.csv", head=T)

data$symptoms <- factor(data$symptoms)
data$time <- factor(data$time)
data$index <-factor(data$index)

m1<-clmm2(symptoms ~ index+time, random=UserID, data = data, Hess=TRUE, nAGQ=10)

In my specific case, I was interested whether there was a significant interaction between index and time, so I ran one additional model and compared them:

m2<-clmm2(symptoms ~ index+time, random=UserID, data = data, Hess=TRUE, nAGQ=10)

CLMM2 uses a random intercept model (to the best of my knowledge, the package ordinal does not do random slopes), if you do not with for a random intercept model you can run the models instead using CLM, e.g.:

m3<-clm(symptoms ~ index+time, data = data)
share|improve this answer
aren't m1 and m2 exactly the same, maybe you meant clmm2(symptoms ~ index*time, random=UserID, data = data, Hess=TRUE, nAGQ=10) anova(m1,m2). Do you know by any chance how to express this with the new clmm function? I am not sure about the notation. – toto_tico Nov 15 at 20:41

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.