# STAT3530 -- categorical predictors and interactions

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

# three simulated datasets, two groups each (a and b), n = 60
load(url(paste0(site, 'specpanels.RData')))

# ===========================================================================
# 1. NAME THE SPECIFICATION
# ===========================================================================
# Three candidate models for a continuous predictor x and a two-level group g:
#
#     y ~ x        one line for both groups
#     y ~ x + g    parallel lines at different levels
#     y ~ x * g    different levels AND different slopes
#
# For each panel: plot it, decide which model the plot calls for, fit that
# model, and then test the claim you read off the plot.

## example: panel 1

# group a in grey, group b in red
plot(y ~ x, data = panel1, pch = 16, col = c('grey40', 'red')[panel1$g])

# the slopes plainly differ, so specify different slopes
fit.int <- lm(y ~ x * g, data = panel1)

# intercept and slope for group a, then the two differences
coef(fit.int)

# add a line for each group
abline(a = ..., b = ..., col = 'grey40', lwd = 2)
abline(a = ..., b = ..., col = 'red', lwd = 2)

# now test what you read off the plot. H0 is a common slope, which is the
# model without the product term
fit.add <- lm(y ~ x + g, data = panel1)
anova(fit.add, fit.int)

## your turn: panel 2

# plot it, then write down the model the plot calls for
plot(...)

# fit that model


# test the slope difference against a common slope
anova(..., ...)

# the test does not back up what you read off the plot. what do you fit
# instead, and is the conclusion that the slopes are equal?


## your turn: panel 3

# plot, specify, fit, test



# before you interpret anything, check where each group's x values are
tapply(panel3$x, panel3$g, range)

# the coefficient on g is the gap between groups at x = 0. what is wrong with
# reporting it here, and what would you do to the predictor to fix it?


# ===========================================================================
# 2. REAL DATA: PENGUINS
# ===========================================================================
# 333 penguins from Palmer Station, Antarctica. Body mass in grams, flipper
# and bill dimensions in millimeters, three species, and sex.
load(url(paste0(site, 'penguins.RData')))

cols <- c('grey40', 'red', 'blue')
plot(body_mass ~ flipper_len, data = penguins, pch = 16,
     col = cols[penguins$species])
legend('topleft', legend = levels(penguins$species), pch = 16, col = cols)

## example: does the flipper-mass relationship depend on species?

# the Gentoo points climb more steeply, so specify different slopes
fit.int <- lm(body_mass ~ flipper_len * species, data = penguins)

# slope for the reference species, then the two differences in slope
coef(fit.int)

# three levels means two product terms, so the slope difference costs 2 df and
# needs an F test rather than a look at the coefficient table
fit.add <- lm(body_mass ~ flipper_len + species, data = penguins)
anova(fit.add, fit.int)

## your turn: does the relationship depend on sex instead?



## your turn: fit bill_dep ~ bill_len, then add species. the slope changes
## sign. what was the marginal fit measuring?



## your turn: Adelie is the reference level because it sorts first. which of
## these change when you make Gentoo the reference -- the coefficients, the
## fitted values, R^2, the F test for the slope difference?

levels(penguins$species)
penguins$species2 <- relevel(penguins$species, ref = 'Gentoo')



# ===========================================================================
# HOW THE PANELS WERE BUILT -- read after the debrief
# ===========================================================================
# All three have n = 60, 30 per group, and the same noise level. All three
# look like the slopes differ. Only panel 1's do.
#
#   panel 1   y ~ x * g. b starts 8 higher with slope 0.5 - 1.1 = -0.6, so the
#             lines cross mid-range. Because they cross, the coefficient on g
#             describes the gap at x = 0 only, and there is no single "group
#             effect" to report.
#   panel 2   y ~ x + g. Slopes 1.0 and 1.3. The difference is real and n = 60
#             cannot resolve it. Failing to reject is not equality.
#   panel 3   y ~ x. IDENTICAL truth for both groups, y = 5 + x. Group b is
#             observed only on x in [6.5, 10], so its slope comes from a third
#             of the range and its intercept is a long extrapolation back to
#             x = 0. The apparent slope difference is noise.