In practice, we often have specific hypotheses about the population of interest. We can test these hypotheses using our sample data!!
By the end of this lesson, you should be able to:
This is a discovery activity, so no readings/videos today.
Let’s return to the fish dataset. Recall that rivers contain small concentrations of mercury which can accumulate in fish. Scientists studied this phenomenon among 171 largemouth bass in the Wacamaw and Lumber rivers of North Carolina, recording the following:
| variable | meaning |
|---|---|
| River | Lumber or Wacamaw |
| Station | Station number where the fish was caught (0, 1, …, 15) |
| Length | Fish’s length (in centimeters) |
| Weight | Fish’s weight (in grams) |
| Concen | Fish’s mercury concentration (in parts per million; ppm) |
# Load the data & packages
library(tidyverse)
fish <- read_csv("https://mac-stat.github.io/data/Mercury.csv")
head(fish)# A tibble: 6 × 5
River Station Length Weight Concen
<chr> <dbl> <dbl> <dbl> <dbl>
1 Lumber 0 47 1616 1.6
2 Lumber 0 48.7 1862 1.5
3 Lumber 0 55.7 2855 1.7
4 Lumber 0 45.2 1199 0.73
5 Lumber 0 44.7 1320 0.56
6 Lumber 0 43.8 1225 0.51GOAL
Use this sample data to make inferences about the relationship of mercury concentration (Concen) with the size of a fish, as measured by its Length. This relationship can be described by the population model:
E[Concen | Length] = \beta_0 + \beta_1 Length
Our sample model is below:
E[Concen | Length] = \hat{\beta}_0 + \hat{\beta}_1 Length = -1.132 + 0.058 Length
fish_mod_1 <- lm(Concen ~ Length, data = fish)
summary(fish_mod_1)
Call:
lm(formula = Concen ~ Length, data = fish)
Residuals:
Min 1Q Median 3Q Max
-1.1499 -0.3436 -0.1022 0.3123 1.9100
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.131645 0.213615 -5.298 3.62e-07 ***
Length 0.058127 0.005228 11.119 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5805 on 169 degrees of freedom
Multiple R-squared: 0.4225, Adjusted R-squared: 0.4191
F-statistic: 123.6 on 1 and 169 DF, p-value: < 2.2e-16Notation:
Interpret \hat{\beta}_1.
Using the 68-95-99.7 rule, construct an approximate 95% confidence interval for \beta_1. That is, calculate an interval estimate of \beta_1.
Compare your CI to an exact 95% CI:
___(___, level = 0.95)Our work above provides point and interval estimates of \beta_1. But we also want to test hypotheses related to \beta_1.
Our fundamental hypothesis here is that, among the fish population, mercury levels are associated with fish length – there’s a relationship between the two.
How can we state this hypothesis using notation? Choose 1.
NOTE: We’ll refer to 0 as the null value to which we’re comparing our population parameter of interest.
Since we don’t have access to the entire fish population, we must evaluate our hypothesis using our sample data. We have to ask: does our sample data provide enough evidence to support the hypothesis that, among the fish population, not just our sample, mercury levels are associated with fish length?
Explain why we can’t use our simple estimate \hat{\beta}_1 to evaluate this hypothesis, but we can use the 95% CI for \beta_1. (This is true in general, not just this fish example!)
So, does your 95% CI from Example 1 support the hypothesis that there’s a relationship between mercury levels and fish length? Explain.
Similarly, do the confidence bands below support the hypothesis?
fish %>%
ggplot(aes(y = Concen, x = Length)) +
geom_smooth(method = "lm") 
FORMAL HYPOTHESIS TESTS
Though CIs can help us evaluate hypotheses, hypothesis tests provide a more formal approach. The general idea:
Or, in context:
NOTE
The goal of these exercises is to explore concepts in hypothesis testing. There’s some new code that helps with this exploration, but the code is not the point!! Throughout, focus on the goal of the code and what its output tells us, NOT the details of the code.
Our results will start to use more scientific notation. Examples:
3e-02 = 0.03 (move the decimal 2 places to the left)3e02 = 300 (move the decimal 2 places to the right)4.12e-07 = 0.000000412 (move the decimal 7 places to the left)4.12e03 = 4120 (move the decimal 3 places to the right)< 2e-16 = the result is less than 0.0000000000000002 (it’s really really small!)To begin testing our hypothesis that there is an association between mercury concentration and fish length (\beta_1 \ne 0), let’s explore how our sample data might behave IF in fact there were no association (\beta_1 = 0). We’ll first do this using simulation, and then through theory.
Suppose that when each individual fish in our sample was collected, the scientists wrote its Concen and Length on a piece of paper.
They then ripped the paper in half, with the observed Concen on one half and the observed Length on the other.
They then accidentally dropped the pieces of paper! In the resulting mess, they couldn’t tell which Length observation went with each Concen observation.
To cover up their mistake, they just randomly matched up the Length and Concen pieces of paper. Essentially, in this shuffling process, the scientists broke the bond of the original Length and Concen pairs.
Question: Intuitively, if the researchers plotted the new pairs of Concen and Length values, what do you think their scatterplot would look like? What about their sample model?
# First check out the first 6 fish
fish %>%
select(Concen, Length) %>%
head()Now, randomly shuffle the Length values! That is, break the original Concen and Length pairs, and just randomly assign Length:
# Run this chunk a few times!!
# Shuffle the Length values
fish %>%
select(Concen, Length) %>%
mutate(Length = sample(Length, size = length(Length), replace = FALSE)) %>%
head()Let’s check out the sample models that arise from the shuffled Concen and Length pairs. Summarize your observations! Was your intuition in Part a correct?
# Shuffle the Length values
# Run this chunk a few times!!
# Then plot the resulting sample data and model
fish %>%
select(Concen, Length) %>%
mutate(Length = sample(Length, size = length(Length), replace = FALSE)) %>%
ggplot(aes(x = Length, y = Concen)) +
geom_point() +
geom_smooth(method = "lm")Concen and Length in the broader population of fish. To get a sense of the range of possible outcomes in this scenario, let’s simulate a bunch of shuffled samples. Specifically, take 500 different shuffled samples and use each to estimate the model of Concen by Length (gray lines). Summarize your observations!! If there were truly no relationship between Concen and Length, how would we expect the sample data to behave?set.seed(1)
shuffled_models <- map_df(1:500, function(i){
fish %>%
select(Concen, Length) %>%
sample_n(size = length(Length), replace = TRUE) %>%
mutate(Length = sample(Length, size = length(Length), replace = FALSE)) %>%
lm(Concen ~ Length, data = .) %>%
coef()
}
)
head(shuffled_models)# Plot the 500 shuffled models
fish %>%
ggplot(aes(x = Length, y = Concen)) +
geom_abline(data = shuffled_models,
aes(intercept = `(Intercept)`, slope = Length),
color = "gray", size = 0.25) +
geom_smooth(method = "lm", se = FALSE, size = 0) # Ignore this line. It's a clunky workaroundConcen and Length, i.e. \beta_1 were truly 0, how would we expect the sample slopes to behave?shuffled_models %>%
ggplot(aes(x = Length)) +
geom_density()Now that we’ve simulated how sample data might behave if there were truly no relationship between Concen and Length (i.e. if \beta_1 were 0), let’s consider the next important step in a hypothesis test:
Evaluate how compatible our observed sample data are with this assumption that \beta_1 = 0!
fish %>%
ggplot(aes(x = Length, y = Concen)) +
geom_abline(data = shuffled_models,
aes(intercept = `(Intercept)`, slope = Length),
color = "gray", size = 0.25) +
geom_smooth(method = "lm", se = FALSE)shuffled_models %>%
ggplot(aes(x = Length)) +
geom_density() +
geom_vline(xintercept = 0.05813, color = "blue")Now that we’ve built up some intuition using simulation, let’s formalize these concepts with some theory.
You shouldn’t be surprised to know that sampling distributions are at the heart of hypothesis testing! The distribution of slopes from our shuffled sample models uses simulation to approximate the sampling distribution of slope estimates \hat{\beta}_1 we’d expect if the population slope \beta_1 were actually 0 (the null value):
shuffled_models %>%
ggplot(aes(x = Length)) +
geom_density()We can also approximate the sampling distribution using mathematical theory, specifically the Central Limit Theorem. If \beta_1 were truly 0, then the distribution of possible sample estimates \hat{\beta}_1 would be Normally distributed around 0:
\hat{\beta}_1 \sim N(0, SE(\hat{\beta}_1)^2)
where in Example 1 we approximated the standard error to be 0.005. Thus a plot of the CLT is below. Confirm that this is similar to the simulated sampling distribution above:
# Don't worry about this code!!!
# Save your plot
clt_plot <- data.frame(x = 0 + c(-4:4)*0.005) %>%
mutate(y = dnorm(x, sd = 0.005)) %>%
ggplot(aes(x = x)) +
stat_function(fun = dnorm, args = list(mean = 0, sd = 0.005)) +
geom_segment(aes(x = x, xend = x, y = 0, yend = y), linetype = "dashed") +
scale_x_continuous(breaks = c(-4:4)*0.005) +
xlab("")
clt_plotOur sample slope estimate of \hat{\beta}_1 = 0.05813 is shown on the sampling distribution of sample slopes \hat{\beta}_1 that we’d expect if the population slope \beta_1 were actually 0 (the null value):
clt_plot +
geom_vline(xintercept = 0.05813, color = "blue")To measure how compatible (or incompatible) our estimate is with the null value, we can measure the distance between the two.
In the context of this particular analysis, is that a big or a small number? How can you tell? That is, what info are you using to make this assessment?
Would your answer be the same if you observed this estimate in another context?
\frac{\hat{\beta}_1 - \text{null value}}{SE(\hat{\beta}_1))} = \text{ number of SE that $\hat{\beta}_1$ falls from the null value}
This is called the test statistic. Calculate and interpret the test statistic for our example.
coef(summary(fish_mod_1))Our sample data is not consistent with the null value of \beta_1 = 0, i.e. the idea that there’s no significant association between mercury concentration and length.
Our sample data is consistent with the null value of \beta_1 = 0.
Again, our sample slope estimate of \hat{\beta}_1 = 0.05813 is shown on the sampling distribution of sample slopes \hat{\beta}_1 that we’d expect if the population slope \beta_1 were actually 0 (the null value):
clt_plot +
geom_vline(xintercept = 0.05813, color = "blue")Use the 68-95-99.7 Rule with the plot above to approximate the p-value:
Pr(>|t|) column. What is it?coef(summary(fish_mod_1))It’s very unlikely that we’d have observed such a steep increase in Concen with Length among our sample fish “by chance”, i.e. if in fact there were no relationship between mercury concentration and length in the broader fish population.
Given what we observed in our sample data, it’s very unlikely that mercury concentration is associated with length.
Given what we observed in our sample data, it’s very unlikely that mercury concentration and length are unrelated.
PAUSE: Check the solutions to this question! p-values are commonly misinterpreted.
The fish in our sample come from two different rivers: Lumber and Wacamaw. Research question: Is there evidence that the mercury concentration in fish differs between the 2 rivers?
To answer this question, we’ll control for fish Length – we don’t want to mistakenly detect a difference in mercury levels just because we happened to sample bigger fish in one river and smaller fish in the other. The relevant population model is:
E[Concen | River, Length] = \beta_0 + \beta_1 RiverWacamaw + \beta_2 Length
To address the above research question, which model coefficient / population parameter is of primary interest: \beta_0, \beta_1, or \beta_2? And what’s the null value of this coefficient? That is, what would this coefficient be if there were truly no difference in mercury concentration between the 2 rivers (when controlling for length)?
For the parameter of interest, obtain and report a sample estimate and its corresponding standard error:
fish_mod_2 <- lm(___, data = fish)
coef(summary(fish_mod_2))Adjust the code below to sketch the appropriate sampling distribution. Represent our sample estimate by a blue line.
# Put OUR sample estimate here
# Round to 3 digits
est <- ___
# Put the corresponding standard error here
# Round to 3 digits
se <- ___
data.frame(x = 0 + c(-4:4)*se) %>%
mutate(y = dnorm(x, sd = se)) %>%
ggplot(aes(x = x)) +
stat_function(fun = dnorm, args = list(mean = 0, sd = se)) +
geom_segment(aes(x = x, xend = x, y = 0, yend = y), linetype = "dashed") +
scale_x_continuous(breaks = c(-4:4)*se) +
geom_vline(xintercept = est, color = "blue") +
labs(y = "density", x = "possible estimates IF the null value were true") Remember the next step: We must evaluate how compatible our sample estimate is with the null value. We’ll take 3 approaches: using a visual assessment, test statistic, and p-value.
Visually, in Part c of the previous exercise, does our sample estimate seem compatible with the null value? Mainly, is our estimate consistent with what we’d expect to observe if there were truly no difference in the mercury concentration between the 2 rivers (when controlling for fish length)?
Let’s focus on the test statistic:
coef(summary(fish_mod_2))Use the 68-95-99.7 Rule with the plot in Part c of the previous exercise to approximate the p-value:
Now obtain and report the exact p-value from the model summary table.
Do your test statistic and p-value indicate that your estimate is or is not compatible with the null value? Explain.
Putting this all together, what do you conclude? When controlling for length…
# Load the data & packages
library(tidyverse)
fish <- read_csv("https://mac-stat.github.io/data/Mercury.csv")
head(fish)# A tibble: 6 × 5
River Station Length Weight Concen
<chr> <dbl> <dbl> <dbl> <dbl>
1 Lumber 0 47 1616 1.6
2 Lumber 0 48.7 1862 1.5
3 Lumber 0 55.7 2855 1.7
4 Lumber 0 45.2 1199 0.73
5 Lumber 0 44.7 1320 0.56
6 Lumber 0 43.8 1225 0.51fish_mod_1 <- lm(Concen ~ Length, data = fish)
summary(fish_mod_1)
Call:
lm(formula = Concen ~ Length, data = fish)
Residuals:
Min 1Q Median 3Q Max
-1.1499 -0.3436 -0.1022 0.3123 1.9100
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.131645 0.213615 -5.298 3.62e-07 ***
Length 0.058127 0.005228 11.119 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5805 on 169 degrees of freedom
Multiple R-squared: 0.4225, Adjusted R-squared: 0.4191
F-statistic: 123.6 on 1 and 169 DF, p-value: < 2.2e-16We estimate that for every 1cm increase in length, the expected mercury concentration increases by 0.058ppm.
0.058 \pm 2*0.005 = (0.058 - 0.010, 0.058 + 0.010) = (0.048, 0.068)
confint(fish_mod_1, level = 0.95) 2.5 % 97.5 %
(Intercept) -1.55334248 -0.70994835
Length 0.04780769 0.06844728The simple estimate \hat{\beta}_1 has error that we need to account for when making inferences about the population. The CI accounts for this error!
Yes, the interval doesn’t include 0. Thus 0 is not a plausible value for \beta_1.
Yes, a 0-slope line cannot be drawn within these confidence bands.
fish %>%
ggplot(aes(y = Concen, x = Length)) +
geom_smooth(method = "lm") 
Intuition.
# First check out the first 6 fish
fish %>%
select(Concen, Length) %>%
head()# A tibble: 6 × 2
Concen Length
<dbl> <dbl>
1 1.6 47
2 1.5 48.7
3 1.7 55.7
4 0.73 45.2
5 0.56 44.7
6 0.51 43.8# Run this chunk a few times!!
# Shuffle the Length values
fish %>%
select(Concen, Length) %>%
mutate(Length = sample(Length, size = length(Length), replace = FALSE)) %>%
head()# A tibble: 6 × 2
Concen Length
<dbl> <dbl>
1 1.6 33.5
2 1.5 34.8
3 1.7 34.9
4 0.73 37.4
5 0.56 45.4
6 0.51 40 # Shuffle the Length values
# Run this chunk a few times!!
# Then plot the resulting sample data and model
fish %>%
select(Concen, Length) %>%
mutate(Length = sample(Length, size = length(Length), replace = FALSE)) %>%
ggplot(aes(x = Length, y = Concen)) +
geom_point() +
geom_smooth(method = "lm")
Concen and Length, we’d expect the sample model to have a slope near to (but not exactly) 0.set.seed(1)
shuffled_models <- map_df(1:500, function(i){
fish %>%
select(Concen, Length) %>%
sample_n(size = length(Length), replace = TRUE) %>%
mutate(Length = sample(Length, size = length(Length), replace = FALSE)) %>%
lm(Concen ~ Length, data = .) %>%
coef()
}
)
head(shuffled_models)# A tibble: 6 × 2
`(Intercept)` Length
<dbl> <dbl>
1 0.917 0.0114
2 0.799 0.00754
3 0.958 0.00374
4 0.827 0.00897
5 0.763 0.0118
6 1.08 -0.000455fish %>%
ggplot(aes(x = Length, y = Concen)) +
geom_abline(data = shuffled_models,
aes(intercept = `(Intercept)`, slope = Length),
color = "gray", size = 0.25) +
geom_smooth(method = "lm", se = FALSE, size = 0) # Ignore this line. It's a clunky workaround
shuffled_models %>%
ggplot(aes(x = Length)) +
geom_density()
fish %>%
ggplot(aes(x = Length, y = Concen)) +
geom_abline(data = shuffled_models,
aes(intercept = `(Intercept)`, slope = Length),
color = "gray", size = 0.25) +
geom_smooth(method = "lm", se = FALSE)
shuffled_models %>%
ggplot(aes(x = Length)) +
geom_density() +
geom_vline(xintercept = 0.05813, color = "blue")
The simulated and CLT / theory-informed sampling distributions are similar!
shuffled_models %>%
ggplot(aes(x = Length)) +
geom_density()
clt_plot <- data.frame(x = 0 + c(-4:4)*0.005) %>%
mutate(y = dnorm(x, sd = 0.005)) %>%
ggplot(aes(x = x)) +
stat_function(fun = dnorm, args = list(mean = 0, sd = 0.005)) +
geom_segment(aes(x = x, xend = x, y = 0, yend = y), linetype = "dashed") +
scale_x_continuous(breaks = c(-4:4)*0.005) +
xlab("")
clt_plot
clt_plot +
geom_vline(xintercept = 0.05813, color = "blue")
How far is our sample slope estimate of \hat{\beta}_1 = 0.05813 from the null value of 0? 0.05813 :)
In the context of this particular analysis, is that a big or a small number? How can you tell? This seems pretty big relative to the standard error.
(0.058127 - 0) / 0.005228[1] 11.1184t value column:coef(summary(fish_mod_1)) Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.13164542 0.213614796 -5.297598 3.617750e-07
Length 0.05812749 0.005227593 11.119359 6.641225e-22clt_plot +
geom_vline(xintercept = 0.05813, color = "blue")
coef(summary(fish_mod_1)) Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.13164542 0.213614796 -5.297598 3.617750e-07
Length 0.05812749 0.005227593 11.119359 6.641225e-22Concen with Length among our sample fish “by chance”, i.e. if in fact there were no relationship between mercury concentration and length in the broader fish population.The null value is \beta_1 = 0
fish_mod_2 <- lm(Concen ~ River + Length, data = fish)
coef(summary(fish_mod_2)) Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.19422929 0.216287220 -5.521497 1.255595e-07
RiverWacamaw 0.14202739 0.089496454 1.586961 1.144018e-01
Length 0.05765684 0.005212709 11.060821 1.033822e-21# Put OUR sample estimate here
# Round to 3 digits
est <- 0.142
# Put the corresponding standard error here
# Round to 3 digits
se <- 0.089
data.frame(x = 0 + c(-4:4)*se) %>%
mutate(y = dnorm(x, sd = se)) %>%
ggplot(aes(x = x)) +
stat_function(fun = dnorm, args = list(mean = 0, sd = se)) +
geom_segment(aes(x = x, xend = x, y = 0, yend = y), linetype = "dashed") +
scale_x_continuous(breaks = c(-4:4)*se) +
geom_vline(xintercept = est, color = "blue") +
labs(y = "density", x = "possible estimates IF the null value were true") 
Yes! Our estimate is relatively close to 0 on this scale (within 2 SE).
Our estimate of \beta_1, 0.142, falls only 1.59 SE from 0:
(0.142 - 0) / 0.089[1] 1.595506coef(summary(fish_mod_2)) Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.19422929 0.216287220 -5.521497 1.255595e-07
RiverWacamaw 0.14202739 0.089496454 1.586961 1.144018e-01
Length 0.05765684 0.005212709 11.060821 1.033822e-21An exact p-value is 0.114:
coef(summary(fish_mod_2)) Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.19422929 0.216287220 -5.521497 1.255595e-07
RiverWacamaw 0.14202739 0.089496454 1.586961 1.144018e-01
Length 0.05765684 0.005212709 11.060821 1.033822e-21It is. It is not very far from 0 (when considering SE) and would not be unlikely to observe if the null value were true.
When controlling for length…
Rough version:
0.142 \pm 2*0.089 = (-0.036, 0.320)
Exact version:
confint(fish_mod_2) 2.5 % 97.5 %
(Intercept) -1.62122031 -0.7672383
RiverWacamaw -0.03465518 0.3187100
Length 0.04736599 0.0679477When controlling for length…