# STAT3530 -- centering the predictor
# Three datasets, the same exercise: fit a simple linear model with the
# predictor as given, then again with the predictor centered, and compare the
# estimates and their estimated variances. Run a line at a time.

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

# --- example 1: ozone and temperature --------------------------------------
load(url(paste0(site, 'ozone.RData')))
str(ozone)

# daily ozone concentration (ppm) and temperature (F) in the Los Angeles basin
plot(O3 ~ temp, data = ozone,
     xlab = 'temperature (F)', ylab = 'ozone (ppm)')

fit <- lm(O3 ~ temp, data = ozone)
coef(fit)
abline(fit, col = 'blue')

# the intercept is the fitted ozone at 0 degrees. Two problems with that:
range(ozone$temp)     # 0 is 25 degrees colder than anything observed
coef(fit)[1]          # and it predicts a negative concentration

# --- centering with scale() ------------------------------------------------
# scale() standardizes (centers and scales) data
scale(ozone$temp)[1:5]                  # centered and divided by sd

# to center, add an argument
scale(ozone$temp, scale = FALSE)[1:5]   # centered only

# same as...
head(ozone$temp - mean(ozone$temp), 5)

# the mean it subtracted is kept as an attribute
attr(scale(ozone$temp, scale = FALSE), 'scaled:center')
mean(ozone$temp)

# refit
fit.c <- lm(O3 ~ scale(temp, scale = F), data = ozone)

# the slope is identical; only the intercept moved
rbind(uncentered = coef(fit), centered = coef(fit.c))

# and the new intercept is exactly the mean response
mean(ozone$O3)

# --- the estimated variance-covariance matrix ------------------------------
# vcov() returns sigma-hat^2 (X'X)^{-1}: variances on the diagonal, the
# covariance between the two estimates off it.
vcov(fit)
vcov(fit.c)

# read off the pieces
sqrt(diag(vcov(fit)))     # standard errors, same as summary(fit)
sqrt(diag(vcov(fit.c)))   # the intercept is estimated far more precisely

# easier to read as a correlation
cov2cor(vcov(fit))
cov2cor(vcov(fit.c))

# check it against the formula by hand
X <- model.matrix(fit.c)
t(X) %*% X                          # diagonal, because 1'x = 0 when x is centered
sigma(fit.c)^2 * solve(t(X) %*% X)  # matches vcov(fit.c)

# --- example 2: how long mammals sleep -------------------------------------
load(url(paste0(site, 'mammals.RData')))
str(mammals)

# sleep is hours per day, lifespan is in years, for 51 species
plot(sleep ~ lifespan, data = mammals,
     xlab = 'lifespan (years)', ylab = 'sleep (hours/day)')

# fit as usual (no centering)
fit2 <- lm(sleep ~ lifespan, data = mammals)
abline(fit2, col = 'blue')

# your turn: refit after centering and compare estimates
# does the centering help with interpretations?

# --- example 3: endemic species in the Galapagos ---------------------------
load(url(paste0(site, 'endemics.RData')))
str(endemics)

# proportion of each island's species that are endemic
endemics$prop <- endemics$endemic / endemics$species
endemics$log.area <- log(endemics$area)

# range of areas; is centering on the log scale needed?
range(endemics$area)

# plot the data, fit a simple linear regression model
plot(prop ~ log.area, data = endemics,
     xlab = 'log area (km^2)', ylab = 'proportion endemic')
fit3 <- lm(prop ~ log.area, data = endemics)
abline(fit3, col = 'blue')

# your turn: refit after centering and compare estimates
# does the centering help with interpretations?