# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Lab 3: Simple Linear Regression                                           ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# This script includes the parameter estimation for a simple linear 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
y <- cars$Acceleration

# Independent Variable
x <- cars$Horsepower

# Visualization: Scatter plot
plot(x, y,
     main = "Scatter Plot",
     xlab = "Horsepower",
     ylab = "Acceleration")

# Pearson Correlation Coefficient
cor(x, y)

# Hypothesis test, assuming x and y follow Normal distr.
# H0: p = 0
# H1: p != 0
cor.test(x, y)

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

# lm is a built in R-function
reg <- lm(y ~ x)
summary(reg)
anova(reg)

# Call the attributes
bhat <- reg$coefficients
ehat <- reg$residuals
yhat <- reg$fitted.values
stdehat <- rstandard(reg)

# Plot the scatterplot and add the regression line
plot(x, y,
     main = "Scatterplot with Fitted Regression Line",
     xlab = "Horsepower",
     ylab = "Acceleration")
lines(x, yhat)

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

# e_i ~ N(0, s ^ 2)
# The assumptions are:
# Homoskedasticity
# Independency
# Normality

# Homoskedasticity & Independency will be checked
# only graphically with a residuals plot
plot(x = yhat,
     y = stdehat,
     main = "Standardized Residuals - Fitted Values",
     xlab = "Fitted Values", ylab = "Standardized Residuals",
     pch = 20)

lines(x = c(0, length(x)), y = c(0, 0), col = "red")

# Normality Check

# Plots: Histogram, Boxplot, Q-Q Plot

# Histogram
hist(stdehat, main = "Histogram of Standardized Residuals")

# Box plot
boxplot(stdehat, main = "Box Plot of Standardized Residuals")

# Q-Q Plot
qqnorm(stdehat, main = "Q-Q Plot of Standardized Residuals")
qqline(stdehat)

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

# Automatically create important diagnostic plots
plot(reg)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Manual Par. Estimation ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

# We will estimate the model coefficients (b0, b1) using OLS (Ordinary Least
# Squares)

n <- length(y)

# Coefficients Estimation
t <- sum((x - mean(x)) ^ 2)
b1 <- sum((x - mean(x)) * (y - mean(y))) / t
b0 <- mean(y) - b1 * mean(x)

# Estimated y
yhat <- b0 + b1 * x

# Residuals (kataloipa) of the model
ehat <- y - yhat

# Variance-Standard Deviation
sigmahat <- sum(ehat ^ 2) / (n - 2)

# Variance-Standard Deviation of the intercept
sigma_b0 <- sigmahat * (1 / n + mean(x) ^ 2 / t)
std_b0 <- sqrt(sigma_b0)

# Variance-Standard Deviation of the slope
sigma_b1 <- sigmahat / t
std_b1 <- sqrt(sigma_b1)

# Regression Plot
plot(x, y,
     main = "Scatter Plot with Fitted Regression Line",
     xlab = "Horsepower",
     ylab = "Acceleration")

lines(x, yhat)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Matrix Par. Estimation ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Could we do this better, using matrices?

# Create Design Matrix
ones <- rep(1, n)
X <- cbind(ones, x)
X
colnames(X) <- c("b0", "b1")
X

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

# Estimated y
mat_yhat <- X %*% betas

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

#Check that the two methods give the same results
all.equal(b0, betas[1, ], check.names = FALSE)
all.equal(b1, betas[2, ], check.names = FALSE)
all.equal(ehat, as.vector(mat_ehat))
all.equal(yhat, as.vector(mat_yhat))

# Variance-Standard Deviation
sigma2_hat <- as.vector(t(mat_ehat) %*% mat_ehat / (n - 2))
sigma_hat <- sqrt(sigma2_hat)

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

# Standardized residuals
H <- X %*% solve(t(X) %*% X) %*% t(X)
stdehat <- ehat / (sigma_hat * sqrt(1 - diag(H)))
