context("ANOVAs: replicating published results")

test_that("purely within ANOVA, return='univ': Maxell & Delaney (2004), Table 12.5 and 12.6, p. 578", {
  ### replicate results from Table 12.6
  data(md_12.1)
  # valus from table:
  f <- c(40.72, 33.77, 45.31)
  ss_num <- c(289920, 285660, 105120)
  ss_error <- c(64080, 76140, 20880)
  num_df <- c(2, 1, 2)
  den_df <- c(18, 9, 18)
  
  md_ez_r <- aov_ez("id", "rt", md_12.1, within = c("angle", "noise"))
  md_car_r <- aov_car(rt ~ 1 + Error(id/angle*noise), md_12.1)
  md_aov_4_r <- aov_4(rt ~ 1 + (angle*noise|id), md_12.1)
  
  expect_that(md_ez_r, is_equivalent_to(md_car_r))
  expect_that(md_ez_r, is_equivalent_to(md_aov_4_r))
  expect_that(round(md_ez_r$anova_table[,"F"], 2), is_equivalent_to(f))
  tmp_univ <- suppressWarnings(summary(md_ez_r$Anova)$univariate.tests)
  if ("Sum Sq" %in% colnames(tmp_univ)) { # to allow both car 3.0 and earlier versions
    ss_test <- tmp_univ[,"Sum Sq"][-1]
  } else {
    ss_test <- tmp_univ[,"SS"][-1]
  }
  expect_that(ss_test, is_equivalent_to(ss_num))  
  expect_that(tmp_univ[,"Error SS"][-1], is_equivalent_to(ss_error))
  expect_that(anova(md_ez_r, correction = "none")[,"num Df"], is_equivalent_to(num_df))
  expect_that(anova(md_ez_r, correction = "none")[,"den Df"], is_equivalent_to(den_df))
})

test_that("Analysis of Singmann & Klauer (2011, Exp. 1)", {
  data(sk2011.1, package = "afex")

  out1 <-  aov_ez("id", "response", sk2011.1[ sk2011.1$what == "affirmation",], within = c("inference", "type"), between = "instruction", anova_table=(es = "pes"), fun_aggregate = mean, return = "afex_aov")
  
  df_num <- rep(1, 7)
  df_den <- rep(38, 7)
  MSE <- c(1072.42, 1007.21, 1007.21, 187.9, 187.9, 498.48, 498.48)
  F <- c(0.13, 13.01, 12.44, 0.06, 3.09, 29.62, 10.73)
  pes <- c(0, 0.26, 0.25, 0, 0.08, 0.44, 0.22)
  p <- c(0.72, 0.0009, 0.001, 0.81, 0.09, 0.001, 0.002)
  
  expect_that(out1$anova_table[["num Df"]], is_equivalent_to(df_num))
  expect_that(out1$anova_table[["den Df"]], is_equivalent_to(df_den))
  expect_that(out1$anova_table[["MSE"]], equals(MSE, tolerance = 0.001))
  expect_that(out1$anova_table[["F"]], equals(F, tolerance = 0.001))
  expect_that(out1$anova_table[["pes"]], equals(pes, tolerance = 0.02))
  expect_that(out1$anova_table[["Pr(>F)"]], equals(p, tolerance = 0.01))
  
})


test_that("Data from O'Brien & Kaiser replicates their paper (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"))

  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"))
  
  expect_equal(out1[["anova_table"]], anova_tab, check.attributes = FALSE)
  
})

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)
  
})

test_that("afex_aov printing", {
  data(sk2011.1, package = "afex")
  out_new <-  aov_ez("id", "response", sk2011.1[ sk2011.1$what == "affirmation",], within = c("inference", "type"), between = "instruction", anova_table=(es = "pes"), fun_aggregate = mean, return = "afex_aov")
  
  expect_output(print(out_new), "Signif. codes")
  expect_output(print(anova(out_new)), "Signif. codes")
  expect_output(print(nice(out_new)), "Anova")
  
  load("afex_aov_16_1.rda")
  expect_output(print(out1), "Anova")
  skip_if_not(packageVersion("car") >= '3.0.13')
  ## only throws error with car version '3.0.13' or larger, the
  expect_error(anova(out1),
               regexp = "summary.Anova.mlm() failed for 'afex_aov' object.",
               fixed = TRUE)
  expect_error(nice(out1),
               regexp = "summary.Anova.mlm() failed for 'afex_aov' object.",
               fixed = TRUE)
})