# STAT3530 -- fitting and checking a model with several predictors
# The first section walks through the steps on the cheddar data. The second
# asks you to repeat them on a new dataset and answer a question. Run a line
# at a time.

site <- 'https://tdruiz-stat3530.share.connect.posit.cloud/_data/'

# ===========================================================================
# DEMONSTRATION: what makes cheddar taste like cheddar?
# ===========================================================================
# Taste scores from a panel, and concentrations of three compounds produced
# as the cheese ripens, for 30 mature cheddars.

load(url(paste0(site, 'cheddar.RData')))

head(cheddar)

# --- 1. look at the data ---------------------------------------------------
# Every pair at once: the response against each predictor, and the predictors
# against each other.
pairs(cheddar)

# The predictors are related to one another, so their coefficients will not
# match what you would get fitting them one at a time.
round(cor(cheddar), 2)

# --- 2. fit the model ------------------------------------------------------
fit <- lm(taste ~ Acetic + H2S + Lactic, data = cheddar)

summary(fit)

# --- 3. check the assumptions ----------------------------------------------
# Residuals versus fitted: want a shapeless band centred on zero.
plot(fitted(fit), resid(fit), xlab = 'fitted value', ylab = 'residual')
abline(h = 0, col = 'red', lty = 2, lwd = 2)

# Normal QQ: want points on the line.
qqnorm(resid(fit)); qqline(resid(fit), col = 'red', lty = 2, lwd = 2)

# One residual plot per predictor: a bend in any of these says that predictor
# is on the wrong scale. Run them one at a time.
plot(cheddar$Acetic, resid(fit), xlab = 'Acetic', ylab = 'residual')
abline(h = 0, col = 'red', lty = 2, lwd = 2)

plot(cheddar$H2S, resid(fit), xlab = 'H2S', ylab = 'residual')
abline(h = 0, col = 'red', lty = 2, lwd = 2)

plot(cheddar$Lactic, resid(fit), xlab = 'Lactic', ylab = 'residual')
abline(h = 0, col = 'red', lty = 2, lwd = 2)

# And the unusual observations.
plot(cooks.distance(fit), type = 'h', xlab = 'observation', ylab = "Cook's distance")

# --- 4. partial residual plots ---------------------------------------------
# A residual plot shows only what is left over. Adding back one predictor's
# fitted contribution shows the shape that predictor's effect should have:
#
#     e_j = e + beta-hat_j * x_j
#
# Again, one at a time.
pr <- resid(fit) + coef(fit)['Acetic'] * cheddar$Acetic
plot(cheddar$Acetic, pr, xlab = 'Acetic', ylab = 'partial residual')
abline(lm(pr ~ cheddar$Acetic), col = 'red', lty = 2, lwd = 2)

pr <- resid(fit) + coef(fit)['H2S'] * cheddar$H2S
plot(cheddar$H2S, pr, xlab = 'H2S', ylab = 'partial residual')
abline(lm(pr ~ cheddar$H2S), col = 'red', lty = 2, lwd = 2)

pr <- resid(fit) + coef(fit)['Lactic'] * cheddar$Lactic
plot(cheddar$Lactic, pr, xlab = 'Lactic', ylab = 'partial residual')
abline(lm(pr ~ cheddar$Lactic), col = 'red', lty = 2, lwd = 2)

# Straight lines here: nothing needs transforming.

# --- 5. answer the question ------------------------------------------------
confint(fit)

# H2S and lactic acid are associated with taste; acetic acid is not, once the
# other two are accounted for. Each unit of H2S is worth between 1.4 and 6.5
# points of mean taste score, holding the other compounds fixed.

# A prediction for a new cheese works exactly as it did with one predictor.
new <- data.frame(Acetic = 5.5, H2S = 6.0, Lactic = 1.4)
predict(fit, new, interval = 'confidence')     # for the mean at these values
predict(fit, new, interval = 'prediction')     # for one new cheese

# ===========================================================================
# YOUR TURN: what explains an occupation's prestige?
# ===========================================================================
# In 1971, a sample of Canadians rated 102 occupations on a prestige scale.
# The data pair those ratings with three facts about each occupation from the
# census: average education in years, average income in dollars, and the
# percentage of people in it who are women.
#
# GOAL: decide which of the three are associated with prestige, and say how
# much each one is worth, with an interval.

load(url(paste0(site, 'prestige.RData')))

head(prestige)

# --- 1. look at the data ---------------------------------------------------
# Use pairs() on the four numeric columns. What looks straight? What does not?



# --- 2. fit the model ------------------------------------------------------
# Regress prestige on education, income, and women.



# --- 3. check the assumptions ----------------------------------------------
# Residuals versus fitted, the QQ plot, one residual plot per predictor, and
# Cook's distance.



# --- 4. partial residual plots ---------------------------------------------
# Is any predictor on the wrong scale? If so, transform it and refit, then
# check the diagnostics again.



# --- 5. answer the question ------------------------------------------------
# Report confidence intervals from your final model and answer the goal above
# in a sentence or two. If you transformed a predictor, say what its
# coefficient means on the original scale.