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

Question

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?

I already read the R newsletter, fox chapter appendix, EZanova, and whatever I could find online.

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.

Answer 1

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)
df1<-read.table("clipboard", header=T) #From copying data in the question above
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)

Answer 2

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.

Answer 3

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.

Answer 4

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.