---
title: "Sampling variation & the Normal model (Notes)"
subtitle: "STAT 155"
author: "Your Name"
format:
  html:
    toc: true
    toc-depth: 2
    embed-resources: true
---


```{r setup}
#| include: false
knitr::opts_chunk$set(
  collapse = TRUE, 
  warning = FALSE,
  message = FALSE,
  error = TRUE,
  fig.height = 2.75, 
  fig.width = 4.25,
  fig.env = 'figure',
  fig.pos = 'h',
  fig.align = 'center')
```


::: {.callout-note title = "Learning goals"}

- Recognize the difference between a **population parameter** and a **sample estimate**.
- Review the Normal probability model, a tool we'll need to turn information in our **sample data** into **inferences** about the broader **population**.
- Explore the *ideas* of randomness, sampling distributions, and standard error through a class experiment. (We'll define these more formally in the next class.)

![](https://ajohns24.github.io/155_Fall_2025/images/concept_map_8.png)

:::

::: {.callout-note title="Additional resources"}

**Required videos**

- [Exploration vs inference](https://youtu.be/_KfOoWQygPk)
- [Normal probability model](https://youtu.be/bGOqwcOMdpo)

**Optional**

- Read: Section 6 Introduction, and Section 6.6 in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/)

:::




\
\
\
\



# Warm-up {-}

Today marks a significant shift in the semester!
Thus far, we’ve focused on **exploratory analysis**: visualizations, numerical summaries, and models that capture the *observed* patterns in our sample data.
Today, we’ll move on to **inference**: using our data to make conclusions about the broader population from which it was taken.

\
\


**Context**

The `elections` dataset was compiled by data posted by [Tony McGovern](https://github.com/tonmcg/US_County_Level_Election_Results_08-24).
This dataset includes the 2012, 2016, 2020, and 2024 presidential election results for the population of all 3000+ U.S. counties (except Alaska). 
(Image: [Wikimedia Commons](https://commons.wikimedia.org/wiki/File:Map_of_USA_with_county_outlines_(black_%26_white).png))

![](https://upload.wikimedia.org/wikipedia/commons/6/68/Map_of_USA_with_county_outlines_%28black_%26_white%29.png)


```{r}
library(tidyverse)
elections <- read_csv("https://ajohns24.github.io/data/election_2024.csv") %>% 
  mutate(repub_24 = round(100*per_gop_2024, 2),
         repub_20 = round(100*per_gop_2020, 2),
         repub_16 = round(100*per_gop_2016, 2), 
         repub_12 = round(100*per_gop_2012, 2)) %>% 
  select(state_name, state_abbr, county_name, county_fips, repub_24, repub_20, repub_16, repub_12)

head(elections)
```


\
\
\
\


**POPULATION PARAMETERS**

Let Y be the Republican candidate's (Trump's) 2024 county-level vote percentage, and X be the Republican candidate's (Mitt Romney's) 2012 vote percentage.
Then our goal will be to model the relationship of Y with X among all U.S. counties (except Alaska), our **population of interest**:

$$
E[Y|X] = \beta_0 + \beta_1 X
$$

We refer to the $\beta_0$ and $\beta_1$ coefficients as **population parameters**.
These are the "true", underlying, fixed parameters that describe the "true" linear relationship between Y and X in the population of U.S. counties.
NOTE: We're putting "true" in quotes here to reflect the fact that some statisticians would view this relationship as constantly in flux, instead of "true" and "fixed", with the elections providing 1 snapshot.

In this example, since we have data on all U.S. counties in the population of interest, we can actually obtain the "true" population parameters:

$$
E[Y|X] = \beta_0 + \beta_1 X = 9.699 + 0.957 X
$$

```{r}
population_model <- lm(repub_24 ~ repub_12, data = elections)
summary(population_model)
```


```{r}
# The blue line is the population model
# The red dashed line shows a 1-to-1 relationship, for comparison
elections %>% 
  ggplot(aes(y = repub_24, x = repub_12)) + 
  geom_point() + 
  geom_smooth(method = "lm") +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed")
```




\
\
\
\


**SAMPLE ESTIMATES**

NOW FORGET THAT YOU KNOW ALL OF THE ABOVE.
Let's pretend that we're working within the *typical* scenario - we don't have access to the entire population of interest, thus we don't actually know the "true" population parameters $\beta_0$ and $\beta_1$.
Instead, suppose we only have data on a *sample* of counties.
We can then use the sample data to *estimate* the population parameters, denoted by putting a "hat" symbol, $\hat{ }$, on top of the betas:

$$
E[Y|X] = \hat{\beta}_0 + \hat{\beta}_1 X
$$





\
\
\
\



**EXAMPLE 1: Intuition**

To explore the nuances of using *sample* data to learn about a *population*, you will each take a random sample of 10 U.S. counties.
You will use these counties to *estimate* the model of Y vs X.

a. Do you think that you will get the *same sample* of 10 counties as your neighbors?



b. If not, do you think that you will get the *same coefficient estimates* as your neighbors?



c. Below is a plot of the actual population model. How do you think the students' sample model will compare? (Let's draw.)

```{r}
elections %>% 
  ggplot(aes(y = repub_24, x = repub_12)) + 
  geom_point() + 
  geom_smooth(method = "lm")
```




\
\
\
\


**EXAMPLE 2: Random sampling**

Use the `sample_n()` function to take a *random sample* of 2 U.S. counties:

```{r}
# Try running the following chunk A FEW TIMES
elections %>% 
  sample_n(size = 2, replace = FALSE)
```

Reflect: 
 
- How do your results compare to your neighbors' results? How, if at all, do they change when you re-run your chunk?

- What is the role of `size = 2`? *HINT*: Remember you can look at function documentation by running `?sample_n` in the console!

- What is the role of `replace = FALSE`? *HINT*: Remember you can look at function documentation by running `?sample_n` in the console!


\
\
\
\


**EXAMPLE 3: Setting the seed**

Now, "set the seed" to 155 and re-try your sampling:

```{r}
# Try running the following FULL chunk A FEW TIMES
set.seed(155)
elections %>% 
  sample_n(size = 2, replace = FALSE)
```

What changed?



\
\
\
\



**EXAMPLE 4: Take your own sample of 50 counties**

The underlying *random number generator* plays a role in the random sample we happen to get.
If we `set.seed(some positive integer)` before taking a random sample, we'll get the same results.
This **reproducibility** is important:    
    
- we get the same results every time we render our qmd
- we can share our work with others & ensure they get our same answers
- it wouldn't be great if you submitted your work to, say, a journal and weren't able to back up / confirm / reproduce your results!    

Follow the chunks below to obtain and use your *own* unique sample.

```{r}
# DON'T SKIP THIS STEP! 
# Set the random number seed to the digits of your own phone number (just the numbers)
set.seed(___)

# Take a sample of 10 counties
my_sample <- elections %>% 
  sample_n(size = 50, replace = FALSE)
my_sample                       
```    

```{r}
# Plot the SAMPLE model of Y vs X among YOUR sample data
# Compare it to the actual POPULATION model (red) 
my_sample %>% 
  ggplot(aes(y = repub_24, x = repub_12)) +
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE) + 
  geom_smooth(data = elections, method = "lm", se = FALSE, color = "red", size = 0.5)
```

```{r}
# Model the relationship among your sample
my_model <- lm(repub_24 ~ repub_12, data = my_sample)
coef(summary(my_model))
```

**REPORT YOUR WORK**

Log your intercept and slope *sample estimates* in [this survey](https://forms.gle/WAXn7NNuw6daf2TM8).










\
\
\
\


# Exercises {-}

## Exercise 1: Using the Normal model {-}

Before getting back to the experiment, let's practice the concepts you explored in today’s video.
First, run the following code which defines a `shaded_normal()` function which is specialized to this activity alone:

```{r}
library(tidyverse)
shaded_normal <- function(mean, sd, a = NULL, b = NULL){
  min_x <- mean - 4*sd
  max_x <- mean + 4*sd
  a <- max(a, min_x)
  b <- min(b, max_x)
  ggplot() + 
    scale_x_continuous(limits = c(min_x, max_x), breaks = c(mean - sd*(0:3), mean + sd*(1:3))) +
    stat_function(fun = dnorm, args = list(mean = mean, sd = sd)) + 
    stat_function(geom = "area", fun = dnorm, args = list(mean = mean, sd = sd), xlim = c(a, b), fill = "blue") + 
    labs(y = "density") + 
    theme_minimal()
}
```


Next, suppose that the speeds of cars on a highway, in miles per hour, can be reasonably represented by the Normal model with a *mean* of 55mph and a *standard deviation* of 5mph from car to car:

$$
X \sim N(55, 5^2)
$$

```{r}
shaded_normal(mean = 55, sd = 5)
```

### Part a {-}

Provide the (approximate) range of the middle 68% of speeds, and shade in the corresponding region on your Normal curve.
NOTE: `a` is the lower end of the range and `b` is the upper end.    

```{r}
shaded_normal(mean = 55, sd = 5, a = ___, b = ___)
```
    

### Part b {-}

The area of the shaded region below reflects the probability that a car's speed exceeds 60mph:

```{r}
shaded_normal(mean = 55, sd = 5, a = 60)
```    

    
Use the 68-95-99.7 rule to *estimate* this probability.   


### Part c {-}

The area of the shaded region below reflects the probability that a car's speed exceeds 67mph:

```{r}
shaded_normal(mean = 55, sd = 5, a = 67)
```

Which of the following is the correct *range* for this probability?
Explain your reasoning.

* less than 0.0015
* between 0.0015 and 0.025
* between 0.025 and 0.16
* greater than 0.16



\
\
\
\






## Exercise 2: Z-scores {-}

Inherently important to all of our calculations above is *how many standard deviations a value "X" is from the mean*.
This distance is called a **Z-score** and can be calculated as follows:    
    
$$
\text{Z-score} = \frac{X - \text{mean}}{\text{sd}}
$$
  
For example, under the $N(55, 5^2)$ speed settings in Exercise 1, if I'm traveling 40 miles an hour, my Z-score is -3.
That is, my speed is 3 standard deviations *below* the average speed:    
    
```{r}
(40 - 55) / 5
```

### Part a {-}

Consider 2 other drivers.
Both drivers are speeding.
Who do you think is speeding *more*, relative to the distributions of speed in their area?
(Just give a quick answer without doing any formal calculations yet.)

- Driver A is traveling at 60mph on the highway where speeds are $N(55, 5^2)$ and the speed limit is 55mph.
- Driver B is traveling at 36mph on a residential road where speeds are $N(30, 3^2)$ and the speed limit is 30mph.
    
    
### Part b {-}

Calculate the Z-scores for Drivers A and B.

```{r}
# Driver A


# Driver B


```

### Part c {-}

Now, based on the Z-scores, who is speeding more?
NOTE: The below plots might provide some insights.   

```{r}
# Driver A
shaded_normal(mean = 55, sd = 5) + 
  geom_vline(xintercept = 60) + 
  lims(x = c(18, 75), y = c(0, 0.15)) + 
  labs(title = "Driver A")

# Driver B
shaded_normal(mean = 30, sd = 3) + 
  geom_vline(xintercept = 36) + 
  lims(x = c(18, 75), y = c(0, 0.15)) + 
  labs(title = "Driver B")

```
    


\
\
\
\





## Exercise 3: Sampling variation {-}

Let's return to the election experiment.
Recall that our **population model** of the Republican vote in 2024 (Y) vs 2012 (X) is:

$$
E[Y|X] = \beta_0 + \beta_1 X = 9.699 + 0.957 X
$$

To *estimate* this population model, each student took a *sample* of 10 counties.
Let's explore how our **sample estimates** of the population model *varied* from student to student:


```{r}
# Import the experiment results
library(gsheet)
results <- gsheet2tbl('https://docs.google.com/spreadsheets/d/1xyLff_FVqN419CyrQfa4CO8uRjYC4oP-S9Sla3FA_lE/edit?usp=sharing')
```

Plot each student's *sample estimate* of the model line (gray).
How do these compare to the assumed *population* model (red)?    

```{r}
elections %>% 
  ggplot(aes(y = repub_24, x = repub_12)) +
  geom_abline(data = results, aes(intercept = sample_intercept, slope = sample_slope), color = "gray") + 
  geom_smooth(method = "lm", color = "red", se = FALSE)
```
    




\
\
\
\





## Exercise 4: Sample intercepts {-}

Let's focus on just the *sample estimates* of the *intercept* parameter:

```{r}
results %>% 
  ggplot(aes(x = sample_intercept)) + 
  geom_density() + 
  geom_vline(xintercept = 9.699, color = "red")
```    

Comment on the shape, center, and spread of these sample estimates and how they relate to the "true" population intercept of 9.699 (red line).
What does this imply about the sample estimation procedure?



\
\
\
\











## Exercise 5: Slope estimates {-}

### Part a

In Exercise 4, we constructed a density plot of the sample estimates of the *intercept* for each student.
Suppose we were to construct a density plot of the sample estimates of the `repub_12` coefficient (i.e. the slopes).
*Intuitively*:

- What *shape* do you think this plot will have?

- Around what value do you think this plot will be *centered*?


### Part b

Check your intuition:        

```{r}
results %>% 
  ggplot(aes(x = sample_slope)) + 
  geom_density() + 
  geom_vline(xintercept = 0.957, color = "red")
```    


### Part c

Thinking back to the 68-95-99.7 rule, visually approximate the *standard deviation* among the sample slopes.




\
\
\
\



## Exercise 6: Standard error {-}

You've likely observed that the *typical* or mean slope estimate is roughly equal to the "true" population slope parameter of 0.957:

```{r}
results %>% 
  summarize(mean(sample_slope))
```

Thus the *standard deviation* of the slope estimates measures how *far* we might expect an estimate to fall from the (assumed) population slope parameter.
That is, it measures the typical or **standard error** in our sample estimates:

```{r}
results %>% 
  summarize(sd(sample_slope))
```

### Part a {-}

Recall *your* sample estimate of the slope.
How far is it from the population slope, 0.957?
    
```{r}

```


### Part b {-}

How many *standard errors* does your estimate fall from the population slope?
That is, what's your **Z-score**?        

```{r}

```
    

### Part c {-}

Reflecting upon your Z-score, do you think your sample estimate was one of the "lucky" ones, or one of the "unlucky" ones?
(How would you define "lucky"?)



\
\
\
\




## Exercise 7: Work on the project! {-}


- Find project info [here](https://moodle.macalester.edu/course/section.php?id=128297)

- I.d. a dataset (options are listed in the "Final project info & description" on Moodle)

- Check out and maybe even start working on Milestone 1!
   
  
    


## Done!

- Finalize your notes: (1) Render your notes to an HTML file; (2) Inspect this HTML in your Viewer -- check that your work translated correctly; and (3) Outside RStudio, navigate to your 'Activities' subfolder within your 'STAT155' folder and locate the HTML file. Open and inspect this in a web browser.
- Clean up your RStudio session: End the rendering process by clicking the 'Stop' button in the 'Background Jobs' tab in the lower left pane.
- Check the wrap-up announcements and solutions in the online course manual, at the bottom of the corresponding chapter.
- Work on homework!


