This lab on Polynomial Regression and Step Functions in R comes from p. 288-292 of "Introduction to Statistical Learning with Applications in R" by Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani. It was re-implemented in Fall 2016 in `tidyverse` format by Amelia McNamara and R. Jordan Crouser at Smith College.


```{r}
library(ISLR)
library(dplyr)
library(ggplot2)
```

# 7.8.1 Polynomial Regression and Step Functions

In this lab, we'll explore how to generate the `Wage` dataset models we saw in class. We first fit the polynomial regression model using the following command:


```{r}
fit = lm(wage ~ poly(age, 4), data = Wage)
coef(summary(fit))
```

This syntax fits a linear model, using the `lm()` function, in order to predict
wage using a fourth-degree polynomial in `age`: `poly(age,4)`. The `poly()` command
allows us to avoid having to write out a long formula with powers
of `age`. The function returns a matrix whose columns are a basis of orthogonal
polynomials, which essentially means that each column is a linear
combination of the variables `age`, `age^2`, `age^3` and `age^4`.

If we prefer, we can also use `poly()` to obtain `age`, `age^2`, `age^3` and `age^4`
directly. We can do this by using the `raw = TRUE` argument to
the `poly()` function. Later we see that this does not affect the model in a
meaningful way -- though the choice of basis clearly affects the coefficient
estimates, it does not affect the fitted values obtained.


```{r}
fit2 = lm(wage~poly(age, 4, raw = TRUE), data = Wage)
coef(summary(fit2))
```

We now create a grid of values for `age` at which we want predictions, and
then call the generic `predict()` function, specifying that we want standard
errors as well.


```{r}
# Get min/max values of age using the range() function
agelims = Wage %>%
    select(age) %>%
    range

# Generate a sequence of age values spanning the range
age_grid = seq(from = min(agelims), to = max(agelims))

# Predict the value of the generated ages,
# returning the standard error using se = TRUE
preds = predict(fit, newdata = list(age = age_grid), se = TRUE)

# Compute error bands (2*SE)
se_bands = cbind("upper" = preds$fit+2*preds$se.fit, 
                 "lower" = preds$fit-2*preds$se.fit)
```

Finally, we plot the data and add the fit from the degree-4 polynomial.


```{r}
ggplot() +
  geom_point(data = Wage, aes(x = age, y = wage)) +
  geom_line(aes(x = age_grid, y = preds$fit), color = "#0000FF") +
  geom_ribbon(aes(x = age_grid, 
                  ymin = se_bands[,"lower"], 
                  ymax = se_bands[,"upper"]), 
              alpha = 0.3) +
  xlim(agelims) +
  labs(title = "Degree-4 Polynomial")
```

We mentioned earlier that whether or not an orthogonal set of basis functions
is produced in the `poly()` function will not affect the model obtained
in a meaningful way. What do we mean by this? The fitted values obtained
in either case are identical (up to a miniscule rounding error caused by building our models on a computer):


```{r}
preds2 = predict(fit2, newdata = list(age = age_grid), se = TRUE)

# Calculate the difference between the two estimates, print out the first few values
head(abs(preds$fit - preds2$fit))
```

# Deciding on a degree

In performing a polynomial regression we must decide on the degree of
the polynomial to use. One way to do this is by using hypothesis tests. We
now fit models ranging from linear to a degree-5 polynomial and seek to
determine the simplest model which is sufficient to explain the relationship
between `wage` and `age`.

We can do this using the `anova()` function, which performs an
analysis of variance (ANOVA, using an F-test) in order to test the null
hypothesis that a model $M_1$ is sufficient to explain the data against the 
alternative hypothesis that a more complex model $M_2$ is required. In order
to use the `anova()` function, $M_1$ and $M_2$ must be **nested models**: the
predictors in $M_1$ must be a subset of the predictors in $M_2$. In this case,
we fit five different models and sequentially compare the simpler model to
the more complex model:


```{r}
fit_1 = lm(wage~age, data = Wage)
fit_2 = lm(wage~poly(age,2), data = Wage)
fit_3 = lm(wage~poly(age,3), data = Wage)
fit_4 = lm(wage~poly(age,4), data = Wage)
fit_5 = lm(wage~poly(age,5), data = Wage)
print(anova(fit_1,fit_2,fit_3,fit_4,fit_5))
```

The $p$-value comparing the linear Model 1 to the quadratic Model 2 is
essentially zero $(<10^{-15})$, indicating that a linear fit is not sufficient. Similarly the $p$-value comparing the quadratic Model 2 to the cubic Model 3
is very low (0.0017), so the quadratic fit is also insufficient. The $p$-value
comparing the cubic and degree-4 polynomials, Model 3 and Model 4, is approximately
0.05 while the degree-5 polynomial Model 5 seems unnecessary because its $p$-value is 0.37. Hence, either a cubic or a quartic polynomial appear to provide a reasonable fit to the data, but lower- or higher-order models are not justified.

In this case, instead of using the `anova()` function, we could also have obtained
these $p$-values more succinctly by exploiting the fact that `poly()` creates
orthogonal polynomials.


```{r}
print(coef(summary(fit_5)))
```

Notice that the p-values are the same, and in fact the square of the
t-statistics are equal to the F-statistics from the `anova()` function; for
example:


```{r}
(-11.983)^2
```

However, the ANOVA method works whether or not we used orthogonal
polynomials; it also works when we have other terms in the model as well.
For example, we can use `anova()` to compare these three models:


```{r}
fit_1 = lm(wage~education+age, data = Wage)
fit_2 = lm(wage~education+poly(age,2), data = Wage)
fit_3 = lm(wage~education+poly(age,3), data = Wage)
print(anova(fit_1,fit_2,fit_3))
```

As an alternative to using hypothesis tests and ANOVA, we could choose
the polynomial degree using cross-validation as we have in previous labs.

Next we consider the task of predicting whether an individual earns more
than \$250,000 per year. We proceed much as before, except that first we
create the appropriate response vector, and then apply the `glm()` function
using `family = "binomial"` in order to fit a polynomial logistic regression
model:


```{r}
fit = glm(I(wage>250)~poly(age,4), data = Wage, family = binomial)
```

Note that we again use the wrapper `I()` to create this binary response
variable on the fly. The expression `wage>250` evaluates to a logical variable
containing `TRUEs` and `FALSEs`, which `glm()` coerces to binary by setting the
`TRUEs` to 1 and the `FALSEs` to 0.

Once again, we make predictions using the `predict()` function:


```{r}
preds = predict(fit, newdata = list(age = age_grid), se = TRUE)
```

However, calculating the confidence intervals is slightly more involved than
in the linear regression case. The default prediction type for a glm() model
is type="link", which is what we use here. This means we get predictions
for the logit: that is, we have fit a model of the form

$$log\left(\frac{Pr(Y = 1|X)}{1 − Pr(Y = 1|X)}\right)= X\beta$$

and the predictions given are of the form $X\hat \beta$. The standard errors given are
also of this form. In order to obtain confidence intervals for $Pr(Y = 1|X)$,
we use the transformation:

$$Pr(Y = 1|X) = \frac{e^{X\beta}}{1 + e^{X\beta}}$$


```{r}
pfit = exp(preds$fit) / (1+exp(preds$fit))

se_bands_logit = cbind("upper" = preds$fit+2*preds$se.fit, 
                       "lower" = preds$fit-2*preds$se.fit)

se_bands = exp(se_bands_logit) / (1+exp(se_bands_logit))
```

We could have directly computed the probabilities by selecting
the `type = "response"` option in the `predict()` function. However, the corresponding confidence intervals would not have been sensible because we would end up with negative probabilities!

Now we're ready to draw the second plot we saw in class:


```{r}
high = Wage %>%
  filter(wage > 250)

low = Wage %>%
  filter(wage <= 250)

ggplot() +
  geom_rug(data = low, aes(x = jitter(age), y = wage), sides = "b", alpha = 0.3) +
  geom_rug(data = high, aes(x = jitter(age), y = wage), sides = "t", alpha = 0.3) +
  geom_line(aes(x = age_grid, y = pfit), color = "#0000FF") +
  geom_ribbon(aes(x = age_grid, 
                  ymin = se_bands[,"lower"], 
                  ymax = se_bands[,"upper"]), 
              alpha = 0.3) +
  xlim(agelims) +
  ylim(c(0,1)) +
  labs(title = "Degree-4 Polynomial",
       x = "Age",
       y = "P(wage > 250)")
```

We have drawn the `age` values corresponding to the observations with `wage`
values above 250 as gray marks on the top of the plot, and those with `wage`
values below 250 are shown as gray marks on the bottom of the plot. We
used the `jitter()` function to jitter the `age` values a bit so that observations
with the same `age` value do not cover each other up. This is often called a
*rug plot*.

# Step functions
In order to fit a step function, we use the `cut()` function:


```{r}
table(cut(Wage$age,4))
fit_step = lm(wage~cut(age,4), data = Wage)
print(coef(summary(fit)))
```

Here `cut()` automatically picked the cutpoints at 33.5, 49, and 64.5 years
of age. We could also have specified our own cutpoints directly using the
`breaks` option. The function `cut()` returns an ordered categorical variable;
the lm() function then creates a set of dummy variables for use in the regression.
The `age<33.5` category is left out, so the intercept coefficient of
\$94,160 can be interpreted as the average salary for those under 33.5 years
of age, and the other coefficients can be interpreted as the average additional
salary for those in the other age groups. 

We can produce predictions
and plots just as we did in the case of the polynomial fit_


```{r}
# Predict the value of the generated ages, returning the standard error using se = TRUE
preds = predict(fit_step, newdata = list(age = age_grid), se = TRUE)

# Compute error bands (2*SE)
se_bands = cbind("upper" = preds$fit+2*preds$se.fit, 
                 "lower" = preds$fit-2*preds$se.fit)

# Plot
ggplot() +
  geom_point(data = Wage, aes(x = age, y = wage)) +
  geom_line(aes(x = age_grid, y = preds$fit), color = "#0000FF") +
  geom_ribbon(aes(x = age_grid, 
                  ymin = se_bands[,"lower"], 
                  ymax = se_bands[,"upper"]), 
              alpha = 0.3) +
  xlim(agelims) +
  labs(title = "Step Function")
```

# Getting credit

To get credit for this lab, post your responses to the following questions:
 - What is one real-world example where you might try polynomial regression? Why?
 - What is one real-world example where you might try using a step function? Why?
 
to [Moodle](https://moodle.smith.edu/mod/quiz/view.php?id=262426).
