# STAT3530 -- inference and diagnostics
# Part 1 shows inference from a fitted model. Part 2 checks model diagnostics.
# Each part has a demonstration followed by a short exercise. Run a line at a time.

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

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

# ===========================================================================
# PART 1: INFERENCE
# ===========================================================================
# The prevend data: RFFT score (executive cognitive function) and age, for
# 208 adults.

fit <- lm(rfft ~ age, data = prevend)

# --- estimates and standard errors -----------------------------------------
# coefficient table shows estimate, standard error, t statistic, p-value.
summary(fit)

# The t statistic is just the estimate divided by its standard error:
coef(fit)[2] / summary(fit)$coefficients[2, 2]

# and the p-value is the area beyond it in both tails of a t distribution
2 * pt(-abs(summary(fit)$coefficients[2, 3]), df = df.residual(fit))

# --- confidence intervals for the coefficients -----------------------------
confint(fit, level = 0.95)

# verify by hand
b  <- coef(fit)[2]
se <- summary(fit)$coefficients[2, 2]
b + c(-1, 1) * qt(0.975, df = df.residual(fit)) * se

# --- prediction at a new value ---------------------------------------------
# point estimate
predict(fit, newdata = data.frame(age = 55))

# confidence interval for the MEAN score among 55 year olds
predict(fit, newdata = data.frame(age = 55), interval = 'confidence')

# prediction interval for ONE 55 year old's score
predict(fit, newdata = data.frame(age = 55), interval = 'prediction')

# Note how much wider the second one is, and why: the first covers where the
# line is, the second also has to cover where a person sits around the line.

# --- your turn -------------------------------------------------------------
# The cats data record body weight (kg) and heart weight (g) for 144 adult
# cats. You met them on the homework. Model heart weight using body weight.
# Data are in the MASS package, which comes with R.

data(cats, package = 'MASS')

# 1. Fit the model and inspect the coefficient summary.
plot(Hwt ~ Bwt, data = cats)
fit <- lm(Hwt ~ Bwt, data = cats)
abline(fit)
summary(fit)

# 2. Is there evidence that heart weight is associated with body weight?
#    (No code necessary for this part -- discuss with a neighbor.)


# 3. Give a 95% confidence interval for the slope.
confint(fit, level = 0.95)

# 4. Predict heart weight for a cat weighing 2.5 kg.
predict(fit, newdata = data.frame(Bwt = c(2.5, 3.0)), interval = "confidence")

# ===========================================================================
# PART 2: DIAGNOSTICS
# ===========================================================================
# Check the model assumptions for the RFFT data.

# --- 1. the data with the fitted line --------------------------------------
# Always look at this first. A curve here is worth more than any test.
plot(rfft ~ age, data = prevend)
abline(fit, col = 'blue', lwd = 2)

# --- 2. residuals against fitted values ------------------------------------
# Checks the mean (linearity) and the spread (constant variance).
# Want: a horizontal band centered around zero
plot(fitted(fit), resid(fit), xlab = 'fitted value', ylab = 'residual')
abline(h = 0, col = 'red', lty = 2, lwd = 2)

# --- 3. the normal QQ plot -------------------------------------------------
# Checks normality. Want: points on the line.
# Read the ends. The middle almost always looks fine.
qqnorm(resid(fit))
qqline(resid(fit), col = 'red', lty = 2, lwd = 2)

# Independence is the one assumption these plots cannot check. The prevend
# rows are one per person, in no meaningful order, so there is nothing to
# plot against. You have to argue for it from how the data were collected.

# Verdict for prevend: nothing here would stop us reporting Part 1.

# --- your turn -------------------------------------------------------------
# Now check the model you fitted in Part 1.

# 5. Check the diagnostic plots for your cats model. 
#    Do any of the model assumptions seem questionable? Consider why or why not.

# rvf plot
plot(fitted(fit), resid(fit))
abline(h = 0, lwd = 2, col = 'red')
lm(resid(fit) ~ fitted(fit)) |> abline(col = 'blue')

# qq
qqnorm(resid(fit))
qqline(resid(fit), col = 'red', lwd = 2)