From: DarrenWeber <Darren.Weber_at_radiology.ucsf.edu>

Date: Sun, 24 Jun 2007 18:54:28 +0000

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

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.
*