# from http://www.maths.bath.ac.uk/~jjf23/LMR/scripts2/prediction.R
# with some comments added by Charlotte Wickham
library(faraway)
data(fat,package="faraway")
lmod <- lm(brozek ~ age + weight + height + neck + chest + abdom + hip + thigh + knee + ankle + biceps + forearm + wrist, data = fat)
x <- model.matrix(lmod)
# prediction for an "average" man
(x0 <- apply(x, 2, median)) #
(y0 <- sum(x0 * coef(lmod))) # x_0^T \hat{\beta}
# prediction intervals for an "average" man
predict(lmod, new = data.frame(t(x0)), se = TRUE)
# se.fit is the standard error in the predicted mean repsonse.
# residual.scale is estimate of sigma
# === Q1 ====
# **Find se.fit using matrix algebra and residual.scale**
# ===========
# Let R do the work of finding a prediction interval
predict(lmod, new = data.frame(t(x0)), interval = "prediction")
# === Q2 ====
# **Using the values from Q1 replicate the calculation of this interval**
# ===========
# Let R do the work of finding the confidence interval on the mean response
predict(lmod, new = data.frame(t(x0)), interval = "confidence")
# === Q3 ====
# **Replicate the calculation of this interval, too.**
# ===========
# Intervals are wider the further we are from the average
# explanatory values.
# Let's look at a prediction for someone who is at the
# (sample) 95th percentile
(x1 <- apply(x, 2, quantile, probs = 0.95))
predict(lmod, new = data.frame(t(x1)), interval = "prediction")
predict(lmod, new = data.frame(t(x1)), interval = "confidence")