```{r echo = FALSE, results = 'hide', message = FALSE}
library(tidyverse)
library(mosaic)

theme_set(theme_bw())                 # change theme for ggplot2
trellis.par.set(theme=theme.mosaic()) # change default color scheme for lattice

knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    tidy    = FALSE)
```

The `plotly` package will let us make some interactive 3D plots

```{r loadpackages, message = FALSE}
library(mosaic)
library(plotly)
```


## Education Spending and SAT Scores

First we load the data.

```{r loadData}
data(SAT)
```

### Controlling for Participation Rate

Let's fit a model to predict SAT using 
1. taking test, and
2. expenditure / pupil

```{r mainEffectsModel}
m_frac_expend <- lm(sat ~ frac + expend, data = SAT)
m_frac_expend
```

```{r mainEffectsCoefficients}
coefficients(m_frac_expend) %>% round(digits = 2)
```

This chunk defines some convenience functions to make plotting easier

```{r convenienceFunctions}
get_prediction_grid <- function(x, z, pred_fun, data)
{
  require(mosaic)
  xmin <- min(data[[x]])
  xmax <- max(data[[x]])
  zmin <- min(data[[z]])
  zmax <- max(data[[z]])
  xmesh <- seq(xmin, xmax, length = 100)
  zmesh <- seq(zmin, zmax, length = 100)
  yhat <- outer(xmesh, zmesh, pred_fun)
  return(list(x1 = xmesh, x2 = zmesh, yhat = t(yhat)))
}

plot_sat_model_3D <- function(model)
{
  f_hat <- makeFun(model)
  surface <- get_prediction_grid("frac", "expend", f_hat, SAT)
  SAT %>%
    plot_ly(x = ~frac, y = ~expend, z = ~sat, size = 0.5, opacity = 0.5) %>%
    add_markers() %>%
    add_surface(x = surface$x1, y = surface$x2, z = surface$yhat)
}
```

### Visualizing the Data and Model in 3D

Let's make a 3D plot of the two variable model.

```{r}
plot_sat_model_3D(m_frac_expend)
```

How do the residuals look?

```{r}
gf_point(residuals(m_frac_expend) ~ fitted(m_frac_expend)) %>%
  gf_smooth()
```

Not great. We probably want a quadratic term for `frac`... but first...


### An Interaction Model

Let's fit a model with an interaction

```{r, message = FALSE}
library(plotly)
m_interaction <- lm(sat ~ frac + expend + frac:expend, data = SAT)
```

```{r}
coefficients(m_interaction) %>% round(2)
```

```{r}
plot_sat_model_3D(m_interaction)
```

```{r}
gf_point(residuals(m_interaction) ~ fitted(m_interaction)) %>%
  gf_smooth()
```

Not necessarily what we needed in this case.

### A Quadratic Control for `frac`

Let's add a quadratic term for `frac` to the no-interaction model

```{r}
m_frac2_expend <- lm(sat ~ frac + I(frac^2) + expend, data = SAT)
```

```{r}
coefficients(m_frac2_expend) %>% round(digits = 2)
```

```{r}
plot_sat_model_3D(m_frac2_expend)
```

```{r}
gf_point(residuals(m_frac2_expend) ~ fitted(m_frac2_expend)) %>%
  gf_smooth()
```

That definitely looks better...


### Quadratic and Interaction Together

What if we have both the quadratic term and the interaction?

```{r}
m_frac2_expend_interaction <- 
  lm(sat ~ frac + expend + I(frac^2) + expend:frac, data = SAT)
```

```{r}
coefficients(m_frac2_expend_interaction) %>% round(digits = 2)
```

```{r}
plot_sat_model_3D(m_frac2_expend_interaction)
```

```{r}
gf_point(residuals(m_frac2_expend_interaction) ~ fitted(m_frac2_expend_interaction)) %>%
  gf_smooth()
```

#### Nested Tests to Answer Targeted Questions

Does expenditure have additional predictive value after controlling for participation rate?

Let's fit a reduced model that doesn't use expenditure.

```{r}
m_frac2 <- lm(sat ~ frac + I(frac^2), data = SAT)
```

Now we can do a nested test to compare this to our model that has expenditure as well as the interaction.

```{r}
anova(m_frac2, m_frac2_expend_interaction)
```

Seems like the evidence is on the borderline of the standard threshold for statistical significance.