# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Lab 4: Multiple Linear Regression                                         ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# This lab includes the basic computations for estimating the parameters of a
# Multiple Linear Regression Model.

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 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
horse <- cars$Horsepower
engine <- cars$Engine
weight <- cars$Weight
price <- cars$Price

# Visualization - Scatter plots
pairs(cbind(accel, horse, engine, weight, price))

# Which of the scatter plots seem linear?
# We want accel to be our y, so we only look the first column.

# accel ~ horse
# accel ~ engine
# accel ~ weight
# accel ~ price

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# The lm function        ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Use lm for the linear model
reg <- lm(accel ~ horse + engine + weight + price)
summary(reg)

# Alternative way to specify the same model
head(cars)
reg <- lm(Acceleration ~ Horsepower + Engine + Weight + Price, data = cars)
summary(reg)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Model Selection        ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Null model
m0 <- lm(Acceleration ~ 1, data = cars)
summary(m0)

# Full model
m1 <- lm(Acceleration ~ Horsepower + Engine + Weight + Price, data = cars)
summary(m1)

# Alternative way to create the full model (the . sign includes all variables,
# the - sign removes a variable).
m1 <- lm(Acceleration ~ . - Year - Fuel - Brand - Type, data = cars)
summary(m1)

# We want to find the "best" model possible. Based on what criterion do we
# choose the "best" model?

# Common Criteria:
# 1. Significance Level
#    All coefficients must be statistically significant.
# 2. MSE (Mean Square Error) or MAE (Mean Absolute Error)
#    Minimization of errors
# 3. AIC (Akaike Information Criterion) or BIC (Bayesian Information Criterion)
#    Information criteria that estimate the distance between the model and the
#    true distribution (Kullback-Leibler Divergence).

# If we have too many models to compare them all, how can we proceed?

# Common Model Selection Procedures:
# 1. Forward Selection
#    Start with simple linear models. Choose the best based on the criterion
#    of choice. Add the second best variable to the model. If the new, multiple
#    model is better than the simple one, proceed by adding another variable.
# 2. Backward Elimination
#    Start with the full linear model (containing all independent variables).
#    Remove one variable based on the criterion of choice. If the new model is
#    superior, proceed by removing another variable.

# Step-wise
# At each step, recheck whether a variable previously added should be removed.
# At each step, recheck whether a variable previously removed should be added.

# Example 1
# - Criterion: Significance Level
# - Strategy: Backward Elimination

summary(m1)

# We notice that the coefficient for price is non-significant (p-value = 0.86).
# We need to remove the variable and repeat.

m2 <- lm(Acceleration ~ Horsepower + Engine + Weight, data = cars)
summary(m2)

# All coefficients are statistically significant. The procedure ends here.
# m2 is the model selected.

# Instead of performing these procedures manually, R offers a set of functions
# that automatically detect the best model according to the criterion and
# strategy specified. Function step() can perform all three aforementioned
# procedures, using the AIC as criterion.

# Example 2
# - Criterion: AIC
# - Strategy: Stepwise Selection starting from the null model.

step(m0,
     scope = list(lower = formula(m0), upper = formula(m1)),
     direction = "both")

# Example 3
# - Criterion: adjusted R squared
# - Strategy: Backward Selection

library(StepReg)

stepwise(formula = formula(m1),
         data = cars,
         strategy = "backward",
         metric = "adjRsq")

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Residuals Diagnostics  ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Suppose model selection highlighted m2 as the best model.
m2 <- lm(Acceleration ~ Horsepower + Engine + Weight, data = cars)
summary(m2)

# Call the attributes
coef <- m2$coefficients
fitted <- m2$fitted.values
res <- m2$residuals

# We need to check the model assumptions. The assumptions are:
  # 1. Independency
  # 2. Homoskedasticity
  # 3. Normality

# 2. Homoskedasticity Assumption

# Plot: Residuals - Fitted Values
plot(x = fitted, y = res,
     main = "Residuals - Fitted Values",
     xlab = "Fitted Values", ylab = "Residuals",
     pch = 20)
lines(x = c(0, length(fitted)), y = c(0, 0), col = "red")

library(lmtest)

# Breusch-Pagan test
  # H0: Homoscedasticity
  # H1: Heteroscedasticity
bptest(m2)

# 3. Normality Assumption

# Plots: Histogram, Boxplot, Q-Q Plot
par(mfrow = c(1, 3))
hist(res, main = "Histogram of Residuals")
boxplot(res, main = "Box Plot of Residuals")
qqnorm(res, main = "Q-Q Plot of Residuals")
qqline(ehat)

# Test: Shapiro-Wilk
  # H0: Normal Distribution
  # H1: Not Normal Distribution
shapiro.test(res)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Model Transformation   ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Let us revisit the initial scatter-plots

# Visualization - Scatter plots
pairs(cbind(accel, horse, engine, weight))
pairs(cars[,2:5])

# None of the three plots seem "ideal".
# Could we transform one of the x?

lhorse <- log(horse)
cars$logHorsepower <- log(cars$Horsepower)

# Take a look at the transformation
par(mfrow = c(1, 2))

plot(horse, accel,
     main = "Scatter Plot",
     xlab = "Horsepower",
     ylab = "Acceleration")

plot(lhorse, accel,
     main = "Scatter Plot",
     xlab = "log of Horsepower",
     ylab = "Acceleration")

# Repeat the procedure
m3 <- lm(Acceleration ~ logHorsepower + Engine + Weight, data = cars)
summary(m3)
resnew <- m3$residuals

# Repeat the heteroscedasticity test
bptest(m3)

# Repeat the normality test
shapiro.test(resnew)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Parameter Estimation   ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# If we did not have lm(), could we code everything we learned ourselves?

# We will estimate the coefficients using matrices. You will notice that the
# method is identical in the simple and the multiple model.

y <- accel
n <- length(y)

# Create Design Matrix
ones <- rep(1, n)
X <- cbind(ones, horse, engine, weight, price)
X

# Coefficients
betas <- solve(t(X) %*% X) %*% t(X) %*% y
betas

# Estimated y
yhat <- X %*% betas

# Residuals
ehat <- y - X %*% betas
ehat

# Variance-Standard Deviation
sigma2_hat <- t(ehat) %*% ehat / (n - ncol(X))
sigma_hat <- sqrt(sigma2_hat)

var_betas <- as.vector(sigma2_hat) * solve(t(X) %*% X)
std_betas <- sqrt(diag(var_betas))
var_betas
std_betas

# Check if a coefficient is statistically significant.

# H0: bi = 0
# H1: bi != 0

# Under Ho, bi.hat ~ N( 0, var(bi) )
# We will use the student t distribution to
# find the p-value (pt function)

t <- -abs(betas / std_betas)
pvalues <- 2 * pt(t, n - ncol(X))

mat <- cbind(t, pvalues)
colnames(mat) <- c("t Statistic", "p-value")
mat
