---
title: "lecture: Alpha Adjustment and ANOVA"
output: html_document
---

```{r}
simulation <- function(E, N, Mean1, Mean2, Type) {
  occ <- c()
  for (i in 1:E){
  family <- replicate(n = N, t <- t.test(rnorm(30, mean = Mean1), rnorm(30, mean = Mean2))$p.value)
  
  # Two-Step
  l1 <- length(family[1][family[1] <= 0.05]) 
  if (l1 > 0) {
    l1 <- l1 + length(family[2:N][family[2:N] <= 0.05/(N-1)])
  }
  
  # Holm
  f_sorted <- sort(family)
  f_counter <- 0
  for (b in 1:length(f_sorted)) {
    if (f_sorted[b] <= 0.05/(length(f_sorted) - b + 1)) {
      f_counter <- f_counter + 1
    } else {
      break
    }
  }
  
  # Bind everything together
  if (Type == "Type-I") {
    occ <- rbind(occ, c(length(family[family <= 0.05]), length(family[family <= 0.05/N]), l1, f_counter))
  }
  else {
    occ <- rbind(occ, c(length(family) - length(family[family <= 0.05]), length(family) - length(family[family <= 0.05/N]), length(family) - l1, length(family) - f_counter))
  }

  }
  return(occ)
}
```


```{r}

E <- 1000

# Given a set of families of true H0; In how many families will a Type-1 error be made?
occ <- simulation(E, 10, 0.0, 0.0, "Type-I")
df <- data.frame(occ)
names(df) <- c("uncorrected", "bonferroni", "twostep", "holm")

# FWER
cat("In how many families did at least one Type-1 error occur? (FWER)\n")
cat("Uncorrected: ", length(df$uncorrected[df$uncorrected > 0]), length(df$uncorrected[df$uncorrected > 0]) / E, "\n")
cat("Bonferroni: ", length(df$bonferroni[df$bonferroni > 0]), length(df$bonferroni[df$bonferroni > 0]) / E, "\n")
cat("Two-Step: ", length(df$twostep[df$twostep > 0]), length(df$twostep[df$twostep > 0]) / E, "\n")
cat("Holm: ", length(df$holm[df$holm > 0]), length(df$twostep[df$holm > 0]) / E, "\n")

# PFER
cat("How many Type-1 errors occured in a family on average? (PFER)\n")
cat("Uncorrected: ", mean(df$uncorrected), ", total: ", mean(df$uncorrected) * E, "\n")
cat("Bonferroni: ", mean(df$bonferroni), ", total: ", mean(df$bonferroni) * E, "\n")
cat("Two-Step: ", mean(df$twostep), ", total: ", mean(df$twostep) * E, "\n")
cat("Holm: ", mean(df$holm), ", total: ", mean(df$holm) * E,"\n")

```


```{r}

# Given a set of families of false H0; In how many families will a Type-2 error be made?

occ <- simulation(E, 10, 0.0, 1.0, "Type-II")
df <- data.frame(occ)
names(df) <- c("uncorrected", "bonferroni", "twostep", "holm")

# FWER
cat("In how many families did at least one Type-2 error occur? (FWER)\n")
cat("Uncorrected: ", length(df$uncorrected[df$uncorrected > 0]), length(df$uncorrected[df$uncorrected > 0]) / E, "\n")
cat("Bonferroni: ", length(df$bonferroni[df$bonferroni > 0]), length(df$bonferroni[df$bonferroni > 0]) / E, "\n")
cat("Two-Step: ", length(df$twostep[df$twostep > 0]), length(df$twostep[df$twostep > 0]) / E, "\n")
cat("Holm: ", length(df$holm[df$holm > 0]), length(df$twostep[df$holm > 0]) / E, "\n")

# PFER
cat("How many Type-2 errors occured in a family on average? (PFER)\n")
cat("Uncorrected: ", mean(df$uncorrected), ", total: ", mean(df$uncorrected) * E, "\n")
cat("Bonferroni: ", mean(df$bonferroni), ", total: ", mean(df$bonferroni) * E, "\n")
cat("Two-Step: ", mean(df$twostep), ", total: ", mean(df$twostep) * E, "\n")
cat("Holm: ", mean(df$holm), ", total: ", mean(df$holm) * E,"\n")
```


```{r}
robot.humor <- c(12, 16, 17, 18, 11)
robot.neutral <- c(13, 12, 15, 18, 9)
robot.strict <- c(20, 21, 17, 23, 24)
treatment <- factor(c(rep("RN", 5), rep("RH", 5), rep("RS", 5)))
data <- data.frame(treatment, c(robot.neutral, robot.humor, robot.strict))
names(data) <- c("Treatment", "Exam.Scores")
data
plot(Exam.Scores ~ Treatment, data = data)
```


```{r}
a <- aov(Exam.Scores ~ Treatment, data = data)
a
```

```{r}
summary(a)
```

```{r}
contrasts(data$Treatment) <- cbind(c(-1/2, 1, -1/2), c(-1, 0, 1))
contrasts(data$Treatment)
a <- aov(Exam.Scores ~ Treatment, data = data)
summary(a, split=list(Treatment=list("Neutral vs. (Humor & Strict)"=1, "Humor vs. Strict"=2)))
```

```{r}
t.test(robot.strict, robot.humor, var.equal = T)
summary(aov(Score ~ Cond, data = data.frame(Cond = c(rep("RH", 5), rep("RS", 5)), Score = c(robot.humor, robot.strict))))
```