#install.packages("MASS")
library(MASS)
data(whiteside)

plot(whiteside$Temp, whiteside$Gas,
     pch  = 16, col = "steelblue", cex = 0.9,
     xlab = "Temperature (°C)",
     ylab = "Gas consumption (m^3)",
     main = "Gas consumption vs Temperature")

# Model with temperature

mod1 <- lm(Gas ~ Temp, data = whiteside)

plot(whiteside$Temp, whiteside$Gas,
     pch  = 16, col = "steelblue", cex = 0.9,
     xlab = "Temperature (°C)",
     ylab = "Gas consumption (m^3)",
     main = "Gas consumption vs Temperature")
abline(coef(mod1), col = "red")

summary(mod1)

plot(mod1, which = 1)


plot(whiteside$Temp, whiteside$Gas,
     pch  = c(16, 17)[as.numeric(whiteside$Insul)],
     col  = c("steelblue", "darkorange")[as.numeric(whiteside$Insul)],
     xlab = "Temperature (°C)", ylab = "Gas consumption (m³)",
     main = "Gas vs Temperature — coloured by Insulation period")
legend("topright", legend = c("Before", "After"), pch    = c(16, 17),
       col    = c("steelblue", "darkorange"), bty    = "n")


# Including Insul 
mod2 <- lm(Gas ~ Temp + Insul, data = whiteside) 
summary(mod2)


b0 <- coef(mod2)["(Intercept)"]
b1 <- coef(mod2)["Temp"]
b2 <- coef(mod2)["InsulAfter"]

plot(whiteside$Temp, whiteside$Gas,
     pch  = c(16, 17)[as.numeric(whiteside$Insul)],
     col  = c("steelblue", "darkorange")[as.numeric(whiteside$Insul)],
     xlab = "Temperature (°C)",
     ylab = "Gas consumption (m³)",
     main = "Multiple regression: two parallel lines")

abline(b0,        b1, col = "steelblue",   lwd = 2)   # Before
abline(b0 + b2,   b1, col = "darkorange",  lwd = 2)   # After

legend("topright",
       legend = c("Before", "After"),
       col    = c("steelblue", "darkorange"),
       pch    = c(16, 17), lwd = 2, bty = "n")


plot(mod2, which = 1)


predict(mod2, newdata = data.frame(Temp = 10, Insul = "After"))-
  predict(mod2, newdata = data.frame(Temp = 5, Insul = "Before"))


# Interaction 
mod3 <- lm(Gas ~ Temp*Insul, data = whiteside) 
summary(mod3)


b0 <- coef(mod3)["(Intercept)"]
b1 <- coef(mod3)["Temp"]
b2 <- coef(mod3)["InsulAfter"]
b3 <- coef(mod3)["Temp:InsulAfte"]

plot(whiteside$Temp, whiteside$Gas,
     pch  = c(16, 17)[as.numeric(whiteside$Insul)],
     col  = c("steelblue", "darkorange")[as.numeric(whiteside$Insul)],
     xlab = "Temperature (°C)",
     ylab = "Gas consumption (m³)",
     main = "Multiple regression: interaction")

abline(b0,        b1, col = "steelblue",   lwd = 2)   # Before
abline(b0 + b2,   b1 + b3, col = "darkorange",  lwd = 2)   # After

legend("topright",
       legend = c("Before", "After"),
       col    = c("steelblue", "darkorange"),
       pch    = c(16, 17), lwd = 2, bty = "n")


predict(mod3, newdata = data.frame(Temp = 2, Insul = "After")) -
  predict(mod3, newdata = data.frame(Temp = 2, Insul = "Before"))


predict(mod3, newdata = data.frame(Temp = 10, Insul = "After")) -
  predict(mod3, newdata = data.frame(Temp = 10, Insul = "Before"))


predict(mod3, newdata = data.frame(Temp = 3, Insul = "After")) -
  predict(mod3, newdata = data.frame(Temp = 2, Insul = "After"))


predict(mod3, newdata = data.frame(Temp = 3, Insul = "Before")) -
  predict(mod3, newdata = data.frame(Temp = 2, Insul = "Before"))


plot(mod3, which = 1)