explore-numerical.qmd

# Exploring numerical data {#sec-explore-numerical}

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

::: {.chapterintro data-latex=""}
This chapter focuses on exploring **numerical** data using summary statistics and visualizations.
The summaries and graphs presented in this chapter are created using statistical software; however, since this might be your first exposure to the concepts, we take our time in this chapter to detail how to create them.
Mastery of the content presented in this chapter will be crucial for understanding the methods and techniques introduced in the rest of the book.
:::

Consider the `loan_amount` variable from the `loan50` dataset, which represents the loan size for each of 50 loans in the dataset.

This variable is numerical since we can sensibly discuss the numerical difference of the size of two loans.
On the other hand, area codes and zip codes are not numerical, but rather they are categorical variables.

Throughout this chapter, we will apply numerical methods using the `loan50` and `county` datasets, which were introduced in @sec-data-basics.
If you'd like to review the variables from either dataset, see Tables @tbl-loan-50-variables and @tbl-county-variables.

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

## Scatterplots for paired data {#sec-scatterplots}

A **scatterplot**\index{plot!scatterplot}\index{scatterplot} provides a case-by-case view of data for two numerical variables.
In @fig-county-multi-unit-homeownership, a scatterplot was used to examine the homeownership rate against the percentage of housing units that are in multi-unit structures (e.g., apartments) in the `county` dataset.
Another scatterplot is shown in @fig-loan50-amount-income, comparing the total income of a borrower `total_income` and the amount they borrowed `loan_amount` for the `loan50` dataset.
In any scatterplot, each point represents a single case.
Since there are `r nrow(loan50)` cases in `loan50`, there are `r nrow(loan50)` points in @fig-loan50-amount-income.

```{r}
#| include: false
terms_chp_05 <- c("scatterplot")
```

```{r}
#| label: fig-loan50-amount-income
#| fig-cap: |
#|   A scatterplot of loan amount versus total income for the `loan50`
#|   dataset.
#| fig-alt: |
#|   A scatterplot with total income on the x-axis and loan amount on the
#|   y-axis. The relationship is moderately positive.
#| fig-asp: 0.48
ggplot(loan50, aes(x = total_income, y = loan_amount)) +
  geom_point(alpha = 0.6, shape = 21, size = 3) +
  labs(x = "Total income", y = "Loan amount") +
  scale_x_continuous(labels = dollar_format(scale = 0.001, suffix = "K")) +
  scale_y_continuous(labels = dollar_format(scale = 0.001, suffix = "K"))
```

Looking at @fig-loan50-amount-income, we see that there are many borrowers with income below \$100,000 on the left side of the graph, while there are a handful of borrowers with income above \$250,000.

```{r}
#| label: fig-median-hh-income-poverty
#| fig-cap: |
#|   A scatterplot of the median household income against the poverty
#|   rate for the `county` dataset. Data are from 2017. A statistical model has
#|   also been fit to the data and is shown as a dashed line.
#| fig-alt: |
#|   A scatterplot with poverty rate on the x-axis and median household
#|   income on the y-axis. The relationship is negative and curvilinear.
#|   The bulk of the points fall with a poverty rate of 10% to 20% and a
#|   median household income of $25K to $75K.
#| fig-asp: 0.48
ggplot(county, aes(x = poverty / 100, y = median_hh_income)) +
  geom_point(
    alpha = 0.3,
    fill = IMSCOL["black", "full"],
    shape = 21,
    size = 3
  ) +
  geom_smooth(linetype = "dashed", color = IMSCOL["red", "full"], se = FALSE) +
  labs(x = "Poverty rate", y = "Median household income") +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  scale_y_continuous(labels = dollar_format(scale = 0.001, suffix = "K"))
```

\vspace{-5mm}

::: {.workedexample data-latex=""}
@fig-median-hh-income-poverty shows a plot of median household income against the poverty rate for `r nrow(county)` counties in the US.
What can be said about the relationship between these variables?

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

The relationship is evidently **nonlinear**\index{nonlinear}, as highlighted by the dashed line.
This is different from previous scatterplots we have seen, which indicate very little, if any, curvature in the trend.
:::

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "nonlinear")
```

\vspace{-5mm}

::: {.guidedpractice data-latex=""}
What do scatterplots reveal about the data, and how are they useful?[^05-explore-numerical-1]
:::

[^05-explore-numerical-1]: Answers may vary.
    Scatterplots are helpful in quickly spotting associations relating variables, whether those associations come in the form of simple trends or whether those relationships are more complex.

\vspace{-5mm}

::: {.guidedpractice data-latex=""}
Describe two variables that would have a horseshoe-shaped association in a scatterplot $(\cap$ or $\frown).$[^05-explore-numerical-2]
:::

[^05-explore-numerical-2]: Consider the case where your vertical axis represents something "good" and your horizontal axis represents something that is only good in moderation.
    Health and water consumption fit this description: we require some water to survive, but consume too much and it becomes toxic and can kill a person.

\clearpage

## Dot plots and the mean {#sec-dotplots}

Sometimes we are interested in the distribution of a single variable.
In these cases, a dot plot provides the most basic of displays.
A **dot plot**\index{plot!dot}\index{dot plot} is a one-variable scatterplot; an example using the interest rate of `r nrow(loan50)` loans is shown in @fig-loan-int-rate-dotplot.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "dot plot")
```

```{r}
#| label: fig-loan-int-rate-dotplot
#| fig-cap: |
#|   A dot plot of interest rate for the `loan50` dataset. The rates
#|   have been rounded to the nearest percent in this plot, and the
#|   distribution's mean is shown as a red triangle.
#| fig-alt: |
#|   A dot plot of interest rate (ranging from about 5% to 25%).
#|   The distribution is right skewed, and the mean is shown at an interest
#|   rate of about 11%.
#| fig-asp: 0.4
loan50_interest_rate_mean <- mean(loan50$interest_rate)

ggplot(loan50, aes(x = interest_rate)) +
  geom_dotplot() +
  labs(x = "Interest rate") +
  scale_x_continuous(labels = label_percent(scale = 1)) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  geom_polygon(
    data = tibble(
      x = c(
        loan50_interest_rate_mean - 1,
        loan50_interest_rate_mean + 1,
        loan50_interest_rate_mean
      ),
      y = c(-0.1, -0.1, 0)
    ),
    aes(x = x, y = y),
    fill = IMSCOL["red", "full"]
  )
```

The **mean**\index{mean}, often called the **average**\index{average} is a common way to measure the center of a **distribution**\index{distribution} of data.
To compute the mean interest rate, we add up all the interest rates and divide by the number of observations.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "mean", "average", "distribution")
```

The sample mean is often labeled $\bar{x}.$ The letter $x$ is being used as a generic placeholder for the variable of interest and the bar over the $x$ communicates we are looking at the average interest rate, which for these 50 loans is `r round(loan50_interest_rate_mean, 2)`%.
It's useful to think of the mean as the balancing point of the distribution, and it's shown as a triangle in @fig-loan-int-rate-dotplot.

::: {.important data-latex=""}
**Mean.**

The sample mean can be calculated as the sum of the observed values divided by the number of observations:

\vspace{-5mm}

$$ \bar{x} = \frac{x_1 + x_2 + \cdots + x_n}{n} $$
:::

\vspace{-5mm}

::: {.guidedpractice data-latex=""}
Examine the equation for the mean.
What does $x_1$ correspond to?
And $x_2$?
Can you infer a general meaning to what $x_i$ might represent?[^05-explore-numerical-3]
:::

[^05-explore-numerical-3]: $x_1$ corresponds to the interest rate for the first loan in the sample, $x_2$ to the second loan's interest rate, and $x_i$ corresponds to the interest rate for the $i^{th}$ loan in the dataset.
    For example, if $i = 4,$ then we are examining $x_4,$ which refers to the fourth observation in the dataset.

\vspace{-5mm}

::: {.guidedpractice data-latex=""}
What was $n$ in this sample of loans?[^05-explore-numerical-4]
:::

[^05-explore-numerical-4]: The sample size was $n = 50.$

The `loan50` dataset represents a sample from a larger population of loans made through Lending Club.
We could compute a mean for the entire population in the same way as the sample mean.
However, the population mean has a special label: $\mu.$ The symbol $\mu$ is the Greek letter *mu* and represents the average of all observations in the population.
Sometimes a subscript, such as $_x,$ is used to represent which variable the population mean refers to, e.g., $\mu_x.$ Oftentimes it is too expensive to measure the population mean precisely, so we often estimate $\mu$ using the sample mean, $\bar{x}.$

\clearpage

::: {.content-visible when-format="html"}
::: {.pronunciation data-latex=""}
The Greek letter $\mu$ is pronounced *mu*, listen to the pronunciation [here](https://youtu.be/PStgY5AcEIw?t=47).
:::
:::

::: {.content-visible when-format="pdf"}
::: {.pronunciation data-latex=""}
The Greek letter $\mu$ is pronounced *mu*.
:::
:::

::: {.workedexample data-latex=""}
Although we do not have an ability to *calculate* the average interest rate across all loans in the populations, we can *estimate* the population value using the sample data.
Based on the sample of 50 loans, what would be a reasonable estimate of $\mu_x,$ the mean interest rate for all loans in the full dataset?

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

The sample mean, `r round(loan50_interest_rate_mean, 2)`, provides a rough estimate of $\mu_x.$ While it is not perfect, this is our single best guess **point estimate**\index{point estimate} of the average interest rate of all the loans in the population under study.
In @sec-foundations-randomization and beyond, we will develop tools to characterize the accuracy of point estimates, like the sample mean.
As you might have guessed, point estimates based on larger samples tend to be more accurate than those based on smaller samples.
:::

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "point estimate")
```

The mean is useful because it allows us to rescale or standardize a metric into something more easily interpretable and comparable.
Suppose we would like to understand if a new drug is more effective at treating asthma attacks than the standard drug.
A trial of 1,500 adults is set up, where 500 receive the new drug, and 1000 receive a standard drug in the control group.
Results of this trial are summarized in @tbl-drug-asthma-results.

```{r}
#| label: tbl-drug-asthma-results
#| tbl-cap: Results of a trial of 1500 adults that suffer from asthma.
#| tbl-pos: H
drug_asthma <- tribble(
  ~x                     , ~`New drug` , ~`Standard drug` ,
  "Number of patients"   ,         500 ,             1000 ,
  "Total asthma attacks" ,         200 ,              300
)

drug_asthma |>
  kbl(
    linesep = "",
    booktabs = TRUE,
    col.names = c("", "New drug", "Standard drug"),
    align = "lcc"
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped"),
    full_width = FALSE
  ) |>
  column_spec(1, width = "12em") |>
  column_spec(2:3, width = "8em")
```

Comparing the raw counts of 200 to 300 asthma attacks would make it appear that the new drug is better, but this is an artifact of the imbalanced group sizes.

Instead, we should look at the average number of asthma attacks per patient in each group:

-   New drug: $200 / 500 = 0.4$ asthma attacks per patient
-   Standard drug: $300 / 1000 = 0.3$ asthma attacks per patient

The standard drug has a lower average number of asthma attacks per patient than the average in the treatment group.

::: {.workedexample data-latex=""}
Come up with another example where the mean is useful for making comparisons.

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

Emilio opened a food truck last year where he sells burritos, and his business has stabilized over the last 3 months.
Over that 3-month period, he has made \$11,000 while working 625 hours.
Emilio's average hourly earnings provides a useful statistic for evaluating whether his venture is, at least from a financial perspective, worth it:

$$ \frac{\$11000}{625\text{ hours}} = \$17.60\text{ per hour} $$

By knowing his average hourly wage, Emilio now has put his earnings into a standard unit that is easier to compare with many other jobs that he might consider.
:::

\clearpage

::: {.workedexample data-latex=""}
Suppose we want to compute the average income per person in the US.
To do so, we might first think to take the mean of the per capita incomes across the 3,142 counties in the `county` dataset.
What would be a better approach?

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

The `county` dataset is special in that each county actually represents many individual people.
If we were to simply average across the `income` variable, we would be treating counties with 5,000 and 5,000,000 residents equally in the calculations.
Instead, we should compute the total income for each county, add up all the counties' totals, and then divide by the number of people in all the counties.
If we completed these steps with the `county` data, we would find that the per capita income for the US is \$30,861.
Had we computed the *simple* mean of per capita income across counties, the result would have been just \$26,093!

This example used what is called a **weighted mean**\index{weighted mean}.
For more information on this topic, check out the following online supplement regarding [weighted means](https://www.openintro.org/go/?id=stat_extra_weighted_mean).
:::

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "weighted mean")
```

## Histograms and shape {#sec-histograms}

Dot plots show the exact value for each observation.
They are useful for small datasets but can become hard to read with larger samples.
Rather than showing the value of each observation, we prefer to think of the value as belonging to a *bin*.
For example, in the `loan50` dataset, we created a table of counts for the number of loans with interest rates between 5.0% and 7.5%, then the number of loans with rates between 7.5% and 10.0%, and so on.
Observations that fall on the boundary of a bin (e.g., 10.00%) are allocated to the lower bin.
The tabulation is shown in @tbl-binnedIntRateAmountTable, and the binned counts are plotted as bars in @fig-loan50IntRateHist into what is called a **histogram**\index{plot!histogram}\index{histogram}.
Note that the histogram resembles a more heavily binned version of the stacked dot plot shown in @fig-loan-int-rate-dotplot.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "histogram")
```

```{r}
#| label: tbl-binnedIntRateAmountTable
#| tbl-cap: Counts for the binned interest rate data.
#| tbl-pos: H
loan50 |>
  mutate(interest_rate_cat = cut(interest_rate, breaks = seq(5, 27.5, 2.5))) |>
  count(interest_rate_cat, name = "Count") |>
  separate(interest_rate_cat, into = c("lower", "upper"), sep = ",") |>
  mutate(
    lower = str_remove(lower, "\\("),
    upper = str_remove(upper, "]"),
    lower = paste0(lower, "%"),
    upper = paste0(upper, "%"),
    lower = str_c("(", lower),
    upper = str_c(upper, "]")
  ) |>
  unite("Interest rate", lower:upper, sep = " - ") |>
  kbl(linesep = "", booktabs = TRUE, align = "lr") |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped"),
    full_width = FALSE
  ) |>
  column_spec(1:2, width = "9em")
```

```{r}
#| label: fig-loan50IntRateHist
#| fig-cap: |
#|   A histogram of interest rate. This distribution is strongly skewed
#|   to the right.
#| fig-alt: |
#|   A histogram of interest rate (ranging from about 5% to 25%).
#|   The distribution is right skewed.
#| fig-asp: 0.4
ggplot(loan50, aes(x = interest_rate)) +
  geom_histogram(breaks = seq(5, 27.5, 2.5)) +
  labs(x = "Interest rate", y = "Count") +
  scale_x_continuous(
    breaks = seq(5, 25, 5),
    labels = label_percent(scale = 1, accuracy = 1)
  )
```

Histograms provide a view of the **data density**\index{data density}.
Higher bars represent where the data are relatively more common.
For instance, there are many more loans with rates between 5% and 10% than loans with rates between 20% and 25% in the dataset.
The bars make it easy to see how the density of the data changes relative to the interest rate.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "data density")
```

Histograms are especially convenient for understanding the shape of the data distribution.
@fig-loan50IntRateHist suggests that most loans have rates under 15%, while only a handful of loans have rates above 20%.
When the distribution of a variable trails off to the right in this way and has a longer right **tail**\index{tail}, the shape is said to be **right skewed**\index{right skewed}.[^05-explore-numerical-5]

[^05-explore-numerical-5]: Other ways to describe data that are right skewed: skewed to the right, skewed to the high end, or skewed to the positive end.

```{r}
#| label: fig-loan50IntRateDensity
#| fig-cap: |
#|   A density plot of interest rate. Again, the distribution is strongly skewed
#|   to the right.
#| fig-alt: |
#|   A density plot of interest rate (ranging from about 5% to 25%).
#|   The distribution is right skewed.
#| fig-asp: 0.4
ggplot(loan50, aes(x = interest_rate)) +
  geom_density(fill = IMSCOL[1, 1], alpha = 0.5) +
  labs(x = "Interest rate", y = "Density") +
  scale_x_continuous(
    breaks = seq(5, 25, 5),
    labels = label_percent(scale = 1, accuracy = 1)
  )
```

@fig-loan50IntRateDensity shows a **density plot**\index{plot!density}\index{density plot} which is a smoothed out histogram.
The technical details for how to draw density plots (precisely how to smooth out the histogram) are beyond the scope of this text, but you will note that the shape, scale, and spread of the observations are displayed similarly in a histogram as in a density plot.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "density plot")
```

Variables with the reverse characteristic -- a long, thinner tail to the left -- are said to be **left skewed**\index{left skewed}.
We also say that such a distribution has a long left tail.
Variables that show roughly equal trailing off in both directions are called **symmetric**\index{symmetric}.

```{r}
#| include: false
terms_chp_05 <- c(
  terms_chp_05,
  "tail",
  "right skewed",
  "left skewed",
  "symmetric"
)
```

::: {.important data-latex=""}
When data trail off in one direction, the distribution has a **long tail**\index{tail}.
If a distribution has a long left tail, it is left skewed.
If a distribution has a long right tail, it is right skewed.
:::

::: {.guidedpractice data-latex=""}
Besides the mean (since it was labeled), what can you see in the dot plot in @fig-loan-int-rate-dotplot that you cannot see in the histogram in @fig-loan50IntRateHist?[^05-explore-numerical-6]
:::

[^05-explore-numerical-6]: The interest rates for individual loans.

In addition to looking at whether a distribution is skewed or symmetric, histograms can be used to identify modes.
A **mode** is represented by a prominent peak in the distribution.
There is only one prominent peak in the histogram of `interest_rate`.

A definition of *mode* sometimes taught in math classes is the value with the most occurrences in the dataset.
However, for many real-world datasets, it is common to have *no* observations with the same value in a dataset, making this definition impractical in data analysis.

\clearpage

@fig-singleBiMultiModalPlots shows histograms that have one, two, or three prominent peaks.
Such distributions are called **unimodal**\index{unimodal}, **bimodal**\index{bimodal}, and **multimodal**\index{multimodal}, respectively.
Any distribution with more than two prominent peaks is called multimodal.
Notice that there was one prominent peak in the unimodal distribution with a second less prominent peak that was not counted since it only differs from its neighboring bins by a few observations.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "unimodal", "bimodal", "multimodal")
```

```{r}
#| label: fig-singleBiMultiModalPlots
#| fig-cap: |
#|   Counting only prominent peaks, the distributions are (left to right)
#|   unimodal, bimodal, and multimodal. Note that the left plot is unimodal
#|   because we are counting prominent peaks, not just any peak.
#| fig-alt: |
#|   Three separate histograms on fabricated data. The first histogram
#|   shows a unimodal distribution, the second histogram shows a bimodal
#|   distribution, and the third histogram shows a multimodal distribution.
#| fig-asp: 0.3
#| out-width: 95%
df_modes <- tibble(
  uni = rchisq(65, 6),
  bi = c(rchisq(25, 5.8), rnorm(40, 20, 2)),
  multi = c(rchisq(25, 3), rnorm(25, 15), rnorm(15, 25, 1.5))
)

p_uni <- ggplot(df_modes, aes(x = uni)) +
  geom_histogram(binwidth = 2) +
  labs(x = NULL, y = NULL) +
  ylim(0, 23) +
  xlim(0, 30)
p_bi <- ggplot(df_modes, aes(x = bi)) +
  geom_histogram(binwidth = 2) +
  labs(x = NULL, y = NULL) +
  ylim(0, 23) +
  xlim(0, 30)
p_multi <- ggplot(df_modes, aes(x = multi)) +
  geom_histogram(binwidth = 2) +
  labs(x = NULL, y = NULL) +
  ylim(0, 23) +
  xlim(0, 30)

p_uni + p_bi + p_multi
```

::: {.workedexample data-latex=""}
@fig-loan50IntRateHist reveals only one prominent mode in the interest rate.
Is the distribution unimodal, bimodal, or multimodal?

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

Remember that *uni* stands for 1 (think *uni*cycles), and *bi* stands for 2 (think *bi*cycles).
:::

::: {.guidedpractice data-latex=""}
Height measurements of young students and adult teachers at an elementary school were taken.
How many modes would you expect in this height dataset?[^05-explore-numerical-8]
:::

[^05-explore-numerical-8]: There might be two height groups visible in the dataset: the children (students) and the adults (teachers).
    That is, the data are probably bimodal.

Looking for modes isn't about finding a clear and correct answer about the number of modes in a distribution, which is why *prominent*\index{prominent} is not rigorously defined in this book.
The most important part of this examination is to better understand your data.

## Variance and standard deviation {#sec-variance-sd}

The mean was introduced as a method to describe the center of a variable, and **variability**\index{variability} in the data is also important.
Here, we introduce two measures of variability: the variance and the standard deviation.
Both of these are very useful in data analysis, even though their formulas are a bit tedious to calculate by hand.
The standard deviation is the easier of the two to comprehend, as it roughly describes how far away the typical observation is from the mean.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "variability")
```

We call the distance of an observation from its mean its **deviation**\index{deviation}.
Below are the deviations for the $1^{st},$ $2^{nd},$ $3^{rd},$ and $50^{th}$ observations in the `interest_rate` variable:

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "deviation")
```

```{r}
#| include: false
loan50_interest_rate_deviation <- round(
  loan50$interest_rate - loan50_interest_rate_mean,
  2
)
loan50_interest_rate_deviation_squared <- round(
  loan50_interest_rate_deviation^2,
  2
)
loan50_interest_rate_var <- var(loan50$interest_rate)
loan50_interest_rate_sd <- sd(loan50$interest_rate)
```

$$
\begin{aligned}
x_1 - \bar{x} &= `r round(loan50$interest_rate[1], 2)` - `r round(loan50_interest_rate_mean, 2)` = `r loan50_interest_rate_deviation[1]` \\
x_2 - \bar{x} &= `r round(loan50$interest_rate[2], 2)` - `r round(loan50_interest_rate_mean, 2)` = `r loan50_interest_rate_deviation[2]` \\
x_3 - \bar{x} &= `r round(loan50$interest_rate[3], 2)` - `r round(loan50_interest_rate_mean, 2)` = `r loan50_interest_rate_deviation[3]` \\
&\vdots \\
x_{50} - \bar{x} &= `r round(loan50$interest_rate[50], 2)` - `r round(loan50_interest_rate_mean, 2)` = `r loan50_interest_rate_deviation[50]` \\
\end{aligned}
$$

\clearpage

If we square these deviations and then take an average, the result is equal to the sample **variance**\index{variance}, denoted by $s^2$:

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "variance")
```

$$
s^2 = \frac{(`r loan50_interest_rate_deviation[1]`)^2 + (`r loan50_interest_rate_deviation[2]`)^2 + (`r loan50_interest_rate_deviation[3]`)^2 + \cdots + (`r loan50_interest_rate_deviation[50]`)^2}{50 - 1} = \frac{`r loan50_interest_rate_deviation_squared[1]` + `r loan50_interest_rate_deviation_squared[2]` + \cdots + `r loan50_interest_rate_deviation_squared[50]`}{49} = `r round(loan50_interest_rate_var, 2)`
$$

We divide by $n - 1,$ rather than dividing by $n,$ when computing a sample's variance.
There's some mathematical nuance here, but the end result is that doing this makes this statistic slightly more reliable and useful.

Notice that squaring the deviations does two things.
First, it makes large values relatively much larger.
Second, it gets rid of any negative signs.

::: {.important data-latex=""}
**Standard deviation.**

The sample standard deviation can be calculated as the square root of the sum of the squared distance of each value from the mean divided by the number of observations minus one:

\vspace{-5mm}

$$s = \sqrt{\frac{\sum_{i=1}^n (x_i - \bar{x})^2}{n-1}}$$
:::

The **standard deviation**\index{standard deviation} is defined as the square root of the variance:

$$s = \sqrt{`r round(loan50_interest_rate_var, 2)`} = `r round(loan50_interest_rate_sd, 2)`$$

While often omitted, a subscript of $_x$ may be added to the variance and standard deviation, i.e., $s_x^2$ and $s_x^{},$ if it is useful as a reminder that these are the variance and standard deviation of the observations represented by $x_1,$ $x_2,$ ..., $x_n.$

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "standard deviation")
```

::: {.important data-latex=""}
**Variance and standard deviation.**

The variance is the average squared distance from the mean.
The standard deviation is the square root of the variance.
The standard deviation is useful when considering how far the data are distributed from the mean.

The standard deviation represents the typical deviation of observations from the mean.
Often about 68% of the data will be within one standard deviation of the mean and about 95% will be within two standard deviations.
However, these percentages are not strict rules.
:::

Like the mean, the population values for variance and standard deviation have special symbols: $\sigma^2$ for the variance and $\sigma$ for the standard deviation.

::: {.content-visible when-format="html"}
::: {.pronunciation data-latex=""}
The Greek letter $\sigma$ is pronounced *sigma*, listen to the pronunciation [here](https://youtu.be/PStgY5AcEIw?t=72).
:::
:::

::: {.content-visible when-format="pdf"}
::: {.pronunciation data-latex=""}
The Greek letter $\sigma$ is pronounced *sigma*.
:::
:::

::: {.guidedpractice data-latex=""}
A good description of the shape of a distribution should include modality and whether the distribution is symmetric or skewed to one side.
Using @fig-severalDiffDistWithSdOf1 as an example, explain why such a description is important.[^05-explore-numerical-9]
:::

[^05-explore-numerical-9]: @fig-severalDiffDistWithSdOf1 shows three distributions that look quite different, but all have the same mean, variance, and standard deviation.
    Using modality, we can distinguish between the first plot (bimodal) and the last two (unimodal).
    Using skewness, we can distinguish between the last plot (right skewed) and the first two.
    While a picture, like a histogram, tells a more complete story, we can use modality and shape (symmetry/skew) to characterize basic information about a distribution.

\clearpage

```{r}
#| label: sdRuleForIntRate
#| fig-cap: |
#|   For the interest rate variable, 34 of the 50 loans (68%) had
#|   interest rates within 1 standard deviation of the mean, and 48 of the 50
#|   loans (96%) had rates within 2 standard deviations. Usually about 68% of
#|   the data are within 1 standard deviation of the mean and 95% within 2
#|   standard deviations, though this is far from a hard rule.
#| fig-alt: |
#|   A histogram of the right-skewed variable, interest rates. Grey
#|   shading indicates the empirical rule cutoffs of one and two standard
#|   deviations. Because the variable is right-skewed, the empirical rule
#|   ranges do not capture the desired amoung of data, 68% and 95% respectively.
#| fig-asp: 0.5
box_sd1 <- tibble(
  x = c(
    loan50_interest_rate_mean - loan50_interest_rate_sd,
    loan50_interest_rate_mean - loan50_interest_rate_sd,
    loan50_interest_rate_mean + loan50_interest_rate_sd,
    loan50_interest_rate_mean + loan50_interest_rate_sd
  ),
  y = c(0, 17, 17, 0)
)

box_sd2 <- tibble(
  x = c(
    loan50_interest_rate_mean - 2 * loan50_interest_rate_sd,
    loan50_interest_rate_mean - 2 * loan50_interest_rate_sd,
    loan50_interest_rate_mean + 2 * loan50_interest_rate_sd,
    loan50_interest_rate_mean + 2 * loan50_interest_rate_sd
  ),
  y = c(0, 17, 17, 0)
)

box_sd3 <- tibble(
  x = c(
    loan50_interest_rate_mean - 3 * loan50_interest_rate_sd,
    loan50_interest_rate_mean - 3 * loan50_interest_rate_sd,
    loan50_interest_rate_mean + 3 * loan50_interest_rate_sd,
    loan50_interest_rate_mean + 3 * loan50_interest_rate_sd
  ),
  y = c(0, 17, 17, 0)
)

ggplot(loan50, aes(x = interest_rate)) +
  labs(x = "Interest rate", y = "Count") +
  geom_histogram(breaks = seq(5, 27.5, 2.5)) +
  geom_polygon(
    data = box_sd1,
    aes(x = x, y = y),
    fill = IMSCOL["gray", "f5"],
    alpha = 0.3
  ) +
  geom_polygon(
    data = box_sd2,
    aes(x = x, y = y),
    fill = IMSCOL["gray", "f5"],
    alpha = 0.3
  ) +
  geom_polygon(
    data = box_sd3,
    aes(x = x, y = y),
    fill = IMSCOL["gray", "f5"],
    alpha = 0.3
  ) +
  scale_x_continuous(
    breaks = seq(-5, 25, 5),
    labels = label_percent(scale = 1, accuracy = 1)
  ) +
  theme(
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank()
  )
```

```{r}
#| label: fig-severalDiffDistWithSdOf1
#| fig-cap: |
#|   Three different population distributions with the same mean (0)
#|   and standard deviation (1).
#| fig-alt: |
#|   Three histograms of fabricated data with empirical rule shading of
#|   one and two standard deviations superimposed. The bimodal and right
#|   skewed histogram demonstrate that the empirical rule is not a good
#|   approximation. The empirical rule fits well to the unimodal bell-shaped
#|   distribution.
#| fig-asp: 0.5
x1 <- c(rep(-0.975, 1000), rep(0.975, 1000))
x1 <- (x1 - mean(x1)) / sd(x1)
x2 <- rnorm(2000)
x2 <- (x2 - mean(x2)) / sd(x2)
x3 <- qchisq(seq(0.26, 0.8, 0.0005), 4)
x3 <- (x3 - mean(x3)) / sd(x3)

dists_mean0_sd1 <- tibble(
  x = c(x1, x2, x3),
  group = c(rep("A", length(x1)), rep("B", length(x2)), rep("C", length(x3)))
)

ggplot(dists_mean0_sd1, aes(x = x)) +
  geom_histogram(binwidth = 1) +
  facet_grid(group ~ ., scales = "free_y") +
  theme(
    # remove y axis
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    # strip facet labels
    strip.background = element_blank(),
    strip.text.y = element_blank()
  ) +
  scale_x_continuous(breaks = seq(-3, 3, 1)) +
  labs(x = NULL) +
  geom_polygon(
    data = tibble(x = c(-1, -1, 1, 1), y = c(0, 1000, 1000, 0)),
    aes(x = x, y = y),
    fill = IMSCOL["gray", "f5"],
    alpha = 0.3
  ) +
  geom_polygon(
    data = tibble(x = c(-2, -2, 2, 2), y = c(0, 1000, 1000, 0)),
    aes(x = x, y = y),
    fill = IMSCOL["gray", "f5"],
    alpha = 0.3
  ) +
  geom_polygon(
    data = tibble(x = c(-3, -3, 3, 3), y = c(0, 1000, 1000, 0)),
    aes(x = x, y = y),
    fill = IMSCOL["gray", "f5"],
    alpha = 0.3
  ) +
  theme(
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank()
  )
```

::: {.workedexample data-latex=""}
Describe the distribution of the `interest_rate` variable using the histogram in @fig-loan50IntRateHist.
The description should incorporate the center, variability, and shape of the distribution, and it should also be placed in context.
Also note any especially unusual cases.

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

The distribution of interest rates is unimodal and skewed to the high end.
Many of the rates fall near the mean at 11.57%, and most fall within one standard deviation (5.05%) of the mean.
There are a few exceptionally large interest rates in the sample that are above 20%.
:::

In practice, the variance and standard deviation are sometimes used as a means to an end, where the "end" is being able to accurately estimate the uncertainty associated with a sample statistic.
For example, in @sec-foundations-mathematical the standard deviation is used in calculations that help us understand how much a sample mean varies from one sample to the next.

\clearpage

## Box plots, quartiles, and the median {#sec-boxplots}

A **box plot**\index{plot!box}\index{box plot} summarizes a dataset using five statistics while also identifying unusual observations.
@fig-loan-int-rate-boxplot-dotplot provides a dot plot and a box plot of the `interest_rate` variable from the `loan50` dataset.^[Box plots were introducted by Mary Eleanor Spear who considered them to be a particular type of bar plot, see page 166 of @Spear1952. Mistakenly, box plots are often attributed to John Tukey who was the first person to call them "box-and-whisker plots."]

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "box plot")
```

```{r}
#| label: fig-loan-int-rate-boxplot-dotplot
#| fig-cap: Distribution of interest rates from the `loan50` dataset.
#| fig-subcap:
#|   - Dot plot
#|   - Box plot
#| fig-alt: |
#|   Two plots showing the same variable, interests rates. The top image is a
#|   dot plot, the bottom image is a box plot. Both images re-establish that
#|   the distribution of interest rates are right skewed.
#| fig-asp: 0.3
ggplot(loan50, aes(x = interest_rate)) +
  geom_dotplot(binwidth = 1, dotsize = 0.6) +
  labs(x = NULL, y = NULL) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  scale_x_continuous(
    labels = label_percent(scale = 1, accuracy = 1),
    limits = c(0, 30)
  )

ggplot(loan50, aes(x = interest_rate)) +
  geom_boxplot(outlier.size = 2.5) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  labs(x = "Interest rate") +
  scale_x_continuous(
    labels = label_percent(scale = 1, accuracy = 1),
    limits = c(0, 30)
  )
```

The dark line inside the box represents the **median**\index{median}, which splits the data in half.
50% of the data fall below this value and 50% fall above it.
Since in the `loan50` dataset there are `r nrow(loan50)` observations (an even number), the median is defined as the average of the two observations closest to the $50^{th}$ percentile.
@tbl-loan50-int-rate-sorted shows all interest rates, arranged in ascending order.
We can see that the $25^{th}$ and the $26^{th}$ values are both `r sort(loan50$interest_rate)[25]`, which corresponds to the thick line in  @fig-loan-int-rate-boxplot-dotplot-2.

```{r}
#| label: tbl-loan50-int-rate-sorted
#| tbl-cap: |
#|   Interest rates from the `loan50` dataset, arranged in ascending order.
#| tbl-pos: H
loan50_interest_rate_sorted <- sort(loan50$interest_rate)
matrix(
  c(
    loan50_interest_rate_sorted[1:10],
    loan50_interest_rate_sorted[11:20],
    loan50_interest_rate_sorted[21:30],
    loan50_interest_rate_sorted[31:40],
    loan50_interest_rate_sorted[41:50]
  ),
  byrow = TRUE,
  ncol = 10
) |>
  data.frame(row.names = c("1", "10", "20", "30", "40")) |>
  kbl(
    linesep = "",
    booktabs = TRUE,
    col.names = as.character(1:10)
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped")
  )
```

When there are an odd number of observations, there will be exactly one observation that splits the data into two halves, and in such a case that observation is the median (no average needed).

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "median")
```

::: {.important data-latex=""}
**Median: the number in the middle.**

If the data are ordered from smallest to largest, the **median** is the observation right in the middle.
If there are an even number of observations, there will be two values in the middle, and the median is taken as their average.
:::

\clearpage

The second step in building a box plot is drawing a rectangle to represent the middle 50% of the data.
The length of the box is called the **interquartile range**\index{interquartile range}, or **IQR**\index{IQR} for short.
It, like the standard deviation, is a measure of \index{variability}variability in data.
The more variable the data, the larger the standard deviation and IQR tend to be.
The two boundaries of the box are called the **first quartile**\index{first quartile} (the $25^{th}$ percentile, i.e., 25% of the data fall below this value) and the **third quartile**\index{third quartile} (the $75^{th}$ percentile, i.e., 75% of the data fall below this value), and these are often labeled $Q_1$ and $Q_3,$ respectively.

```{r}
#| include: false
terms_chp_05 <- c(
  terms_chp_05,
  "interquartile range",
  "IQR",
  "first quartile",
  "third quartile",
  "percentile"
)
```

::: {.important data-latex=""}
**Interquartile range (IQR).**

The IQR interquartile range is the length of the box in a box plot.
It is computed as $IQR = Q_3 - Q_1,$ where $Q_1$ and $Q_3$ are the $25^{th}$ and $75^{th}$ percentiles, respectively.

A $\alpha$ **percentile**\index{percentile} is a number with $\alpha$% of the observations below and $100-\alpha$% of the observations above.
For example, the $90^{th}$ percentile of SAT scores is the value of the SAT score with 90% of students below that value and 10% of students above that value.
:::

::: {.guidedpractice data-latex=""}
What percent of the data fall between $Q_1$ and the median?
What percent is between the median and $Q_3$?[^05-explore-numerical-10]
:::

[^05-explore-numerical-10]: Since $Q_1$ and $Q_3$ capture the middle 50% of the data and the median splits the data in the middle, 25% of the data fall between $Q_1$ and the median, and another 25% falls between the median and $Q_3.$

Extending out from the box, the **whiskers**\index{whisker} attempt to capture the data outside of the box.
The whiskers of a box plot reach to the minimum and the maximum values in the data, unless there are points that are considered unusually high or unusually low, which are identified as potential **outliers**\index{outlier} by the box plot.
These are labeled with a dot on the box plot.
The purpose of labeling the outlying points -- instead of extending the whiskers to the minimum and maximum observed values -- is to help identify any observations that appear to be unusually distant from the rest of the data.
There are a variety of formulas for determining whether a particular data point is considered an outlier, and different statistical software use different formulas.
A commonly used formula is that any observation beyond $1.5\times IQR$ away from the first or the third quartile is considered an outlier.
In a sense, the box is like the body of the box plot and the whiskers are like its arms trying to reach the rest of the data, up to the outliers.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "outlier", "whiskers")
```

::: {.important data-latex=""}
**Outliers are extreme.**

An **outlier** is an observation that appears extreme relative to the rest of the data.
Examining data for outliers serves many useful purposes, including

-   identifying strong skew \index{strong skew} in the distribution,
-   identifying possible data collection or data entry errors, and
-   providing insight into interesting properties of the data.

Keep in mind, however, that some datasets have a naturally long skew and outlying points do **not** represent any sort of problem in the dataset.
:::

::: {.guidedpractice data-latex=""}
Using the box plot in @fig-loan-int-rate-boxplot-dotplot-2, estimate the values of the $Q_1,$ $Q_3,$ and IQR for `interest_rate` in the `loan50` dataset.[^05-explore-numerical-11]
:::

[^05-explore-numerical-11]: These visual estimates will vary a little from one person to the next: $Q_1 \approx$ 8%, $Q_3 \approx$ 14%, IQR $\approx$ 14 - 8 = 6%.

\clearpage

## Robust statistics

How are the **sample statistics** \index{sample statistic} of the `interest_rate` dataset affected by the observation, 26.3%?
What would have happened if this loan had instead been only 15%?
What would happen to these summary statistics \index{summary statistic} if the observation at 26.3% had been even larger, say 35%?
The three conjectured scenarios are plotted alongside the original data in @fig-loan-int-rate-robust-ex, and sample statistics are computed under each scenario in @tbl-robustOrNotTable.

```{r}
#| label: tbl-robustOrNotTable
#| tbl-cap: |
#|   A comparison of how the median, IQR, mean, and standard deviation
#|   change as the value of an extreme observation from the original interest
#|   data changes.
#| tbl-pos: H
loan50_original <- loan50 |>
  select(interest_rate) |>
  mutate(
    mark = if_else(interest_rate == 26.30, TRUE, FALSE)
  )

loan50_15 <- loan50 |>
  select(interest_rate) |>
  mutate(
    mark = if_else(interest_rate == 26.30, TRUE, FALSE),
    interest_rate = if_else(interest_rate == 26.30, 15, interest_rate)
  )

loan50_35 <- loan50 |>
  select(interest_rate) |>
  mutate(
    mark = if_else(interest_rate == 26.30, TRUE, FALSE),
    interest_rate = if_else(interest_rate == 26.30, 35, interest_rate)
  )

loan50_robust_check <- bind_rows(
  loan50_original |> mutate(Scenario = "Original data"),
  loan50_15 |> mutate(Scenario = "Move 26.3% to 15%"),
  loan50_35 |> mutate(Scenario = "Move 26.3% to 35%")
) |>
  mutate(
    Scenario = fct_relevel(
      Scenario,
      "Original data",
      "Move 26.3% to 15%",
      "Move 26.3% to 35%"
    )
  )

loan50_robust_check |>
  group_by(Scenario) |>
  summarise(
    Median = median(interest_rate),
    IQR = IQR(interest_rate),
    Mean = mean(interest_rate),
    SD = sd(interest_rate)
  ) |>
  kbl(linesep = "", booktabs = TRUE, align = "lcccc") |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped")
  ) |>
  add_header_above(c(" " = 1, "Robust" = 2, "Not robust" = 2))
```

```{r}
#| label: fig-loan-int-rate-robust-ex
#| fig-cap: Dot plots of the original interest rate data and two modified datasets.
#| fig-subcap:
#|   - Original data
#|   - Move 26.3% to 15%
#|   - Move 26.3% to 35%
#| fig-alt: |
#|   Three dot plots are displayed where a single point is highlighted.
#|   In the top plot is the original data where the highlighted point is the
#|   largest value, at 26.3%. The second plot shows the distribution of the
#|   interest rate data after moving the highlighted point from 26.3% to 15%. The
#|   third plot shows the distribution of the interest rate data after moving the
#|   highlighted point from 26.3% to 35%. The corresponding change in statistics
#|   are given in @tbl--robustOrNotTable.
#| fig-asp: 0.3
#| fig-width: 5
ggplot(loan50_original, aes(x = interest_rate)) +
  geom_dotplot(binwidth = 1, dotsize = 0.7) +
  scale_x_continuous(
    labels = label_percent(scale = 1, accuracy = 1),
    limits = c(0, 35)
  ) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  gghighlight(mark) +
  labs(x = NULL)

ggplot(loan50_15, aes(x = interest_rate)) +
  geom_dotplot(binwidth = 1, dotsize = 0.7) +
  scale_x_continuous(
    labels = label_percent(scale = 1, accuracy = 1),
    limits = c(0, 35)
  ) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  gghighlight(mark) +
  labs(x = NULL)

ggplot(loan50_35, aes(x = interest_rate)) +
  geom_dotplot(binwidth = 1, dotsize = 0.7) +
  labs(x = "Interest rate") +
  scale_x_continuous(
    labels = label_percent(scale = 1, accuracy = 1),
    limits = c(0, 35)
  ) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) +
  gghighlight(mark) +
  labs(x = NULL)
```

\clearpage

::: {.guidedpractice data-latex=""}
Which is more affected by extreme observations, the mean or median?
Is the standard deviation or IQR more affected by extreme observations?[^05-explore-numerical-12]
:::

[^05-explore-numerical-12]: Mean is affected more than the median.
    Standard deviation is affected more than the IQR.

The median and IQR are called **robust statistics**\index{robust statistic} because extreme observations have little effect on their values: moving the most extreme value generally has little influence on these statistics.
On the other hand, the mean and standard deviation are more heavily influenced by changes in extreme observations, which can be important in some situations.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "robust statistics")
```

::: {.workedexample data-latex=""}
The median and IQR did not change under the three scenarios in @tbl-robustOrNotTable.
Why might this be the case?

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

The median and IQR are only sensitive to numbers near $Q_1,$ the median, and $Q_3.$ Since values in these regions are stable in the three datasets, the median and IQR estimates are also stable.
:::

You might not be surprised that the answer to the question "which is better, the mean or the median?" is: *it depends*. The two statistics measure different things, and so their use is dependent on the context in the analysis. Consider the following scenarios:

* Is it better to measure the average profit per customer or the median profit per customer?
  - If concern is about the overall profit margin of the company, the mean is a better measure to assess what is happening across the company. The company could have a positive median profit per customer and still be unprofitable.
  - If concern is around understanding the profit per typical customer, possibly to understand the growth headroom for the company's profit, the median profit per customer would tell you more about individual customer profits.
  
* If you operate an app and want to know how long it takes for the app to open on your customers' phones, do you want the mean amount of time or the median amount of time?
  - The mean leads to an understanding of the overall amount of time being wasted in opening the app.
  - The median tells you about the typical user experience.
  - However, if the app takes less than 5 milliseconds to launch for 50% of the users but more than 10 seconds to launch for 10% of the users, the median doesn't give the information you need. In that scenario, you might want an upper percentile, like the 95th percentile.

::: {.guidedpractice data-latex=""}
The distribution of loan amounts in the `loan50` dataset is right skewed, with a few large loans lingering out into the right tail.
If you were wanting to understand the typical loan size, should you be more interested in the mean or median?[^05-explore-numerical-13]
:::

[^05-explore-numerical-13]: If we are looking to simply understand what a typical individual loan looks like, the median is probably more useful.
    However, if the goal is to understand something that scales well, such as the total amount of money we might need to have on hand if we were to offer 1,000 loans, then the mean would be more useful.

Regardless of the choice of centrality statistic (either mean or median), for most analyses, it is important to consider more than just the centrality. Other statistics like upper and lower percentiles, IQR, or standard deviation provide information about the variability of the observations. And visualizing the data through a graphical representation will typically provide a wealth of information necessary for understanding the full data picture associated with the research question.

\clearpage

## Transforming data {#sec-transforming-data}

When data are very strongly skewed, we sometimes transform them, so they are easier to model.
@fig-county-unemployed-pop-transform-1 and @fig-county-unemployed-pop-transform-2 show right-skewed distributions: distribution of the percentage of unemployed people and the distribution of the population in all counties in the United States.
The distribution of population is more strongly skewed than the distribution of percentage unemployed, hence the log transformation results in a much bigger change in the shape of the distribution.

::: {.workedexample data-latex=""}
Consider the histogram of county populations shown in  @fig-county-unemployed-pop-transform-3, which shows extreme skew.
What characteristics of the plot keep it from being useful?

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

Nearly all of the data fall into the left-most bin, and the extreme skew obscures many of the potentially interesting details at the low values.
:::

```{r}
#| label: fig-county-unemployed-pop-transform
#| fig-cap: |
#|   Histograms of percentage unenmployed, population, and their log transformed
#|   versions in all US counties. For the plots of transformed variables, the
#|   x-value corresponds to the power of 10, e.g., 1 on the x-axis corresponds
#|   to $10^1 =$ 10 and 5 on the x-axis corresponds to $10^5=$ 100,000. Data
#|   are from 2017.
#| fig-subcap:
#|   - Percentage unemployed
#|   - log$_{10}$-transformed percentage unemployed
#|   - Population
#|   - log$_{10}$-transformed populations
#| fig-alt: |
#|   Two pairs of histograms are shown. In the first pair of histograms,
#|   percentage of unemployed in all US counties is shown next to the log-10
#|   transformed percentage of unemployed. The raw percent unemployed shows a
#|   right-skewed distribution, and the log-10 transformed data is symmetric and
#|   bell-shaped. In the second pair of histograms, population of each US county
#|   is shown next to the log-10 transformed population. The raw population is
#|   extremely right-skewed, and the log-10 transformed data is reasonably
#|   symmetric and bell-shaped.
#| out-width: 100%
#| layout-ncol: 2
#| fig-width: 5
#| fig-asp: 0.5
ggplot(county, aes(x = unemployment_rate)) +
  geom_histogram(binwidth = 1) +
  scale_x_continuous(labels = label_percent(scale = 1)) +
  labs(
    x = "Percentage unemployed",
    y = "Count"
  )

ggplot(county, aes(x = log(unemployment_rate, base = 10))) +
  geom_histogram(binwidth = 0.1) +
  labs(
    x = expression(log[10] * "(Percentage unemployed)"),
    y = "Count"
  )

ggplot(county_complete, aes(x = pop_2019)) +
  geom_histogram(binwidth = 500000) +
  scale_x_continuous(
    labels = label_number(scale = 0.000001, suffix = "M", accuracy = 1)
  ) +
  labs(
    x = "Population (M = millions)",
    y = "Count"
  )

ggplot(county_complete, aes(x = log(pop_2019, base = 10))) +
  geom_histogram(binwidth = 0.5) +
  labs(
    x = expression(log[10] * "(Population)"),
    y = "Count"
  )
```

There are some standard transformations that may be useful for strongly right skewed data where much of the data is positive but clustered near zero.
A **transformation**\index{transformation} is a rescaling of the data using a function.
For instance, a plot of the logarithm (base 10) of unemployment rates and county populations results in the new histograms in @fig-county-unemployed-pop-transform-2.
The transformed data are symmetric, and any potential outliers appear much less extreme than in the original dataset.
By reigning in the outliers and extreme skew, transformations often make it easier to build statistical models for the data.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "transformation")
```

Transformations can also be applied to one or both variables in a scatterplot.
A scatterplot of the population change from 2010 to 2017 against the population in 2010 is shown in @fig-county-pop-change-transform-1.
It's difficult to decipher any interesting patterns because the population variable is so strongly skewed.
However, if we apply a log$_{10}$ transformation to the population variable, as shown in @fig-county-pop-change-transform-2, a positive association between the variables is revealed.
In fact, we may be interested in fitting a trend line to the data when we explore methods around fitting regression lines in @sec-model-slr.

\clearpage

```{r}
#| label: fig-county-pop-change-transform
#| fig-cap: |
#|   Scatterplots of population change and log$_{10}$-transformed population
#|   change vs. population before change.
#| fig-subcap:
#|   - Population change
#|   - log$_{10}$-transformed population change
#| fig-alt: |
#|   A pair of scatterplots is given. In both plots the y-axis is
#|   population change. In the first plot, the x-axis is given by the raw variable
#|   population before change, and most of the points are quite small, obscuring
#|   any  x-y relationship. In the second plot, the x-axis is given by the
#|   transformed variable log-10 of the population before change. The relationship
#|   between log-10 population before change and population change is seen to be
#|   moderately strong, positive, and linear.
#| layout-ncol: 2
#| fig-width: 4
#| fig-asp: 0.8
ggplot(county, aes(y = pop_change, x = pop2010)) +
  geom_point(alpha = 0.5) +
  scale_x_continuous(
    labels = label_number(scale = 0.000001, suffix = "M", accuracy = 1)
  ) +
  labs(
    x = "Population before change\n(M = millions)",
    y = "Population change"
  )

ggplot(county, aes(y = pop_change, x = log(pop2010, base = 10))) +
  geom_point(alpha = 0.5) +
  labs(
    x = expression(log[10] * "(Population before change)"),
    y = "Population change"
  )
```

Transformations other than the logarithm can be useful, too.
For instance, the square root $(\sqrt{\text{original observation}})$ and inverse $\bigg ( \frac{1}{\text{original observation}} \bigg )$ are commonly used by data scientists.
Common goals in transforming data are to see the data structure differently, reduce skew, assist in modeling, or straighten a nonlinear relationship in a scatterplot.

## Mapping data

The `county` dataset offers many numerical variables that we could plot using dot plots, scatterplots, or box plots, but they can miss the true nature of the data as geographic.
When we encounter geographic data, we should create an **intensity map**\index{plot!intensity map}\index{intensity map}, where colors are used to show higher and lower values of a variable.
@fig-county-intensity-maps shows intensity maps for poverty rate in percent (`poverty`), unemployment rate in percent (`unemployment_rate`), homeownership rate in percent (`homeownership`), and median household income in \$1000s (`median_hh_income`).
The color key indicates which colors correspond to which values.
The intensity maps are not generally very helpful for getting precise values in any given county, but they are very helpful for seeing geographic trends and generating interesting research questions or hypotheses.

```{r}
#| include: false
terms_chp_05 <- c(terms_chp_05, "intensity map")
```

::: {.workedexample data-latex=""}
What interesting features are evident in the poverty and unemployment rate intensity maps in @fig-county-intensity-maps-1 and @fig-county-intensity-maps-2?

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

Poverty rates are evidently higher in a few locations.
Notably, the deep south shows higher poverty rates, as does much of Arizona and New Mexico.
High poverty rates are evident in the Mississippi flood plains a little north of New Orleans and in a large section of Kentucky.

The unemployment rate follows similar trends, and we can see correspondence between the two variables.
In fact, it makes sense for higher rates of unemployment to be closely related to poverty rates.
One observation that stands out when comparing the two maps: the poverty rate is much higher than the unemployment rate, meaning while many people may be working, they are not making enough to break out of poverty.
:::

::: {.guidedpractice data-latex=""}
What interesting features are evident in the median household income intensity map in @fig-county-intensity-maps-4?[^05-explore-numerical-14]
:::

[^05-explore-numerical-14]: Answers will vary.
    There is some correspondence between high earning and metropolitan areas, where we can see darker spots (higher median household income), though there are several exceptions.
    You might look for large cities you are familiar with and try to spot them on the map as dark spots.

\clearpage

```{r}
#| label: county-map-prep
dfips <- maps::county.fips |>
  as_tibble() |>
  extract(polyname, c("region", "subregion"), "^([^,]+),([^,]+)$")

map_county <- map_data("county") |>
  as_tibble() |>
  left_join(dfips) |>
  mutate(
    fips = case_when(
      subregion == "okaloosa" & region == "florida" ~ 12091L,
      subregion == "st martin" & region == "louisiana" ~ 22099L,
      subregion == "currituck" & region == "north carolina" ~ 37053L,
      # Oglala Lakota Count, see https://en.wikipedia.org/wiki/Oglala_Lakota_County,_South_Dakota
      subregion == "shannon" & region == "south dakota" ~ 46113L,
      subregion == "galveston" & region == "texas" ~ 48167L,
      subregion == "accomack" & region == "virginia" ~ 51001L,
      subregion == "pierce" & region == "washington" ~ 53053L,
      subregion == "san juan" & region == "washington" ~ 53055L,
      TRUE ~ fips
    )
  )

county_for_map <- county_complete |>
  select(
    fips,
    name,
    state,
    poverty_2017,
    unemployment_rate_2017,
    homeownership_2010,
    median_household_income_2017
  )

map_county <- map_county |>
  left_join(county_for_map, by = "fips")
```

```{r}
#| label: fig-county-intensity-maps
#| fig-cap: Intensity maps of US counties
#| fig-subcap:
#|   - Poverty rate
#|   - Unemployment rate
#|   - Homeownership rate
#|   - Median household income, in thousands of USD
#| fig-alt: |
#|   Four intensity map (also called heatmaps) of the United States.
#|   Each county is colored by the value of a variables. In the first plot, the
#|   county is colored by the poverty rate which is higher in the south and parts
#|   of the midwest. In the second plot, the county is colored by the unemployment
#|   rate which seems higher in the southwest and parts of the post-industrial
#|   Midwest. In the third plot, the county is colored by home ownership rate
#|   which is lowest on the west coast. In the fourth and last plot, the county
#|   is colored by median household income which is highest in the Bay Area
#|   and along the New England coast.
#| fig-asp: 0.42
ggplot(map_county, aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = poverty_2017)) +
  scale_fill_viridis_c(option = "D", labels = label_percent(scale = 1)) +
  labs(x = NULL, y = NULL, fill = "Poverty") +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  coord_quickmap()

ggplot(map_county, aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = unemployment_rate_2017)) +
  scale_fill_viridis_c(
    option = "A",
    labels = label_percent(scale = 1, accuracy = 1)
  ) +
  labs(x = NULL, y = NULL, fill = "Unemployment\nrate") +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  coord_quickmap()

ggplot(map_county, aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = homeownership_2010)) +
  scale_fill_viridis_c(option = "B", labels = label_percent(scale = 1)) +
  labs(x = NULL, y = NULL, fill = "Homeownership\nrate") +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  coord_quickmap()

ggplot(map_county, aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = median_household_income_2017)) +
  scale_fill_viridis_c(
    option = "C",
    labels = label_dollar(scale = 0.001, suffix = "K")
  ) +
  labs(x = NULL, y = NULL, fill = "Median\nhousehold\nincome") +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  coord_quickmap()
```

\clearpage

## Chapter review {#sec-chp5-review}

### Summary

Fluently working with numerical variables is an important skill for data analysts.
In this chapter we have introduced different visualizations and numerical summaries applied to numeric variables.
The graphical visualizations are even more descriptive when two variables are presented simultaneously.
We presented scatterplots, dot plots, histograms, and box plots.
Numerical variables can be summarized using the mean, median, quartiles, standard deviation, and variance.

### Terms

The terms introduced in this chapter are presented in @tbl-terms-chp-05.
If you're not sure what some of these terms mean, we recommend you go back in the text and review their definitions.
You should be able to easily spot them as **bolded text**.

```{r}
#| label: tbl-terms-chp-05
#| tbl-cap: Terms introduced in this chapter.
#| tbl-pos: H
make_terms_table(terms_chp_05)
```

\clearpage

## Exercises {#sec-chp5-exercises}

Answers to odd-numbered exercises can be found in [Appendix -@sec-exercise-solutions-05].

::: {.exercises data-latex=""}
{{< include exercises/_05-ex-explore-numerical.qmd >}}
:::