01-Motivation.Rmd

# Introduction and Example Datasets

```{r include=FALSE}
library(ggplot2)
library(dplyr)
library(stringr)
library(readr)
library(palmerpenguins)

knitr::opts_chunk$set(fig.align="center")
```

Regression is one of the most flexible and widely-used tools for inferential data analysis. This book introduces the statistical method of linear regression, starting with simple linear regression and then expanding to multiple linear regression.


:::: {.examplebox}

<!-- ![](img/lter_penguins.png) -->
<!-- (Artwork by \@allison_horst)  -->
<!-- ^[https://github.com/allisonhorst/palmerpenguins]) -->

```{example penguin-flip-mass-question} 
At the Palmer research station in Antarctica^[https://pal.lternet.edu], researchers made measurements on three different penguin species: Adélie, Chinstrap, and Gentoo. The `palmerpenguin` R package (@palmerpenguins) contains these measurements, which we will use in many of the examples. There are 342 measurements of penguins' flipper length and body mass:
```

```{r g-penguin-flip-mass, echo=FALSE, fig.cap="Flipper length and body mass in the Palmer Penguin dataset.", out.width="75%"}
g_penguin_flip_mass <- penguins %>%
  rename(Species=species) %>%
  filter(!is.na(body_mass_g)) %>%
  ggplot(aes(x=flipper_length_mm,
           y=body_mass_g)) +
  theme_bw() + 
  geom_point(aes(col=Species)) +
    xlab("Flipper Length (mm)") +
    ylab("Body Mass (g)")
g_penguin_flip_mass
```

From Figure \@ref(fig:g-penguin-flip-mass), there clearly seems to be a relationship between flipper length and body mass of the penguins. Can we quantify the size and strength of this relationship?

Yes! Using linear regression we can find the best fitting line through the data:

```{r g-penguin-flip-mass-lm, echo=FALSE, message=FALSE, fig.cap="Flipper length and body mass in the Palmer Penguin dataset, with best fitting line.", out.width="70%"}
g_penguin_flip_mass_lm <- g_penguin_flip_mass +
  geom_smooth(se=FALSE, method="lm", col="black")
g_penguin_flip_mass_lm
```
The slope of this line tells us about the size and direction of the relationship between flipper length and average body mass. We can also estimate the amount of uncertainty in the intercept and slope of the line.

::::

## What is regression?

Finding the best fitting line in Figure \@ref(fig:g-penguin-flip-mass-lm) is an example of **simple linear regression**, which is a statistical method for explaining variability in an one quantity (e.g. penguin body mass) using variation in a different quantity (e.g. penguin flipper length).

More generally, regression is a statistical method for explaining variability in an **outcome variable**  (other names include: dependent variable, response variable) using variation in one or more **predictor variables** (other names include: independent variable, covariate, explanatory variable, regressor).

In linear regression, the outcome variable is always a continuous variable, meaning that it takes on any numerical values. The values might be only over a specific range (e.g. body mass is necessarily greater than zero) and have a limited precision (e.g. body mass is only measured out to a certain number of decimal places). If the outcome variable is binary (e.g. sick/healthy) or categorical (e.g. species), then a different type of regression model is needed. This book will focus primarily on linear regression, but we will cover other types of regression in Chapters \@ref(logistic) and \@ref(logisticinference).

On the other hand, the predictor variables can be almost any type, whether continuous or categorical. In Simple Linear Regression (Chapters \@ref(slr) through \@ref(slrscaling)), there is only one continuous predictor variable. But in its general form, Multiple Linear Regression (Chapter \@ref(mlr)) can include an arbitrary number of parameters. The predictors can be any combination of variable types, including interactions (Chapter \@ref(indinter)) and non-linear transformations (Chapter \@ref(transformations)).

## Regression Goals
We can describe the goals of regression in two ways: the scientific goals and the practical (mathematical) goals.

### Scientific Goals {#sec:scigoals}

The **scientific goals** of regression are: description, inference, and prediction.

The **descriptive** goal of regression is focused on showing the relationship between $x$ and $y$. Emphasis is on quantifying and visualizing the relationship, rather than on drawing specific conclusions. Descriptive goals are inherently exploratory, and don't always require a specific question. Sometimes, we have data on an entire population, in which case the goal of the analysis is about describing relationships rather than making inference for a larger population.

One of the fundamental goals of statistics, and one that sets it apart from other data science fields, is **inference**. Inference takes a descriptive analysis one step further, and uses information about uncertainty to quantify the *strength* of a relationship. In short, inference is used to answer the question: *Is there a relationship between $x$ and $y$, beyond what we would just expect by chance?*
Fundamental tools in statistical inference are confidence intervals and hypothesis testing. Inferential goals are also sometimes called *association* goals, since the objective is to learn about an association (or lack thereof) between variables.

A third, and fundamentally different, goal of regression analysis is **prediction**, when we want to predict values of $y$ for new observations using their value of $x$. Prediction is the primary goal in the area of statistical (machine) learning, and regression is only one of many tools that can be used to make predictions for a dataset. This book covers prediction in Chapter \@ref(modelselectionprediction). 

:::: {.examplebox}
```{example }
Consider the following questions using the penguin data in Figure \@ref(fig:g-penguin-flip-mass-lm):
```

1. **Descriptive**: What is the shape of the relationship between flipper length and body mass?
2. **Inferential**: Is there a relationship between flipper length and average body mass in penguins?
3. **Prediction**: What is the predicted body mass for a penguin with 200mm flippers?
4. **Inferential**: What is the estimated differences in body mass between penguins who differ in flipper length by 50mm?
    
::::


### Mathematical Goals

The **mathematical goals** goals of regression are what we need to accomplish to achieve our scientific goal (regardless of whether that is description, inference, or prediction). Visually, the goal of linear regression is to find the best fitting line for the $y$'s. In later chapters, we will see how to generalize this to any nonlinear function. Mathematically, we do this by estimating **parameters** and their **uncertainty**.

## Example Datasets

Throughout this book, we will make repeated use of several example datasets, which we now describe.

### Palmer Penguins (Part 2)

:::: {.examplebox}

```{example penguin-exposition-billdepth-sex}
Other values in the Palmer Penguins dataset include penguin sex and  measurements of bill depth. The following histograms show the distributions of bill depth by sex:
```  
  
```{r g-penguin-bill-sex, echo=F, fig.cap="Distribution of penguin bill length by sex.", out.width="70%"}
penguins %>%
  filter(!is.na(sex)) %>%
  mutate(Sex=str_to_title(sex)) %>%
ggplot(aes(x=bill_depth_mm,
                        fill=Sex)) + 
  theme_bw() + 
  geom_histogram(alpha=0.5, position="identity", bins=25) +
    xlab("Bill Depth (mm)")
```

In Figure \@ref(fig:g-penguin-bill-sex), it appears that there might be a difference in the average bill length by sex. We can use linear regression to quantify this difference.

:::: 


:::: {.examplebox}
```{example penguin-exposition-billdepth-species}
The following histograms show the distributions of bill depth by species:
```

```{r g-penguin-bill-species, echo=F, message=FALSE, fig.cap="Distribution of penguin bill length by species.", out.width="70%"}
penguins %>%
  filter(!is.na(bill_depth_mm)) %>%
    rename(Species=species) %>%
ggplot(aes(x=bill_depth_mm,
                        fill=Species)) + 
  theme_bw() + 
  geom_histogram(alpha=0.5, position="identity", bins=25) +
    xlab("Bill Depth (mm)")
```

**Q:** It looks like there might be a difference in the distributions of bill length by species. Can we quantify this?

**A:** Yes! Using linear regression, we can estimate the difference in mean bill length between the three penguin species.

:::: 


:::: {.examplebox}
```{example penguin-exposition-billdepth-speciessex}
Figure \@ref(fig:g-penguin-bill-species-sex) shows the distribution of bill length by species and sex.
```

```{r g-penguin-bill-species-sex, echo=F, message=FALSE, fig.cap="Distribution of penguin bill length by species."}
penguins %>%
  filter(!is.na(sex)) %>%
ggplot(aes(x=bill_depth_mm)) + 
  theme_bw() + 
  facet_grid(species~sex) + 
  geom_histogram(alpha=0.5, position="identity", bins=25) +
    xlab("Bill Depth (mm)")
```

**Q:**  It looks like there might be a difference in the distributions of bill length by *both* species and sex. Can we quantify this?

**A:** Yes! Using linear regression, we can estimate the mean bill length for all possible combinations of sex and species.

:::: 

### Baseball Hits

More so than almost any sport, baseball is full of statistics. While traditional statistics such as batting average and earned run average have long been collected, the tracking systems in modern baseball stadiums allow a plethora of information to be collected for each batted ball. 

```{r echo=FALSE, include=FALSE}
bb <- read_csv("data/rockies_hits_2019.csv")
bb <- bb %>%
  mutate(events=factor(events, levels=c("single", "double", "triple", "home_run"))) %>%
  mutate(Outcome=factor(hit, levels=c(T, F), labels=c("Hit", "Out")))
```

:::: {.examplebox}
```{example bb-exposition-lspeed-hitdist}
Figure \@ref(fig:bb-scatter-lspeed-hitdist) shows all of the balls in play from the Colorado Rockies in the 2019 season. Along the horizontal axis is the launch speed, which is the speed at which the ball left the bat when hit. The vertical axis shows the distance the ball traveled. The color of each point shows whether the batted ball resulted in a hit or an out.
```
```{r bb-scatter-lspeed-hitdist, echo=FALSE, fig.cap="Distance travelled and launch speed for balls in play from Colorado Rockies in 2019."}
g_bb_lspeed_hitdist <- bb %>%
  filter(!is.na(launch_speed)) %>%
  ggplot() + theme_bw() + 
  geom_point(aes(x=launch_speed, y=hit_distance, col=Outcome), alpha=0.3) +
  xlab("Launch Speed (mph)") +
  ylab("Distance Travelled (ft)")
g_bb_lspeed_hitdist
```

There are several questions we could ask using the data from Figure \@ref(fig:bb-scatter-lspeed-hitdist):

* Is there a relationship between launch speed and hit distance? If yes, can that relationship be quantified?
* For a ball hit at 90 mph, what is the predicted distance traveled?
* Using launch speed, can we predict whether or not a batted ball will result in a hit or an out?

:::: 

:::: {.examplebox}
```{example bb-exposition-langle-hitdist}
Another variable that is measured is the launch angle of balls in play, which is the angle (in degrees) that the ball leaves the bat. A launch angle of 0 means the ball was hit straight out; a positive angle means it was hit up in to the air and a negative angle means it was hit down into the ground.
```
```{r bb-scatter-langle-hitdist, echo=FALSE, fig.cap="Distance travelled and launch angle for balls in play from Colorado Rockies in 2019."}
g_bb_langle_hitdist <- bb %>%
  filter(!is.na(launch_speed)) %>%
  ggplot() + theme_bw() + 
  geom_point(aes(x=launch_angle, y=hit_distance, col=Outcome), alpha=0.3) +
  xlab("Launch Angle (degrees)") +
  ylab("Distance Travelled (ft)")
g_bb_langle_hitdist
```

There are several questions we could ask using the data from Figure \@ref(fig:bb-scatter-langle-hitdist):

* Is there a relationship between launch angle and hit distance? If yes, can that relationship be quantified?
* For a ball hit at 45 degrees, what is the predicted distance traveled?
* Using launch angle, can we predict whether or not a batted ball will result in a hit or an out?

:::: 


### Housing Price {#housingprices}

Many factors impact housing prices, perhaps most importantly the economic demographics of the region. Zillow, the online real estate company, maintains a public database of housing-related data at <https://www.zillow.com/research/data/>. Combining this data with demographic information from the Census can let us analyze different housing trends.

:::: {.examplebox}
```{example housing-exposition-income-priace}
Figure \@ref(fig:g-housing-inc-price) shows data from Zillow on the median price for a single-family residence in January 2018 in different metropolitan areas across the U.S. On the horizontal axis is the 2018 median household income for the area, taken from the U.S. Census' American Community Survey.
```

```{r include=FALSE}
housing_inc <- read_csv("data/housing_income.csv")
housing_inc <- subset(housing_inc, !is.na(median_income))
```

```{r g-housing-inc-price, echo=FALSE, fig.cap="Median single-famly residence prices and annual income in U.S. metropolitan areas."}
gg_housing_inc <- ggplot(housing_inc, aes(x=median_income,
                 y=median_sale_price)) + theme_bw() + 
  geom_point() +
  ylab("Median House Price ($)") +
  xlab("Median Household Income ($)")
gg_housing_inc
```

In Figure \@ref(fig:g-housing-inc-price), there seems to be a relationship between median house price and median annual income. One question we might ask is, how robust is the relationship to removing the outlier point in the top right? The data for Figure \@ref(fig:g-housing-inc-price) can be downloaded here: [housing_income.csv](data/housing_income.csv).

:::: 


### Bike Share Programs {#bikedata}

:::: {.examplebox}
A bike sharing program recorded meteorological factors as part of efforts to understand factors related to bike usage on weekdays. A subset of this data is plotted in Figure \@ref(fig:bike-temp-usage)

```{r include=FALSE}
bike <- read_csv("data/bike_sharing.csv",
                 col_types = list(year=col_character()))
```

```{r bike-temp-usage, echo=FALSE, fig.cap="Bike sharing data."}
bike %>%
  rename(Year=year) %>%
ggplot() + theme_bw() + 
    geom_point(aes(x=temp, y=registered_users, col=Year)) +
  xlab("Temperature (Celsius)") +
  ylab("Number of Active Registered Users")
```


It looks like there is a positive relationship between temperature and the number of active bike users. There also appears to be a year effect, with more users in 2012 compared to 2011. Can we use temperature, year, and other factors such as humidity and windspeed to predict the number of users active for a given day?

:::: 


### Car fuel efficiency {#mpgdata}


:::: {.examplebox}
The `ggplot2` package contains the `mpg` dataset on fuel efficiency for vehicles from 1999-2008. Figure \@ref(fig:mpg-displ-cyl) shows the relationship between displacement and city miles per gallon.

```{r mpg-displ-cyl, echo=FALSE, fig.cap="Fuel efficiency data from `mpg` dataset."}
mpg %>%
  mutate(cyl=as.character(cyl)) %>%
ggplot() + theme_bw() + 
  geom_point(aes(x=displ,
                 y=cty,
                 col=cyl)) +
    labs(color="Number of\nCylinders",
         x="Engine Displacement (Liters)",
         y="City Miles Per Gallon")
```
:::: 

<!-- ### Photosynthesis Data -->

<!-- (to be added) -->

<!-- ### Car Mileage Data -->

<!-- (to be added) -->