# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Confidence Intervals & Hypothesis Testing                                 ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# This script shows how to calculate confidence intervals (CI) and perform
#mhypothesis testing (HT) for simple & multiple linear regression models.

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 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

n <- length(accel)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Confidence Intervals   ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# We will use the following model (the selected one from the previous lab).
m2 <- lm(Acceleration ~ Horsepower + Engine + Weight, data = cars)
sm2 <- summary(m2)

# The confidence interval for a parameter takes the form
# CI = [ bhat - t_{n - k - 1, a / 2} * s.e.(bhat),
#        bhat + t_{n - k - 1, a / 2} * s.e.(bhat) ]

bhat <- sm2$coefficients[, 1]
bhatse <- sm2$coefficients[, 2]

alpha <- 0.05
width <- qt(1 - alpha / 2, df = n - length(bhat), lower.tail = TRUE) * bhatse
lower <- bhat - width
upper <- bhat + width
CI <- cbind(lower, upper)
CI

# The confint function gives the lower and upper limits for each parameter for
# the 95% confidence interval (default).

confint(m2)
confint(m2, level = 0.95)

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Prediction Intervals   ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# For a given value of x we can calculate the predicted y and produce intervals
# for it. These intervals are either confidence or prediction intervals.

newcars <- data.frame(Horsepower = c(100, 155, 190),
                      Engine = c(3180, 3800, 4372),
                      Weight = c(1310, 1672, 1921))

predict(m2, newdata = newcars, level = 0.95, interval = "confidence")
predict(m2, newdata = newcars, level = 0.95, interval = "prediction")

# We can also compute the intervals manually

sigmahat <- sm2$sigma
X <- cbind(rep(1, nrow(m2$model)), as.matrix(m2$model[-1]))
xnew <- c(1, as.numeric(newcars[1, ]))
ynew <- xnew %*% bhat

# Confidence Interval
sigmaconf <- sigmahat * sqrt(t(xnew) %*% solve(t(X) %*% X) %*% xnew)
confwidth <- qt(1 - alpha / 2, df = n - ncol(X), lower.tail = TRUE) * sigmaconf
confCI <- c(ynew - confwidth, ynew + confwidth)
ynew
confCI

# Prediction Interval
sigmapred <- sigmahat * sqrt(1 + t(xnew) %*% solve(t(X) %*% X) %*% xnew)
predwidth <- qt(1 - alpha / 2, df = n - ncol(X), lower.tail = TRUE) * sigmapred
predCI <- c(ynew - predwidth, ynew + predwidth)
ynew
predCI

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Hypothesis Testing     ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~

# We will see tests of the following forms:
  # H0: bi = c
  # H1: bi != c or bi < c or bi > c

# If c = 0 the test is a test of statistical significance for bi

summary(m2)

# Nevertheless, we give here one more example, for a two-sided test this time.
  # H0: b3 = 0.005
  # H1: b3 != 0.005

c <- 0.005
tstat <- (bhat[4] - c) / bhatse[4]
tstat

# The test here is two-sided, so we need to compare our observed tstat value
# with t*, the value for the t distribution, with n - d degrees of
# freedom that satisfies:
# P(|T| > t*) = a <=> P(T > t* or T < - t*) = a
# where T ~ t_(n-d).

a <- 0.05
qt_a <- qt(1 - a / 2, n - ncol(X))
qt_a <- qt(a / 2, n - ncol(X), lower.tail = FALSE)
qt_a

# We notice that t_stat > qt_a, so we reject the null hypothesis H0 in favor of
# the alternative hypothesis H1 at significance level a = 0.05.

# An other way to decide whether we have strong evidence to reject the null
# hypothesis, or not is by calculating the p-value.

# p-value = P(|T|  >|t_stat|) = P(T > |t_stat| or T < - |t_stat|),
# where T ~ t_(n-d) (for a two-sided test).
# If p-value < a we reject H0.

p_value <- 2 * pt(abs(tstat), n - ncol(X), lower.tail = FALSE)
p_value

# p_value <- pt(abs(tstat), n - ncol(X), lower.tail = FALSE)
#            + pt(-abs(tstat), n - ncol(X), lower.tail = TRUE)
