--- title: "Hypothesis Testing- Details" subtitle: "Notes and in-class exercises" format: html: embed-resources: true toc: true --- You can download the .qmd file for this activity [here](../activity_templates/22-hypothesis-testing-details.qmd) and open in R-studio. The rendered version is posted in the [course website](https://mutasim221b.github.io/Mac-STAT-155-Fall-25/) (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. ```{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') ``` # Notes - You can download a template file to work with [here](../activity_templates/23_hypothesis_testing_details.qmd). - **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: - Apply the procedure for a formal hypothesis test - Articulate how we can formalize a research question as a testable, statistical hypothesis ## Readings and videos Please complete the following reading or videos **before** class: - Reading: Section 7.3 (stop when you get to Section 7.3.4) in the [STAT 155 Notes](https://mac-stat.github.io/Stat155Notes/) - Video 1: [Introduction to Statistical Inference](https://voicethread.com/share/15687646/) - Video 2: [Hypothesis Testing Framework](https://macalester.voicethread.com/myvoice/thread/15687700) - Video 3: [Hypothesis Testing Procedure](https://voicethread.com/share/15713570/) # Warm-up {-} **GOAL** Understand the *population* model of hiking `time` (in hours) by the highest `elevation` of the hike (in 1000s of feet): E(time | elevation) = $\beta_0$ + $\beta_1$ elevation We can make *inferences* about this relationship using data on a sample of hikes in the `peaks` data: ```{r} # Load packages & data library(tidyverse) peaks <- read.csv("https://mac-stat.github.io/data/high_peaks.csv") %>% mutate(elevation = elevation/1000) # Plot the relationship peaks %>% ggplot(aes(x = elevation, y = time)) + geom_point() # Model the relationship hike_model <- lm(time ~ elevation, data = peaks) coef(summary(hike_model)) # Get CIs confint(hike_model) ``` The *sample model* is E(time | elevation) = $\hat{\beta_0}$ + $\hat{\beta_1}$ elevation = 11.21 - 0.13 elevation Thus: - We **estimate** that, for every additional 1000 feet in elevation, the expected / average hiking time decreases by 0.13 hours (~8 minutes). - We expect that this estimate might be off by / have an **error** of 1.18 hours per 1000 feet. - We're **95% confident** that, in the *broader population of hikes*, an additional 1000ft in elevation is associated with anywhere from a 2.50 hour *decrease* to a 2.24 hour *increase* in expected / average hiking time. Thus we do NOT have *statistically significant* evidence that hiking time is associated with elevation (a change of 0 is in this interval, hence is a plausible value). ## Example 1: Starting a formal hypothesis test $H_0$: hiking time is *not* associated with elevation $H_a$: hiking time *is* associated with elevation Translate this to notation: $H_0$: $\beta_1 = 0$ $H_a$: $\beta_1 \ne 0$ NOTE: We will evaluate these hypotheses using a 0.05 **significance level**. It's important to choose our significance level *before* starting our test to avoid post-hoc "tweaks" that suit our narrative. ## Example 2: What would we expect under $H_0$? Suppose $H_0$ were true, i.e. $\beta_1$ were 0. By the Central Limit Theorem, we'd expect the sampling distribution of possible $\hat{\beta}_1$ sample estimates to be Normally distributed around 0 with the standard error obtained from our model summary table: $$ \hat{\beta}_1 \sim N(\beta_1, s.e.(\hat{\beta}_1)^2) \;\;\; \Rightarrow \;\;\; \hat{\beta}_1 \sim N(0, 1.18^2) $$ Based on the corresponding plot below, is *our* sample estimate of $\hat{\beta}_1 = -0.13$ consistent with $H_0$? ```{r fig.width = 8} # IGNORE THE SYNTAX!!! # Put OUR sample estimate here est <- -0.13 # Put the corresponding standard error here se <- 1.18 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 H0 were true") ``` > Yes! It's less than 1 standard error away from 0. ## Example 3: Test statistic ```{r} coef(summary(hike_model)) ``` a. Calculate the test statistic for our hypothesis test. ```{r eval = TRUE} # By hand (-0.1269 - 0) / 1.1755 # This is also reported in the t value column coef(summary(hike_model)) ``` b. Interpret the test statistic: Our sample estimate of $\beta_1$ is ___ standard errors ___ the null value of ___. > Our sample estimate of $\beta_1$ is 0.108 standard errors below the null value of 0. c. So is our sample data consistent with $H_0$? > Yes -- it's "close to" 0 (0.108 s.e. is small). ## Example 4: p-value a. Use the 68-95-99.7 Rule with the sketch from Example 2 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 > bigger than 0.32 Our estimate is less than 1 s.e. from 0. b. Report the *exact* p-value from the model `summary()`. ```{r} coef(summary(hike_model)) ``` > 0.915 (as reported in the `Pr(>|t|)` column) c. How can we interpret this p-value? Choose all that apply! - a. Given our sample data, it's likely that hiking time is associated with elevation (i.e. that $H_a$ is true). $$P(H_a | data) = \text{prob of $H_a$ being true given our data}$$ - b. Given our sample data, it's likely that hiking time is not associated with elevation (i.e. that $H_0$ is true). $$P(H_0 | data) = \text{prob of $H_0$ being true given our data}$$ - c. If in fact there were no relationship between hiking time and elevation in the broader population of hikes (i.e. $H_0$ were true), it's likely that we'd have gotten our observed decrease of hiking time with elevation "by chance". $$P(data | H_0) = \text{prob of observing our data given $H_0$ is true}$$ - d. 0.915 is the probability of observing a test statistic as extreme or more extreme than the once we obtained (|-0.107|), assuming $H_0$ is true > c and d d. So is our sample data consistent with $H_0$? > yes -- there's a high chance we would've gotten a sample slope this far from 0 if $H_0$ were true. In other words, we are likely to see a test statistic more extreme than the one we obtained under $H_0$ is true. ## Example 5: Conclusion Throughout the test, we've been comparing *our* data to the null hypothesis $H_0$. Thus at the conclusion of a hypothesis test, either... - Our results are "statistically significant". - We have enough evidence to *reject* $H_0$ in favor of $H_a$ - NOTE: this does *not* mean that we "accept" $H_a$, we just "have evidence" for it. - Our results are *not* "statistically significant". - We do *not* have enough evidence to reject $H_0$, i.e. we *fail to reject* $H_0$. - NOTE: This does *not* mean that we "accept" $H_0$, we just "don't have enough evidence" to to reject it. a. So, at the 0.05 significance level, what's our conclusion? Is there a **statistically significant** association between hiking time and elevation? > p-value > 0.05. Thus we fail to reject $H_0$. We do not have statistically significant evidence of an association between hiking time and elevation. b. Does this conclusion agree with the one we made using the confidence interval for $\beta_1$, (-2.50, 2.24)? > yes c. Does this conclusion agree with what we'd conclude from the confidence bands below? ```{r} peaks %>% ggplot(aes(x = elevation, y = time)) + geom_point() + geom_smooth(method = "lm") ``` > Yes, we can draw an example of a 0-slope line that falls within the bands. \ \ \ \ ::: {.callout-note title = "CIs, test statistics, and p-values"} We've observed 3 equivalent approaches to determining whether or not to reject the null hypothesis that $\beta$ equals some null value. Assuming a 95% confidence level: - If the **95%** CI for $\beta$ does not include the null value, reject $H_0$! - If the p-value < **0.05**, reject $H_0$! - If the |test statistic| > **2**, reject $H_0$! ::: \ \ \ \ ::: {.callout-note title = "Hypothesis t-Tests for Model coefficients"} Consider some population model $$ E[Y | X_1, X_2, X_3] = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 $$ In our model `summary()` table, the reported test statistic (`t value`) and p-value (`Pr(>|t|)`) correspond to the following "**t-test**": $H_0$: $\beta = 0$ $H_a$: $\beta \ne 0$ The *meaning* or *intepretation* of these hypothesis depends upon the meaning of $\beta$, which itself depends upon: - whether the corresponding $X$ predictor is quantitative or categorical - what other predictors are in the model ::: # Exercises ## Exercise 1 **Research Question**: Can we predict whether or not a mushroom is poisonous based on the shape of its cap? For this exercise, we will look at data from various species of gilled mushrooms in the Agaricus and Lepiota Family. We have information on whether a mushroom is `poisonous` (TRUE if it is, FALSE if it's edible), the shape of its cap (`cap_shape`, a categorical variable with 6 categories), the texture of its cap surface (`cap_surface`, a categorical variable with 4 categories), and the size of its gills (`gill_size`, a categorical variable with two categories) ```{r} # Load the data & packages library(tidyverse) mushrooms <- read_csv("https://Mac-STAT.github.io/data/mushrooms.csv") mushrooms <- mushrooms %>% mutate(cap_shape = relevel(as.factor(cap_shape), ref="flat")) %>% dplyr::select(poisonous, cap_shape) head(mushrooms) ``` ### Part a One of the most poisonous species of mushrooms is the *Amanita phalloides* or ["Death Cap" mushroom](http://www.bccdc.ca/health-info/prevention-public-health/death-cap-mushrooms), which typically has a flat cap shape when mature. Based on this anecdote, we hypothesize that species of mushrooms with flat caps in general *associated* with the likelihood of being poisonous. First, let's translate this question to an appropriate null and alternative hypothesis that we can compare with a formal hypothesis test. Remember that `poisonous` is a binary outcome, so we need to frame our null and alternative hypotheses in terms of odds (i.e., Odds(poisionous | flat cap) = P(poisonous|flat cap)/P(edible | flat cap)). > $H_0$: Odds(poisonous | flat cap) = 1 > > $H_a$: Odds(poisonous | flat cap) ≠ 1 ### Part b - Fit a logistic regression model to investigate whether `cap_shape` is associated with a mushroom being `poisonous`. (Note that in the setup code above, we have forced the reference category for the `cap_shape` predictor to be `flat`; without this, the reference category by default would be set as `bell`, which is the first category when sorted alphabetically). ```{r eval = TRUE} mushroom_mod1 <- glm(poisonous ~ cap_shape, data=mushrooms, family="binomial") coef(mushroom_mod1) ``` ### Part c Provide an appropriate interpretation of the intercept coefficient **on the odds scale**. Based on this interpretation, do you believe mushrooms with flat caps are more likely to be poisonous, or more likely to be edible? ```{r eval = TRUE} exp(coef(mushroom_mod1)) ``` > The odds of a flat-capped mushroom being poisonous are 0.975:1--that is, mushrooms with flat caps are very slightly *less* likely to be poisonous than they are edible. ### Part d Let's look at the full model summary: ```{r} summary(mushroom_mod1) ``` Report and interpret the test statistic for the intercept term (our coefficient of interest): > > The test statistic is -0.712—this means that the coefficient estimate of interest is 0.712 standard errors away from (specifically, below) the null value of 0 (note that this is on the log-odds scale). ### Part e - Report and interpret the p-value for the intercept term. - Based on this p-value and a significance level of 0.05, do we have evidence that mushrooms with flat caps are more likely to be poisonous than edible? > The p-value for the intercept term is 0.4762. > Interpretation: If the null hypothesis were true (i.e., the odds of being poisonous were 1), the probability of seeing a test statistic as or more extreme than |-0.712| is 0.4762. Because the p-value is greater than our significance level of 0.05, we have no evidence to suggest that a flat-capped mushroom is more or less likely to be poisonous. ### Part f Now suppose we are interested in whether the odds of being poisonous are different for mushrooms with other cap shapes. By hand, calculate the odds of being poisonous for mushrooms with knobbed caps, conical caps, and sunken caps (remember that the non-exponentiated coefficients represent a difference in log-odds compared to the reference category): > odds(poisonous | knobbed cap) = > odds(poisonous | conical cap) = > odds(poisonous | sunken cap) = ```{r} # knobbed #conical #sunken ``` ### Part g Based on these odds, which of the 4 mushroom cap shapes we've investigated (flat, knobbed, conical, and sunken) do you believe is the best indicator that it's edible? Which cap shape do you expect is most likely to be poisonous? > **Your answer** ### Part h Let's get the full model summary again: ```{r} summary(mushroom_mod1) ``` Now report and interpret the p-values for the coefficients corresponding to `cap_shapeknobbed`, `cap_shapeconical`, and `cap_shapesunken`: > cap_shapeknobbed: Our null hypothesis is that the odds ratio between flat-capped and knob-capped mushrooms is 1 (i.e., the odds of a knob-capped mushroom being poisonous are equal to the odds of a flat-capped mushroom being poisonous). If we assume the null hypothesis is true, then the probability of seeing a test statistic as or more extreme than |11.60| is (<2e-16). Because the p-value is far below our significance level of 0.05, we take this as strong evidence that knob-capped mushrooms are much more likely to be poisonous than flat-capped mushrooms. > ** Your Response on cap_shapeconical ** > ** Your Response on cap_shapesunken ** **Notes on (h)** For interpreting p-value as mentioned in (part h), check if the odds ratio (OR) differ from 1 at all. > $H_0$: OR = 1 > > $H_a$: OR ≠ 1 If it does, then look at the sign of the coefficient to decide whether OR>1 or OR<1. If positive coefficient then OR>1 and hence use 'more likely' in p-value interpretation. Whereas, if negative coefficient then OR<1 and hence use 'less likely' in p-value interpretation. ### Part i Based on the model summary output in **part h** above, if you were given a plate of mushrooms with different cap shapes and had to pick one to eat, which one would you choose? Which cap shape would you absolutely avoid at all costs? Are your decisions guided by the coefficient estimates, the p-values, or both? > **Your answer** ### Part j Let's look at the data a slightly different way, using a 6x2 table of counts: ```{r} mushrooms %>% mutate(cap_shape=as.factor(cap_shape), poisonous=as.factor(poisonous)) %>% dplyr::count(cap_shape, poisonous, .drop=FALSE) %>% pivot_wider(names_from=poisonous, values_from=n, names_prefix="Poisonous = ") ``` Now, if you were given a plate of mushrooms with different cap shapes and had to pick one shape to eat and one to absolutely avoid, would you choose the same shapes? Why or why not? > **Your answer** # Additional Problem ## Exercise 2 For this exercise, let's return to the fish dataset from a previous activity. ```{r} fish <- read_csv("https://Mac-STAT.github.io/data/Mercury.csv") head(fish) ``` **Research question:** We believe the length of a fish (measured in centimeters) is causally associated with its mercury concentration (measured in parts per million [ppm]). We suspect that the river a fish is sampled from may be a confounder, since differences in the river environment may causally influence both the average length of fish (e.g. due to differences in water temperature or food availability) as well as mercury concentration (e.g. due to differences between the two rivers in mercury pollution levels). ### Part a Fit a linear regression model that can be used to answer our research question. ```{r} mod_fish1 <- lm(Concen ~ Length + River, data=fish) summary(mod_fish1) ``` ### Part b Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw` coefficient. Assume we have specified a significance level of 0.05. > **Response** ### Part c Suppose we now want to determine if the causal effect of fish length on mercury concentration *differs* according to the river a fish was sampled from. First, modify the code chunk below to visualize the 3-way relationship between the `Concen`, `Length`, and `River` variables. ```{r} fish %>% ggplot(aes(x = Length, y = Concen, colour = River)) + geom_point()+ geom_smooth(method="lm", se=F) ``` Next, fit an appropriate linear regression model with an interaction term to investigate this question. ```{r} mod_fish2 <- lm(Concen ~ Length * River, data=fish) summary(mod_fish2) ``` ### Part d Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw:Length` interaction term in this revised model (`mod_fish2`). Assume we've set a significance level of 0.05. > **Response** ### Part e Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw` coefficient in this revised model (`mod_fish2`). (again, you can assume we've set a significance level of 0.05). ### Part f (CHALLENGE) Suppose another researcher runs the same model we fit in part c above (`mod_fish2`), but they claim that a more appropriate alternative hypothesis should be Beta_1 < 0, (and not Beta_1 ≠ 0, as is assumed by default when running a regression model). Because of this, they reported a smaller p-value for the coefficient, and claim that the Wacamaw River has a lower baseline mercury concentration (i.e., when `Length = 0cm`). What is the p-value they would have reported for the `RiverWacamaw` coefficient in `mod_fish2`? > **Response** What is a potential ethical problem with the other researcher's claim that the alternative hypothesis should be Beta_1 < 0? > **Response** ### Part g (CHALLENGE) You point out to the other researcher that the intercept and `RiverWacamaw` coefficients are both negative, so whatever difference in mercury concentration between the two rivers your model predicts "at baseline" is not useful or meaningful--you cannot have a fish that is 0cm long, nor a mercury concentration <0ppm. You propose that a more appropriate model should transform the `Length` variable in some way to make the intercept more interpretable. Create a new variable named `Length_adj` with this transformation and use it to re-fit the model: ```{r} mod_fish3 <- lm(Concen ~ Length_adj*River, data=fish) summary(mod_fish3) ``` Compare the output of this model to that of `mod_fish2`. What happened to the estimate, test statistic, and p-value for the `RiverWacamaw` coefficient? How does this affect your conclusion? How about the other researcher's conclusion? > **Response** \ \ \ \ # Solutions ```{r eval = TRUE, echo = FALSE} # Load the data & packages library(tidyverse) mushrooms <- read_csv("https://Mac-STAT.github.io/data/mushrooms.csv") mushrooms <- mushrooms %>% mutate(cap_shape = relevel(as.factor(cap_shape), ref="flat")) %>% dplyr::select(poisonous, cap_shape) fish <- read_csv("https://Mac-STAT.github.io/data/Mercury.csv") ``` ## Exercise 1 ### Part a One of the most poisonous species of mushrooms is the *Amanita phalloides* or ["Death Cap" mushroom](http://www.bccdc.ca/health-info/prevention-public-health/death-cap-mushrooms), which typically has a flat cap shape when mature. Based on this anecdote, we hypothesize that species of mushrooms with flat caps in general may be *more* likely to be poisonous than edible. First, let's translate this question to an appropriate null and alternative hypothesis that we can compare with a formal hypothesis test. Remember that `poisonous` is a binary outcome, so we need to frame our null and alternative hypotheses in terms of odds (i.e., Odds(poisonous | flat cap) = P(poisonous|flat cap)/P(edible | flat cap)). > $H_0$: Odds(poisonous | flat cap) = 1 > > $H_a$: Odds(poisonous | flat cap) ≠ 1 ### Part b - Fit a logistic regression model to investigate whether `cap_shape` is associated with a mushroom being `poisonous`. (Note that in the setup code chunk above, we have forced the reference category for the `cap_shape` predictor to be `flat`; otherwise, the reference category by default would be set as `bell`, which is the first category when sorted alphabetically). ```{r eval = TRUE} mushroom_mod1 <- glm(poisonous ~ cap_shape, data=mushrooms, family="binomial") coef(mushroom_mod1) ``` ### Part c Provide an appropriate interpretation of the intercept coefficient **on the odds scale**. Based on this interpretation, do you believe mushrooms with flat caps are more likely to be poisonous, or more likely to be edible? ```{r eval = TRUE} exp(coef(mushroom_mod1)) ``` > The odds of a flat-capped mushroom being poisonous are 0.975:1--that is, mushrooms with flat caps are very slightly *less* likely to be poisonous than they are edible. ### Part d Let's look at the full model summary: ```{r eval = TRUE} summary(mushroom_mod1) ``` Report and interpret the test statistic for the intercept term (our coefficient of interest): > The test statistic is -0.712—this means that the coefficient estimate of interest is 0.712 standard errors away from (specifically, below) the null value of 0 (note that this is on the log-odds scale). ### Part e - Report and interpret the p-value for the intercept term. - Based on this p-value and a significance level of 0.05, do we have evidence that mushrooms with flat caps are more likely to be poisonous than edible? > The p-value for the intercept term is 0.4762. > Interpretation: If the null hypothesis were true (i.e., the odds of being poisonous were 1), the probability of seeing a test statistic as or more extreme than |-0.712| is 0.4762. Because the p-value is greater than our significance level of 0.05, we have no evidence to suggest that a flat-capped mushroom is more or less likely to be poisonous. ### Part f Now suppose we are interested in whether the odds of being poisonous are different for mushrooms with other cap shapes. Calculate the odds of being poisonous for mushrooms with knobbed caps, conical caps, and sunken caps (remember that the non-exponentiated coefficients represent a difference in log-odds compared to the reference category): ```{r eval = TRUE} # knobbed exp(-0.025+0.992) #conical exp(-0.025+14.59) #sunken exp(-0.025-14.54) ``` ### Part g Based on these odds, which of the 4 mushroom cap shapes we've investigated (flat, knobbed, conical, and sunken) do you believe is the best indicator that it's edible? Which cap shape do you expect is most likely to be poisonous? > Using only the coefficient estimates, it appears that mushrooms with a sunken cap shape appear to be most likely to be edible, as the odds they are poisonous are approximately $4.7 \times 10^{-7}$ to 1. Mushrooms with conical caps appear to be most likely to be poisonous (odds of being poisonous are >2 million to 1). ### Part h Let's get the full model summary again: ```{r eval = TRUE} summary(mushroom_mod1) ``` Now report and interpret the p-values for the coefficients corresponding to `cap_shapeknobbed`, `cap_shapeconical`, and `cap_shapesunken`: > `cap_shapeknobbed`: Our null hypothesis is that the odds ratio between flat-capped and knob-capped mushrooms is 1 (i.e., the odds of a knob-capped mushroom being poisonous are equal to the odds of a flat-capped mushroom being poisonous). If we assume the null hypothesis is true, then the probability of seeing a test statistic as or more extreme than |11.60| is 3.91e-31. Because the p-value is far below our significance level of 0.05, we take this as strong evidence that knob-capped mushrooms are much more likely to be poisonous than flat-capped mushrooms. > `cap_shapeconical`: Our null hypothesis is that the odds ratio between flat-capped and cone-capped mushrooms is 1 (i.e., the odds of a cone-capped mushroom being poisonous are equal to the odds of a flat-capped mushroom being poisonous). If we assume the null hypothesis is true, then the probability of seeing a test statistic as or more extreme than |0.03| is 0.974 (i.e., we are very likely to see a test statistic more extreme than |-0.03| if in fact there were no difference between flat-capped and cone-cap mushrooms in odds of being poisonous). Because the p-value is far above our significance level of 0.05, we do not have evidence that the odds of a cone-capped mushroom being poisonous differ from the odds of a flat-capped mushroom being poisonous. > `cap_shapesunken`: Our null hypothesis is that the odds ratio between flat-capped and sunken-cap mushrooms is 1 (i.e., the odds of a sunken-cap mushroom being poisonous are equal to the odds of a flat-capped mushroom being poisonous). If we assume the null hypothesis is true, then the probability of seeing a test statistic as or more extreme than |-0.09| is 0.926 (i.e., we are very likely to see a test statistic more extreme than |-0.09| if in fact there were no difference between flat-capped and sunken-cap mushrooms in odds of being poisonous). Because the p-value is far above our significance level of 0.05, we do not have evidence that the odds of a sunken-capped mushroom being poisonous differ from the odds of a flat-capped mushroom being poisonous. ### Part i Based on the model summary output in **part h** above, if you were given a plate of mushrooms with different cap shapes and had to pick one to eat, which one would you choose? Which cap shape would you absolutely avoid at all costs? Are your decisions guided by the coefficient estimates, the p-values, or both? > Answers may vary--if only considering coefficient estimates, then cone-shaped caps are most likely to be poisonous and sunken-shaped caps are most likely to be edible. But if we only look at p-values, then knob-shaped caps have the strongest evidence that they are more likely to be poisonous, and bell-shaped caps have the strongest evidence that they are more likely to be edible. ### Part j Let's look at the data a slightly different way, using a 6x2 table of counts: ```{r eval = TRUE} mushrooms %>% mutate(cap_shape=as.factor(cap_shape), poisonous=as.factor(poisonous)) %>% dplyr::count(cap_shape, poisonous, .drop=FALSE) %>% pivot_wider(names_from=poisonous, values_from=n, names_prefix="Poisonous = ") ``` Now, if you were given a plate of mushrooms with different cap shapes and had to pick one shape to eat and one to absolutely avoid, would you choose the same shapes? Why or why not? > Personally, I'd stick with the sunken-shaped caps for eating. Even though our model suggests there's no evidence to believe they are less likely to be poisonous, 0 out of 32 in the sample are poisonous, which seems like the least risky choice. However, I'd tend to avoid the knob-capped mushrooms more than the cone-capped mushrooms—even though the latter are all poisonous in the sample, there were only 4 observations, so it's possible that due to sampling variation, the odds of being poisonous for cone-capped mushrooms is lower than that of knob-capped mushrooms (where we have many more observations). # Additional Probelem ## Exercise 2 For this exercise, let's return to the fish dataset from the previous activity (Activity 22). ```{r eval = TRUE} fish <- read_csv("https://Mac-STAT.github.io/data/Mercury.csv") head(fish) ``` **Research question:** We believe the length of a fish (measured in centimeters) is causally associated with its mercury concentration (measured in parts per million [ppm]). We suspect that the river a fish is sampled from may be a confounder, since differences in the river environment may causally influence both the average length of fish (e.g. due to differences in water temperature or food availability) as well as mercury concentration (e.g. due to differences between the two rivers in mercury pollution levels). ### Part a Fit a linear regression model that can be used to answer our research question. ```{r eval = TRUE} mod_fish1 <- lm(Concen ~ Length + River, data=fish) summary(mod_fish1) ``` ### Part b Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw` coefficient. Assume we have specified a significance level of 0.05. > coefficient: Holding fish length constant, we estimate the average mercury concentration among fish in the Wacamaw River to be 0.14ppm higher than fish in the Lumber River. > Test statistic: The estimate we observe (0.14) is 1.587 standard errors higher than the null value of a 0ppm difference in mercury concentration. > p-value: Assuming the null hypothesis is true and there is no actual difference in mercury concentration among the two fish populations (adjusting for fish length), the probability of observing a test statistic as or more extreme than |1.587| is 0.114. Because 0.114 > 0.05, we do not have sufficient evidence to reject the null hypothesis, and conclude that the average mercury concentration does not differ between the two rivers. ### Part c Suppose we now want to determine if the causal effect of fish length on mercury concentration *differs* according to the river a fish was sampled from. First, modify the code chunk below to visualize the 3-way relationship between the `Concen`, `Length`, and `River` variables. ```{r eval = TRUE} fish %>% ggplot(aes(x = Length, y = Concen, colour = River)) + geom_point()+ geom_smooth(method="lm", se=F) ``` Next, fit an appropriate linear regression model with an interaction term to investigate this question. ```{r eval = TRUE} mod_fish2 <- lm(Concen ~ Length * River, data=fish) summary(mod_fish2) ``` ### Part d Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw:Length` interaction term in this revised model (`mod_fish2`). Assume we've set a significance level of 0.05. > coefficient: First, we interpret the `Length` coefficient--that is, among fish in the Lumber river, we expect that a 1cm increase in length is associated with a 0.043ppm increase in mercury concentration. The interaction coefficient tells us the expected change in that relationship when considering fish in the Wacamaw River instead: we expect an *additional* 0.024ppm increase in mercury concentration associated with a 1cm increase in length (i.e., in the Wacamaw River, we expect mercury concentration to increase by 0.067ppm per 1cm increase in fish length). > Test statistic: The estimate we observe (0.024326) is 2.321 standard errors higher than the null value of zero. > p-value: Assuming the null hypothesis is true and there is no difference in the relationship between fish length and mercury between the 2 rivers, the probability of observing a test statistic as or more extreme than |2.321| is 0.02. Because 0.02 < 0.05, we take this as evidence to reject the null hypothesis, and conclude that the effect of fish length on mercury concentration does differ slightly between the two rivers. ### Part e Interpret the coefficient estimate, test statistic, and p-value for the `RiverWacamaw` coefficient in this revised model (`mod_fish2`). (again, you can assume we've set a significance level of 0.05). > coefficient: Visully, the `RiverWacamaw` coefficient represents the difference in the y-intercepts for the best fit lines for the Lumber and Wacamaw Rivers in the part d plot. Interpretation: Among fish that are 0cm long, average fish mercury concentrations are 0.826ppm lower in the Wacamaw River than in the Lumber River. > Test statistic: The estimate we observe (-0.826291) is 1.937 standard errors lower than the null value of 0. > p-value: Assuming the null hypothesis is true, the probability of observing a test statistic as or more extreme than |-1.937| is 0.0544. Because 0.0544 > 0.05, we do not have evidence to reject the null. ### Part f (CHALLENGE) Suppose another researcher runs the same model we fit in part c above (`mod_fish2`), but they claim that a more appropriate alternative hypothesis should be Beta_1 < 0, (and not Beta_1 ≠ 0, as is assumed by default when running a regression model). Because of this, they reported a smaller p-value for the coefficient, and claim that the Wacamaw River has a lower baseline mercury concentration (i.e., when `Length = 0cm`). What is the p-value they would have reported for the `RiverWacamaw` coefficient in `mod_fish2`? > 0.0544/2 = 0.0272 (we divide the "two-tailed" p-value in half to obtain the p-value for a "one-tailed" test) What is a potential ethical problem with the other researcher's claim that the alternative hypothesis should be Beta_1 < 0? > It is possible that the researchers had a particular reason or incentive to publish evidence in support of their hypothesis (some potential reasons are that scientific journals are generally less interested in publishing results, a financial conflict of interests, or favoring a pet hypothesis). They could have first looked at the results of a "two-tailed" statistical test and since the p-value is very close to the traditional significance threshold of 0.05, come up with a post-hoc rationalization to perform a hypothesis test resulting in a "statistically significant" p-value. This unethical practice is known in the field as "p-hacking." ### Part g (CHALLENGE) You point out to the other researcher that the intercept and `RiverWacamaw` coefficients are both negative, so whatever difference in mercury concentration between the two rivers your model predicts "at baseline" is not useful or meaningful--you cannot have a fish that is 0cm long, nor a mercury concentration <0ppm. You propose that a more appropriate model should transform the `Length` variable in some way to make the intercept more interpretable. Create a new variable named `Length_adj` with this transformation and use it to re-fit the model: ```{r eval = TRUE} fish <- fish %>% mutate(Length_adj=Length-min(Length)) mod_fish3 <- lm(Concen ~ Length_adj*River, data=fish) summary(mod_fish3) ``` Compare the output of this model to that of `mod_fish2`. What happened to the estimate, standard error, test statistic, and p-value for the `RiverWacamaw` coefficient? How does this affect your conclusion? How about the other researcher's conclusion? > The `RiverWacamaw` coefficient increased (and became closer to 0). The standard error decreased, but the test statistic decreased in magnitude and the p-value increased. > What happened with this transformation is that the vertical axis got shifted so that the new "zero" was at 25.2cm (the minimum fish length in the data). At this point there is a smaller difference between the Lumber and Wacamaw River lines. However, as we saw above, there does seem to be a true modest difference in the slopes of these lines so there are larger differences between the 2 rivers at larger fish lengths.