A non-parametric repeated-measures multi-way Anova in R? 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:


*

*http://cran.r-project.org/web/packages/repolr/

*http://cran.r-project.org/web/packages/ordinal/
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") ?
 A: 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. 
A: 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 gold...it's been done before...at 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
A: 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)
anova(m1,m2)

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)

