# STAT3530 -- diagnose and fix
# Twelve simulated datasets, each with one predictor. For some of them the
# simple linear model is fine as it is; for others it is not. For each
# dataset: diagnose, try a fix, and check that the fix worked. Run a line at
# a time.

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

load(url(paste0(site, 'diagnose.RData')))    # loads d1, d2, ..., d12

# ===========================================================================
# THE TOOLKIT
# ===========================================================================
# The same loop works for every dataset. Here it is on d1.

# --- 1. look at the data ---------------------------------------------------
plot(y ~ x, data = d1)

# --- 2. fit the simple linear model and check it ---------------------------
fit <- lm(y ~ x, data = d1)
abline(fit, col = 'blue', lwd = 2)

# residuals vs fitted: want a horizontal 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)

# and read the coefficient table: is there a relationship at all, and does
# the intercept mean anything?
summary(fit)

# --- 3. only for data recorded in time order -------------------------------
# each residual against the one before it: want a shapeless blob.
# d1 is not recorded in time order, so this plot means nothing for d1 --
# it is here only to show the code.
e <- resid(fit)
plot(head(e, -1), tail(e, -1), xlab = 'residual i', ylab = 'residual i + 1')
abline(h = 0, v = 0, col = 'grey')

# --- 4. try a fix, then go back to step 2 ----------------------------------
# A fix can change the response, the predictor, or the model. Write the change
# straight into the formula, for example:
#
#   lm(log(y) ~ x, data = d1)
#   lm(y ~ I(x - 5), data = d1)
#
# Careful: arithmetic on the RIGHT of ~ must be wrapped in I(). Without it,
# lm(y ~ x^2) runs without complaint but quietly fits the plain line.
#
# When you change the response, redraw the scatterplot on the new scale too,
# e.g. plot(log(y) ~ x, data = d1).

# --- comparing candidate fits ----------------------------------------------
# Two fits with the SAME response can be compared by their residual standard
# deviation: smaller is better.
summary(fit)$sigma

# This does not work once the response changes. y and log(y) are on different
# scales, so their residual standard deviations cannot be compared. To choose
# between transformations of y, compare the diagnostic plots instead.

# --- if you get stuck ------------------------------------------------------
# - The log is one member of a family of transformations of the response:
#   powers y^p. The further p falls below 1, the harder large values are
#   pulled in -- powers between 0 and 1 (say y^0.25 or y^0.75) are gentler
#   than the log, and negative powers (say y^-0.5 or y^-2) are harsher.
#   If the log overcorrects, so that the curve now bends the other way, try a
#   power between 0 and 1. If it undercorrects, so that the curve still bends
#   the same way, try a negative power. Negative powers reverse the order of
#   y, so expect the slope to change sign.
#
# - Logs and powers cannot straighten a relationship that changes direction,
#   rising and then falling or falling and then rising. For those, think about
#   what happens to x when you centre it and then square it.

# ===========================================================================
# YOUR TURN
# ===========================================================================
# For each dataset, note what is wrong (or "nothing"), which plot or output
# told you, what you tried, and whether the diagnostics were clean afterwards.
#
# d7 has one value per year, recorded in time order. Its predictor is called
# year rather than x.
#
# Finish d1 first, since you have already started it.

# --- d1 --------------------------------------------------------------------


# --- d2 --------------------------------------------------------------------


# --- d3 --------------------------------------------------------------------


# --- d4 --------------------------------------------------------------------


# --- d5 --------------------------------------------------------------------


# --- d6 --------------------------------------------------------------------


# --- d7 --------------------------------------------------------------------


# --- d8 --------------------------------------------------------------------


# --- d9 --------------------------------------------------------------------


# --- d10 -------------------------------------------------------------------


# --- d11 -------------------------------------------------------------------


# --- d12 -------------------------------------------------------------------


# ===========================================================================
# AFTERWARDS
# ===========================================================================
# 1. Which datasets could you resolve, and how?
#
# 2. Which are you unsure of, and what did you try?
#
# 3. Which could you not resolve?