# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# One sample and Independent Samples t test                                 ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# This script includes the basic analysis for one sample and independent samples
# t tests.

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 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-Benz')
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)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# One Sample t test      ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Question: Assuming acceleration is a random variable following a Normal
# distribution with unknown expectation mu and std. deviation s, test whether
# the expected acceleration (mu) is equal to 15.

# First, we need to check the normality assumption for acceleration.

# Shapiro-Wilk test
# H0: acceleration is normally distributed.
# H1: acceleration is NOT normally distributed.

shapiro.test(accel)

# Result
# p-value = 0.2599, we do not reject H0 in the a = 0.05 level.
# If H0 was rejected, we should use other (usually non-parametric) tests, such
# as the Wilcoxon signed rank test.

# Option 1: Use the t.test function.

t.test(accel, mu = 15)

# The default test is two-sided, but we can choose less or greater for the H1.
# mu is the value of H0 for which we test.
# The test also produces a confidence interval (95% by default).

# Option 2: Conduct a one-sample t test manually.

m <- mean(accel) # 15.62118
s <- sd(accel)
m0 <- 15
t <- sqrt(n) * (m - m0) / s
t

# Test
#   H0: mean = 15
#   H1: mean != 15

2 * pt(abs(t), df = n - 1, lower.tail = FALSE) # P(T < -|t|) + P(T > |t|)
t.test(accel, mu = 15)

# Test
#   H0: mean = 15
#   H1: mean > 15

pt(t, df = n - 1, lower.tail = FALSE) # P(T > t)
t.test(accel, mu = 15, alternative = "greater")

# Test
#   H0: mean = 15
#   H1: mean < 15

pt(t, df = n - 1, lower.tail = TRUE) # P(T < t)
t.test(accel, mu = 15, alternative = "less", conf.level = 0.90)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Ind. Samples t test    ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Question: Assuming acceleration of the std and the turbo cars are two random
# variables following Normal distributions with unknown expectations mu1, mu2,
# and std. deviation s, test whether the two types of cars have the same
# expected acceleration (mu1 = mu2).

# We want to compare the acceleration between two independent groups of cars,
# divided by their Type (std vs turbo).

# 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"]

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

# We need to check the normality assumption for both categories.

shapiro.test(accel_std)
shapiro.test(accel_turbo)

# Error Bar:
  # The dot represents the mean of each category.
  # The blue lines are confidence intervals (95% by default).

# 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")

# Question: Doesn't this look like a boxplot? What's their difference?
  # What does the center point / line indicate in each graph?
  # What do the interval limits indicate in each graph?

# Option 1: Use the t.test function.

t.test(accel_std, accel_turbo, var.equal = TRUE)

# Attention: We assume equal variances.

# Option 2: Conduct a one-sample t test manually.

m0 <- mean(accel_std)
m1 <- mean(accel_turbo)
var0 <- var(accel_std)
var1 <- var(accel_turbo)
var <- ((n0 - 1) * var0 + (n1 - 1) * var1) / (n0 + n1 - 2)
t <- (m0 - m1) / sqrt(var * (1 / n0 + 1 / n1))
t

# Test
#   H0: m0 = m1         (m0 is for std and m1 is for turbo)
#   H1: m0 != m1

2 * pt(abs(t), df = n - 2, lower.tail = FALSE) # P(T < -|t|) + P(T > |t|)
t.test(accel_std, accel_turbo, var.equal = TRUE)

# Test
#   H0: m0 = m1
#   H1: m0 > m1

pt(t, df = n - 2, lower.tail = FALSE) # P(T > t)
t.test(accel_std, accel_turbo, alternative = "greater", var.equal = TRUE)

# Test
#   H0: m0 = m1
#   H1: m0 < m1

pt(t, df = n - 2, lower.tail = TRUE) # P(T < t)
t.test(accel_std, accel_turbo, alternative = "less",
       conf.level = 0.90, var.equal = TRUE)
