# STAT3530 -- stepwise selection

# install.packages('StepReg')
library(StepReg)

# load data
site <- 'https://tdruiz-stat3530.share.connect.posit.cloud/_data/'
load(url(paste0(site, 'pollution.RData')))

# log scale for the three pollution variables
pollution$lso2 <- log(pollution$so2)
pollution$lhc  <- log(pollution$hc)
pollution$lnox <- log(pollution$nox)
poll <- pollution[, setdiff(names(pollution), c('so2', 'hc', 'nox'))]

# 15 predictors, so 2^15 possible models
ncol(poll) - 1
2^(ncol(poll) - 1)

# ===========================================================================
# warm-up: forward selection by hand
# ===========================================================================

# fit an intercept only model
fit.null <- lm(mortality ~ 1, data = poll)

# and the largest model we are willing to consider
fit.full <- lm(mortality ~ ., data = poll)

## candidate 1: precip

# update() allows adjusting the model formula
update(fit.null, . ~ . + precip, data = poll)

# test the added term
update(fit.null, . ~ . + precip, data = poll) |>
  anova(fit.null)

## candidate 2: lnox

# now you try
update(fit.null, ...) |>
  anova(fit.null)

## now try every term

# add1() fits every single-term addition and reports the F test for each
add1(fit.null, scope = formula(fit.full), test = 'F')

# which p value is smallest?
add1(fit.null, scope = formula(fit.full), test = 'F') |>
  dplyr::slice_min(`Pr(>F)`)

# nonwhite wins
fit.fwd1 <- update(fit.null, . ~ . + nonwhite)

## one more forward step

# now you try (find the smallest p value by inspection)
add1(...)

# store the next step
fit.fwd2 <- ...

# ===========================================================================
# forward selection
# ===========================================================================
# Stepwise selection is add1() repeated: add the winner, re-run add1() on
# what is left, stop when the best candidate is no longer significant. 

# read the documentation
?stepwise

# forward selection
fwd <- stepwise(mortality ~ ., data = poll, type = 'linear',    # model space
                strategy = 'forward',                           # algorithm
                metric = 'SL',                                  # criterion
                sle = 0.05)                                     # alpha

# inspect the path (check this matches your answer from before)
fwd$overview$forward$SL

# retrieve selected model
fit.fwd <- fwd$forward$SL
summary(fit.fwd)

# check your understanding: how do you adjust alpha to obtain a smaller model?
stepwise(mortality ~ ., data = poll, type = 'linear',
         strategy = 'forward', metric = 'SL', sle = 0.05)

# ===========================================================================
# backward elimination
# ===========================================================================

# candidate 1: precip
update(fit.full, . ~ . - precip) |> 
  anova(fit.full)

# candidate 2: humidity
update(fit.full, . ~ . - humidity) |> 
  anova(fit.full)

# a single backward step (verify that the two examples above match first two rows)
drop1(fit.full, scope = formula(fit.full), test = 'F')

# backward elimination via stepwise()
bwd <- stepwise(mortality ~ ., data = poll, type = 'linear',    # model space
                strategy = 'backward',                          # algorithm
                metric = 'SL',                                  # criterion
                sls = 0.05)                                     # alpha

# the path shows REMOVED predictors
bwd$overview$backward$SL

# retrieve selected model to see RETAINED predictors
fit.bwd <- bwd$backward$SL
summary(fit.bwd)

# check your understanding: how do you adjust alpha to obtain a smaller model?
stepwise(mortality ~ ., data = poll, type = 'linear',
         strategy = 'backward', metric = 'SL', sls = 0.05)

# ===========================================================================
# AIC and BIC
# ===========================================================================

# forward and backward selection differ on lso2
formula(fit.fwd)
formula(fit.bwd)

# compare AIC: smaller is better
c(forward = AIC(fit.fwd), backward = AIC(fit.bwd))

# compare BIC: smaller is better
c(forward = BIC(fit.fwd), backward = BIC(fit.bwd))

# stepwise() also allows a stopping rule by minimum AIC/BIC
# this picks a much larger model!
stepwise(mortality ~ ., data = poll, type = 'linear',
         strategy = 'backward', metric = 'AIC')

# ===========================================================================
# all subsets
# ===========================================================================

# why not just fit every possible model? this is called "all subsets". it's 
# more computationally expensive, but feasible with, e.g., 15 predictors

# choose by smallest aic
allsub <- stepwise(mortality ~ ., data = poll, type = 'linear',
                    strategy = 'subset', metric = 'AIC')
fit.all.aic <- allsub$subset$AIC
formula(fit.all.aic)

# your turn: choose by smallest *bic*
allsub <- stepwise(...)
fit.all.bic <- ...
formula(fit.all.bic)

# ===========================================================================
# YOUR TURN: DIABETES
# ===========================================================================
# 442 patients. The response is a quantitative measure of disease progression
# one year after baseline; the ten predictors are age, sex, body mass index,
# mean arterial pressure (map), and six blood serum measurements.

# install.packages('lars')
library(lars)
data(diabetes)

# rearrange diabetes data into a dataframe
diab <- data.frame(y = diabetes$y, unclass(diabetes$x))
str(diab)

# select a model: you choose algorithm and tuning parameters
step.diab <- stepwise(...)
fit.diab <- ...