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?

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

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)

Attribution
Source : Link , Question Author : Adam SA , Answer Author : Matt Albrecht

Leave a Comment