[R] ANOVA non-sphericity test and corrections (eg, Greenhouse-Geisser)

From: DarrenWeber <Darren.Weber_at_radiology.ucsf.edu>
Date: Sun, 24 Jun 2007 18:54:28 +0000

I'm an experimental psychologist and when I run ANOVA analysis in SPSS, I normally ask for a test of non-sphericity (Box's M-test). I also ask for output of the corrections for non-sphericity, such as Greenhouse-Geisser and Huhn-Feldt. These tests and correction factors are commonly used in the journals for experimental and other psychology reports. I have been switching from SPSS to R for over a year now, but I realize now that I don't have the non-sphericity test and correction factors.

The backgroud to this question is that I'm doing a psychophysics experiment and the data analysis uses an aov model that looks like this:

roiAOV <- aov( roi ~ (Cue*Hemisphere) + Error(Subject/ (Cue*Hemisphere)), data=roiDataframe)

where Cue is 2 levels (an arrow on a screen that points left or right) and Hemisphere is 2 levels (brain activity in the left or right cerebral hemisphere). There are 8 subjects and they all have one observation for all levels of Cue * Hemisphere (ie, a within-subjects design).

I need some functions for calculation of Box's M and the Greenhouse-  correction factor. Below is some code example that I have for the latter, but I don't know how to adapt it so it can take the aov object.

# BEGIN CODE BLOCK
# see pp. 45-47 of
# http://www.psych.upenn.edu/~baron/rpsych.pdf

# ---
# create some data

x0 <- rnorm(30)
x2 <- rnorm(30)
x4 <- rnorm(30)

data <- cbind(x0,x2,x4)

# n is the number of 'subjects' or rows; while
# k is the number of levels in the factor (columns),
# which is the number of repeated measures
n <- dim(data)[1]
k <- dim(data)[2]

# if k <= 2, the epsilon correction is not required.
# When the data are perfectly spherical, epsilon = 1.
# The minimum value possible for epsilon is

epsi = 1 / (k - 1)

# if epsi == 1 and k == 2, quit now...

# ---
# calculate variance-covariance matrix
# diagonal entries are variance,
# off-diagonal are covariance

S <- var(data)

# ---
# Now estimate Epsilon

D <- k^2 * ( mean(diag(S)) - mean(S) )^2

N1 <- sum(S^2)
N2 <- 2 * k * sum( apply(S,1,mean)^2 )
N3 <- k^2 * mean(S)^2

epsi <- D / ( (k-1) * (N1 - N2 + N3) )

# e is used to modify the degrees of freedom for the
# F test, so we have (k-1) becomes epsi(k-1) and also
# (k - 1)(n - 1) becomes epsi(k - 1)(n - 1). The new
# p-value for the F statistic is found with the
# pf() function. For example, if we have
# df1 = k - 1 = 3 - 1 = 2
# df2 = (k - 1)(n - 1) = (3 - 1)(10 - 1) = 18
# F = 40.719
# then the adjusted 2-tailed p-value is given by:
#

Fvalue <- 40.719
Pepsi <- 2 * (1 - pf(Fvalue, df1=epsi*(k-1), df2=epsi*(k-1)*(n-1) ) )

# Huynh-Feldt correction
# The Greenhouse-Geisser epsilon tends to underestimate
# epsilon when epsilon is greater than 0.70 (Stevens, 1990).
# An estimated e=0.96 may be actually 1. Huynh-Feldt
# correction is less conservative. The Huynh-Feldt
# epsilon is calculated from the Greenhouse-Geisser epsilon,

epsiHF <- (n * (k-1) * epsi - 2) / ((k-1) * ((n-1) - (k-1)*epsi))

PepsiHF <- 2 * (1 - pf(Fvalue, df1=epsiHF*(k-1), df2=epsiHF*(k-1)*(n-1) ) )

# MANOVA (and multivariate tests) may be better if the
# Greenhouse-Geisser and the Huynh-Feldt corrections do not agree,
# which may happen when epsilon drops below 0.70. When epsilon

# drops below 0.40, both the G-G and H-F corrections may indicate
# that the violation of sphericity is affecting the adjusted p-values.
# MANOVA is not always appropriate, though. MANOVA usually requires
# a larger sample size. Maxwell and Delaney (1990, p. 602) suggest
# a rough rule of thumb that the sample size n should be greater
# than k+10.
# END CODE BLOCK



R-help_at_stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 24 Jun 2007 - 23:29:21 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 25 Jun 2007 - 07:32:20 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.