inf-model-applications.qmd

# Applications: Model and infer {#sec-inf-model-applications}

```{r}
#| include: false
source("_common.R")
```

## Case study: Mario Kart {#sec-case-study-mario-cart}

In this case study, we consider Ebay auctions of a video game called *Mario Kart* for the Nintendo Wii.
The outcome variable of interest is the total price of an auction, which is the highest bid plus the shipping cost.
We will try to determine how total price is related to each characteristic in an auction while simultaneously controlling for other variables.
For instance, all other characteristics held constant, are longer auctions associated with higher or lower prices?
And, on average, how much more do buyers tend to pay for additional Wii wheels (plastic steering wheels that attach to the Wii controller) in auctions?
Multiple regression will help us answer these and other questions.

::: {.data data-latex=""}
The [`mariokart`](http://openintrostat.github.io/openintro/reference/mariokart.html) data can be found in the [**openintro**](http://openintrostat.github.io/openintro) R package.
:::

The `mariokart` dataset includes results from 141 auctions.
Four observations from this dataset are shown in @tbl-mariokart-data-frame, and descriptions for each variable are shown in @tbl-mariokart-var-def.
Notice that the condition and stock photo variables are **indicator variables**\index{indicator variable}\index{variable!indicator}, similar to `bankruptcy` in the `loans` dataset from @sec-inf-model-mlr.

```{r}
#| label: tbl-mariokart-data-frame
#| tbl-cap: Top four rows of the `mariokart` dataset.
#| tbl-pos: H
mariokart <- mariokart |>
  mutate(price = total_pr, cond_new = forcats::fct_rev(cond)) |>
  filter(price < 100)

mariokart |>
  select(price, cond_new, stock_photo, duration, wheels) |>
  slice_head(n = 4) |>
  kbl(linesep = "", booktabs = TRUE, row.names = FALSE, align = "ccccc") |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"), 
    latex_options = c("striped"),
    full_width = FALSE
  ) |>
  column_spec(1:5, width = "7em")
```

```{r}
#| label: tbl-mariokart-var-def
#| tbl-cap: Variables and their descriptions for the `mariokart` dataset.
#| tbl-pos: H
mariokart_var_def <- tribble(
  ~variable, ~description,
  "price", "Final auction price plus shipping costs, in US dollars.",
  "cond_new", "Indicator variable for if the game is new (1) or used (0).",
  "stock_photo", "Indicator variable for if the auction's main photo is a stock photo.",
  "duration", "The length of the auction, in days, taking values from 1 to 10.",
  "wheels", "The number of Wii wheels included with the auction. A Wii wheel is an optional steering wheel accessory that holds the Wii controller."
)

mariokart_var_def |>
  kbl(linesep = "", booktabs = TRUE, col.names = c("Variable", "Description")) |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped"), 
    full_width = TRUE
  ) |>
  column_spec(1, monospace = TRUE) |>
  column_spec(2, width = "30em")
```

\clearpage

### Mathematical approach to linear models

In @tbl-mariokart-model-output we fit a mathematical linear regression model with the game's condition as a predictor of auction price.

$$E[\texttt{price}] = \beta_0 + \beta_1\times \texttt{cond\_new}$$

Results of the model are summarized below:

```{r}
#| label: tbl-mariokart-model-output
#| tbl-cap: Summary of a linear model for predicting `price` based on `cond_new`.
#| tbl-pos: H
mariokart <- mariokart |>
  mutate(
    cond_new = ifelse(cond_new == "used", 0 , 1),
    stock_photo = ifelse(stock_photo == "yes", 1, 0)
  )

lm(price ~ cond_new, data = mariokart) |> 
  tidy() |>
  mutate(p.value = ifelse(p.value < 0.0001, "<0.0001", round(p.value, 4))) |>
  kbl(linesep = "", booktabs = TRUE, digits = 2, align = "lrrrr") |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped")
  ) |>
  column_spec(1, width = "10em", monospace = TRUE) |>
  column_spec(2:5, width = "5em")
```

::: {.guidedpractice data-latex=""}
Write down the equation for the model, note whether the slope is statistically different from zero, and interpret the coefficient.[^27-inf-model-applications-1]
:::

[^27-inf-model-applications-1]: The equation for the line may be written as $\widehat{\texttt{price}} = 47.15 + 10.90\times \texttt{cond\_new}$.
    Examining the regression output in @tbl-mariokart-model-output we can see that the p-value for `cond_new` is very close to zero, indicating there is strong evidence that the coefficient is different from zero when using this one-variable model.
    The variable `cond_new` is a two-level categorical variable that takes value 1 when the game is new and value 0 when the game is used.
    This means the 10.90 model coefficient predicts a price of an extra \$10.90 for those games that are new versus those that are used.

Sometimes there are underlying structures or relationships between predictor variables.
For instance, new games sold on Ebay tend to come with more Wii wheels, which may have led to higher prices for those auctions.
We would like to fit a model that includes all potentially important variables simultaneously, which would help us evaluate the relationship between a predictor variable and the outcome while controlling for the potential influence of other variables.

We want to construct a model that accounts for not only the game condition but simultaneously accounts for three other variables:

$$
E[\texttt{price}] = \beta_0 + \beta_1 \times \texttt{cond\_new} + \beta_2\times \texttt{stock\_photo} + \beta_3 \times \texttt{duration} + \beta_4 \times  \texttt{wheels}
$$

@tbl-mariokart-full-model-output summarizes the full model.
Using the output, we identify the point estimates of each coefficient and the corresponding impact (measured with information on the standard error used to compute the p-value).

```{r}
#| label: tbl-mariokart-full-model-output
#| tbl-cap: |
#|   Summary of a linear model for predicting `price` based on `cond_new`,
#|   `stock_photo`, `duration`, and `wheels`.
#| tbl-pos: H
lm(price ~ cond_new + stock_photo + duration + wheels, data = mariokart) |> 
  tidy() |>
  mutate(p.value = ifelse(p.value < 0.0001, "<0.0001", round(p.value, 4))) |>
  kbl(linesep = "", booktabs = TRUE, digits = 2, align = "lrrrr") |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped")
  ) |>
  column_spec(1, width = "10em", monospace = TRUE) |>
  column_spec(2:5, width = "5em")
```

::: {.guidedpractice data-latex=""}
Write out the model's equation using the point estimates from @tbl-mariokart-full-model-output.
How many predictors are there in the model?
How many coefficients are estimated?[^27-inf-model-applications-2]
:::

[^27-inf-model-applications-2]: $\widehat{\texttt{price}} = 36.21 + 5.13 \times \texttt{cond\_new} + 1.08 \times \texttt{stock\_photo} - 0.03 \times \texttt{duration} + 7.29 \times \texttt{wheels},$ with 4 predictors but 5 coefficients (including the intercept).

::: {.guidedpractice data-latex=""}
What does $\beta_4,$ the coefficient of variable $x_4$ (Wii wheels), represent?
What is the point estimate of $\beta_4?$[^27-inf-model-applications-3]
:::

[^27-inf-model-applications-3]: In the population of all auctions, it is the average difference in auction price for each additional Wii wheel included when holding the other variables constant.
    The point estimate is $b_4 = 7.29$

::: {.guidedpractice data-latex=""}
Compute the residual of the first observation in @tbl-mariokart-data-frame using the equation identified in @tbl-mariokart-full-model-output.[^27-inf-model-applications-4]
:::

[^27-inf-model-applications-4]: $e_i = y_i - \hat{y_i} = 51.55 - 49.62 = 1.93$.

::: {.workedexample data-latex=""}
In @tbl-mariokart-model-output, we estimated a coefficient for `cond_new` in of $b_1 = 10.90$ with a standard error of $SE_{b_1} = 1.26$ when using simple linear regression.
Why might there be a difference between that estimate and the one in the multiple regression setting?

------------------------------------------------------------------------

If we examined the data carefully, we would see that there is **multicollinearity**\index{multicollinearity} among some predictors.
For instance, when we estimated the connection of the outcome `price` and predictor `cond_new` using simple linear regression, we were unable to control for other variables like the number of Wii wheels included in the auction.
That model was biased by the confounding variable `wheels`.
When we use both variables, this particular underlying and unintentional bias is reduced or eliminated (though bias from other confounding variables may still remain).
:::

### Computational approach to linear models

Previously, using a mathematical model, we investigated the coefficients associated with `cond_new` when predicting `price` in a linear model.

```{r}
#| label: fig-mariokart-rand-dist
#| fig-cap: |
#|   Estimated slopes from linear models (`price` regressed on `cond_new`) 
#|   built on 1,000 randomized datasets. Each dataset was permuted under the 
#|   null hypothesis.
#| fig-alt: |
#|   Estimated slopes from linear models (of price regressed on cond_new) 
#|   built on 1,000 randomized datasets. Each dataset was permuted under the 
#|   null hypothesis. The randomized slopes range from about -5 to +5, and the
#|   observed slope value is 10.9 which is not near the distribution of 
#|   randomized slopes.
#| fig-asp: 0.5
set.seed(47)
mariokart |>
  specify(price ~ cond_new) |>
  hypothesize(null = "independence") |>
  generate(reps = 1000, type = "permute") |>
  calculate(stat = "slope") |>
  ggplot(aes(x = stat)) +
  geom_histogram(binwidth = 0.5, fill = IMSCOL["green", "full"]) +
  labs(
    title = "Slopes from models fit to 1,000 permutations of the data",
    x = "Slope of model fit to permuted data",
    y = "Count"
    ) +
  geom_vline(xintercept = 10.9, color = IMSCOL["red","full"], size = 1) +
  theme(plot.title.position = "plot")
```


::: {.workedexample data-latex=""}
In @fig-mariokart-rand-dist, the red line (the observed slope) is far from the bulk of the histogram.
Explain why the randomly permuted datasets produce slopes that are quite different from the observed slope.

------------------------------------------------------------------------

The null hypothesis is that, in the population, there is no linear relationship between the `price` and the `cond_new` of the *Mario Kart* games.
When the data are randomly permuted, prices are randomly assigned to a condition (new or used), so that the null hypothesis is forced to be true, i.e., permutation is done under the assumption that no relationship between the two variables exists.
In the actual study, the new *Mario Kart* games do actually cost more (on average) than the used games!
So the slope describing the actual observed relationship is not one that is likely to have happened in a randomly dataset permuted under the assumption that the null hypothesis is true.
:::

::: {.guidedpractice data-latex=""}
Using the histogram in @fig-mariokart-rand-dist, find the p-value and conclude the hypothesis test in the context of the problem.[^27-inf-model-applications-5]
:::

[^27-inf-model-applications-5]: The observed slope is 10.9 which is nowhere near the range of values for the permuted slopes (roughly -5 to +5).
    Because the observed slope is not a plausible value under the null distribution, the p-value is essentially zero.
    We reject the null hypothesis and claim that there is a relationship between whether the game is new (or not) and the average predicted price of the game.

::: {.guidedpractice data-latex=""}
Is the conclusion based on the histogram of randomized slopes consistent with the conclusion obtained using the mathematical model?
Explain.[^27-inf-model-applications-6]
:::

[^27-inf-model-applications-6]: The p-value in @tbl-mariokart-model-output is also essentially zero, so the null hypothesis is also rejected when the mathematical model approach is taken.
    Often, the mathematical and computational approaches to inference will give quite similar answers.

Although knowing there is a relationship between the condition of the game and its price, we might be more interested in the difference in price, here given by the slope of the linear regression line.
That is, $\beta_1$ represents the population value for the difference in price between new *Mario Kart* games and used games.

```{r}
#| label: fig-mariokart-boot-dist
#| fig-cap: |
#|   Estimated slopes from linear models (`price` regressed on `cond_new`) 
#|   built on 1,000 bootstrapped datasets. Each bootstrap sample was
#|   taken from the original Mario Kart auction data.
#| fig-alt: |
#|   Estimated slopes from linear models (of price regressed on cond new) 
#|   built on 1,000 bootstrapped datasets. Each bootstrap sample was 
#|   taken from the original Mario Kart auction data. The histogram of the
#|   boostrapped slopes is bell-shaped, symmetric, and centered around the 
#|   observed value of 10.9.
#| fig-asp: 0.5
set.seed(47)
mariokart |>
  specify(price ~ cond_new) |>
  generate(reps = 1000, type = "bootstrap") |>
  calculate(stat = "slope") |>
  ggplot(aes(x = stat)) +
  geom_histogram(binwidth = 0.5, fill = IMSCOL["green", "full"]) +
  labs(
    title = "Slopes from models fit to 1,000 bootstrap samples of the data",
    x = "Slope of model fit to bootstrapped data",
    y = "Count"
    ) +
  theme(plot.title.position = "plot")
```


::: {.workedexample data-latex=""}
@fig-mariokart-boot-dist displays the slope estimates taken from bootstrap samples of the original data.
Using the histogram, estimate the standard error of the slope.
Is your estimate similar to the value of the standard error of the slope provided in the output of the mathematical linear model?

------------------------------------------------------------------------

The slopes seem to vary from approximately 8 to 14.
Using the empirical rule, we know that if a variable has a bell-shaped distribution, most of the observations will be with 2 standard errors of the center.
Therefore, a rough approximation of the standard error is 1.5.
The standard error given in @tbl-mariokart-model-output is 1.26 which is not too different from the value computed using the bootstrap approach.
:::

::: {.guidedpractice data-latex=""}
Use @fig-mariokart-boot-dist to create a 90% standard error bootstrap confidence interval for the true slope.
Interpret the interval in context.[^27-inf-model-applications-7]
:::

[^27-inf-model-applications-7]: Using the bootstrap SE method, we know the normal percentile is $z^\star = 1.645$, which gives a CI of $b_1 \pm 1.645 \cdot SE \rightarrow 10.9 \pm 1.645 \cdot 1.5 \rightarrow (8.43, 13.37).$ For games that are new, the average price is higher by between \$8.43 and \$13.37 than games that are used, with 90% confidence.

::: {.guidedpractice data-latex=""}
Use @fig-mariokart-boot-dist to create a 90% bootstrap percentile confidence interval for the true slope.
Interpret the interval in context.[^27-inf-model-applications-8]
:::

[^27-inf-model-applications-8]: Because there were 1,000 bootstrap resamples, we look for the cutoffs which provide 50 bootstrap slopes on the left, 900 in the middle, and 50 on the right.
    Looking at the bootstrap histogram, the rough 90% confidence interval is \$9 to \$13.10.
    For games that are new, the average price is higher by between \$9.00 and \$13.10 than games that are used, with 90% confidence.

### Cross-validation

In @sec-model-mlr, models were compared using $R^2_{adj}.$ In @sec-inf-model-mlr, however, a computational approach was introduced to compare models by removing chunks of data one at a time and assessing how well the variables predicted the observations that had been held out.

@fig-mariokart-cv-residuals was created by cross-validating models with the same variables as in @tbl-mariokart-model-output and @tbl-mariokart-full-model-output.
We applied 3-fold cross-validation, so 1/3 of the data was removed while 2/3 of the observations were used to build each model (first on `cond_new` only and then on `cond_new`, `stock_photo`, `duration`, and `wheels`).
Note that each time 1/3 of the data is removed, the resulting model will produce slightly different model coefficients.

The points in @fig-mariokart-cv-residuals represent the prediction (x-axis) and residual (y-axis) for each observation run through the cross-validated model.
In other words, the model is built (using the other 2/3) without the observation (which is in the 1/3) being used.
The residuals give us a sense for how well the model will do at predicting observations which were not a part of the original dataset, e.g., future studies.

```{r}
#| label: fig-mariokart-cv-residuals
#| fig-cap: |
#|   Cross-validation predictions and errors from linear models built on 
#|   two different sets of variables.
#| fig-subcap: 
#|   - "`price` vs. `cond_new`"
#|   - "`price` vs. `cond_new`, `stock_photo`, `duration`, and `wheels`"
#| fig-alt: |
#|   Cross-validation predictions and errors from linear models built on 
#|   two different sets of variables. Left regressed price on cond new; 
#|   right regressed price on cond new, stock photo, duration, and wheels. For 
#|   each of the two models, the x-axis plots the predicted value and the 
#|   y-axis plots the prediction error. More variables lead to smaller 
#|   prediction errors.
#| out-width: 100%
#| layout-ncol: 2
#| fig-width: 5
#| fig-asp: 0.7
set.seed(4)
mario_CV1 <- train(price ~ cond_new,
  data = mariokart,
  na.action = na.omit,
  method = "lm",
  trControl = trainControl(
    method = "cv", number = 3,
    savePredictions = TRUE
  )
)

mario_CV2 <- train(price ~ cond_new + stock_photo + duration + wheels,
  data = mariokart,
  na.action = na.omit,
  method = "lm",
  trControl = trainControl(
    method = "cv", number = 3,
    savePredictions = TRUE
  )
)

SSE_CV1 <- mario_CV1$pred |>
  mutate(resid = obs - pred) |>
  select(resid) |>
  pull() %>%
  `^`(2) |>
  sum() |>
  round(2)

SSE_CV2 <- mario_CV2$pred |>
  mutate(resid = obs - pred) |>
  select(resid) |>
  pull() %>%
  `^`(2) |>
  sum() |>
  round(2)

mario_CV1$pred |>
  mutate(resid = obs - pred) |>
  ggplot() +
  geom_point(aes(x = pred, y = resid, color = Resample, shape = Resample)) +
  geom_rug(aes(y = resid)) +
  geom_hline(yintercept = 0) +
  ylim(c(-15, 25)) +
  xlim(c(20, 65)) +
  labs(
    color = "CV Fold", shape = "CV Fold",
    x = "Predicted", y = "Prediction error (residual)"
  ) +
  annotate(geom = "text", x = 55, y = 25, label = paste("CV SSE =", SSE_CV1)) +
  labs(title = "Model with one predictor") +
  theme(
    legend.position = c(0.2, 0.7),
    legend.background = element_rect(color = "gray", fill = "white")
  ) +
  scale_color_openintro("three")

mario_CV2$pred |>
  mutate(resid = obs - pred) |>
  ggplot() +
  geom_point(aes(x = pred, y = resid, color = Resample, shape = Resample)) +
  geom_rug(aes(y = resid)) +
  geom_hline(yintercept = 0) +
  ylim(c(-15, 25)) +
  xlim(c(20, 65)) +
  labs(
    color = "CV Fold", shape = "CV Fold",
    x = "Predicted", y = "Prediction error (residual)"
  ) +
  annotate(geom = "text", x = 55, y = 25, label = paste("CV SSE =", SSE_CV2)) +
  labs(title = "Model with four predictors") +
  scale_color_openintro("three")
```

::: {.guidedpractice data-latex=""}
In @fig-mariokart-cv-residuals-2, note the point at roughly predicted = 50 and prediction error = 10.
Estimate the observed and predicted value for that observation.[^27-inf-model-applications-9]
:::

[^27-inf-model-applications-9]: The predicted value is roughly $\widehat{\texttt{price}} = \$50.$ The observed value is roughly $\texttt{price}_i = \$60$ riders (using $e_i = y_i - \hat{y}_i).$

::: {.guidedpractice data-latex=""}
In @fig-mariokart-cv-residuals-2, for the same point at roughly predicted = 50 and prediction error = 10, describe which cross-validation fold(s) were used to build its prediction model.[^27-inf-model-applications-10]
:::

[^27-inf-model-applications-10]: The point appears to be in fold 2, so folds 1 and 3 were used to build the prediction model.

::: {.guidedpractice data-latex=""}
By noting the spread of the cross-validated prediction errors (on the y-axis) in @fig-mariokart-cv-residuals, which model should be chosen for a final report on these data?[^27-inf-model-applications-11]
:::

[^27-inf-model-applications-11]: The cross-validated residuals on `cond_new` vary roughly from -15 to 15, while the cross-validated residuals on the four predictor model vary less, roughly from -10 to 10.
    Given the smaller residuals from the four predictor model, it seems as though the larger model is better.

::: {.guidedpractice data-latex=""}
Using the summary statistic cross-validation sum of squared errors (CV SSE), which model should be chosen for a final report on these data?[^27-inf-model-applications-12]
:::

[^27-inf-model-applications-12]: The CV SSE is smaller (by a factor of almost two!) for the model with four predictors.
    Using a single valued criterion (CV SSE) allows us to make a decision to choose the model with four predictors.

\clearpage

## Interactive R tutorials {#inferential-modeling-tutorials}

Navigate the concepts you've learned in this part in R using the following self-paced tutorials.
All you need is your browser to get started!

::: {.alltutorials data-latex=""}
[Tutorial 6: Inferential modeling](https://openintrostat.github.io/ims-tutorials/06-model-infer/)

::: {.content-hidden unless-format="pdf"}
https://openintrostat.github.io/ims-tutorials/06-model-infer
:::
:::

::: {.singletutorial data-latex=""}
[Tutorial 6 - Lesson 1: Inference in regression](https://openintro.shinyapps.io/ims-06-model-infer-01/)

::: {.content-hidden unless-format="pdf"}
https://openintro.shinyapps.io/ims-06-model-infer-01
:::
:::

::: {.singletutorial data-latex=""}
[Tutorial 6 - Lesson 2: Randomization test for slope](https://openintro.shinyapps.io/ims-06-model-infer-02/)

::: {.content-hidden unless-format="pdf"}
https://openintro.shinyapps.io/ims-06-model-infer-02
:::
:::

::: {.singletutorial data-latex=""}
[Tutorial 6 - Lesson 3: t-test for slope](https://openintro.shinyapps.io/ims-06-model-infer-03/)

::: {.content-hidden unless-format="pdf"}
https://openintro.shinyapps.io/ims-06-model-infer-03
:::
:::

::: {.singletutorial data-latex=""}
[Tutorial 6 - Lesson 4: Checking technical conditions for slope inference](https://openintro.shinyapps.io/ims-06-model-infer-04/)

::: {.content-hidden unless-format="pdf"}
https://openintro.shinyapps.io/ims-06-model-infer-04
:::
:::

::: {.singletutorial data-latex=""}
[Tutorial 6 - Lesson 5: Inference beyond the simple linear regression model](https://openintro.shinyapps.io/ims-06-model-infer-05/)

::: {.content-hidden unless-format="pdf"}
https://openintro.shinyapps.io/ims-06-model-infer-05
:::
:::

::: {.content-hidden unless-format="pdf"}
You can also access the full list of tutorials supporting this book at\
<https://openintrostat.github.io/ims-tutorials>.
:::

::: {.content-visible when-format="html"}
You can also access the full list of tutorials supporting this book [here](https://openintrostat.github.io/ims-tutorials).
:::

## R labs

Further apply the concepts you've learned in this part in R with computational labs that walk you through a data analysis case study.

::: {.singlelab data-latex=""}
[Multiple linear regression - Grading the professor](https://www.openintro.org/go?id=ims-r-lab-model-infer)

::: {.content-hidden unless-format="pdf"}
https://www.openintro.org/go?id=ims-r-lab-model-infer
:::

:::

::: {.content-hidden unless-format="pdf"}
You can also access the full list of labs supporting this book at\
<https://www.openintro.org/go?id=ims-r-labs>.
:::

::: {.content-visible when-format="html"}
You can also access the full list of labs supporting this book [here](https://www.openintro.org/go?id=ims-r-labs).
:::