Hypothesis Testing- Discovery

Notes and in-class exercises

You can download the .qmd file for this activity here and open in R-studio. The rendered version is posted in the course website (Activities tab). I often experiment with the class activities (and see it in live!) and make updates, but I always post the final version before class starts. To be sure you have the most up-to-date copy, please download it once you’ve settled in before class begins.

Notes

  • You can download a template file to work with here.
  • File organization: Save this file in the “Activities” subfolder of your “STAT155” folder.

Learning goals

By the end of this lesson, you should be able to:

  • Understand how standard errors and confidence intervals enable us to make statistical inferences
  • Articulate how we can formalize a research question as a testable, statistical hypothesis

Readings and videos

This is a discovery activity, so no assigned readings/videos today.

Exercises

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.51

Exercise 1

Research question: Is there evidence that the mercury concentration in fish (Concen) differs according to the River they were sampled from?

part a: fit the model

Fit a simple linear regression model that would address our research question

mod_fish <- lm(Concen ~ River, data=fish)
summary(mod_fish)
## 
## Call:
## lm(formula = Concen ~ River, data = fish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1664 -0.5681 -0.1764  0.4219  2.4219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.07808    0.08866  12.160   <2e-16 ***
## RiverWacamaw  0.19835    0.11712   1.694   0.0922 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7575 on 169 degrees of freedom
## Multiple R-squared:  0.01669,    Adjusted R-squared:  0.01087 
## F-statistic: 2.868 on 1 and 169 DF,  p-value: 0.09218

Interpret the intercept from this model.

Our model estimates an average mercury concentration of 1.078ppm among fish in the Lumber River.

part b: construct a CI

Using the 68-95-99.7 rule, construct an approximate 95% confidence interval for the intercept term, and provide an appropriate interpretation.

# approximate 95% confidence interval

1.078 + 2*0.089
## [1] 1.256
1.078 - 2*0.089
## [1] 0.9

Preferred interpretation: It is plausible that the true mean mercury concentration among fish in the Lumber River is between 0.90ppm and 1.25ppm.

(technical addendum to this interpretation): …specifically, we expect that if we take many different samples and obtain a set of corresponding parameter estimates and confidence intervals, we expect that 95% of the resulting intervals will contain the true mean mercury concentration of the entire Lumber River fish population. We hope that our interval is one of the lucky 95% and not one of the unlucky 5% that don’t contain the true population parameter.

Not as preferred interpretation, but still okay for this course: We are 95% confident that the mean mercury concentration among fish in the Lumber River is between 0.90ppm and 1.25ppm.

Compare your CI to an exact 95% confidence interval for the model coefficients:

confint(mod_fish, level=0.95)
##                    2.5 %    97.5 %
## (Intercept)   0.90305825 1.2531061
## RiverWacamaw -0.03285077 0.4295435

part c: what can we conclude from multiple samples?

Suppose we take 200 different samples of fish from the Lumber River. Based on these results, in how many of those samples would you expect to observe mean mercury concentration greater than 1.25ppm?

We don’t/can’t actually know! This depends on the true population parameter and the accuracy of our sampling distribution.

What we can say is that if our sampling distribution model is accurate, then we should expect that about 10 out of 200 samples (5% of them) will produce confidence intervals that don’t contain the population parameter. We should expect that half of these–so 5 samples–are overestimates and the other half are underestimates.

part d: intuition for constructing & interpreting test statistics

Suppose previous environmental studies have found little evidence of mercury pollution in other rivers in the area, so perhaps our “default” assumption is that fish from the Lumber river should have an expected mercury concentration of 0ppm. How many standard errors is our sample estimate (1.078ppm) away from this expectation? What are three possible conclusions?

If we assume that 0ppm is the “true” mercury concentration, then our estimate of Beta_0 = 1.078ppm with a standard error of 0.08866 means that our estimate is (1.07808-0)/0.08866 = 12.16 standard errors away from what we should expect.

Possible conclusions:

  1. Our working assumption that 0ppm should be the “true” mercury concentration in the Lumber river fish population was wrong! The confidence interval we constructed above suggests that a true value of 0ppm is extremely implausible.
  1. Perhaps ~0ppm is actually the true average mercury concentration in the population, we just got extremely, outrageously unlucky with our sample.
  1. Perhaps there was a measurement/data entry error, and the units are actually parts per billion, not million.

part e: do individual observations contradict our conclusions?

Now suppose we sample a single fish from the Lumber River and find it has a mercury concentration of 2.5ppm. Are you surprised by this result? Why or why not? (Hint: create a code chunk that calculates the mean, standard deviation, and maximum of the Concen variable in each river in our original sample)

fish %>% 
  group_by(River) %>% 
  summarise(mean=mean(Concen), 
            sd=sd(Concen), 
            max=max(Concen))
## # A tibble: 2 × 4
##   River    mean    sd   max
##   <chr>   <dbl> <dbl> <dbl>
## 1 Lumber   1.08 0.649   3.5
## 2 Wacamaw  1.28 0.829   3.6

Observing a single fish with a mercury concentration of 2.5ppm is actually not that surprising! 2.5ppm is a little more than 2 standard deviations away from the mean mercury concentration in our sample of fish from the Lumber River (1.08+2*0.64=2.43), but there are certainly fish in the sample with even higher mercury concentrations (max=3.5ppm), so this isn’t outside the bounds of what we’d expect.

Class Notes

Additional Reading (after class): FORMAL HYPOTHESIS TESTS

Though CIs can help us evaluate hypotheses, hypothesis tests provide a more formal approach. The general idea:

  1. Assume that the hypothesis is not true.
  2. Determine what sample observations would be expected under this assumption.
  3. Evaluate how compatible our observed sample data are with this assumption (are they consistent with what we’d expect?).

Or, in context:

  1. Assume there were truly no association between mercury concentration and fish length, i.e. \(\beta_1 = 0\) (the “null value”).
  2. Determine how our fish sample might have “behaved” if there were truly no association.
  3. Evaluate how compatible our fish sample observations are with the null value (are they consistent with what we’d expect if there were no association?).

Exercise 2:

Let’s look at the model summary output again:

summary(mod_fish)
## 
## Call:
## lm(formula = Concen ~ River, data = fish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1664 -0.5681 -0.1764  0.4219  2.4219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.07808    0.08866  12.160   <2e-16 ***
## RiverWacamaw  0.19835    0.11712   1.694   0.0922 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7575 on 169 degrees of freedom
## Multiple R-squared:  0.01669,    Adjusted R-squared:  0.01087 
## F-statistic: 2.868 on 1 and 169 DF,  p-value: 0.09218

part a: interpret model coefficient

Now, let’s interpret the RiverWacamaw coefficient. Based only on the coefficient (don’t think about the standard error yet), what can we say about the difference in mercury concentration among fish in the two rivers?

Response

part b: construct a CI

Using the 68-95-99.7 rule, construct an approximate 95% confidence interval for the RiverWacamaw coefficient, and provide an appropriate interpretation.

Response

part c: interpreting the CI

Do you believe it plausible that the mean mercury concentration of the fish population in the Wacamaw River is approximately the same as that of the fish population in the Lumber River? How would you confirm this? What assumptions are you making?

Response

part d: effect of sample size on our conclusions

Suppose we sample 10 times as many fish from the Wacamaw River, and get a similar coefficient estimate (0.2). Thinking back to the Central Limit Theorem, what should happen to the standard error of the RiverWacamaw coefficient? How small of a standard error would we need to more conclusively say that there is an actual difference in mean mercury concentrations of the Lumber River and Wacamaw River fish populations?

Response

part e: reconciling parameter estimates and uncertainty

Suppose the true population coefficient for the RiverWacamawparameter is 0.02 (i.e. the average mercury concentration is 0.02ppm higher for the Wacamaw River fish population compared to that of the Lumber River). Is this meaningful?

Response

part f (CHALLENGE)

Using the model summary output, report the mean mercury concentration for our sample of fish from the Wacamaw River:

summary(mod_fish)
## 
## Call:
## lm(formula = Concen ~ River, data = fish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1664 -0.5681 -0.1764  0.4219  2.4219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.07808    0.08866  12.160   <2e-16 ***
## RiverWacamaw  0.19835    0.11712   1.694   0.0922 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7575 on 169 degrees of freedom
## Multiple R-squared:  0.01669,    Adjusted R-squared:  0.01087 
## F-statistic: 2.868 on 1 and 169 DF,  p-value: 0.09218

Response:

Which of the following values do you think is the standard error of the sample mean for the Wacamaw River?

  • 0.11712
  • 0.08866
  • 0.11712 + 0.08866 = 0.20578
  • 0.11712 - 0.08866 = 0.02846
  • something else

To answer this question, look at the code chunk below, which fits the same model, but uses the Wacamaw River as our reference category instead of the Lumber River:

mod_fish2 <- lm(Concen ~ River, data=fish %>% mutate(River=ifelse(River == "Wacamaw", paste0("_", River), River)))
summary(mod_fish2)
## 
## Call:
## lm(formula = Concen ~ River, data = fish %>% mutate(River = ifelse(River == 
##     "Wacamaw", paste0("_", River), River)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1664 -0.5681 -0.1764  0.4219  2.4219 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.27643    0.07652  16.681   <2e-16 ***
## RiverLumber -0.19835    0.11712  -1.694   0.0922 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7575 on 169 degrees of freedom
## Multiple R-squared:  0.01669,    Adjusted R-squared:  0.01087 
## F-statistic: 2.868 on 1 and 169 DF,  p-value: 0.09218

Compare this to the output for mod_fish. What do you notice about the standard errors of the intercepts (i.e., the standard errors of the means for each river) compared to the standard errors of the RiverWacamaw and RiverLumber coefficients (i.e., the standard errors of the differences between the means)?

Response:

Additional Exercise

NOTE

  • The goal of these additional 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!)

Exercise 1: Simulating data under the null value

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.

  1. Read carefully:

  • 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?

  1. We can simulate this idea in RStudio! Throughout, focus on concepts over code.
# 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

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()
## # A tibble: 6 × 2
##   Concen Length
##    <dbl>  <dbl>
## 1   1.6    33  
## 2   1.5    48  
## 3   1.7    40.4
## 4   0.73   27.3
## 5   0.56   35.5
## 6   0.51   31

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")

  1. These shuffled samples give us a sense of the sample models we’d expect IF \(\beta_1\) were actually 0, i.e. if there were no relationship between 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 <- mosaic::do(500)*(
  fish %>% 
  select(Concen, Length) %>% 
  sample_n(size = length(Length), replace = TRUE) %>% 
  mutate(Length = sample(Length, size = length(Length), replace = FALSE)) %>%
  with(lm(Concen ~ Length))
)
head(shuffled_models)
##   Intercept     Length     sigma  r.squared        F numdf dendf .row .index
## 1 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      1
## 2 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      2
## 3 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      3
## 4 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      4
## 5 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      5
## 6 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      6
# 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 workaround

  1. Finally, focus on just the slope estimates \(\hat{\beta}_1\) behind the lines above. Summarize your observations!! If there were truly no relationship between Concen 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()

Exercise 2: Comparing our sample results to the null value (intuition)

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\)!

  1. Let’s start with a visual assessment. Is our sample model (blue line) consistent / compatible with the “null models” simulated under the assumption the \(\beta_1 = 0\) (gray lines)?
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)

  1. Consider another visual assessment, focusing just on the slopes of the lines in the plot above, i.e. our sample estimates \(\hat{\beta}_1\). Is our sample slope of \(\hat{\beta}_1 = 0.05813\) (blue line) consistent / compatible with the collection of slopes from the “null models” (density plot)?
shuffled_models %>% 
  ggplot(aes(x = Length)) + 
  geom_density() + 
  geom_vline(xintercept = 0.05813, color = "blue")

  1. In Part b, you evaluated the compatibility of our sample slope \(\hat{\beta}_1 = 0.05813\) model with slopes of the “null models” using visual cues alone. But for a formal hypothesis test, we need some numerical measures of compatibility. Brainstorm some ideas!

Exercise 3: CLT

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, s.e.(\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_plot

Exercise 4: Comparing our sample results to the null value (test statistic)

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")

To measure how compatible (or incompatible) our estimate is with the null value, we can measure the distance between the two.

  1. Our sample slope estimate of \(\hat{\beta}_1 = 0.05813\) is 0.05813 units (ppm/cm) from the null value of 0.

  • 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?

  1. In practice, we standardize our distance calculations so that they have the same meaning no matter the context of the analysis. We can do this by calculating a z-score:

\[ \frac{\hat{\beta}_1 - \text{null value}}{s.e.(\hat{\beta}_1))} = \text{ number of s.e. that $\hat{\beta}_1$ falls from the null value} \]

This is called the test statistic. Calculate and interpret the test statistic for our example.

  1. Within rounding, the same test statistic appears in our model summary table! Where is it?!
coef(summary(fish_mod_1))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'fish_mod_1' not found
  1. What can we conclude from your test statistic calculation and interpretation:

  • 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\).

Exercise 5: Comparing our sample results to the null value (p-value)

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")

  1. Another way to evaluate the compatibility of our estimate with the null value of \(\beta_1 = 0\) is to ask how likely we are to have gotten this estimate if the null value were indeed true! That is, what’s the probability that we would have gotten a sample slope that’s at least 0.058 above or below 0 IF in fact there were no association between mercury concentration and length? Equivalently, what’s the probability that we would have gotten a test statistic so far from 0 IF in fact \(\beta_1\) were 0? Use the 68-95-99.7 Rule with the plot above to approximate the p-value:
  • less than 0.003
  • between 0.003 and 0.05
  • between 0.05 and 0.32
  • bigger than 0.32
  1. The calculation above is called a p-value. A more accurate p-value is reported in the model summary table in the Pr(>|t|) column. What is it?
coef(summary(fish_mod_1))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'fish_mod_1' not found
  1. How can we interpret the p-value? IMPORTANT: Think about the steps that went into calculating the p-value!!

  • 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.

Exercise 6: River test Part I

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 \]

  1. 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)?

  2. 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))
## Error in parse(text = input): <text>:1:19: unexpected input
## 1: fish_mod_2 <- lm(__
##                       ^
  1. Let’s use this data to evaluate our research question. To do so, remember the next step: Assume that the null value from Part a is correct and determine what sample estimates we would expect to get in this scenario.

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") 
## Error in parse(text = input): <text>:3:9: unexpected input
## 2: # Round to 3 digits
## 3: est <- __
##            ^

Exercise 7: River test Part II

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.

  1. 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)?

  2. Let’s focus on the test statistic:

  • Obtain the test statistic from the model summary table below.
  • Show how this was calculated from the sample estimate and its standard error.
  • Interpret the test statistic.
coef(summary(fish_mod_2))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'fish_mod_2' not found
  1. Recall that the p-value calculates the probability that, IF in fact the null value were correct (i.e. there was actually no difference in mercury concentration by river when controlling for length), we would have gotten a sample estimate that’s as far from the null value as ours, either above or below. Use the 68-95-99.7 Rule with the plot in Part c of the previous exercise to approximate the p-value:
  • less than 0.003
  • between 0.003 and 0.05
  • between 0.05 and 0.32
  • bigger than 0.32

Now obtain and report the exact p-value from the model summary table.

  1. Do your test statistic and p-value indicate that your estimate is or is not compatible with the null value? Explain.

  2. Putting this all together, what do you conclude? When controlling for length…

  • we have evidence of a statistically significant difference in mercury concentration by river.
  • we do not have sufficient evidence to conclude that there’s a statistically significant difference in mercury concentration by river.

Exercise 8: River CI

  1. Finally, construct a 95% CI for the population parameter of interest in our river analysis.
  1. What do you conclude from this CI? When controlling for length…
  • we have evidence of a statistically significant difference in mercury concentration by river.
  • we do not have sufficient evidence to conclude that there’s a statistically significant difference in mercury concentration by river.
  1. Your conclusions from the hypothesis test (test statistic and p-value) should agree with those from your CI! If this is not the case, review your work.





Solutions

Exercise 1

Research question: Is there evidence that the mercury concentration in fish (Concen) differs according to the River they were sampled from?

part a: fit the model

Fit a simple linear regression model that would address our research question

mod_fish <- lm(Concen ~ River, data=fish)
summary(mod_fish)
## 
## Call:
## lm(formula = Concen ~ River, data = fish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1664 -0.5681 -0.1764  0.4219  2.4219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.07808    0.08866  12.160   <2e-16 ***
## RiverWacamaw  0.19835    0.11712   1.694   0.0922 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7575 on 169 degrees of freedom
## Multiple R-squared:  0.01669,    Adjusted R-squared:  0.01087 
## F-statistic: 2.868 on 1 and 169 DF,  p-value: 0.09218

Interpret the intercept from this model.

Our model estimates an average mercury concentration of 1.078ppm among fish in the Lumber River.

part b: construct a CI

Using the 68-95-99.7 rule, construct an approximate 95% confidence interval for the intercept term, and provide an appropriate interpretation.

1.078 +/- 2*0.089 –> [0.90, 1.256]

Preferred interpretation: It is plausible that the true mean mercury concentration among fish in the Lumber River is between 0.90ppm and 1.25ppm.

(technical addendum to this interpretation): …specifically, we expect that if we take many different samples and obtain a set of corresponding parameter estimates and confidence intervals, we expect that 95% of the resulting intervals will contain the true mean mercury concentration of the entire Lumber River fish population. We hope that our interval is one of the lucky 95% and not one of the unlucky 5% that don’t contain the true population parameter.

Not as preferred interpretation, but still okay for this course: We are 95% confident that the mean mercury concentration among fish in the Lumber River is between 0.90ppm and 1.25ppm.

Compare your CI to an exact 95% confidence interval for the model coefficients:

confint(mod_fish, level=0.95)
##                    2.5 %    97.5 %
## (Intercept)   0.90305825 1.2531061
## RiverWacamaw -0.03285077 0.4295435

part c: what can we conclude from multiple samples?

Suppose we take 200 different samples of fish from the Lumber River. In how many of those samples would you expect to observe an estimated mean mercury concentration greater than 1.25ppm?

We don’t/can’t actually know! This depends on the true population parameter and the accuracy of our sampling distribution.

What we can say is that if our sampling distribution model is accurate, then we should expect that about 10 out of 200 samples (5% of them) will produce confidence intervals that don’t contain the population parameter. We should expect that half of these–so 5 samples–are overestimates and the other half are underestimates.

part d: intuition for constructing & interpreting test statistics

Suppose previous environmental studies have found little evidence of mercury pollution in other rivers in the area, so perhaps our “default” assumption is that fish from the Lumber river should have an expected mercury concentration of 0ppm. How many standard errors is our sample estimate (1.078ppm) away from this expectation? What are three possible conclusions?

If we assume that 0ppm is the “true” mercury concentration, then our estimate of Beta_0 = 1.078ppm with a standard error of 0.08866 means that our estimate is (1.07808-0)/0.08866 = 12.16 standard errors away from what we should expect.

Possible conclusions:

  1. Our working assumption that 0ppm should be the “true” mercury concentration in the Lumber river fish population was wrong! The confidence interval we constructed above suggests that a true value of 0ppm is extremely implausible.
  1. Perhaps ~0ppm is actually the true average mercury concentration in the population, we just got extremely, outrageously unlucky with our sample.
  1. Perhaps there was a measurement/data entry error, and the units are actually parts per billion, not million.

part e: do individual observations contradict our conclusions?

Now suppose we sample a single fish from the Lumber River and find it has a mercury concentration of 2.5ppm. Are you surprised by this result? Why or why not? (Hint: create a code chunk that calculates the mean, standard deviation, and maximum of the Concen variable in each river in our original sample)

fish %>% 
  group_by(River) %>% 
  summarise(mean=mean(Concen), 
            sd=sd(Concen), 
            max=max(Concen))
## # A tibble: 2 × 4
##   River    mean    sd   max
##   <chr>   <dbl> <dbl> <dbl>
## 1 Lumber   1.08 0.649   3.5
## 2 Wacamaw  1.28 0.829   3.6

Observing a single fish with a mercury concentration of 2.5ppm is actually not that surprising! 2.5ppm is a little more than 2 standard deviations away from the mean mercury concentration in our sample of fish from the Lumber River (1.08+2*0.64=2.43), but there are certainly fish in the sample with even higher mercury concentrations (max=3.5ppm), so this isn’t outside the bounds of what we’d expect.

Exercise 1

Let’s look at the model summary output again:

summary(mod_fish)
## 
## Call:
## lm(formula = Concen ~ River, data = fish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1664 -0.5681 -0.1764  0.4219  2.4219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.07808    0.08866  12.160   <2e-16 ***
## RiverWacamaw  0.19835    0.11712   1.694   0.0922 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7575 on 169 degrees of freedom
## Multiple R-squared:  0.01669,    Adjusted R-squared:  0.01087 
## F-statistic: 2.868 on 1 and 169 DF,  p-value: 0.09218

part a: interpret model coefficient

Now, let’s interpret the RiverWacamaw coefficient. Based only on the coefficient (don’t think about the standard error yet), what can we say about the difference in mercury concentration among fish in the two rivers?

The RiverWacamaw coefficient is 0.19835, meaning that the mean mercury concentration among fish in the Wacamaw River is, on average, about 0.20ppm higher than that of fish in the Lumber River.

part b: construct a CI

Using the 68-95-99.7 rule, construct an approximate 95% confidence interval for the RiverWacamaw coefficient, and provide an appropriate interpretation.

0.20 +/- 2*0.11 –> [-0.02, 0.42]

Preferred interpretation: It is plausible that the true difference in mean mercury concentration among fish in the Wacamaw River compared to the Lumber River is between -0.02 ppm and 0.42ppm.

Not as preferred interpretation: We are 95% confident that the mean mercury concentration among fish in the Wacamaw River somewhere between 0.02ppm less than that of fish in the Lumber River and 0.42ppm more than that of fish in the Lumber River.

part c: interpreting the CI

Do you believe it plausible that the mean mercury concentration of the fish population in the Wacamaw River is approximately the same as that of the fish population in the Lumber River? How would you confirm this? What assumptions are you making?

Answers may vary–this is certainly plausible, since our 95% CI contains 0 (i.e., there is no difference in means between the two rivers). However, we might also argue that there is SOME evidence of a difference, since most of the CI is > 0.

part d: effect of sample size on our conclusions

Suppose we sample 10 times as many fish from the Wacamaw River, and get a similar coefficient estimate (0.2). Thinking back to the Central Limit Theorem, what should happen to the standard error of the RiverWacamaw coefficient? How small of a standard error would we need to more conclusively say that there is an actual difference in mean mercury concentrations of the Lumber River and Wacamaw River fish populations?

A larger sample should result in a smaller standard error of the RiverWacamaw coefficient. If the standard error is smaller than 0.1 (say 0.098), then a 95% confidence interval would be [0.004, 0.396]. Since this interval doesn’t include 0, we could conclude that fish in the Wacamaw River, on average, have a higher mercury concentration than fish in the Lumber River. More importantly, the lower standard error of the coefficient allows us to say there is evidence that this difference should be observable across new samples.

part e: reconciling parameter estimates and uncertainty

Suppose the true population coefficient for the RiverWacamawparameter is 0.02 (i.e. the average mercury concentration is 0.02ppm higher for the Wacamaw River fish population compared to that of the Lumber River). Is this meaningful?

This will depend on context–a priori, this difference appears to be negligible, and we could potentially chalk it up to uncontrolled confounders (e.g., perhaps fish in one river tend to be older/bigger and therefore have slightly higher mercury concentrations, even if there is no underlying difference in mercury pollution). We also might consider: what is considered a “harmful” mercury concentration, and are fish in either river near that threshold? Has this changed over time, and by how much?

part f: (CHALLENGE)

Using the model summary output, report the mean mercury concentration for our sample of fish from the Wacamaw River:

summary(mod_fish)
## 
## Call:
## lm(formula = Concen ~ River, data = fish)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1664 -0.5681 -0.1764  0.4219  2.4219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.07808    0.08866  12.160   <2e-16 ***
## RiverWacamaw  0.19835    0.11712   1.694   0.0922 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7575 on 169 degrees of freedom
## Multiple R-squared:  0.01669,    Adjusted R-squared:  0.01087 
## F-statistic: 2.868 on 1 and 169 DF,  p-value: 0.09218

1.07808 + 0.19835 = 1.27643ppm

Which of the following values do you think is the standard error of the sample mean for the Wacamaw River?

To answer this question, look at the code chunk below, which fits the same model, but uses the Wacamaw River as our reference category instead of the Lumber River:

mod_fish2 <- lm(Concen ~ River, data=fish %>% mutate(River=ifelse(River == "Wacamaw", paste0("_", River), River)))
summary(mod_fish2)
## 
## Call:
## lm(formula = Concen ~ River, data = fish %>% mutate(River = ifelse(River == 
##     "Wacamaw", paste0("_", River), River)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1664 -0.5681 -0.1764  0.4219  2.4219 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.27643    0.07652  16.681   <2e-16 ***
## RiverLumber -0.19835    0.11712  -1.694   0.0922 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7575 on 169 degrees of freedom
## Multiple R-squared:  0.01669,    Adjusted R-squared:  0.01087 
## F-statistic: 2.868 on 1 and 169 DF,  p-value: 0.09218

Compare this to the output for mod_fish. What do you notice about the standard errors of the intercepts (i.e., the standard errors of the means for each river) compared to the standard errors of the RiverWacamaw and RiverLumber coefficients (i.e., the standard errors of the differences between the means)?

SEs for the means are 0.08866 and 0.07652 for the Lumber and Wacamaw rivers, respectively. The SE for the difference is the same in both models (0.11712), which is greater than the SE of either mean. Because standard errors quantify uncertainty in a given parameter estimate, this tells us that the uncertainty of the estimated difference between two means is greater than the uncertainty in our estimate of either mean by itself.

Solution of Additional Exercise

Exercise 1: Simulating data under the null value

Part a

Intuition.

Part b

# 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    50  
## 2   1.5    47.2
## 3   1.7    30  
## 4   0.73   46.8
## 5   0.56   34.9
## 6   0.51   34.7

# 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")

Part c

If there were truly no relationship between Concen and Length, we’d expect the sample model to have a slope near to (but not exactly) 0.

set.seed(1)
shuffled_models <- mosaic::do(500)*(
  fish %>% 
  select(Concen, Length) %>% 
  sample_n(size = length(Length), replace = TRUE) %>% 
  mutate(Length = sample(Length, size = length(Length), replace = FALSE)) %>%
  with(lm(Concen ~ Length))
)
head(shuffled_models)
##   Intercept     Length     sigma  r.squared        F numdf dendf .row .index
## 1 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      1
## 2 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      2
## 3 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      3
## 4 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      4
## 5 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      5
## 6 0.9165014 0.01136718 0.8865891 0.01204811 2.060961     1   169    1      6

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 workaround

Part d

We’d expect slopes to be Normally distributed around 0 (the null value).

shuffled_models %>% 
  ggplot(aes(x = Length)) + 
  geom_density()

Exercise 2: Comparing our sample results to the null value (intuition)

Part a

No! Its slope is much bigger than for the sample models simulated using the null value!

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)

Part b

No! Our slope is much bigger than for the sample models simulated using the null value!

shuffled_models %>% 
  ggplot(aes(x = Length)) + 
  geom_density() + 
  geom_vline(xintercept = 0.05813, color = "blue")

Part c

Will vary.

Exercise 3: CLT

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

Exercise 4: Comparing our sample results to the null value (test statistic)

clt_plot + 
  geom_vline(xintercept = 0.05813, color = "blue")

Part a

  • 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.

Part b

Our sample slope of 0.058 falls more than 11 standard errors above the null value of 0. That’s far!!!

(0.058127 - 0) / 0.005228
## [1] 11.1184

Part c

In the t value column:

coef(summary(fish_mod_1))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'fish_mod_1' not found

Part d

  • 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.

Exercise 5: Comparing our sample results to the null value (p-value)

clt_plot + 
  geom_vline(xintercept = 0.05813, color = "blue")

Part a

Our estimate is more than 3 s.e. away from 0. Since 99.7% of estimates fall within 3 s.e. the probability of this happening is…

  • less than 0.003 (1 - 0.997)

Part b

< 2e-16 (which is very very close to 0)

coef(summary(fish_mod_1))
## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'fish_mod_1' not found

Part c

  • 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.

Exercise 6: River test Part I

Part a

\(\beta_1\), the RiverWacamaw coefficient

The null value is \(\beta_1 = 0\)

Part b

\(\hat{\beta}_1 = 0.142\) with a standard error of 0.089:

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

Part c

# 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")

Exercise 7: River test Part II

Part a

Yes! Our estimate is relatively close to 0 on this scale (within 2 s.e.).

Part b

Our estimate of \(\beta_1\), 0.142, falls only 1.59 s.e. from 0:

(0.142 - 0) / 0.089
## [1] 1.595506
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

Part c

Our estimate is somewhere between 1 and 2 s.e. from 0:

  • between 0.05 and 0.32

An 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-21

Part d

It is. It is not very far from 0 (when considering s.e.) and would not be unlikely to observe if the null value were true.

Part e

When controlling for length…

  • we do not have sufficient evidence to conclude that there’s a statistically significant difference in mercury concentration by river.

Exercise 8: River CI

Part a

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.0679477

Part b

When controlling for length…

  • we do not have sufficient evidence to conclude that there’s a statistically significant difference in mercury concentration by river.