Added anova.afex_aov() and nice() option to adjust F-tests for multiple comparisons

R/methods.afex_aov.R

@@ -2,6 +2,11 @@
 #' 
 #' Methods defined for objects returned from the ANOVA functions \code{\link{aov_car}} et al. of class \code{afex_aov} containing both the ANOVA fitted via \code{car::Anova} and base R's \code{aov}.
 #' 
+#' @param object,x object of class \code{afex_aov} as returned from \code{\link{aov_car}} and related functions.
+#' @param p.adjust.method \code{character} indicating if p-values for individual effects should be adjusted for multiple comparisons (see \link[stats]{p.adjust} and details).
+#' @param ... further arguments passed through, see description of return value for details.
+#' @param trms,xlev,grid same as for \code{\link{lsm.basis}}.
+#' 
 #' @return
 #' \describe{
 #'   \item{\code{anova}}{Returns an ANOVA table of class \code{c("anova", "data.frame")}. Information such as effect size (\code{es}) or df-correction are calculated each time this method is called.}
@@ -11,9 +16,12 @@
 #'   
 #' }
 #'
-#' @param object,x object of class \code{afex_aov} as returned from \code{\link{aov_car}} and related functions.
-#' @param ... further arguments passed through, see description of return value for details.
-#' @param trms,xlev,grid same as for \code{\link{lsm.basis}}.
+#' @details 
+#' Exploratory ANOVA, for which no detailed hypotheses have been specified a priori, harbor a multiple comparison problem (Cramer et al., 2015). To avoid an inflation of familywise Type I error rate, results need to be corrected for multiple comparisons using \code{p.adjust.method}.
+#' \code{p.adjust.method} defaults to the method specified in the call to \code{\link{aov_car}} in \code{anova_table}. If no method was specified and \code{p.adjust.method = NULL} p-values are not adjusted.
+#' 
+#' @references 
+#' Cramer, A. O. J., van Ravenzwaaij, D., Matzke, D., Steingroever, H., Wetzels, R., Grasman, R. P. P. P., ... Wagenmakers, E.-J. (2015). Hidden multiplicity in exploratory multiway ANOVA: Prevalence and remedies.  \emph{Psychonomic Bulletin & Review}, 1–8. doi:\href{http://doi.org/10.3758/s13423-015-0913-5}{10.3758/s13423-015-0913-5}
 #' 
 #' @name afex_aov-methods
 NULL
@@ -23,7 +31,7 @@ NULL
 #' @rdname afex_aov-methods
 #' @inheritParams nice
 #' @export
-anova.afex_aov <- function(object, es = afex_options("es_aov"), observed = NULL, correction = afex_options("correction_aov"), MSE = TRUE, intercept = FALSE, ...) {
+anova.afex_aov <- function(object, es = afex_options("es_aov"), observed = NULL, correction = afex_options("correction_aov"), MSE = TRUE, intercept = FALSE, p.adjust.method = NULL, ...) {
   # internal functions:
   # check arguments
   es <- match.arg(es, c("none", "ges", "pes"), several.ok = TRUE)
@@ -76,7 +84,7 @@ anova.afex_aov <- function(object, es = afex_options("es_aov"), observed = NULL,
       }
       obs_SSn1 <- sum(tmp2$SS*obs)
       obs_SSn2 <- tmp2$SS*obs
-    }else{
+    } else {
       obs_SSn1 <- 0
       obs_SSn2 <- 0
     }
@@ -88,6 +96,10 @@ anova.afex_aov <- function(object, es = afex_options("es_aov"), observed = NULL,
   #browser()
   if (!MSE) anova_table$MSE <- NULL 
   if (!intercept) if (row.names(anova_table)[1] == "(Intercept)")  anova_table <- anova_table[-1,, drop = FALSE]
+  # Correct for multiple comparisons
+  if(is.null(p.adjust.method)) p.adjust.method <- ifelse(is.null(attr(object$anova_table, "p.adjust.method")), "none", attr(object$anova_table, "p.adjust.method"))
+  anova_table[,"Pr(>F)"] <- p.adjust(anova_table[,"Pr(>F)"], method = p.adjust.method)
+  attr(anova_table, "p.adjust.method") <- p.adjust.method
   anova_table
 }
 
@@ -106,6 +118,7 @@ print.afex_aov <- function(x, ...) {
 #' @export
 summary.afex_aov <- function(object, ...) {
   if (class(object$Anova)[1] == "Anova.mlm") {
+    if(attr(object$anova_table, "p.adjust.method") != "none") message("Note, results are NOT adjusted for multiple comparisons as requested (p.adjust.method = '", attr(object$anova_table, "p.adjust.method"), "') because the desired method of sphericity correction is unknown. For adjusted p-values print the object (to see object$anova_table), or call one of anova.afex_aov() or nice().")
     return(summary(object$Anova, multivariate = FALSE))
   } else if (class(object$Anova)[1] == "anova") {
     return(object$anova_table)

R/nice.R

@@ -7,6 +7,7 @@
 #' @param es Effect Size to be reported. The default is given by \code{afex_options("es_aov")}, which is initially set to \code{"ges"} (i.e., reporting generalized eta-squared, see details). Also supported is partial eta-squared (\code{"pes"}) or \code{"none"}.
 #' @param observed character vector referring to the observed (i.e., non manipulated) variables/effects in the design. Important for calculation of generalized eta-squared (ignored if \code{es} is not \code{"ges"}), see details.
 #' @param correction Character. Which sphericity correction of the degrees of freedom should be reported for the within-subject factors.  The default is given by \code{afex_options("correction_aov")}, which is initially set to \code{"GG"} corresponding to the Greenhouse-Geisser correction. Possible values are \code{"GG"}, \code{"HF"} (i.e., Hyunh-Feldt correction), and \code{"none"} (i.e., no correction).
+#' @param p.adjust.method \code{character} indicating if p-values for individual effects should be adjusted for multiple comparisons (see \link[stats]{p.adjust} and details).
 #' @param sig.symbols Character. What should be the symbols designating significance? When entering an vector with \code{length(sig.symbol) < 4} only those elements of the default (\code{c(" +", " *", " **", " ***")}) will be replaced. \code{sig.symbols = ""} will display the stars but not the \code{+}, \code{sig.symbols = rep("", 4)} will display no symbols.
 #' @param MSE logical. Should the column containing the Mean Sqaured Error (MSE) be displayed? Default is \code{TRUE}.
 #' @param intercept logical. Should intercept (if present) be included in the ANOVA table? Default is \code{FALSE} which hides the intercept.
@@ -21,13 +22,17 @@
 #' Conversion functions to other formats (such as HTML, ODF, or Word) can be found at the \href{http://cran.r-project.org/web/views/ReproducibleResearch.html}{Reproducible Research Task View}.
 #'
 #' The default reports generalized eta squared (Olejnik & Algina, 2003), the "recommended effect size for repeated measured designs" (Bakeman, 2005). Note that it is important that all measured variables (as opposed to experimentally manipulated variables), such as e.g., age, gender, weight, ..., must be declared via \code{observed} to obtain the correct effect size estimate. Partial eta squared (\code{"pes"}) does not require this.
-#'
+#' 
+#' Exploratory ANOVA, for which no detailed hypotheses have been specified a priori, harbor a multiple comparison problem (Cramer et al., 2015). To avoid an inflation of familywise Type I error rate, results need to be corrected for multiple comparisons using \code{p.adjust.method}.
+#' \code{p.adjust.method} defaults to the method specified in the call to \code{\link{aov_car}} in \code{anova_table}. If no method was specified and \code{p.adjust.method = NULL} p-values are not adjusted.
+#' 
 #' @seealso \code{\link{aov_ez}} and \code{\link{aov_car}} are the convenience functions to create the object appropriate for \code{nice_anova}.
 #'
 #' @author The code for calculating generalized eta-squared was written by Mike Lawrence.\cr Everything else was written by Henrik Singmann.
 #'
 #' @references Bakeman, R. (2005). Recommended effect size statistics for repeated measures designs. \emph{Behavior Research Methods}, 37(3), 379-384. doi:10.3758/BF03192707
-
+#'
+#' Cramer, A. O. J., van Ravenzwaaij, D., Matzke, D., Steingroever, H., Wetzels, R., Grasman, R. P. P. P., ... Wagenmakers, E.-J. (2015). Hidden multiplicity in exploratory multiway ANOVA: Prevalence and remedies.  \emph{Psychonomic Bulletin & Review}, 1–8. doi:\href{http://doi.org/10.3758/s13423-015-0913-5}{10.3758/s13423-015-0913-5}
 #'
 #' Olejnik, S., & Algina, J. (2003). Generalized Eta and Omega Squared Statistics: Measures of Effect Size for Some Common Research Designs. \emph{Psychological Methods}, 8(4), 434-447. doi:10.1037/1082-989X.8.4.434
 #' 
@@ -76,8 +81,9 @@ nice <- function(object, ...) UseMethod("nice", object)
 #' @rdname nice
 #' @method nice afex_aov
 #' @export
-nice.afex_aov <- function(object, es = afex_options("es_aov"), observed = NULL, correction = afex_options("correction_aov"), MSE = TRUE, intercept = FALSE, sig.symbols = c(" +", " *", " **", " ***"), ...) { 
-  anova_table <- as.data.frame(anova(object, es = es, observed = observed, correction = correction, MSE = MSE, intercept = intercept))
+nice.afex_aov <- function(object, es = afex_options("es_aov"), observed = NULL, correction = afex_options("correction_aov"), MSE = TRUE, intercept = FALSE, sig.symbols = c(" +", " *", " **", " ***"), p.adjust.method = NULL, ...) { 
+  if(is.null(p.adjust.method)) p.adjust.method <- ifelse(is.null(attr(object$anova_table, "p.adjust.method")), "none", attr(object$anova_table, "p.adjust.method"))
+  anova_table <- as.data.frame(anova(object, es = es, observed = observed, correction = correction, MSE = MSE, intercept = intercept, p.adjust.method = p.adjust.method))
   nice.anova(anova_table, MSE = MSE, intercept = intercept, sig.symbols = sig.symbols)
 }
 

man/afex_aov-methods.Rd

@@ -11,7 +11,7 @@
 \usage{
 \method{anova}{afex_aov}(object, es = afex_options("es_aov"),
   observed = NULL, correction = afex_options("correction_aov"),
-  MSE = TRUE, intercept = FALSE, ...)
+  MSE = TRUE, intercept = FALSE, p.adjust.method = NULL, ...)
 
 \method{print}{afex_aov}(x, ...)
 
@@ -34,6 +34,8 @@
 
 \item{intercept}{logical. Should intercept (if present) be included in the ANOVA table? Default is \code{FALSE} which hides the intercept.}
 
+\item{p.adjust.method}{\code{character} indicating if p-values for individual effects should be adjusted for multiple comparisons (see \link[stats]{p.adjust} and details).}
+
 \item{...}{further arguments passed through, see description of return value for details.}
 
 \item{trms,xlev,grid}{same as for \code{\link{lsm.basis}}.}
@@ -50,4 +52,11 @@
 \description{
 Methods defined for objects returned from the ANOVA functions \code{\link{aov_car}} et al. of class \code{afex_aov} containing both the ANOVA fitted via \code{car::Anova} and base R's \code{aov}.
 }
+\details{
+Exploratory ANOVA, for which no detailed hypotheses have been specified a priori, harbor a multiple comparison problem (Cramer et al., 2015). To avoid an inflation of familywise Type I error rate, results need to be corrected for multiple comparisons using \code{p.adjust.method}.
+\code{p.adjust.method} defaults to the method specified in the call to \code{\link{aov_car}} in \code{anova_table}. If no method was specified and \code{p.adjust.method = NULL} p-values are not adjusted.
+}
+\references{
+Cramer, A. O. J., van Ravenzwaaij, D., Matzke, D., Steingroever, H., Wetzels, R., Grasman, R. P. P. P., ... Wagenmakers, E.-J. (2015). Hidden multiplicity in exploratory multiway ANOVA: Prevalence and remedies.  \emph{Psychonomic Bulletin & Review}, 1–8. doi:\href{http://doi.org/10.3758/s13423-015-0913-5}{10.3758/s13423-015-0913-5}
+}
 

man/nice.Rd

@@ -13,7 +13,7 @@ nice(object, ...)
 \method{nice}{afex_aov}(object, es = afex_options("es_aov"),
   observed = NULL, correction = afex_options("correction_aov"),
   MSE = TRUE, intercept = FALSE, sig.symbols = c(" +", " *", " **",
-  " ***"), ...)
+  " ***"), p.adjust.method = NULL, ...)
 
 \method{nice}{anova}(object, MSE = TRUE, intercept = FALSE,
   sig.symbols = c(" +", " *", " **", " ***"), ...)
@@ -36,6 +36,8 @@ nice(object, ...)
 \item{intercept}{logical. Should intercept (if present) be included in the ANOVA table? Default is \code{FALSE} which hides the intercept.}
 
 \item{sig.symbols}{Character. What should be the symbols designating significance? When entering an vector with \code{length(sig.symbol) < 4} only those elements of the default (\code{c(" +", " *", " **", " ***")}) will be replaced. \code{sig.symbols = ""} will display the stars but not the \code{+}, \code{sig.symbols = rep("", 4)} will display no symbols.}
+
+\item{p.adjust.method}{\code{character} indicating if p-values for individual effects should be adjusted for multiple comparisons (see \link[stats]{p.adjust} and details).}
 }
 \value{
 A \code{data.frame} with the ANOVA table consisting of characters. The columns that are always present are: \code{Effect}, \code{df} (degrees of freedom), \code{F}, and \code{p}.
@@ -51,6 +53,9 @@ The returned \code{data.frame} is print-ready when adding to a document with pro
 Conversion functions to other formats (such as HTML, ODF, or Word) can be found at the \href{http://cran.r-project.org/web/views/ReproducibleResearch.html}{Reproducible Research Task View}.
 
 The default reports generalized eta squared (Olejnik & Algina, 2003), the "recommended effect size for repeated measured designs" (Bakeman, 2005). Note that it is important that all measured variables (as opposed to experimentally manipulated variables), such as e.g., age, gender, weight, ..., must be declared via \code{observed} to obtain the correct effect size estimate. Partial eta squared (\code{"pes"}) does not require this.
+
+Exploratory ANOVA, for which no detailed hypotheses have been specified a priori, harbor a multiple comparison problem (Cramer et al., 2015). To avoid an inflation of familywise Type I error rate, results need to be corrected for multiple comparisons using \code{p.adjust.method}.
+\code{p.adjust.method} defaults to the method specified in the call to \code{\link{aov_car}} in \code{anova_table}. If no method was specified and \code{p.adjust.method = NULL} p-values are not adjusted.
 }
 \examples{
 ## example from Olejnik & Algina (2003)
@@ -89,6 +94,8 @@ The code for calculating generalized eta-squared was written by Mike Lawrence.\c
 \references{
 Bakeman, R. (2005). Recommended effect size statistics for repeated measures designs. \emph{Behavior Research Methods}, 37(3), 379-384. doi:10.3758/BF03192707
 
+Cramer, A. O. J., van Ravenzwaaij, D., Matzke, D., Steingroever, H., Wetzels, R., Grasman, R. P. P. P., ... Wagenmakers, E.-J. (2015). Hidden multiplicity in exploratory multiway ANOVA: Prevalence and remedies.  \emph{Psychonomic Bulletin & Review}, 1–8. doi:\href{http://doi.org/10.3758/s13423-015-0913-5}{10.3758/s13423-015-0913-5}
+
 Olejnik, S., & Algina, J. (2003). Generalized Eta and Omega Squared Statistics: Measures of Effect Size for Some Common Research Designs. \emph{Psychological Methods}, 8(4), 434-447. doi:10.1037/1082-989X.8.4.434
 }
 \seealso{

tests/testthat/test-aov_car-basic.car-basic.R

@@ -84,3 +84,41 @@ test_that("Data from O'Brien & Kaiser replicates their paper (p. 328, Table 8, c
 
 })
 
+test_that("Data from O'Brien & Kaiser adjusted for familywise error rate (p. 328, Table 8, column 'average'", {
+  data(obk.long, package = "afex")
+  out1 <- aov_car(value ~ treatment * gender + Error(id/(phase*hour)), data = obk.long, observed = "gender", return = "afex_aov", anova_table = list(correction = "none", p.adjust.method = "bonferroni"))
+
+  expect_that(unname(unlist(out1[["anova_table"]]["treatment", c("num Df", "den Df", "F")])), equals(c(2, 10, 3.94), tolerance = 0.001))
+  expect_that(unname(unlist(out1[["anova_table"]]["gender", c("num Df", "den Df", "F")])), equals(c(1, 10, 3.66), tolerance = 0.001))
+  expect_that(round(unname(unlist(out1[["anova_table"]]["treatment:gender", c("num Df", "den Df", "F")])), 2), equals(c(2, 10, 2.86), tolerance = 0.001))
+
+  ## check against own results:
+  anova_tab <- structure(list(`num Df` = c(2, 1, 2, 2, 4, 2, 4, 4, 8, 4, 8, 
+  8, 16, 8, 16), `den Df` = c(10, 10, 10, 20, 20, 20, 20, 40, 40, 
+  40, 40, 80, 80, 80, 80), MSE = c(22.8055555555555, 22.8055555555555, 
+  22.8055555555555, 4.01388888888889, 4.01388888888889, 4.01388888888889, 
+  4.01388888888889, 1.5625, 1.5625, 1.5625, 1.5625, 1.20208333333333, 
+  1.20208333333333, 1.20208333333333, 1.20208333333333), F = c(3.940494501098, 
+  3.65912050065102, 2.85547267441343, 16.1329196993199, 4.85098375975551, 
+  0.282782484190432, 0.636602429722426, 16.6856704980843, 0.0933333333333336, 
+  0.450268199233716, 0.620437956204379, 1.17990398215104, 0.345292160558641, 
+  0.931293452060798, 0.735935938468544), ges = c(0.198248507309966, 
+  0.114806410630587, 0.179183259116394, 0.151232705544895, 0.0967823866181358, 
+  0.00312317714869712, 0.0140618480455475, 0.12547183572154, 0.00160250371109459, 
+  0.0038716854273722, 0.010669821220833, 0.0153706689696344, 0.00905399063368842, 
+  0.012321395080303, 0.0194734697889242), `Pr(>F)` = c(0.0547069269265198, 
+  0.0848002538616402, 0.104469234023772, 6.73163655770545e-05, 
+  0.00672273209545241, 0.756647338927411, 0.642369488905348, 4.02664339633774e-08, 
+  0.999244623719389, 0.771559070589063, 0.755484449904079, 0.32158661418337, 
+  0.990124565656718, 0.495611922963992, 0.749561639456282)), .Names = c("num Df", 
+  "den Df", "MSE", "F", "ges", "Pr(>F)"), heading = c("Anova Table (Type 3 tests)\n", 
+  "Response: value"), row.names = c("treatment", "gender", "treatment:gender", 
+  "phase", "treatment:phase", "gender:phase", "treatment:gender:phase", 
+  "hour", "treatment:hour", "gender:hour", "treatment:gender:hour", 
+  "phase:hour", "treatment:phase:hour", "gender:phase:hour", "treatment:gender:phase:hour"
+  ), class = c("data.frame"))
+  anova_tab$`Pr(>F)` <- p.adjust(anova_tab$`Pr(>F)`, method = "bonferroni")
+
+  expect_equal(out1[["anova_table"]], anova_tab, check.attributes = FALSE)
+
+})