# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# One - Way and Two - Way ANOVA                                             ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# This script includes the basic analysis One - Way and Two - Way ANOVA.

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Import the Data        ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Read the data (the first row contains variable names)
cars <- read.table(file.path("data", "cars.txt"), header = TRUE)

# Check the first observations of the data
head(cars)

# Dependent Variable
accel <- cars$Acceleration

# Independent Variables
type <- cars$Type
brand <- cars$Brand

n <- length(accel)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Factors               ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# The variables Type and Brand are categorical
# Type
#   0: std
#   1: turbo
# Brand
#   1: BMW
#   2: Datsun
#   3: Mercedes-Benz

# We have to explicitly inform R about the nature of these variables, since
# 0, 1, 2, and 3 are by default understood as numbers, not coded categories.

# The class() function informs us about the class of the variables.

class(type)
class(brand)

# We can turn the variables into factors:

type <- factor(type, levels = c(0, 1), labels = c("std", "turbo"))

brand <- factor(brand)
levels(brand)
levels(brand) <- c('BMW', 'Datsun', 'Mercedes')
levels(brand)

class(type)
class(brand)

# If a categorical variable is ordered, we can code:
# factor(type, levels = c('std' = "0", 'turbo' = "1"), ordered = TRUE)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Example 1: Type        ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# We will test and compare the mean acceleration times of the two types of
# engines by applying one-way ANOVA method as there is only one factor or
# criterion to classify the sample observations.

# This example shows that when the factor has only two categories (types has
# std and turbo), t-test and ANOVA produce the same results. In fact, One - Way
# ANOVA is a generalization of t-test in that direction, as it can be used for
# factors with more than two categories.

# Error Bar
# install.packages("gplots")
library(gplots)
plotmeans(accel ~ type, p = 0.95, connect = FALSE, n.label = FALSE,
          main = "Error Bar for Acceleration by Type",
          ylab = "Acceleration")

# ANOVA
#   H0: mu_std  = mu_turbo
#   H1: mu_std != mu_turbo

# Option 1: We can use the built-in function aov.
aov_type <- aov(accel ~ type)
summary(aov_type)

# Residual Diagnostics (Normality Assumption)
resid <- aov_type$residuals
qqnorm(resid)
qqline(resid)
shapiro.test(resid)

# Option 2: We can also perform the One - Way ANOVA manually.

# First, we break our sample in two pieces:

# Acceleration for the cars of type 1 - turbo
accel_turbo <- accel[type == 'turbo']

# Acceleration for the cars of type 0 - std
accel_std <- accel[type == 'std']

# Number of factor's levels
m <- 2

# Calculate the overall and the group means.

y <- accel
y1 <- accel_turbo
y0 <- accel_std

n1 <- length(accel_turbo)
n0 <- n - n1

y_bar <- mean(accel)
y1_bar <- mean(accel_turbo)
y0_bar <- mean(accel_std)

SSE <- sum((y1 - y1_bar) ^ 2) + sum((y0 - y0_bar) ^ 2)
SSF <- n0 * (y0_bar - y_bar) ^ 2 + n1 * (y1_bar - y_bar) ^ 2
SST <- sum((y - y_bar) ^ 2)
SST
SSE + SSF

F0 <- (SSF / (m - 1)) / (SSE / (n - m))
F0 # 48.21515
# F0 ~ F_(m-1),(n-m)

# How do we respond to a hypothesis test?

# Option A: Critical Value

a <- 0.05
qf_a <- qf(a, m - 1, n - m, lower.tail = FALSE)
qf_a # 3.865475

# We notice that F0 > qf_a, so we reject the null hypothesis H0 at the a = 0.05
# significance level.

# Option B: p-value

pvalue <- pf(F0, m - 1, n - m, lower.tail = FALSE)
pvalue # 1.610205e-11

# p_value = 1.610205e-11 < 0.05, so we reject H0 not only at the a = 0.05,
# significance level, but also at a = 0.01, or even a = 0.001.

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Example 2: Brand       ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# We can also test and compare the acceleration times of the three different
# brands, by applying one - way ANOVA method as there is, again, only one factor
# or criterion to classify the sample observations. The null and alternative
# hypothesis to be used are given as:

# H0: mu_BMW = mu_Datsun = mu_Mercedes
# H1: not all means are equal

# Error Bar
plotmeans(accel ~ brand, p = 0.95, connect = FALSE,
          xlab = "Brands", ylab = "Acceleration",
          main = "Mean Plot with 95% CI", n.label = FALSE)

# ANOVA
aov_brand <- aov(accel ~ brand)
summary(aov_brand)

# The one-way ANOVA test does not inform which group has a different mean.
# Instead, we can perform a Tukey test with the function TukeyHSD().

# Pairwise comparison
TukeyHSD(aov_brand)

# We are interested in the last column, titled "p adj", where we can see which
# difference is statistically significant, or, equivalently, in the two middle
# columns, titled "lwr" and "upr", where we can see the lower and upper bounds
# of each difference's confidence interval. Here we notice that the only not
# statistically significant difference is that between the means of the brands
# Mercedes-Benz and Datsun (p-value > 0.05 or the confidence interval contains
# zero). Specifically:

# We reject mu_Datsun = mu_BMW (p adj < 0.05)
# We reject mu_Mercedes = mu_BMW (p adj < 0.05)
# We do not reject mu_Mercedes = mu_BMW (p adj > 0.05)

# Residual Diagnostics (Normality Assumption)
resid <- aov_brand$residuals
qqnorm(resid)
qqline(resid)
shapiro.test(resid)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Two-Way ANOVA          ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# We could use ANOVA with more than one factor. We want to look for differences
# in the acceleration means when both Type and Brand are taken into account
# (2 X 3 = 6 categories). This model is named Two-Way ANOVA.

# We can use the built-in function aov.

# We first include in our analysis the main effect of each variable and
# interaction effect.
aov1 <- aov(y ~ type + brand + type:brand)
summary(aov1)

# or simply
aov1 <- aov(y ~ type * brand)
summary(aov1)

# We notice that the interaction term, "type:brand" is not statistically
# significant so we continue our analysis without it.
aov2 <- aov(y ~ type + brand)
summary(aov2)

# Residual Diagnostics (Normality Assumption)
resid <- aov2$residuals
qqnorm(resid)
qqline(resid)
shapiro.test(resid)

# Two-way Interaction Plot
interaction.plot(brand, type, y, type = "b", col = c(1:3),
                 leg.bty = "n", lwd = 2, pch = c(18, 24),
                 xlab = "Brands",
                 ylab = "Mean Acceleration",
                 main = "Interaction Plot")
