options(repr.plot.width=3, repr.plot.height=3)

plot <- function(...) {
    par(mar=c(4,4,0.1,0.1))
    graphics::plot(...)
}


library(faraway)
data(stat500)

head(stat500)

plot(final ~ midterm, data=stat500)

model_a1 <- lm(final ~ 1 + midterm, data=stat500)
print(model_a1)

plot(final ~ midterm, data=stat500)
abline(model_a1$coefficients[1], model_a1$coefficients[2])

head( model.matrix(~ 1 + midterm, data=stat500) )

head( model.matrix(~ 1, data=stat500) )

model_a0 <- lm(final ~ 1, data=stat500)
model_a0

print( anova(model_a0, model_a1) )

cor.test(stat500$final, stat500$midterm)

summary(model_a1)

stat500$midterm_squared <- stat500$midterm ^ 2

model_a2 <- lm(final ~ 1 + midterm + midterm_squared, data=stat500)
model_a2

print( anova(model_a1, model_a2) )

library(faraway)
data(coagulation)
# For brevity's sake, I'm going to use a subset of the data
coagulation <- coagulation[c(1,2,5,6,11,12,17,18,19),]

coagulation
plot(coag ~ diet, data=coagulation)

print(coagulation$diet)

model.matrix(~ diet, data=coagulation)

model.matrix(~ 0 + diet, data=coagulation)

model_b0 <- lm(coag ~ 1, data=coagulation)
model_b1 <- lm(coag ~ diet, data=coagulation)

print( model_b0 )
print( model_b1 )

print( anova(model_b0, model_b1) )

library(MASS)
coefficient_estimator <- ginv(model.matrix(~ diet, data=coagulation))
round(coefficient_estimator, 3)

library(faraway)
data(pvc)

head(pvc)
summary(pvc)

model_c_interaction <- lm(psize ~ operator*resin, data=pvc)
# Shorthand for
#   psize ~ 1 + operator + resin + operator:resin

model_c_operator_resin <- lm(psize ~ operator + resin, data=pvc)
model_c_operator <- lm(psize ~ operator, data=pvc)
model_c_resin <- lm(psize ~ resin, data=pvc)
model_c_null <- lm(psize ~ 1, data=pvc)

model_c_interaction

print( anova( model_c_operator_resin, model_c_interaction) )

print( anova( model_c_operator, model_c_operator_resin) )

print( anova( model_c_resin, model_c_operator_resin) )