foundations-mathematical.qmd

# Inference with mathematical models {#sec-foundations-mathematical}

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

```{r}
#| include: false
opts_knit$set(global.par = TRUE)
par(mar = c(2, 0, 0, 0))
```

::: {.chapterintro data-latex=""}
In [Chapter -@sec-foundations-randomization] and [Chapter -@sec-foundations-bootstrapping] questions about population parameters were addressed using computational techniques.
With randomization tests, the data were permuted assuming the null hypothesis.
With bootstrapping, the data were resampled in order to measure the variability.
In many cases (indeed, with sample proportions), the variability of the statistic can be described by the computational method (as in previous chapters) or by a mathematical formula (as in this chapter).

The normal distribution is presented here to describe the variability associated with sample proportions which are taken from either repeated samples or repeated experiments.
The normal distribution is quite powerful in that it describes the variability of many different statistics, and we will encounter the normal distribution throughout the remainder of the book.

For now, however, focus is on the parallels between how data can provide insight about a research question either through computational methods or through mathematical models.
:::

## Central Limit Theorem {#sec-CLTsection}

In recent chapters, we have encountered four case studies.
While they differ in the settings, in their outcomes, and in the technique we have used to analyze the data, they all have something in common: the general shape of the distribution of the statistics (called the **sampling distribution**).\index{sampling distribution} You may have noticed that the distributions were symmetric and bell-shaped.

```{r}
#| include: false
terms_chp_13 <- c("sampling distribution", "null distribution")
```

::: {.important data-latex=""}
**Sampling distribution.**

A sampling distribution is the distribution of all possible values of a *sample statistic* from samples of a given sample size from a given population.
We can think about the sample distribution as describing how sample statistics (e.g., the sample proportion $\hat{p}$ or the sample mean $\bar{x}$) varies from one study to another.
A sampling distribution is contrasted with a data distribution which shows the variability of the *observed* data values.
The data distribution can be visualized from the observations themselves.
However, because a sampling distribution describes sample statistics computed from many studies, it cannot be visualized directly from a single dataset.
Instead, we use either computational or mathematical structures to estimate the sampling distribution and hence to describe the expected variability of the sample statistic in repeated studies.
:::

@fig-FourCaseStudies shows the null distributions in each of the four case studies where we ran 10,000 simulations.
Note that the **null distribution**\index{null distribution} is the sampling distribution of the statistic created under the setting where the null hypothesis is true.
Therefore, the null distribution will always be centered at the value of the parameter given by the null hypothesis.
In the case of the opportunity cost study, which originally had just 1,000 simulations, we have included an additional 9,000 simulations.

```{r}
#| include: false
opportunity_cost_sim_dist <- opportunity_cost |>
  specify(response = decision, explanatory = group, success = "buy video") |>
  hypothesize(null = "independence") |>
  generate(reps = 10000, type = "permute") |>
  calculate(stat = "diff in props", order = c("treatment", "control"))

sex_discrimination_sim_dist <- sex_discrimination |>
  specify(response = decision, explanatory = sex, success = "promoted") |>
  hypothesize(null = "independence") |>
  generate(reps = 10000, type = "permute") |>
  calculate(stat = "diff in props", order = c("male", "female"))

medical_consultant_sim_dist <- tibble(stat = rbinom(10000, 62, 0.1) / 62)

tappers_listeners_sim_dist <- tibble(stat = rbinom(10000, 120, 0.5) / 120)
```

```{r}
#| label: fig-FourCaseStudies
#| fig-cap: |
#|   The null distribution for each of the four case studies presented previously.
#|   Note that the center of each distribution is given by the value of the parameter
#|   set in the null hypothesis.
#| fig-alt: |
#|   Four null hypothesis histograms for examples in the book including 
#|   opportunity cost, sex discrimination, medical consultant, and tappers 
#|   and listeners. Each histogram is centered at the null value of the 
#|   proportion which is 0, 0, 0.1, and 0.5, respectively.
#| fig-asp: 0.7
# opportunity cost -------------------------------------------------------------

p_oc <- ggplot(opportunity_cost_sim_dist, aes(x = stat)) +
  geom_histogram(binwidth = 0.013, fill = IMSCOL["blue", "full"]) +
  scale_x_continuous(breaks = seq(-0.3, 0.3, 0.1), labels = label_number(accuracy = 0.1), limits = c(-0.35, 0.35)) +
  theme(
    axis.text.y = element_blank(),
    panel.grid = element_blank()
  ) +
  labs(x = expression(H[0]:~ p[T] - p[C] == 0), y = NULL) +
  annotate("rect", xmin = -0.13, xmax = 0.13, ymin = 150, ymax = 450, fill = "white", alpha = 0.5) +
  annotate("text", x = 0, y = 300, label = "Opportunity cost", size = 6)

# sex discrimination ----------------------------------------------------------

p_gd <- ggplot(sex_discrimination_sim_dist, aes(x = stat)) +
  geom_histogram(binwidth = 0.04, fill = IMSCOL["green", "full"]) +
  scale_x_continuous(breaks = seq(-0.6, 0.6, 0.2), labels = label_number(accuracy = 0.1), limits = c(-0.6, 0.6)) +
  theme(
    axis.text.y = element_blank(),
    panel.grid = element_blank()
  ) +
  labs(x = expression(H[0]:~ p[M] - p[F] == 0), y = NULL) +
  annotate("rect", xmin = -0.18, xmax = 0.18, ymin = 250, ymax = 750, fill = "white", alpha = 0.5) +
  annotate("text", x = 0, y = 500, label = "Sex discrimination", size = 6)

# medical consultant -----------------------------------------------------------

p_mc <- ggplot(medical_consultant_sim_dist, aes(x = stat)) +
  geom_histogram(binwidth = 0.016, fill = IMSCOL["red", "full"]) +
  scale_x_continuous(breaks = seq(-0.1, 0.3, 0.05), labels = label_number(accuracy = 0.01), limits = c(-0.1, 0.30)) +
  theme(
    axis.text.y = element_blank(),
    panel.grid = element_blank()
  ) +
  labs(x = expression(H[0]:~ p == 0.1), y = NULL) +
  annotate("rect", xmin = 0, xmax = 0.2, ymin = 200, ymax = 600, fill = "white", alpha = 0.5) +
  annotate("text", x = 0.10, y = 400, label = "Medical consultant", size = 6)

# tappers and listeners --------------------------------------------------------

p_tl <- ggplot(tappers_listeners_sim_dist, aes(x = stat)) +
  geom_histogram(binwidth = 0.016, fill = IMSCOL["yellow", "full"]) +
  scale_x_continuous(breaks = seq(0.30, 0.70, 0.05), labels = label_number(accuracy = 0.01), limits = c(0.3, 0.7)) +
  theme(
    axis.text.y = element_blank(),
    panel.grid = element_blank()
  ) +
  labs(x = expression(H[0]:~ p == 0.5), y = NULL) +
  annotate("rect", xmin = 0.3, xmax = 0.7, ymin = 150, ymax = 450, fill = "white", alpha = 0.5) +
  annotate("text", x = 0.5, y = 300, label = "Tappers + listeners", size = 6)

(p_oc + p_gd) /
  (p_mc + p_tl)
```

::: {.guidedpractice data-latex=""}
Describe the shape of the distributions and note anything that you find interesting.[^13-foundations-mathematical-1]
:::

[^13-foundations-mathematical-1]: In general, the distributions are reasonably symmetric.
    The case study for the medical consultant is the only distribution with any evident skew (the distribution is skewed right).

The case study for the medical consultant is the only distribution with any evident skew.
As we observed in [Chapter -@sec-data-hello], it's common for distributions to be skewed or contain outliers.
However, the null distributions we have so far encountered have all looked somewhat similar and, for the most part, symmetric.
They all resemble a bell-shaped curve.
The bell-shaped curve similarity is not a coincidence, but rather, is guaranteed by mathematical theory.

::: {.important data-latex=""}
**Central Limit Theorem for proportions.**\index{Central Limit Theorem}

If we look at a proportion (or difference in proportions) and the scenario satisfies certain conditions, then the sample proportion (or difference in proportions) will appear to follow a bell-shaped curve called the *normal distribution*.
:::

```{r}
#| include: false
terms_chp_13 <- c(terms_chp_13, "Central Limit Theorem")
```

An example of a perfect normal distribution is shown in @fig-simpleNormal.
Imagine laying a normal curve over each of the four null distributions in @fig-FourCaseStudies.
While the mean (center) and standard deviation (width or spread) may change for each plot, the general shape remains roughly intact.

```{r}
#| label: fig-simpleNormal
#| fig-cap: A normal curve.
#| fig-alt: A normal curve centered at 0 with a standard deviation of one.
#| fig-asp: 0.4
#| out-width: 60%
#| fig-width: 5
par(mar = c(0, 0, 0, 0))
normTail(axes = FALSE)
```

Mathematical theory guarantees that if repeated samples are taken a sample proportion or a difference in sample proportions will follow something that resembles a normal distribution when certain conditions are met.
(Note: we typically only take **one** sample, but the mathematical model lets us know what to expect if we *had* taken repeated samples.) These conditions fall into two general categories describing the independence between observations and the need to take a sufficiently large sample size.

1.  Observations in the sample are **independent**.
    Independence is guaranteed when we take a random sample from a population.
    Independence can also be guaranteed if we randomly divide individuals into treatment and control groups.

2.  The sample is **large enough**.
    The sample size cannot be too small.
    What qualifies as "small" differs from one context to the next, and we'll provide suitable guidelines for proportions in [Chapter -@sec-inference-one-prop].

So far we have had no need for the normal distribution.
We've been able to answer our questions somewhat easily using simulation techniques.
However, soon this will change.
Simulating data can be non-trivial.
For example, some of the scenarios encountered in [Chapter -@sec-model-mlr] where we introduced regression models with multiple predictors would require complex simulations in order to make inferential conclusions.
Instead, the normal distribution and other distributions like it offer a general framework for statistical inference that applies to a very large number of settings.

::: {.important data-latex=""}
**Technical Conditions.**

In order for the normal approximation to describe the sampling distribution of the sample proportion as it varies from sample to sample, two conditions must hold.
If these conditions do not hold, it is unwise to use the normal distribution (and related concepts like Z scores, probabilities from the normal curve, etc.) for inferential analyses.

1.  **Independent observations**
2.  **Large enough sample:** For proportions, at least 10 expected successes and 10 expected failures in the sample.
:::

## Normal Distribution {#sec-normalDist}

Among all the distributions we see in statistics, one is overwhelmingly the most common.
The symmetric, unimodal, bell curve is ubiquitous throughout statistics.
It is so common that people know it as a variety of names including the **normal curve**\index{normal curve}, **normal model**\index{normal model}, or **normal distribution**\index{normal distribution}.[^13-foundations-mathematical-2]
Under certain conditions, sample proportions, sample means, and sample differences can be modeled using the normal distribution.
Additionally, some variables such as SAT scores and heights of US adult males closely follow the normal distribution.

[^13-foundations-mathematical-2]: It is also introduced as the Gaussian distribution after Frederic Gauss, the first person to formalize its mathematical expression.

```{r}
#| include: false
terms_chp_13 <- c(terms_chp_13, "normal curve", "normal model", "normal distribution")
```

::: {.important data-latex=""}
**Normal distribution facts.**

Distributions of many variables are nearly normal, but none are exactly normal.
Thus, the normal distribution, while not perfect for any single problem, is very useful for a variety of problems.
:::

\vspace{-2mm}

In this section, we will discuss the normal distribution in the context of data to become familiar with normal distribution techniques.
In the following sections and beyond, we'll move our discussion to focus on applying the normal distribution and other related distributions to model point estimates for hypothesis tests and for constructing confidence intervals.

### Normal distribution model

The normal distribution always describes a symmetric, unimodal, bell-shaped curve.
However, normal curves can look different depending on the details of the model.
Specifically, the normal model can be adjusted using two parameters: mean and standard deviation.
As you can probably guess, changing the mean shifts the bell curve to the left or right, while changing the standard deviation stretches or constricts the curve.
@fig-twoSampleNormals-1 shows the normal distribution with mean $0$ and standard deviation $1$ (which is commonly referred to as the **standard normal distribution**\index{standard normal distribution}).
A normal distribution with mean $19$ and standard deviation $4$ is shown on the right.
@fig-twoSampleNormalsStacked shows the same two normal distributions on the same axis.

\vspace{-3mm}

```{r}
#| label: fig-twoSampleNormals
#| fig-cap: |
#|   Two normal distributions with different centers and spreads.
#| fig-subcap:
#|   - Mean = 0, SD = 1
#|   - Mean = 19, SD = 4
#| fig-alt: |
#|   Two normal curves. The first one has a mean at zero with a standard
#|   deviation of one and is called the standard normal distribution. The 
#|   second one has a mean at 19 with a standard deviation of four.
#| fig-asp: 0.5
#| out-width: 100%
#| layout-ncol: 2
#| fig-width: 5
normals <- tibble(
  x = c(rnorm(100000, mean = 0, sd = 1), rnorm(100000, mean = 19, sd = 4)),
  group = c(rep(1, 100000), rep(2, 100000))
)

ggplot(normals |> filter(group == 1), aes(x = x)) +
  geom_histogram(aes(y = ..density..), alpha = 0.5) +
  geom_function(fun = dnorm, args = list(mean = 0, sd = 1), color = IMSCOL["blue", "full"], size = 1) +
  labs(y = NULL, x = NULL) +
  scale_x_continuous(breaks = -3:3) +
  theme(
    axis.text.y = element_blank(),
    panel.grid = element_blank()
  )

ggplot(normals |> filter(group == 2), aes(x = x)) +
  geom_histogram(aes(y = ..density..), alpha = 0.5, fill = IMSCOL["green", "full"]) +
  geom_function(
    fun = dnorm, args = list(mean = 19, sd = 4), color = IMSCOL["green", "full"], size = 1,
    linetype = "dashed"
  ) +
  labs(y = NULL, x = NULL) +
  scale_x_continuous(breaks = 19 + 4 * (-3:3)) +
  theme(
    axis.text.y = element_blank(),
    panel.grid = element_blank()
  )
```

```{r}
#| include: false
terms_chp_13 <- c(terms_chp_13, "standard normal distribution")
```

```{r}
#| label: fig-twoSampleNormalsStacked
#| fig-cap: |
#|   The two normal models shown in @fig-twoSampleNormals but plotted together
#|   on the same scale.
#| fig-alt: |
#|   Two normal curves plotted on the same axis. One is a
#|   standard normal curve with a mean at zero and a standard
#|   deviation of one. The other has a mean of 19 with a standard
#|   deviation of four.
#| fig-asp: 0.25
#| out-width: 100%

ggplot(normals, aes(x = x)) +
  geom_function(fun = dnorm, args = list(mean = 0, sd = 1), color = IMSCOL["blue", "full"], size = 1) +
  geom_function(fun = dnorm, args = list(mean = 19, sd = 4), color = IMSCOL["green", "full"], size = 1, linetype = "dashed") +
  labs(y = NULL, x = NULL) +
  theme(
    axis.text.y = element_blank(),
    panel.grid = element_blank()
  )
```

\clearpage

If a normal distribution has mean $\mu$ and standard deviation $\sigma,$ we may write the distribution as $N(\mu, \sigma).$ The two distributions in @fig-twoSampleNormalsStacked can be written as

$$ N(\mu = 0, \sigma = 1)\quad\text{and}\quad N(\mu = 19, \sigma = 4) $$

Because the mean and standard deviation describe a normal distribution exactly, they are called the distribution's **parameters**\index{parameter}.

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

::: {.workedexample data-latex=""}
Write down the short-hand for a normal distribution with the following parameters.

a.  mean 5 and standard deviation 3
b.  mean -100 and standard deviation 10
c.  mean 2 and standard deviation 9

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

a.  $N(\mu = 5,\sigma = 3)$
b.  $N(\mu = -100, \sigma = 10)$
c.  $N(\mu = 2, \sigma = 9)$
:::

### Standardizing with Z scores

::: {.guidedpractice data-latex=""}
SAT scores follow a nearly normal distribution with a mean of 1500 points and a standard deviation of 300 points.
ACT scores also follow a nearly normal distribution with mean of 21 points and a standard deviation of 5 points.
Suppose Nel scored 1800 points on their SAT and Sian scored 24 points on their ACT.
Who performed better?[^13-foundations-mathematical-3]
:::

[^13-foundations-mathematical-3]: We use the standard deviation as a guide.
    Nel is 1 standard deviation above average on the SAT: $1500 + 300 = 1800.$ Sian is 0.6 standard deviations above the mean on the ACT: $21+0.6 \times 5 = 24.$ In @fig-satActNormals, we can see that Nel did better compared to other test takers than Sian did, so their score was better.

```{r}
#| label: fig-satActNormals
#| fig-cap: |
#|   Nel's and Sian's scores shown with the distributions of SAT and ACT scores.
#| fig-alt: |
#|   Nel's SAT score is plotted on a normal curve with a
#|   mean of 1500 and a standard deviation of 300. Sian's ACT
#|   score is plotted on a normal curve with a mean of 21 and a
#|   standard deviation of 5. Nel's score is relatively higher
#|   than Sian's score.
#| fig-asp: 0.5
sat_mean <- 1500
sat_sd <- 300
act_mean <- 21
act_sd <- 5

tests <- tibble(
  score = c(rnorm(100000, mean = sat_mean, sd = sat_sd), rnorm(100000, mean = act_mean, sd = act_sd)),
  test  = c(rep("SAT", 100000), rep("ACT", 100000))
)

p_sat <- ggplot(tests |> filter(test == "SAT"), aes(x = score)) +
  geom_function(fun = dnorm, args = list(mean = sat_mean, sd = sat_sd), color = IMSCOL["black", "full"]) +
  labs(y = NULL, x = NULL) +
  theme(
    axis.text.y = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank()
  ) +
  scale_x_continuous(breaks = sat_mean + sat_sd * (-3:3)) +
  annotate("segment",
    x = 1800, xend = 1800, y = 0, yend = 0.001,
    color = IMSCOL["blue", "full"], linetype = "dashed", size = 1
  ) +
  annotate("segment", x = sat_mean - 4 * sat_sd, xend = sat_mean + 4 * sat_sd, y = 0, yend = 0, color = IMSCOL["black", "full"]) +
  annotate("text", x = 1800, y = 0.0011, label = "Nel", color = IMSCOL["blue", "full"]) +
  coord_cartesian(xlim = sat_mean + c(-1, 1) * 4 * sat_sd)

p_act <- ggplot(tests |> filter(test == "ACT"), aes(x = score)) +
  geom_function(
    fun = dnorm, args = list(mean = act_mean, sd = act_sd),
    color = IMSCOL["black", "full"]
  ) +
  labs(y = NULL, x = NULL) +
  theme(
    axis.text.y = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank()
  ) +
  scale_x_continuous(breaks = act_mean + act_sd * (-3:3)) +
  annotate("segment",
    x = 24, xend = 24, y = 0, yend = 0.075,
    color = IMSCOL["green", "full"], linetype = "dashed", size = 1
  ) +
  annotate("segment", x = act_mean - 4 * act_sd, xend = act_mean + 4 * act_sd, y = 0, yend = 0, color = IMSCOL["black", "full"]) +
  annotate("text", x = 24, y = 0.08, label = "Sian", color = IMSCOL["green", "full"]) +
  coord_cartesian(xlim = act_mean + c(-1, 1) * 4 * act_sd)

p_sat /
  p_act
```

\clearpage

The solution to the previous example relies on a standardization technique called a Z score, a method most commonly employed for nearly normal observations (but that may be used with any distribution).
The **Z score**\index{Z score} of an observation is defined as the number of standard deviations it falls above or below the mean.
If the observation is one standard deviation above the mean, its Z score is 1.
If it is 1.5 standard deviations *below* the mean, then its Z score is -1.5.
If $x$ is an observation from a distribution $N(\mu, \sigma),$ we define the Z score mathematically as

```{r}
#| include: false
terms_chp_13 <- c(terms_chp_13, "Z score")
```

$$ Z = \frac{x-\mu}{\sigma} $$

Using $\mu_{SAT}=1500,$ $\sigma_{SAT}=300,$ and $x_{Nel}=1800,$ we find Nel's Z score:

$$ Z_{Nel} = \frac{x_{Nel} - \mu_{SAT}}{\sigma_{SAT}} = \frac{1800-1500}{300} = 1 $$

::: {.important data-latex=""}
**The Z score.**

The Z score of an observation is the number of standard deviations it falls above or below the mean.
We compute the Z score for an observation $x$ that follows a distribution with mean $\mu$ and standard deviation $\sigma$ using

$$Z = \frac{x-\mu}{\sigma}$$

If the observation $x$ comes from a *normal* distribution centered at $\mu$ with standard deviation of $\sigma$, then the Z score will be distributed according to a *normal* distribution with a center of 0 and a standard deviation of 1.
That is, the normality remains when transforming from $x$ to $Z$ with a shift in both the center as well as the spread.
:::

::: {.guidedpractice data-latex=""}
Use Sian's ACT score, 24, along with the ACT mean and standard deviation to compute their Z score.[^13-foundations-mathematical-4]
:::

[^13-foundations-mathematical-4]: $Z_{Sian} = \frac{x_{Sian} - \mu_{ACT}}{\sigma_{ACT}} = \frac{24 - 21}{5} = 0.6$

Observations above the mean always have positive Z scores while those below the mean have negative Z scores.
If an observation is equal to the mean (e.g., SAT score of 1500), then the Z score is $0.$

::: {.workedexample data-latex=""}
Let $X$ represent a random variable from $N(\mu=3, \sigma=2),$ and suppose we observe $x=5.19.$ Find the Z score of $x.$ Then, use the Z score to determine how many standard deviations above or below the mean $x$ falls.

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

Its Z score is given by $Z = \frac{x-\mu}{\sigma} = \frac{5.19 - 3}{2} = 2.19/2 = 1.095.$ The observation $x$ is 1.095 standard deviations *above* the mean.
We know it must be above the mean since $Z$ is positive.
:::

::: {.guidedpractice data-latex=""}
Head lengths of brushtail possums follow a nearly normal distribution with mean 92.6 mm and standard deviation 3.6 mm.
Compute the Z scores for possums with head lengths of 95.4 mm and 85.8 mm.[^13-foundations-mathematical-5]
:::

[^13-foundations-mathematical-5]: For $x_1=95.4$ mm: $Z_1 = \frac{x_1 - \mu}{\sigma} = \frac{95.4 - 92.6}{3.6} = 0.78.$ For $x_2=85.8$ mm: $Z_2 = \frac{85.8 - 92.6}{3.6} = -1.89.$

We can use Z scores to roughly identify which observations are more unusual than others.
One observation $x_1$ is said to be more unusual than another observation $x_2$ if the absolute value of its Z score is larger than the absolute value of the other observation's Z score: $|Z_1| > |Z_2|.$ This technique is especially insightful when a distribution is symmetric.

::: {.guidedpractice data-latex=""}
Which of the two brushtail possum observations in the previous guided practice is more *unusual*?[^13-foundations-mathematical-6]
:::

[^13-foundations-mathematical-6]: Because the *absolute value* of Z score for the second observation is larger than that of the first, the second observation has a more unusual head length.

\vspace{-5mm}

### Normal probability calculations

::: {.workedexample data-latex=""}
Nel from the SAT Guided Practice earned a score of 1800 on their SAT with a corresponding $Z=1.$ They would like to know what percentile they fall in among all SAT test-takers.

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

Nel's **percentile**\index{percentile} is the percentage of people who earned a lower SAT score than Nel.
We shade the area representing those individuals in @fig-satBelow1800.
The total area under the normal curve is always equal to 1, and the proportion of people who scored below Nel on the SAT is equal to the *area* shaded in @fig-satBelow1800: 0.8413.
In other words, Nel is in the $84^{th}$ percentile of SAT takers.
:::

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

```{r}
#| label: fig-satBelow1800
#| fig-cap: |
#|   The normal model for SAT scores, shading the area of those individuals who
#|   scored below Nel.
#| fig-alt: |
#|   Nel's SAT score of 1800 is plotted on a normal curve
#|   with mean of 1500 and standard deviation of 300. The scores
#|   lower than Nel's are shaded in blue and represent roughly
#|   84% of the distribution.
#| fig-asp: 0.3
#| out-width: 60%
#| fig-width: 5
par(mar = c(2, 0, 0, 0))
normTail(m = 1500, s = 300, L = 1800, col = IMSCOL["blue", "full"])
```

We can use the normal model to find percentiles or probabilities.
A **normal probability table**\index{normal probability table}, which lists Z scores and corresponding percentiles, can be used to identify a percentile based on the Z score (and vice versa).
Statistical software can also be used.

```{r}
#| include: false
terms_chp_13 <- c(terms_chp_13, "normal probability table")
```

Normal probabilities are most commonly found using statistical software which we will show here using R.
We use the software to identify the percentile corresponding to any particular Z score.
For instance, the percentile of $Z=0.43$ is 0.6664, or the $66.64^{th}$ percentile.
The `pnorm()` function is available in default R and will provide the percentile associated with any cutoff on a normal curve.
The `normTail()` function is available in the [**openintro**](http://openintrostat.github.io/openintro/) R package and will draw the associated normal distribution curve.

```{r}
#| echo: true
#| fig-asp: 0.3
#| out-width: 60%
#| fig-width: 5
pnorm(0.43, mean = 0, sd = 1)
openintro::normTail(m = 0, s = 1, L = 0.43)
```

We can also find the Z score associated with a percentile.
For example, to identify Z for the $80^{th}$ percentile, we use `qnorm()` which identifies the **quantile** for a given percentage.
The quantile represents the cutoff value.
(To remember the function `qnorm()` as providing a cutoff, notice that both `qnorm()` and "cutoff" start with the sound "kuh".
To remember the `pnorm()` function as providing a probability from a given cutoff, notice that both `pnorm()` and probability start with the sound "puh".) We determine the Z score for the $80^{th}$ percentile using `qnorm()`: 0.84.

```{r}
#| echo: true
#| fig-asp: 0.3
#| out-width: 60%
#| fig-width: 5
qnorm(0.80, mean = 0, sd = 1)
openintro::normTail(m = 0, s = 1, L = 0.842)
```

::: {.guidedpractice data-latex=""}
Determine the proportion of SAT test takers who scored better than Nel on the SAT.[^13-foundations-mathematical-7]
:::

[^13-foundations-mathematical-7]: If 84% had lower scores than Nel, the number of people who had better scores must be 16%.
    (Generally ties are ignored when the normal model, or any other continuous distribution, is used.)

### Normal probability examples

Cumulative SAT scores are approximated well by a normal model, $N(\mu=1500, \sigma=300).$

::: {.workedexample data-latex=""}
Shannon is a randomly selected SAT taker, and nothing is known about Shannon's SAT aptitude.
What is the probability that Shannon scores at least 1630 on their SATs?

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

First, always draw and label a picture of the normal distribution.
(Drawings need not be exact to be useful.) We are interested in the chance they score above 1630, so we shade the upper tail.
See the normal curve below.

The $x$-axis identifies the mean and the values at 2 standard deviations above and below the mean.
The simplest way to find the shaded area under the curve makes use of the Z score of the cutoff value.
With $\mu=1500,$ $\sigma=300,$ and the cutoff value $x=1630,$ the Z score is computed as

$$ Z = \frac{x - \mu}{\sigma} = \frac{1630 - 1500}{300} = \frac{130}{300} = 0.43 $$

We use software to find the percentile of $Z=0.43,$ which yields 0.6664.
However, the percentile describes those who had a Z score *lower* than 0.43.
To find the area *above* $Z=0.43,$ we compute one minus the area of the lower tail, as seen below.

The probability Shannon scores at least 1630 on the SAT is 0.3336.
This calculation is visualized in @fig-subtractingArea.
:::

```{r}
#| label: fig-subtractingArea
#| fig-cap: |
#|   Visual calculation of the probability that Shannon scores at 
#|   least 1630 on the SAT.
#| fig-alt: |
#|   Three normal curves visualizing that probabilities can
#|   be calculated using arithmetic on areas. The area above
#|   1630 can be calculated by taking the total area of one and
#|   subtracting the area below 1630 (which is 0.6664). The
#|   resulting value is 0.3336.
#| fig-asp: 0.45
#| fig-width: 7
par_og <- par(no.readonly = TRUE) # save original par
par(mar = c(0, 0, 0, 0), mfrow = c(2, 1))

normTail(m = 1500, s = 300, U = 1630, col = IMSCOL["blue", "full"])

X <- seq(-3.2, 3.2, 0.01)
Y <- dnorm(X)
plot(X, Y, type = "l", axes = F, xlim = c(-3.4, 16 + 3.4), ylim = c(0, 0.652))
lines(X, rep(0, length(X)))
these <- which(X <= 8)
polygon(c(X[these[1]], X[these], X[rev(these)[1]]), c(0, Y[these], 0), col = IMSCOL["blue", "full"])
lines(X, Y)
# abline(h=0)
lines(c(0, 0), dnorm(0) * c(0.01, 0.99), col = COL[6], lty = 3)
lines(c(3, 8 - 3), c(0.2, 0.2), lwd = 3)
text(0, 0.45, format(c(1, 0.0001), scientific = FALSE, digits = 4)[1], cex = 1)
lines(X + 8, Y, type = "l", xlim = c(-3.4, 3.4))
lines(X + 8, rep(0, length(X)))
these <- which(X <= 0.43)
polygon(c(X[these[1]], X[these], X[rev(these)[1]]) + 8, c(0, Y[these], 0), col = IMSCOL["blue", "full"])
lines(X + 8, Y)
lines(c(0, 0), dnorm(0) * c(0.01, 0.99), col = COL[6], lty = 3)
lines(8 + c(3, 8 - 3), c(0.23, 0.23), lwd = 3)
lines(8 + c(3, 8 - 3), c(0.17, 0.17), lwd = 3)
lines(c(3.72, 4.28), rep(0.45, 2), lwd = 2)
text(8, 0.45, format(0.6664, scientific = FALSE, digits = 4)[1], cex = 1)
lines(X + 8 + 8, Y, type = "l", xlim = c(-3.4, 3.4))
lines(X + 8 + 8, rep(0, length(X)))
these <- which(X > 0.43)
polygon(c(X[these[1]], X[these], X[rev(these)[1]]) + 8 + 8, c(0, Y[these], 0), col = IMSCOL["blue", "full"])
lines(X + 8 + 8, Y)
lines(c(0, 0), dnorm(0) * c(0.01, 0.99), col = COL[6], lty = 3)
text(12, 0.45, "=", cex = 1.5)
text(16, 0.45, format(0.3336, scientific = FALSE, digits = 4)[1], cex = 1)
par(par_og) # restore original par
```

::: {.tip data-latex=""}
**Always draw a picture first, and find the Z score second.**

For any normal probability situation, *always always always* draw and label the normal curve and shade the area of interest first.
The picture will provide an estimate of the probability.

After drawing a figure to represent the situation, identify the Z score for the observation of interest.
:::

::: {.guidedpractice data-latex=""}
If the probability of Shannon scoring at least 1630 is 0.3336, then what is the probability they score less than 1630?
Draw the normal curve representing this exercise, shading the lower region instead of the upper one.[^13-foundations-mathematical-8]
:::

[^13-foundations-mathematical-8]: We found the probability to be 0.6664.
    A picture for this exercise is represented by the shaded area below "0.6664".

::: {.workedexample data-latex=""}
Edward earned a 1400 on their SAT.
What is their percentile?

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

First, a picture is needed.
Edward's percentile is the proportion of people who do not get as high as a 1400.
These are the scores to the left of 1400, as shown below.

```{r}
#| label: edward-percentile
#| fig-asp: 0.4
#| out-width: 60%
#| fig-align: center
#| fig-alt: |
#|   A normal curve with a mean of 1500 and a standard
#|   deviation of 300. Values less than 1400 are shaded blue in
#|   the lower part of the curve. The shaded part is roughly
#|   37 percent of the graph.
#| fig-width: 5
par(mar = c(2, 0, 0, 0))
normTail(m = 1500, s = 300, L = 1400, col = IMSCOL["blue", "full"])
```

The mean $\mu=1500,$ the standard deviation $\sigma=300,$ and the cutoff for the tail area $x=1400$ are used to compute the Z score:

$$ Z = \frac{x - \mu}{\sigma} = \frac{1400 - 1500}{300} = -0.33$$

Statistical software can be used to find the proportion of the $N(0,1)$ curve to the left of $-0.33$ which is 0.3707.
Edward is at the $37^{th}$ percentile.
:::

::: {.workedexample data-latex=""}
Use the results of the previous example to compute the proportion of SAT takers who did better than Edward.
Also draw a new picture.

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

If Edward did better than 37% of SAT takers, then about 63% must have done better than them, as shown below.

```{r}
#| label: edward-percentile-better
#| fig-asp: 0.4
#| out-width: 60%
#| fig-align: center
#| fig-width: 5
#| fig-alt: |
#|   A normal curve with a mean of 1500 and a standard deviation of 300.
#|   Values greater than 1400 are shaded blue in the upper part of the curve.
#|   The shaded part is roughly 63 percent of the graph.
par(mar = c(2, 0, 0, 0))
normTail(m = 1500, s = 300, U = 1400, col = IMSCOL["blue", "full"])
```
:::

::: {.tip data-latex=""}
**Areas to the right.**

Most statistical software, as well as normal probability tables in most books, give the area to the left.
If you would like the area to the right, first find the area to the left and then subtract the amount from one.
:::

::: {.guidedpractice data-latex=""}
Stuart earned an SAT score of 2100.
Draw a picture for each part.
(a) What is their percentile?
(b) What percent of SAT takers did better than Stuart?[^13-foundations-mathematical-9]
:::

[^13-foundations-mathematical-9]: Numerical answers: (a) 0.9772.
    (b) 0.0228.

Based on a sample of 100 men,[^13-foundations-mathematical-10] the heights of adults who identify as male, between the ages 20 and 62 in the US is nearly normal with mean 70.0" and standard deviation 3.3".

[^13-foundations-mathematical-10]: This sample was taken from the USDA Food Commodity Intake Database.

::: {.workedexample data-latex=""}
Kamron is 5'7" (67 inches) and Adrian is 6'4" (76 inches).
(a) What is Kamron's height percentile?
(b) What is Adrian's height percentile?
Also draw one picture for each part.

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

Numerical answers, calculated using statistical software (e.g., `pnorm()` in R): (a) `r round(pnorm(67, mean = 70, sd = 3.3), 4) * 100`th percentile.
(b) `r round(pnorm(76, mean = 70, sd = 3.3), 4) * 100`th percentile.
:::

The last several problems have focused on finding the probability or percentile for a particular observation.
What if you would like to know the observation corresponding to a particular percentile?

::: {.workedexample data-latex=""}
Yousef's height is at the $40^{th}$ percentile.
How tall are they?

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

As always, first draw the picture.

```{r}
#| fig-asp: 0.4
#| out-width: 60%
#| fig-align: center
#| fig-width: 5
par(mar = c(2, 0, 0, 0))
normTail(70, 3.3, L = qnorm(0.4, 70, 3.3), col = IMSCOL["blue", "full"])
text(67, 0.03, "40%\n(0.40)", cex = 1, col = IMSCOL["black", "full"])
```

In this case, the lower tail probability is known (0.40), which can be shaded on the diagram.
We want to find the observation that corresponds to the known probability of 0.4.
We can find the observation in two different ways: using the height curve seen above or using the Z score associated with the standard normal curve centered at zero with a standard deviation of one.

If you have access to software (like R, code seen below) that allows you to specify the mean and standard deviation of the normal curve, you can calculate the observed value on the curve (i.e., Yousef's height) directly.

```{r}
qnorm(0.4, mean = 70, sd = 3.3)
```

Yousef is 69.2 inches tall.
That is, Yousef is about 5'9" (this is notation for 5-feet, 9-inches).

Without access to flexible software, you will need the information given by a standard normal curve (a normal curve centered at zero with a standard deviation of one).
First, determine the Z score associated with the $40^{th}$ percentile.

Because the percentile is below 50%, we know $Z$ will be negative.
Statistical software provides the $Z$ value to be $-0.25.$

```{r}
#| echo: true
qnorm(0.4, mean = 0, sd = 1)
```

Knowing $Z_{Yousef}=-0.25$ and the population parameters $\mu=70$ and $\sigma=3.3$ inches, the Z score formula can be set up to determine Yousef's unknown height, labeled $x_{Yousef}$:

$$ -0.253 = Z_{Yousef} = \frac{x_{Yousef} - \mu}{\sigma} = \frac{x_{Yousef} - 70}{3.3} $$

Solving for $x_{Yousef}$ yields the height 69.2 inches.
Again, Yousef is about 5'9".
:::

::: {.workedexample data-latex=""}
What is the adult male height at the $82^{nd}$ percentile?

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

In order to practice using Z scores, we will use the standard normal curve to solve the problem.

Again, we draw the figure first.

```{r}
#| label: height82Perc
#| fig-asp: 0.4
#| out-width: 60%
#| fig-align: center
#| fig-width: 5
par(mar = c(2, 0, 0, 0))
normTail(70, 3.3, L = qnorm(0.82, 70, 3.3), col = IMSCOL["blue", "full"])
text(70, 0.04, "82%\n(0.82)", cex = 1, col = IMSCOL["black", "full"])
text(74.3, 0.02, "18%\n(0.18)", cex = 1, col = IMSCOL["black", "full"])
```

And calculate the Z value associated with the $82^{nd}$ percentile:

```{r}
#| echo: true
qnorm(0.82, mean = 0, sd = 1)
```

Next, we want to find the Z score at the $82^{nd}$ percentile, which will be a positive value (because the percentile is bigger than 50%).
Using `qnorm()`, the $82^{nd}$ percentile corresponds to $Z=0.92.$ Finally, the height $x$ is found using the Z score formula with the known mean $\mu,$ standard deviation $\sigma,$ and Z score $Z=0.92$:

$$ 0.92 = Z = \frac{x-\mu}{\sigma} = \frac{x - 70}{3.3} $$

This yields 73.04 inches or about 6'1" as the height at the $82^{nd}$ percentile.
:::

::: {.guidedpractice data-latex=""}
Using Z scores, answer the following questions.

(a) What is the $95^{th}$ percentile for SAT scores?\
(b) What is the $97.5^{th}$ percentile of the male heights? As always with normal probability problems, first draw a picture.[^13-foundations-mathematical-11]
:::

[^13-foundations-mathematical-11]: Remember: draw a picture first, then find the Z score.
    (We leave the pictures to you.) The Z score can be found by using the percentiles and the normal probability table.
    (a) We look for 0.95 in the probability portion (middle part) of the normal probability table, which leads us to row 1.6 and (about) column 0.05, i.e., $Z_{95}=1.65.$ Knowing $Z_{95}=1.65,$ $\mu = 1500,$ and $\sigma = 300,$ we setup the Z score formula: $1.65 = \frac{x_{95} - 1500}{300}.$ We solve for $x_{95}$: $x_{95} = 1995.$ (b) Similarly, we find $Z_{97.5} = 1.96,$ again setup the Z score formula for the heights, and calculate $x_{97.5} = 76.5.$

::: {.guidedpractice data-latex=""}
Using Z scores, answer the following questions.

(a) What is the probability that a randomly selected male adult is at least 6'2" (74 inches)?\
(b) What is the probability that a male adult is shorter than 5'9" (69 inches)?[^13-foundations-mathematical-12]
:::

[^13-foundations-mathematical-12]: Numerical answers: (a) 0.1131.
    (b) 0.3821.

::: {.workedexample data-latex=""}
What is the probability that a randomly selected adult male is between 5'9" and 6'2"?

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

These heights correspond to 69 inches and 74 inches.
First, draw the figure.
The area of interest is no longer an upper or lower tail.

```{r}
#| fig-asp: 0.4
#| out-width: 60%
#| fig-align: center
#| fig-width: 5
par(mar = c(2, 0, 0, 0))
normTail(70, 3.3, M = c(69, 74), col = IMSCOL["blue", "full"])
```

The total area under the curve is 1.
If we find the area of the two tails that are not shaded (from the previous Guided Practice, these areas are $0.3821$ and $0.1131$), then we can find the middle area:

```{r}
#| fig-asp: 0.25
#| fig-align: center
#| out-width: 100%
X <- seq(-3.2, 3.2, 0.01)
Y <- dnorm(X)
par(mar = c(2, 0, 0, 0))
plot(X, Y, type = "l", axes = F, xlim = c(-3.4, 24 + 3.4), ylim = c(0, 0.6), xlab = NA, ylab = NA)
lines(X, rep(0, length(X)))
these <- which(X <= 8)
polygon(c(X[these[1]], X[these], X[rev(these)[1]]), c(0, Y[these], 0), col = IMSCOL["blue", "full"])
lines(X, Y)
lines(c(3, 8 - 3), c(0.2, 0.2), lwd = 3)
lines(X + 8, Y, type = "l")
lines(X + 8, rep(0, length(X)))
these <- which(X < -0.303)
polygon(c(X[these[1]], X[these], X[rev(these)[1]]) + 8, c(0, Y[these], 0), col = IMSCOL["blue", "full"])
lines(X + 8, Y)
lines(8 + c(3, 8 - 3), c(0.2, 0.2), lwd = 3)
lines(X + 16, Y, type = "l")
lines(X + 16, rep(0, length(X)))
these <- which(X > 1.212)
polygon(c(X[these[1]], X[these], X[rev(these)[1]]) + 16, c(0, Y[these], 0), col = IMSCOL["blue", "full"])
lines(X + 16, Y)
lines(16 + c(3, 8 - 3), c(0.23, 0.23), lwd = 3)
lines(16 + c(3, 8 - 3), c(0.17, 0.17), lwd = 3)
lines(X + 24, Y, type = "l", xlim = c(-3.4, 3.4))
lines(X + 24, rep(0, length(X)))
these <- which(X > -0.303 & X < 1.212)
polygon(c(X[these[1]], X[these], X[rev(these)[1]]) + 24, c(0, Y[these], 0), col = IMSCOL["blue", "full"])
lines(X + 24, Y)
text(0, 0.53, "1.0000")
text(8, 0.53, "0.3821")
text(16, 0.53, "0.1131")
text(24, 0.53, "0.5048")
```

That is, the probability of being between 5'9" and 6'2" is 0.5048.
:::

::: {.guidedpractice data-latex=""}
Find the percent of SAT takers who earn between 1500 and 2000.[^13-foundations-mathematical-13]
:::

[^13-foundations-mathematical-13]: This is an abbreviated solution.
    (Be sure to draw a figure!) First find the percent who get below 1500 and the percent that get above 2000: $Z_{1500} = 0.00 \to 0.5000$ (area below), $Z_{2000} = 1.67 \to 0.0475$ (area above).
    Final answer: $1.0000-0.5000 - 0.0475 = 0.4525.$

::: {.guidedpractice data-latex=""}
What percent of adult males are between 5'5" and 5'7"?[^13-foundations-mathematical-14]
:::

[^13-foundations-mathematical-14]: 5'5" is 65 inches.
    5'7" is 67 inches.
    Numerical solution: $1.000 - 0.0649 - 0.8183 = 0.1168,$ i.e., 11.68%.

## Quantifying the variability of a statistic

As seen in later chapters, it turns out that many of the statistics used to summarize data (e.g., the sample proportion, the sample mean, differences in two sample proportions, differences in two sample means, the sample slope from a linear model, etc.) vary according to the normal distribution seen above.
The mathematical models are derived from the normal theory, but even the computational methods (and the intuitive thinking behind both approaches) use the general bell-shaped variability seen in most of the distributions constructed so far.

### 68-95-99.7 rule

Here, we present a useful general rule for the probability of falling within 1, 2, and 3 standard deviations of the mean in the normal distribution.
The rule will be useful in a wide range of practical settings, especially when trying to make a quick estimate without a calculator or Z table.

```{r}
#| label: fig-er6895997
#| fig-cap: |
#|   Probabilities for falling within 1, 2, and 3 standard deviations of the mean
#|   in a normal distribution.
#| fig-alt: |
#|   A normal curve showing the area within one standard
#|   deviation of the mean (which is 0.68), the area within two
#|   standard deviations of the mean (which is 0.95), and the are
#|   within three standard deviations of the mean (which is 0.997).
#| fig-asp: 0.3
#| out-width: 60%
par(mar = c(2, 0, 0, 0))
X <- seq(-4, 4, 0.01)
Y <- dnorm(X)
plot(X, Y, type = "n", axes = F, ylim = c(0, 0.4), xlim = c(-3.2, 3.2), xlab = NA, ylab = NA)
abline(h = 0, col = IMSCOL["black", "full"])
axis(1, at = -3:3, label = expression(
  mu - 3 * sigma, mu - 2 * sigma, mu - sigma, mu,
  mu + sigma, mu + 2 * sigma, mu + 3 * sigma
))
ii <- c(1, 2, 3)
jj <- c(1, 1, 1)
for (i in 3:1) {
  these <- (X >= i - 1 & X <= i)
  polygon(c(i - 1, X[these], i), c(0, Y[these], 0), col = IMSCOL[ii[i], jj[i]], border = IMSCOL[ii[i], jj[i]])
  these <- (X >= -i & X <= -i + 1)
  polygon(c(-i, X[these], -i + 1), c(0, Y[these], 0), col = IMSCOL[ii[i], jj[i]], border = IMSCOL[ii[i], jj[i]])
}
# ===> label 99.7 <===#
arrows(-3, 0.03, 3, 0.03, code = 3, col = "#444444", length = 0.15)
text(0, 0.02, "99.7%", pos = 3)
# ===> label 95 <===#
arrows(-2, 0.13, 2, 0.13, code = 3, col = "#444444", length = 0.15)
text(0, 0.12, "95%", pos = 3)
# ===> label 68 <===#
arrows(-1, 0.23, 1, 0.23, code = 3, col = "#444444", length = 0.15)
text(0, 0.22, "68%", pos = 3)
lines(X, Y, col = "#888888")
abline(h = 0, col = "#AAAAAA")
```

::: {.guidedpractice data-latex=""}
Use `pnorm()` (or a Z table) to confirm that about 68%, 95%, and 99.7% of observations fall within 1, 2, and 3, standard deviations of the mean in the normal distribution, respectively.
For instance, first find the area that falls between $Z=-1$ and $Z=1,$ which should have an area of about 0.68.
Similarly there should be an area of about 0.95 between $Z=-2$ and $Z=2.$[^13-foundations-mathematical-15]
:::

[^13-foundations-mathematical-15]: First draw the pictures.
    To find the area between $Z=-1$ and $Z=1,$ use `pnorm()` or the normal probability table to determine the areas below $Z=-1$ and above $Z=1.$ Next verify the area between $Z=-1$ and $Z=1$ is about 0.68.
    Repeat this for $Z=-2$ to $Z=2$ and for $Z=-3$ to $Z=3.$

It is possible for a normal random variable to fall 4, 5, or even more standard deviations from the mean.
However, these occurrences are very rare if the data are nearly normal.
The probability of being further than 4 standard deviations from the mean is about 1-in-30,000.
For 5 and 6 standard deviations, it is about 1-in-3.5 million and 1-in-1 billion, respectively.

::: {.guidedpractice data-latex=""}
SAT scores closely follow the normal model with mean $\mu = 1500$ and standard deviation $\sigma = 300.$ About what percent of test takers score 900 to 2100?
What percent score between 1500 and 2100 ?[^13-foundations-mathematical-16]
:::

[^13-foundations-mathematical-16]: 900 and 2100 represent two standard deviations above and below the mean, which means about 95% of test takers will score between 900 and 2100.
    Since the normal model is symmetric, then half of the test takers from part (a) ($\frac{95\%}{2} = 47.5\%$ of all test takers) will score 900 to 1500 while 47.5% score between 1500 and 2100.

### Standard error {#sec-se}

Point estimates vary from sample to sample, and we quantify this variability with what is called the **standard error (SE)**\index{standard error}.
The standard error is equal to the standard deviation associated with the statistic.
So, for example, to quantify the variability of a point estimate from one sample to the next, the variability is called the standard error of the point estimate.
Almost always, the standard error is itself an estimate, calculated from the sample of data.

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

The way we determine the standard error varies from one situation to the next.
However, typically it is determined using a formula based on the Central Limit Theorem.

### Margin of error {#sec-moe}

Very related to the standard error is the **margin of error**\index{margin of error}.
The margin of error describes how far away observations are from their mean.\
For example, to describe where most (i.e., 95%) observations lie, we say that the margin of error is approximately $2 \times SE$.
That is, 95% of the observations are within two margins of error of the mean.

::: {.important data-latex=""}
**Margin of error for sample proportions.**

The distance given by $z^\star \times SE$ is called the **margin of error**\index{margin of error}.

$z^\star$ is the cutoff value found on the normal distribution.
The most common value of $z^\star$ is 1.96 (often approximated to be 2) indicating that the margin of error describes the variability associated with 95% of the sampled statistics.
:::

Notice that if the spread of the observations goes from some lower bound to some upper bound, a rough approximation of the SE is to divide the range by 4.
That is, if you notice the sample proportions go from 0.1 to 0.4, the SE can be approximated to be 0.075.

```{r}
#| include: false
terms_chp_13 <- c(terms_chp_13, "margin of error")
```

## Case Study (test): Opportunity cost {#sec-caseopp}

The approach for using the normal model in the context of inference is very similar to the practice of applying the model to individual observations that are nearly normal.
We will replace null distributions we previously obtained using the randomization or simulation techniques and verify the results once again using the normal model.
When the sample size is sufficiently large, the normal approximation generally provides us with the same conclusions as the simulation model.

### Observed data

In @sec-caseStudyOpportunityCost we were introduced to the opportunity cost study, which found that students became thriftier when they were reminded that not spending money now means the money can be spent on other things in the future.
Let's re-analyze the data in the context of the normal distribution and compare the results.

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

### Variability of the statistic

@fig-OpportunityCostDiffs-w-normal summarizes the null distribution as determined using the randomization method.
The best fitting normal distribution for the null distribution has a mean of 0.
We can calculate the standard error of this distribution by borrowing a formula that we will become familiar with in [Chapter -@sec-inference-two-props], but for now let's just take the value $SE = 0.078$ as a given.
Recall that the point estimate of the difference was 0.20, as shown in @fig-OpportunityCostDiffs-w-normal.
Next, we'll use the normal distribution approach to compute the p-value.

```{r}
#| label: fig-OpportunityCostDiffs-w-normal
#| fig-cap: |
#|   Null distribution of differences with an overlaid normal curve for the
#|   opportunity cost study. 10,000 simulations were run for this figure.
#| fig-alt: |
#|   Normal curve overlaid on the histogram of the simulated
#|   differences in sample proportions in the opportunity cost study. The
#|   theoretical normal distribution fits the histogram
#|   quite well. The observed difference in sample proportions of
#|   0.2 is in the tail of both the normal curve and the histogram.
#| fig-asp: 0.45
ggplot(opportunity_cost_sim_dist, aes(x = stat)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.05) +
  geom_function(fun = dnorm, args = list(mean = 0, sd = sd(opportunity_cost_sim_dist$stat)), color = IMSCOL["black", "full"], size = 1) +
  scale_x_continuous(breaks = seq(-0.3, 0.3, 0.1), labels = label_number(accuracy = 0.1)) +
  coord_cartesian(xlim = c(-0.35, 0.35)) +
  theme(
    axis.text.y = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank()
  ) +
  labs(
    x = "When the null hypothesis is true,\ndifference in randomized proportions of students who do not buy the video",
    y = NULL
  ) +
  annotate("segment", x = 0.2, y = 0, xend = 0.2, yend = 2, linetype = "dashed") +
  annotate("text", x = 0.2, y = 2.8, label = "observed\ndifference")
```

### Observed statistic vs. null statistics

As we learned in @sec-normalDist, it is helpful to draw and shade a picture of the normal distribution so we know precisely what we want to calculate.
Here we want to find the area of the tail beyond 0.2, representing the p-value.

```{r}
#| label: OpportunityCostDiffs_normal_only
#| fig-asp: 0.3
#| out-width: 60%
#| fig-align: center
#| fig-width: 5
par(mar = c(2, 0, 0, 0))
normTail(m = 0, s = 0.078, U = 0.2, col = IMSCOL["blue", "full"])
```

Next, we can calculate the Z score using the observed difference, 0.20, and the two model parameters.
The standard error, $SE = 0.078,$ is the equivalent of the model's standard deviation.

$$Z = \frac{\text{observed difference} - 0}{SE} = \frac{0.20 - 0}{0.078} = 2.56$$

We can either use statistical software or look up $Z = 2.56$ in the normal probability table to determine the right tail area: 0.0052, which is about the same as what we got for the right tail using the randomization approach (0.006).
Using this area as the p-value, we see that the p-value is less than 0.05, we conclude that the treatment did indeed impact students' spending.

::: {.important data-latex=""}
**Z score in a hypothesis test.**

In the context of a hypothesis test, the Z score for a point estimate is

$$Z = \frac{\text{point estimate} - \text{null value}}{SE}$$

The standard error in this case is the equivalent of the standard deviation of the point estimate, and the null value comes from the claim made in the null hypothesis.
:::

We have confirmed that the randomization approach we used earlier and the normal distribution approach provide almost identical p-values and conclusions in the opportunity cost case study.
Next, let's turn our attention to the medical consultant case study.

## Case study (test): Medical consultant {#sec-casemed}

### Observed data

In @sec-case-study-med-consult we learned about a medical consultant who reported that only 3 of their 62 clients who underwent a liver transplant had complications, which is less than the more common complication rate of 0.10.
In that work, we did not model a null scenario, but we will discuss a simulation method for a one proportion null distribution in @sec-one-prop-null-boot, such a distribution is provided in @fig-MedConsNullSim-w-normal.
We have added the best-fitting normal curve to the figure, which has a mean of 0.10.
Borrowing a formula that we'll encounter in [Chapter -@sec-inference-one-prop], the standard error of this distribution was also computed: $SE = 0.038.$

### Variability of the statistic

Before we begin, we want to point out a simple detail that is easy to overlook: the null distribution we generated from the simulation is slightly skewed, and the histogram is not particularly smooth.
In fact, the normal distribution only sort-of fits this model.

```{r}
#| label: fig-MedConsNullSim-w-normal
#| fig-cap: |
#|   The null distribution for the sample proportion, created from 10,000
#|   simulated studies from the medical consultant, along with the best-fitting
#|   normal model.
#| fig-alt: |
#|   Normal curve overlaid on the histogram of the simulated
#|   sample proportions in the medical consultant study. The
#|   theoretical normal distribution does not fit the histogram
#|   as well as in @fig-OpportunityCostDiffs-w-normal due to the
#|   skewness of the distribution of the sample proportions.
#|   The observed sample proportion of 0.0484 is towards the tail
#|   in both the normal curve and the histogram, but it is not
#|   an extreme value.
#| fig-asp: 0.5   
ggplot(medical_consultant_sim_dist, aes(x = stat)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.0167, fill = IMSCOL["blue", "full"]) +
  geom_function(fun = dnorm, args = list(mean = mean(medical_consultant_sim_dist$stat), sd = sd(medical_consultant_sim_dist$stat)), color = IMSCOL["black", "full"], size = 1) +
  scale_x_continuous(breaks = seq(0, 0.3, 0.05), labels = label_number(accuracy = 0.01)) +
  theme(
    axis.text.y = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank()
  ) +
  labs(
    x = expression(hat(p)[sim]),
    y = NULL
  ) +
  annotate("segment", x = 0.05, y = 0, xend = 0.05, yend = 8, linetype = "dashed") +
  annotate("text", x = 0.05, y = 9, label = "observed\nvalue")
```

### Observed statistic vs. null statistics

As always, we'll draw a picture before finding the normal probabilities.
Below is a normal distribution centered at 0.10 with a standard error of 0.038.

```{r}
#| label: MedConsNullSim-normal-only
#| fig-asp: 0.3
#| out-width: 60%
#| fig-align: center
#| fig-width: 5
par(mar = c(2, 0, 0, 0))
normTail(m = 0.1, s = 0.038, L = 0.048, col = IMSCOL["blue", "full"])
```

Next, we can calculate the Z score using the observed complication rate, $\hat{p} = 0.048$ along with the mean and standard deviation of the normal model.
Here again, we use the standard error for the standard deviation.

$$Z = \frac{\hat{p} - p_0}{SE_{\hat{p}}} = \frac{0.048 - 0.10}{0.038} = -1.37$$

Identifying $Z = -1.37$ using statistical software or in the normal probability table, we can determine that the left tail area is 0.0853 which is the estimated p-value for the hypothesis test.
There is a small problem: the p-value of 0.0853 is almost 30% smaller than the simulation p-value of 0.1222 which will be calculated in @sec-one-prop-null-boot.

The discrepancy is explained by the normal model's poor representation of the null distribution in @fig-MedConsNullSim-w-normal.
As noted earlier, the null distribution from the simulations is not very smooth, and the distribution itself is slightly skewed.
That's the bad news.
The good news is that we can foresee these problems using some simple checks.
We'll learn more about these checks in the following chapters.

In @sec-CLTsection we noted that the two common requirements to apply the Central Limit Theorem are (1) the observations in the sample must be independent, and (2) the sample must be sufficiently large.
The guidelines for this particular situation -- which we will learn in [Chapter -@sec-inference-one-prop] -- would have alerted us that the normal model was a poor approximation.

### Conditions for applying the normal model

The success story in this section was the application of the normal model in the context of the opportunity cost data.
However, the biggest lesson comes from the less successful attempt to use the normal approximation in the medical consultant case study.

Statistical techniques are like a carpenter's tools.
When used responsibly, they can produce amazing and precise results.
However, if the tools are applied irresponsibly or under inappropriate conditions, they will produce unreliable results.
For this reason, with every statistical method that we introduce in future chapters, we will carefully outline conditions when the method can reasonably be used.
These conditions should be checked in each application of the technique.

After covering the introductory topics in this course, advanced study may lead to working with complex models which, for example, bring together many variables with different variability structure.
Working with data that come from normal populations makes higher-order models easier to estimate and interpret.
There are times when simulation, randomization, or bootstrapping are unwieldy in either structure or computational demand.
Normality can often lead to excellent approximations of the data using straightforward modeling techniques.

## Case study (interval): Stents {#sec-casestent}

A point estimate is our best guess for the value of the parameter, so it makes sense to build the confidence interval around that value.
The standard error, which is a measure of the uncertainty associated with the point estimate, provides a guide for how large we should make the confidence interval.
The 68-95-99.7 rule tells us that, in general, 95% of observations are within 2 standard errors of the mean.
Here, we use the value 1.96 to be slightly more precise.

::: {.important data-latex=""}
**Constructing a 95% confidence interval.**

When the sampling distribution of a point estimate can reasonably be modeled as normal, the point estimate we observe will be within 1.96 standard errors of the true value of interest about 95% of the time.
Thus, a **95% confidence interval**\index{95\% confidence interval}\index{confidence interval} for such a point estimate can be constructed:

\vspace{-5mm}

$$\text{point estimate} \pm 1.96 \times SE$$

We can be **95% confident**\index{95\% confident} this interval captures the true value.
:::

\vspace{-3mm}

```{r}
#| include: false
terms_chp_13 <- c(terms_chp_13, "95% confidence interval", "95% confident")
```

::: {.guidedpractice data-latex=""}
Compute the area between -1.96 and 1.96 for a normal distribution with mean 0 and standard deviation 1.[^13-foundations-mathematical-17]
:::

[^13-foundations-mathematical-17]: We will leave it to you to draw a picture.
    The Z scores are $Z_{left} = -1.96$ and $Z_{right} = 1.96.$ The area between these two Z scores is $0.9750 - 0.0250 = 0.9500.$ This is where "1.96" comes from in the 95% confidence interval formula.

::: {.workedexample data-latex=""}
The point estimate in the opportunity cost study was that 20% fewer students would buy a video if they were reminded that money not spent now could be spent later on something else.
This point estimate can reasonably be modeled with a normal distribution with a standard error of $SE = 0.078.$ Construct a 95% confidence interval for the point estimate.

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

Since we're told the point estimate can be modeled with a normal distribution:

\vspace{-5mm}

$$\text{point estimate} \pm 1.96 \times SE = 0.20 \pm 1.96 \times 0.078 = (0.047, 0.353)$$

We are 95% confident that the video purchase rate resulting from the treatment is between 4.7% and 35.3% lower than in the control group.
Since this confidence interval does not contain 0, it is consistent with our earlier hypothesis test where we rejected the notion of "no difference".

Note that we have used SE = 0.078 from the last section.
However, it would more generally be appropriate to recompute the SE slightly differently for this confidence interval using sample proportions.
Don't worry about this detail for now since the two resulting standard errors are, in this case, almost identical.
:::

### Observed data

Consider an experiment that examined whether implanting a stent in the brain of a patient at risk for a stroke helps reduce the risk of a stroke.
The results from the first 30 days of this study, which included 451 patients, are summarized in @tbl-stentStudyResultsCIsection.
These results are surprising!
The point estimate suggests that patients who received stents may have a *higher* risk of stroke: $p_{trmt} - p_{ctrl} = 0.090.$

```{r}
#| label: tbl-stentStudyResultsCIsection
#| tbl-cap: Descriptive statistics for 30-day results for the stent study.
#| tbl-pos: H
stent30 |>
  count(group, outcome) |>
  pivot_wider(names_from = outcome, values_from = n, values_fill = 0) |>
  janitor::adorn_totals(where = c("row", "col")) |>
  kbl(
    linesep = "", booktabs = TRUE,
    col.names = c("Group", "No event", "Stroke", "Total")
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped"), full_width = FALSE
  ) |>
  column_spec(1:4, width = "7em")
```

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

### Variability of the statistic

::: {.workedexample data-latex=""}
Consider the stent study and results.
The conditions necessary to ensure the point estimate $p_{trmt} - p_{ctrl} = 0.090$ is nearly normal have been verified for you, and the estimate's standard error is $SE = 0.028.$ Construct a 95% confidence interval for the change in 30-day stroke rates from usage of the stent.

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

The conditions for applying the normal model have already been verified, so we can proceed to the construction of the confidence interval:

$$\text{point estimate} \pm 1.96 \times SE = 0.090 \pm 1.96 \times 0.028 = (0.035, 0.145)$$

We are 95% confident that implanting a stent in a stroke patient's brain increased the risk of stroke within 30 days by a rate of 0.035 to 0.145.
This confidence interval can also be used in a way analogous to a hypothesis test: since the interval does not contain 0 (is completely above 0), it means the data provide convincing evidence that the stent used in the study changed the risk of stroke within 30 days.
:::

As with hypothesis tests, confidence intervals are imperfect.
About 1-in-20 properly constructed 95% confidence intervals will fail to capture the parameter of interest, simply due to natural variability in the observed data.
@fig-95PercentConfidenceInterval shows 25 confidence intervals for a proportion that were constructed from 25 different datasets that all came from the same population where the true proportion was $p = 0.3.$ However, 1 of these 25 confidence intervals happened not to include the true value.
The interval which does not capture $p=0.3$ is not due to bad science.
Instead, it is due to natural variability, and we should *expect* some of our intervals to miss the parameter of interest.
Indeed, over a lifetime of creating 95% intervals, you should expect 5% of your reported intervals to miss the parameter of interest (unfortunately, you will not ever know which of your reported intervals captured the parameter and which missed the parameter).

```{r}
#| label: fig-95PercentConfidenceInterval
#| fig-cap: |
#|   Twenty-five samples of size $n=300$ were collected from a population with
#|   $p = 0.30.$ For each sample, a confidence interval was created to try to capture
#|   the true proportion $p.$ However, 1 of these 25 intervals did not capture $p = 0.30.$
#| fig-asp: 0.5
#| fig-alt: |
#|   A series of 25 horizontal lines are drawn, representing each of 25 different
#|   samples. Each vertical line starts at the value of the lower bound of the confidence
#|   interval and ends at the value of the upper bound of the confidence interval which
#|   was created from that particular sample. In the center of the line is a solid dot
#|   at the observed proportion of successes for that particular sample. A dashed vertical
#|   line runs through the horizontal lines at p = 0.3 (which is the true value of the
#|   population proportion). 24 of the 25 horizontal lines cross the vertical line at
#|   0.3, but one of the horizontal lines is completely lower than 0.3. The line that
#|   does not cross 0.3 is colored red because the confidence interval from that particular
#|   sample would not have captured the true population proportion.

set.seed(4747)
num.sim <- 25
n <- 300
p <- 0.3
p.hat <- rbinom(num.sim, n, p) / n
SE <- sqrt(p.hat * (1 - p.hat) / n)
par(mar = c(2, 1, 1, 1), mgp = c(2.7, 0.9, 0))
xR <- p + 4 * c(-1, 1) * max(SE)
yR <- c(0, 25)
plot(xR, yR,
  type = "n", xlab = "",
  ylab = "", axes = FALSE
)
abline(v = p, lty = 2, col = IMSCOL["gray", "full"])
axis(1,
  at = p + c(-0.05, 0, 0.05),
  c("0.25", "p = 0.30", "0.35"), cex.axis = 1
)
for (i in 1:num.sim) {
  ci <- p.hat[i] + 2 * c(-1, 1) * SE[i]
  if (abs(p.hat[i] - p) > 2 * SE[i]) {
    col <- IMSCOL["red", "full"]
    points(p.hat[i], i, cex = 1.4, col = col)
  } else {
    col <- IMSCOL["blue", "full"]
  }
  points(p.hat[i], i, pch = 20, cex = 1.2, col = col)
  lines(ci, rep(i, 2), col = col)
}
```

::: {.guidedpractice data-latex=""}
In @fig-95PercentConfidenceInterval, one interval does not contain the true proportion, $p = 0.3.$ Does this imply that there was a problem with the datasets that were selected?[^13-foundations-mathematical-18]
:::

[^13-foundations-mathematical-18]: No.
    Just as some observations occur more than 1.96 standard deviations from the mean, some point estimates will be more than 1.96 standard errors from the parameter.
    A confidence interval only provides a plausible range of values for a parameter.
    While we might say other values are implausible based on the data, this does not mean they are impossible.

### Interpreting confidence intervals

A careful eye might have observed the somewhat awkward language used to describe confidence intervals.

::: {.important data-latex=""}
**Correct confidence interval interpretation.**

We are XX% confident that the population parameter is between *lower* and *upper* (where *lower* and *upper* are both numerical values).

**Incorrect** language might try to describe the confidence interval as capturing the population parameter with a certain probability.

This is one of the most common errors: while it might be useful to think of it as a probability, the confidence level only quantifies how plausible it is that the parameter is in the interval.
:::

Another especially important consideration of confidence intervals is that they *only try to capture the population parameter*.
Our intervals say nothing about the confidence of capturing individual observations, a proportion of the observations, or about capturing point estimates.
Confidence intervals provide an interval estimate for and attempt to capture **population parameters**.

\index{confidence interval}

\clearpage

## Chapter review {#sec-chp13-review}

### Summary

We can summarise the process of using the normal model as follows:

-   **Frame the research question.** The mathematical model can be applied to both the hypothesis testing and the confidence interval framework. Make sure that your research question is being addressed by the most appropriate inference procedure.
-   **Collect data with an observational study or experiment.** To address the research question, collect data on the variables of interest. Note that your data may be a random sample from a population or may be part of a randomized experiment.
-   **Model the randomness of the statistic.** In many cases, the normal distribution will be an excellent model for the randomness associated with the statistic of interest. The Central Limit Theorem tells us that if the sample size is large enough, sample averages (which can be calculated as either a proportion or a sample mean) will be approximately normally distributed when describing how the statistics change from sample to sample.
-   **Calculate the variability of the statistic.** Using formulas, come up with the standard deviation (or more typically, an estimate of the standard deviation called the standard error) of the statistic. The SE of the statistic will give information on how far the observed statistic is from the null hypothesized value (if performing a hypothesis test) or from the unknown population parameter (if creating a confidence interval).
-   **Use the normal distribution to quantify the variability.** The normal distribution will provide a probability which measures how likely it is for your observed and hypothesized (or observed and unknown) parameter to differ by the amount measured. The unusualness (or not) of the discrepancy will form the conclusion to the research question.
-   **Form a conclusion.** Using the p-value or the confidence interval from the analysis, report on the research question of interest. Also, be sure to write the conclusion in plain language so casual readers can understand the results.

@tbl-chp13-summary is another look at the mathematical model approach to inference.

```{r}
#| label: tbl-chp13-summary
#| tbl-cap: Summary of mathematical models as an inferential statistical method.
#| tbl-pos: H
inference_method_summary_table |>
  filter(question != "What are the technical conditions?") |>
  select(question, mathematical) |>
  kbl(
    linesep = "\\addlinespace", booktabs = TRUE,
    col.names = c("Question", "Answer")
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "condensed"),
    latex_options = c("striped"), full_width = TRUE
  ) |>
  column_spec(1, width = "15em")
```

\clearpage

### Terms

The terms introduced in this chapter are presented in @tbl-terms-chp-13.
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-13
#| tbl-cap: Terms introduced in this chapter.
#| tbl-pos: H
make_terms_table(terms_chp_13)
```

\clearpage

## Exercises {#sec-chp13-exercises}

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

::: {.exercises data-latex=""}
{{< include exercises/_13-ex-foundations-mathematical.qmd >}}
:::