# How to get Sphericity in R for a nested within subject design?

Thank you for reading. I am trying to get sphericity values for a purely within subject design. I have been unable to use ezANOVA, or Anova(). Anova works if I add a between subject factor, but I have been unable to get sphericity for a purely within subject design. Any advice?

My original ANOVA

anova(aov(resp ~ sucrose*citral, random =~1 | subject, data = p12bl, subset = exps==1))
anova(aov(resp ~ sucrose*citral, random =~1 | subject/sucrose*citral, data = p12bl, subset = exps==1))

> str(subset(p12bl, exps==1))
'data.frame':   192 obs. of  12 variables:
$exps : int 1 1 1 1 1 1 1 1 1 1 ...$ Order    : int  1 1 1 1 1 1 1 1 1 1 ...
$threshold: Factor w/ 2 levels " Suprathreshold",..: 1 1 1 1 1 1 1 1 1 1 ...$ SET      : Factor w/ 2 levels " A"," B": 1 1 1 1 1 1 1 1 1 1 ...
$subject : Factor w/ 16 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...$ stim     : chr  "S1C1" "S1C1" "S1C1" "S1C1" ...
$resp : num 6.01 5.63 0 2.57 6.81 ...$ id       : int  1 2 3 4 5 6 7 8 9 10 ...
$X1 : Factor w/ 1 level "S": 1 1 1 1 1 1 1 1 1 1 ...$ sucrose  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
$X3 : Factor w/ 1 level "C": 1 1 1 1 1 1 1 1 1 1 ...$ citral   : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...

subset(p12b,exps==1)
exps Order       threshold SET observ S1C1 S1C2 S1C3 S1C4 S2C1 S2C2 S2C3 S2C4 S3C1 S3C2 S3C3 S3C4 S4C1 S4C2 S4C3 S4C4
1     1     1  Suprathreshold   A      1  6.0  7.1  7.5  8.6 15.0 15.4 15.0 13.1 16.9   13 13.1 16.5   24   16   21   20
2     1     1  Suprathreshold   A      2  5.6  0.8  4.0  5.6  5.6 11.3 12.9 14.5 18.5   15 12.9 14.5   24   26   29   28
3     1     1  Suprathreshold   A      3  0.0  0.0  1.7  0.0  5.0  8.4  8.4  5.0 11.7   20 18.5 16.8   29   37   37   30
4     1     1  Suprathreshold   A      4  2.6  3.3  9.1 16.3  5.4 10.0  9.6 16.8 13.5   12 22.2 23.1   19   20   22   23
5     1     1  Suprathreshold   A      5  6.8  5.3 15.4 14.5 11.5  8.3 14.5 14.2  8.9   17 11.2 15.1   24   23   19   19
6     1     1  Suprathreshold   A      6  2.6  2.8  2.6  5.2 13.4 15.6 13.7 13.0 13.7   15 16.0 18.9   22   24   25   25
7     1     1  Suprathreshold   A      7  1.3  5.8 10.2  9.8 11.9 12.3 17.7 16.7 11.4   19 19.2 21.1   16   19   18   19
8     1     1  Suprathreshold   A      8  2.0  5.6  3.9  2.0  4.9  5.2  7.5  4.9 20.2   21  8.2  9.5   30   26   32   45
9     1     1  Suprathreshold   A      9  9.4 11.3 11.7 12.1 14.7 13.8 12.6 14.9 15.2   15 15.9 13.9   17   18   15   18
10    1     1  Suprathreshold   A     10  4.5 17.8 18.5 21.6  5.8 10.9 17.0 20.2  6.6   10 17.8 18.7   12   12   16   19
11    1     1  Suprathreshold   A     11  9.8 13.0 16.1 18.0 10.5 11.6 15.4 17.3 10.1   14 15.2 16.7   13   15   15   17
12    1     1  Suprathreshold   A     12  9.6 10.4 13.3 11.3 12.1 12.6 13.6 13.6 14.9   16 15.1 16.3   16   18   18   17


Sample output

ezANOVA( data = subset(p12bl, exps==1)  , dv= .(resp), sid = .(observ), within = .(sucrose,citral), between = NULL, collapse_within = FALSE)
Note: model has only an intercept; equivalent type-III tests substituted.
\$ANOVA
Effect DFn DFd  SSn  SSd     F       p p<.05   pes
1        sucrose   3  33 4953 3263 16.70 9.0e-07     * 0.603
2         citral   3  33  410  553  8.16 3.3e-04     * 0.426
3 sucrose:citral   9  99   56  791  0.77 6.4e-01       0.066


Warning messages: 1: You have removed one or more Ss from the analysis. Refactoring "observ" for ANOVA. 2: Too few Ss for Anova(), reverting to aov(). See "Warning" section of the help on ezANOVA.

• Do you want to look at sphericity or a factor to correct F-values because of violations in sphericity? – John Aug 25 '10 at 16:50
• @John I need to see if I need the correction – Adam SA Aug 25 '10 at 18:55

Try:

library(ez)
ezANOVA(data=subset(p12bl, exps==1),
within=.(sucrose, citral),
wid=.(subject),
dv=.(resp)
)


and the equivalent aov command, minus sphericity etc, is:

aov(resp ~ sucrose*citral + Error(subject/(sucrose*citral)),
data= subset(p12bl, exps==1))


Here's the equivalent using Anova from car directly:

library(car)
sucrose<-factor(rep(c(1:4), each=4))
citral<-factor(rep(c(1:4), 4))
idata<-data.frame(sucrose,citral)

mod<-lm(cbind(S1C1, S1C2, S1C3, S1C4, S2C1, S2C2, S2C3, S2C4,
S3C1, S3C2, S3C3, S3C4, S4C1, S4C2, S4C3, S4C4)~1, data=df1)
av.mod<-Anova(mod, idata=idata, idesign=~sucrose*citral)
summary(av.mod)

• Thank you for taking the time to answer my question. however, your response doesn't answer my question. I already mentioned in my post that I already came across the ezANOVA function. I tried it without success. I am also already familiar with the aov notation you metnion. – Adam SA Aug 25 '10 at 15:55
• Do you get any error messages when using ezANOVA? Whenever I use it, I get a nice table of sphericity tests and GG/HF corrections. It does not work for ANCOVA models, for this you will need to use the lme4 package, and there are other things to account for variance. – Matt Albrecht Aug 25 '10 at 16:05
• It would probably help if you put the output of "str(p12bl)" or "str(subset(p12bl, exp2==1))" here. – Matt Albrecht Aug 25 '10 at 16:16
• @Matt thank you for the very helpful edit. I was trying other things after the ~ except 1 - as in your answer. Thanks again! – Adam SA Aug 28 '10 at 16:41
• Ah, it seems that there's a bug in the current version of ez (v1.6.1). I've tracked it down to the fact that in converting to the formats that Anova() likes, ez was failing to turn something into a factor. I'm not sure whether this was always a bug, or whether Anova() recently changed in a way that broke ez. Regardless, I'm preparing ez v2.0, and this is one of the bug fixes. Stay tuned... – Mike Lawrence Aug 28 '10 at 17:16

Did you try the car package, from John Fox? It includes the function Anova() which is very useful when working with experimental designs. It should give you corrected p-value following Greenhouse-Geisser correction and Huynh-Feldt correction. I can post a quick R example if you wonder how to use it.

Also, there is a nice tutorial on the use of R with repeated measurements and mixed-effects model for psychology experiments and questionnaires; see Section 6.10 about sphericity.

As a sidenote, the Mauchly's Test of Sphericity is available in mauchly.test(), but it doesn't work with aov object if I remembered correctly. The R Newsletter from October 2007 includes a brief description of this topic.

• I was just about to recommend this. The ez package relies on car for some analysis. – John Aug 25 '10 at 18:47
• Thank you, the R newsletter is where I picked the mlm direction from. alas, I am still figuring how to choose idata idesign for a nested within subject design. I will try to figure out Anova(). – Adam SA Aug 25 '10 at 19:14
• @ John Anova() is supposed to take my aov() object as a model, but for some reason I am having problems making it work. – Adam SA Aug 26 '10 at 13:21
• @Adam In the on-line help for Anova(), you will find a working example with within- and between-subjects factors; the main caveats of Anova() is that you need to pass idata/idesign parameters in addition to your lm object. Also see this thread on the r-help mailing list, j.mp/9lphud. – chl Aug 26 '10 at 13:39
• Here is the correct link, j.mp/c6nDvn. – chl Aug 27 '10 at 15:30

I generally recommend avoiding these types of sphericity tests altogether by using modern mixed modeling methods. If you are not working with few subjects this will give you a great deal of flexibility in modeling an appropriate covariance structure, freeing you from the strict assumption of sphericity when necessary. I infer from the str output that you have 16 subjects with 12 observations each (I assume balance b/c you areusing clasical method-of-moments tools) which should be enough data to fit a mixed model with structured covariance matrices via (restricted) maximum likelihood.

Without being close to your data I can't offer specific model recommendations, but a place to start in R would be to replace aov in your model specifications with lme (after library(nlme)). The reason this will work is that you have mistakenly provided an nlme-style random argument to aov (when, as @Matt Albrecht pointed out, an Error term would have been approriate). In nlme, with the random argument set to ~ 1|<your grouping structure> and no correlation or weight arguments, you are specifying a random intercept for each group, implying the response covariance within groups is ZGZ' + R = 1G1'+ \sigma^2 I ==> compound symmetry with between-group variance off-diagonal and between-group variance + within-group variance on the diagonal ==> a spherical structure. From there you can begin to explore (e.g. using the built-in graphical methods), model, and test (e.g. comparing information criteria or using LRTs for nested models) the various forms of non-sphericity. Some of the tools for the modeling component are:

• Using the weights argument to model non-constant variance (diagonals) within or between groups (e.g. error variance changes between sucrose levels).
• Using the correlation argument to model non-constant covariance (off-diagonals) within groups (e.g. a structure within a group where residual errors that are closer together in time (e.g. AR1 structure) or space (e.g. Spherical structure) are more similar).
• Modeling random slopes by adding terms to the LHS of the | in the random formula.

Though the process can be complex with many potential pitfalls, I believe it will lead you to think more about the data generating mechanism, and when combined with careful graphical checks (I recommed lattice b/c nlme has excellent lattice-based plotting methods -- but ggplot works well too) you are likely to have not only a better scientific understanding of the process, but also less biased and more efficient estimators with which to draw inferences.

ez has now been updated to version 2.0. Among other improvements, the bug that caused it to fail to work for this example has been fixed.

• Comments are not for extended discussion; this conversation has been moved to chat. – whuber Jan 3 '17 at 15:20