26  Nested Models & F-Tests

Settling In





Recap

ImportantStatistical superpowers

We can do more than 1 thing at a time! Specifically, we can test hypotheses about multiple population parameters in a single test.



ImportantLearning goals

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

  • Determine if one model is nested within another
  • Determine which null and alternative hypotheses require an f-test
  • Determine which f-tests require the use of the anova function in R vs. the overall f-test given in regular regression output
  • Interpret the results of an f-test in context



ImportantAdditional resources

Required video

Optional



Hypothesis Testing for Model Selection

Nested Models

“Nested models” are models where the covariates of one model are a subset of another model. In other words, you can get to the smaller, nested model by setting some of the \beta’s to 0.



As an example, consider the following models for estimating the association between forced expiratory volume (FEV) and smoking status:

Model 1:

E[FEV \mid smoke] = \beta_0 + \beta_1 smoke

Model 2:

E[FEV \mid smoke, age] = \beta_0 + \beta_1 smoke + \beta_2 age

Here, Model 1 is “nested” inside Model 2, since the covariates included in Model 1 (only smoke) are a subset of those in Model 2 (both smoke and age).





An example of non-nested models are…

Model 3:

E[FEV \mid smoke, height] = \beta_0 + \beta_1 smoke + \beta_2 height

Model 4:

E[FEV \mid smoke, sex] = \beta_0 + \beta_1 smoke + \beta_2 sex

Here, even though Model 3 and Model 4 both contain smoke as explanatory variables, neither is nested in the other, since sex is not a part of Model 3, and height is not a part of Model 4.




F-test

NoteF-test test statistic

The F-test test statistic is not calculated like the t-test test statistic, thus cannot be interpreted the same way! Check out the theory box below to learn how this statistic is calculated.

NoteTheory Box: Calculating the F-test test statistic

Consider 2 models which we’ll estimate using n data points:

  • Model 1 (Big model)
    • E[Y|X_1,...,X_k,...,X_p] = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k + \cdots + \beta_p X_p
    • “Residual Sum of Squares”: RSS_1 = \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 where \hat{Y}_i are the predictions calculated from Model 1
  • Model 2 (Small model)
    • E[Y|X_1,...,X_k] = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k
    • “Residual Sum of Squares”: RSS_2 = \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 where \hat{Y}_i are the predictions calculated from Model 2

And consider the F-test:

H_0: \beta_{k+1} = \beta_{k+2} = ... = \beta_p = 0

H_a: at least one of \beta_{k+1}, \beta_{k+2}, …, \beta_p is non-0

Then the F-test test statistic measures the additional residual error when using Model 2 vs Model 1, while accounting for the complexity of each model. The greater the test statistic, the greater the residuals when we simplify the model (go from Model 1 to Model 2), thus the more evidence we have to reject H_0:

F = \frac{RSS_2 - RSS_1}{RSS_1} \frac{n - (p + 1)}{p - k}



Hypothesis Test Overview

Name of Test Hypotheses Relevant Code/Output

Test for individual coefficients

  • t-test (linear)

  • z-test (logistic)

H_0: \beta_j = 0

H_a: \beta_j \ne 0

coef(summary())

Look for:

lm(): t value, Pr(>|t|)

glm(): z value, Pr(>|z|)

Nested Test (group of coef)

  • F-test (linear)

  • Chi Square-test (logistic)

H_0: some set of \beta_j’s = 0

H_a: \text{at least one of those } \beta_j \ne 0

anova(small,large)

Look for:

lm(): Pr(>F)

glm(): Pr(>|Chi|)

Overall Nested Test

  • F-test (linear)

H_0: \beta_1=\beta_2=\cdots=\beta_k = 0

H_a: \text{at least one } \beta_j \ne 0

summary()

Look for:

lm(): F-statistic: X on n and m DF, p-value:





Warm-up

Example 1: Nested Models

  1. Which of the following models are nested in the model E[A \mid B, C, D] = \beta_0 + \beta_1 D + \beta_2 B + \beta_3 C + \beta_4 B * C?
  • Model 1: E[A \mid B] = \beta_0 + \beta_1 B
  • Model 2: E[A \mid B, D] = \beta_0 + \beta_1 B + \beta_2 D
  • Model 3: E[B \mid C] = \beta_0 + \beta_1 C
  • Model 4: E[A \mid B, C, D] = \beta_0 + \beta_1 B + \beta_2 C + \beta_3 D
  • Model 5: E[A \mid B, C, D] = \beta_0 + \beta_1 C + \beta_2 B + \beta_3 D + \beta_4 B * D
  • Model 6: E[A \mid D] = \beta_0 + \beta_1 D



  1. Consider the following models involving variables A, B, C, and D:
  • Model 1: E[A \mid B] = \beta_0 + \beta_1 B
  • Model 2: E[A \mid B, C] = \beta_0 + \beta_1 B + \beta_2 C
  • Model 3: E[A \mid B, C] = \beta_0 + \beta_1 B + \beta_2 C + \beta_3 BC
  • Model 4: E[A \mid C, D] = \beta_0 + \beta_1 C + \beta_2 D
  • Model 5: E[B \mid A] = \beta_0 + \beta_1 A
  • Model 6: E[B \mid A, C] = \beta_0 + \beta_1 A + \beta_2 C + \beta_3 AC

Determine for each of the following statements whether that statement is True or False.

  • Model 1 is nested in Model 2
  • Model 1 is nested in Model 3
  • Model 1 is nested in Model 4
  • Model 2 is nested in Model 3
  • Model 3 is nested in Model 2
  • Model 2 is nested in Model 6
  1. What is one (numeric) way to compare the quality of nested models? Explain how you would determine which model is “better” based on this metric.

Example 2: t-tests & “overall” F-tests

To explore the relationship of casual bikeshare ridership with the day of week (weekend or weekday), feels like temperature (F), and actual temperature (F), let’s explore the following population model:

E[rides | weekend, temp_feel, temp_actual] = \beta_0 + \beta_1 weekendTRUE + \beta_2 temp_feel + \beta_3 temp_actual

We estimate this below using bike_model_1:

bikes <- read_csv("https://mac-stat.github.io/data/bikeshare.csv") %>% 
  rename(rides = riders_casual)

bike_model_1 <- lm(rides ~ weekend + temp_feel + temp_actual, bikes)
summary(bike_model_1)

Call:
lm(formula = rides ~ weekend + temp_feel + temp_actual, data = bikes)

Residuals:
     Min       1Q   Median       3Q      Max 
-1592.49  -252.52   -17.41   181.65  1973.12 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1370.647     86.438 -15.857   <2e-16 ***
weekendTRUE   813.556     36.312  22.405   <2e-16 ***
temp_feel      11.997      8.712   1.377   0.1689    
temp_actual    15.884      9.459   1.679   0.0935 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 443.8 on 727 degrees of freedom
Multiple R-squared:  0.584, Adjusted R-squared:  0.5823 
F-statistic: 340.2 on 3 and 727 DF,  p-value: < 2.2e-16
  1. The “overall” F-test is reported in the bottom line of the summary(). What do you conclude from this test? NOTE: The F-test test statistic is not calculated in the same way as the t-test test statistic (it does not give the number of SE’s away from 0).

  2. What do you conclude from the t-tests in the temp_feel and temp_actual rows of the summary() table?

  3. Putting these together, what do you think? Which variables might be “good” predictors of ridership?

Example 3: F-tests for nested models

Let’s take the temperature variables out of the model:

# Refit bike_model_1 (just for ease of comparison)
bike_model_1 <- lm(rides ~ weekend + temp_feel + temp_actual, bikes)

# Fit a new model
bike_model_2 <- lm(rides ~ weekend, bikes)
  1. Note that bike_model_2 is nested in bike_model_1. Thus we can compare them with the following F-test for nested models. State the hypotheses, p-value, and your conclusion.
# Put the smaller model first!!!
anova(bike_model_2, bike_model_1)
Analysis of Variance Table

Model 1: rides ~ weekend
Model 2: rides ~ weekend + temp_feel + temp_actual
  Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
1    729 253858215                                  
2    727 143177975  2 110680240 280.99 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. A series of models and tests can provide more insight than one model or test alone! What did we learn about the relationship of rides with temp_feel and temp_actual from the above examples combined? Why do you think this happened? What would you do next?

Exercises

DIRECTIONS

  • Throughout this activity, test hypotheses at the 0.05 significance level.
  • Make all conclusions and interpretations in context.

CONTEXT

The MacGrades.csv dataset contains a sub-sample (to help preserve anonymity) of every grade assigned to a former Macalester graduating class. For each of the 6414 rows of data, the following information is provided (with a few missing values):

  • sessionID: A section ID number
  • sid: A student ID number
  • grade: The grade obtained, as a numerical value (i.e. an A is a 4, an A- is a 3.67, etc.)
  • dept: A department identifier (these have been made ambiguous to maintain anonymity)
  • level: The course level (e.g. 100-, 200-, 300-, and 600-)
  • sem: A semester identifier
  • enroll: The section enrollment
  • iid: An instructor identifier (these have been made ambiguous to maintain anonymity)
# Load packages & data
library(tidyverse)
MacGrades <- read_csv("https://mac-stat.github.io/data/MacGrades.csv")%>% 
  mutate(level = factor(level)) # make level a factor variable
head(MacGrades)

Exercise 1: Explore

NOTE: This exercise, since it’s exploratory in nature, can suck up a lot of time if you let them! For the sake of getting to the rest of the activity, please spend no more than ~5 minutes on this.

  1. Hypothesize two relationships between the variables in the dataset (pick any two relationships you want!). Your response should be written in a paragraph form.

Response Put your response here

  1. Explore the relationship between course grades and other variables in the data. Make two visualizations, and describe any patterns you observe.

Exercise 2: F-tests for grade vs level

Suppose we are interested in the relationship of student grade with the course level (categorical).

  1. Using grade as your outcome variable, fit a linear regression model to investigate this question. Comment on the nature of the relationship between course level and student grades (this should not be a coefficient interpretation, but instead a description of a general trend, or lack thereof).

  2. State the null and alternative hypotheses associated with the research question in part a.

H_0:

H_a:

  1. What type of test do we need here: a t-test for a single model coefficient, the overall F-test, or a nested F-test?

  2. What is the p-value associated with this hypothesis test? Do we have enough evidence to reject the null hypothesis, using a significance threshold of 0.05?

Exercise 3: F-tests for grade vs enrollment

Suppose we are interested in the relationship between course enrollment and student grades.

  1. Again, use grade as your outcome variable, and fit a linear regression model to investigate this question.

  2. State the null and alternative hypotheses associated with the research question in part a.

H_0:

H_a:

  1. What type of test do we need here: a t-test for a single model coefficient, the overall F-test, or a nested F-test?

  2. What is the p-value associated with this hypothesis test? Do we have enough evidence to reject the null hypothesis, using a significance threshold of 0.05?

Exercise 4: More F-tests

Suppose we are now interested in the association between course grade and enrollment for classes of the same level, i.e. when controlling for class level.

  1. Write a model statement in the form E[Y | X] = ... that will produce a statistical model that will allow us to answer our scientific question. Replace Y and X, where appropriate, with response and predictor variables.

E[Y | X] = ___

Which coefficient(s) in your model is the one that is relevant to your research question?

  1. What are the relevant null and alternative hypotheses that address the scientific question in part (a)?

  2. Fit the model you wrote in part (a), calculate a p-value, and report the results of the hypothesis test in part (b).

Reflection

F-tests are useful when the null hypothesis you wish to test is such that more than one covariate is simultaneously equal to a specific number (typically zero). What scenarios, outside of those shown in this example, can you think of where a relevant scientific hypothesis you want to test involves more than one covariate being simultaneously equal to zero?

Extra Practice

Exercise 5: Repeat

Repeat Exercise 4, supposing we are instead interested in the association between course grade and course level for classes of the same enrollment.

Exercise 6: Guess the p-value

Consider 2 models of grades, 1 of which you’ve used before (but maybe named something else):

grades_model_A <- lm(grade ~ enroll + level, MacGrades)
summary(grades_model_A)

grades_model_B <- lm(grade ~ level, MacGrades)
summary(grades_model_B)
  1. Using your notation from Exercise 4, state the hypotheses that would be tested by running anova(grades_model_B, grades_model_A). DO NOT YET RUN THIS CODE!

  2. WITHOUT RUNNING THE anova() CODE: What will be the p-value reported by anova(grades_model_B, grades_model_A)? What’s your reasoning?

  3. Check your intuition:

anova(grades_model_B, grades_model_A)

Exercise 7: But is it a “good” model?

In the above exercises, you should have concluded that, when controlling for the other, both course enrollments and course level are significantly associated with grades. But are they “good” predictors? Let’s explore grades_model_A in more depth.

  1. The relationships of grade with enroll and level are statistically significant, but are they practically significant / contextually meaningful?

  2. Grades vary from student to student within and across courses. What percentage of this variation is explained by course enrollments and level?

  3. Putting this all together, what’s your overall conclusion about the relationship here?

Exercise 8: Reference categories

Our final research question pertains to whether or not there is a relationship between course grade and department. Again, use course grade as the outcome variable in your linear regression model.

  1. State the null and alternative hypotheses in colloquial language associated with the relevant hypothesis test.

H_0:

H_a:$

  1. Fit a linear regression model, and conduct your hypothesis testing procedure to answer the research question posed in this Exercise. State your conclusions accordingly (you do not need to interpret any regression coefficients, just state and interpret the results of your hypothesis test!).

  2. Are any of the individual department p-values significant? What do these p-values tell us, and why is this not contradictory to your answer in part (b)?





Wrap Up

  • Finish and study the activity!
  • Project Milestone 2 due Friday night
  • PS 8 due next Monday
  • Reminders:
    • Remember the important basics: complete activities, ask for help, study on a regular basis, start assignments early.
    • If you miss class on a project day, it’s your responsibility to communicate with your group, find out what you missed, and make equivalent contributions outside of class. If missing class on project days becomes a pattern, it will impact your project grade.
    • Your final project grade will reflect your individual engagement and contributions to the project. You are not expected to be perfect. Rather, you should put effort, attention, and reflection into the following areas:
      • communication
        • contributing to group discussions
        • including all other group members in discussion
        • creating a space where others felt comfortable sharing ideas and mistakes
        • communicating with your group outside class
      • report contributions
        • contributing to the content / analysis in the report
        • contributing to the writing of the report
        • contributing to the editing of the report





Solutions

Warm Up

Example 1: Nested Models

Solution
  1. Models 1, 2, 4 and 6.

  • Model 1 is nested in Model 2 TRUE
  • Model 1 is nested in Model 3 TRUE
  • Model 1 is nested in Model 4 FALSE
  • Model 2 is nested in Model 3 TRUE
  • Model 3 is nested in Model 2 FALSE
  • Model 2 is nested in Model 6 FALSE
  1. You could compare the Adjusted R^2 values from each model, and note that the one with a higher adjusted R^2 is better by this metric. Multiple R^2 would not be a good metric, because the larger model (within the nesting structure) will always have a higher R^2 value.

Example 2: t-tests & “overall” F-tests

Solution
bikes <- read_csv("https://mac-stat.github.io/data/bikeshare.csv") %>% 
  rename(rides = riders_casual)

bike_model_1 <- lm(rides ~ weekend + temp_feel + temp_actual, bikes)
summary(bike_model_1)

Call:
lm(formula = rides ~ weekend + temp_feel + temp_actual, data = bikes)

Residuals:
     Min       1Q   Median       3Q      Max 
-1592.49  -252.52   -17.41   181.65  1973.12 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1370.647     86.438 -15.857   <2e-16 ***
weekendTRUE   813.556     36.312  22.405   <2e-16 ***
temp_feel      11.997      8.712   1.377   0.1689    
temp_actual    15.884      9.459   1.679   0.0935 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 443.8 on 727 degrees of freedom
Multiple R-squared:  0.584, Adjusted R-squared:  0.5823 
F-statistic: 340.2 on 3 and 727 DF,  p-value: < 2.2e-16
  1. H_0: \beta_1 = \beta_2 = \beta_3 = 0
    H_a: at least 1 of \beta_1, \beta_2, \beta_3 is not 0

    We reject H_0 (p-value < 0.05). We have statistically significant evidence that at least 1 of weekend, temp_feel, or temp_actual are associated with rides.

  2. .

    • H_0: \beta_2 = 0
      H_a: \beta_2 \ne 0

      We fail to reject H_0 (p-value > 0.05). We do not have statistically significant evidence that when controlling for weekend and temp_actual, there’s an association between rides and temp_feel. That is, if we already know weekend and temp_actual, then temp_feel doesn’t provide significantly more information about rides.

    • Similarly, we do not have statistically significant evidence that when controlling for weekend and temp_feel, there’s an association between rides and temp_actual.

  3. It might seem like temperature variables aren’t good predictors of ridership, but weekend status is. As we’ll learn, that’s not correct!!

Example 3: F-tests for nested models

Solution

Let’s take the temperature variables out of the model:

# Refit bike_model_1 (just for ease of comparison)
bike_model_1 <- lm(rides ~ weekend + temp_feel + temp_actual, bikes)

# Fit a new model
bike_model_2 <- lm(rides ~ weekend, bikes)
  1. H_0: \beta_2 = \beta_3 = 0
    H_a: at least 1 of \beta_2 and \beta_3 are non-0

    We reject H_0 (p-value < 0.05). We have statistically significant evidence that when controlling for weekend status, rides is associated with at least one of temp_feel or temp_actual. That is, even if we already know weekend status, at least one of temp_actual or temp_feel provides significantly new information about rides.

  2. temp_feel and temp_actual aren’t significant when the other is in the model, but at least 1 of these is significantly associated with rides. These predictors are likely multicollinear – thus if we already have one of them in our model, we likely don’t need the other. Our next step would be to take 1 out of the model and see what happens to the significance of the other.

Exercises

Exercise 1: Explore

Solution

Answers will vary a lot from student to student!

  1. It is reasonable to assume that course grade varies by department as well as course level and instructor. Certain instructors may grade more strictly (or curve more) than others, and similarly, this can vary across department due to cultural norms within the department. As the level of a course gets higher, I would expect grades to perhaps get lower, since courses with higher numbers are expected to be more difficult. Then again, students perhaps “care” more about such courses, and may put in more effort to get a higher grade. I doubt semester plays a significant role in determining course grades, though it is possible that Fall semester first-years or Spring semester seniors have worse grades, on average. We don’t have course year as a variable in our data, so we would be unable to examine this relationship. As enrollment in a course goes up, I would expect grades to decrease, since professors have less time to dedicate to individual students when course enrollment is high.

  2. Explore the relationship between course grades and other variables in the data. Make at least four visualizations, and describe any patterns you observe.

# Exploratory plots

# course grade vs. enrollment
MacGrades %>%
  ggplot(aes(enroll, grade)) +
  geom_jitter() +
  theme_classic() +
  ggtitle("Course grades by enrollment numbers")

# course grade vs. level
MacGrades %>%
  mutate(level = factor(level)) %>%
  ggplot(aes(y = grade, x = level)) +
  geom_boxplot() +
  ggtitle("Course grades by course level")

# course grade vs. level (treating grade as categorical)
MacGrades %>%
  filter(!is.na(grade)) %>%
  mutate(level = factor(level),
         grade = factor(grade)) %>%
  ggplot(aes(x = level, fill = grade)) +
  geom_bar(position = "fill") +
  scale_fill_viridis_d(option = "H") +
  theme_classic() +
  ggtitle("Course grades by course level")

# course grade vs. semester
MacGrades %>%
  filter(!is.na(grade)) %>%
  mutate(grade = factor(grade)) %>%
  ggplot(aes(x = sem, fill = grade)) +
  geom_bar(position = "fill") +
  scale_fill_viridis_d(option = "H") +
  theme_classic() +
  ggtitle("Course grades by semester") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Let's do something fancy and check out how grades have changed over time... this will require
# some string manipulation

MacGrades$year <- MacGrades$sem %>%
  str_replace("FA", "") %>%
  str_replace("SP", "") %>%
  str_replace("S1", "") %>%
  str_replace("S2", "") %>% 
  str_replace("IT", "") %>% as.numeric()

MacGrades %>%
  ggplot(aes(year, grade)) +
  geom_jitter() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_classic() +
  ggtitle("Course grades by year, with least-squares line")

In general, course grades seem to be associated with enrollment numbers. Specifically, when enrollments are greater than 50, we see very few students receiving a course grade lower than a 2.0, which is different than when enrollments are fewer than 50 students. There does not appear to be a clear relationship between course grade and course level, with the exception of 600-level courses. In these cases, every student received either an A or A-. It does seem like the proportion of students who received lower than a 2.67 is greater for 100-level courses than the other course levels.

Exercise 2: F-tests for grade vs level

Solution
mod <- lm(grade ~ level, data = MacGrades)
summary(mod)

Call:
lm(formula = grade ~ level, data = MacGrades)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5108 -0.3620  0.2088  0.4892  0.6380 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.3123669  0.0173797 190.588  < 2e-16 ***
level       0.0004962  0.0000776   6.394 1.75e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5925 on 5707 degrees of freedom
  (437 observations deleted due to missingness)
Multiple R-squared:  0.007112,  Adjusted R-squared:  0.006938 
F-statistic: 40.88 on 1 and 5707 DF,  p-value: 1.746e-10

We observe that as course level goes up, course grades also tend to increase on average.

  1. .

H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = 0

H_1: \text{At least one of } \beta_1, \beta_2, \beta_3, \beta_4 \neq 0 In words, the null is that there is no relationship between course level and course grades, and the alternative is that there is some relationship (either positive or negative) between course level and course grades.

  1. overall F-test

  2. The p-value associated with this hypothesis test is 9.713 x 10^{-13}. We do have enough evidence to reject the null hypothesis. The difference in grades by level is not statistically significant.

Exercise 3: F-tests for grade vs enrollment

Solution
mod <- lm(grade ~ enroll, data = MacGrades)
summary(mod)

Call:
lm(formula = grade ~ enroll, data = MacGrades)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4529 -0.3871  0.2265  0.5534  0.8448 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.4842350  0.0173790 200.485  < 2e-16 ***
enroll      -0.0031336  0.0006683  -4.689 2.81e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5934 on 5707 degrees of freedom
  (437 observations deleted due to missingness)
Multiple R-squared:  0.003838,  Adjusted R-squared:  0.003664 
F-statistic: 21.99 on 1 and 5707 DF,  p-value: 2.806e-06

H_0: \beta_1 = 0

H_1: \beta_1 \neq 0

  1. t-test (an overall F-test would also work since there’s only one non-intercept coefficient in this model!!!)

  2. The p-value associated with this hypothesis test is 2.806 x 10^{-06}. We do have enough evidence to reject the null hypothesis. Note that this p-value could be obtained from either the overall model fit or from the individual coefficient for enroll (they are the same). They may be ever so slightly different when there are few observations in your dataset, but when there are a lot, they will be exactly identical.

Exercise 4: More F-tests

Solution

E[grade | enroll, level] = \beta_0 + \beta_1 enroll + \beta_2 level200 + \beta_3 level300 + \beta_4 level400 + \beta_5 level600

The relevant coefficient that answers our scientific question is \beta_1, or the coefficient that corresponds to enrollment.

The relevant null and alternative hypotheses are:

H_0: \beta_1 = 0 H_1: \beta_1 \neq 0 We do not need to conduct an F-test to complete this hypothesis testing procedure, since our hypothesis involves only a single regression coefficient.

mod <- lm(grade ~ enroll + level, data = MacGrades)
summary(mod)

Call:
lm(formula = grade ~ enroll + level, data = MacGrades)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5201 -0.3602  0.1930  0.5165  0.7600 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.376e+00  2.656e-02 127.088  < 2e-16 ***
enroll      -2.185e-03  6.896e-04  -3.168  0.00154 ** 
level        4.311e-04  8.021e-05   5.375 7.98e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.592 on 5706 degrees of freedom
  (437 observations deleted due to missingness)
Multiple R-squared:  0.008856,  Adjusted R-squared:  0.008508 
F-statistic: 25.49 on 2 and 5706 DF,  p-value: 9.512e-12

We have statistically significant evidence of a relationship between enrollment and course grade, for courses of the same level (p = 0.001058). We reject the null hypothesis that there is no relationship between enrollment and course grade, adjusting for course level.

Exercise 5: Repeat

Solution

Our model statement is identical to that in Exercise 4, but the relevant coefficients are \beta_2, \beta_3, \beta_4, and \beta_5.

The relevant null and alternative hypotheses are:

H_0: \beta_2 = \beta_3 = \beta_4 = \beta_5 = 0

H_1: \text{At least one of } \beta_2, \beta_3, \beta_4, \beta_5 \neq 0

We do need to conduct an F-test to complete this hypothesis testing procedure, since our hypothesis involves more than one regression coefficient.

# Same model as in Question 10, we just now need to do an F-test!
mod <- lm(grade ~ enroll + level, data = MacGrades)
smaller_mod <- lm(grade ~ enroll, data = MacGrades)

anova(smaller_mod, mod)
Analysis of Variance Table

Model 1: grade ~ enroll
Model 2: grade ~ enroll + level
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   5707 2009.8                                  
2   5706 1999.7  1    10.123 28.887 7.977e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We have statistically significant evidence of a relationship between course level and course grade, for courses of the same enrollment (p = 2.19 \times 10^{-10}). We reject the null hypothesis that there is no relationship between course level and course grade, adjusting for enrollment.

Exercise 6: Guess the p-value

Solution
grades_model_A <- lm(grade ~ enroll + level, MacGrades)
summary(grades_model_A)

Call:
lm(formula = grade ~ enroll + level, data = MacGrades)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5201 -0.3602  0.1930  0.5165  0.7600 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.376e+00  2.656e-02 127.088  < 2e-16 ***
enroll      -2.185e-03  6.896e-04  -3.168  0.00154 ** 
level        4.311e-04  8.021e-05   5.375 7.98e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.592 on 5706 degrees of freedom
  (437 observations deleted due to missingness)
Multiple R-squared:  0.008856,  Adjusted R-squared:  0.008508 
F-statistic: 25.49 on 2 and 5706 DF,  p-value: 9.512e-12
grades_model_B <- lm(grade ~ level, MacGrades)
summary(grades_model_B)

Call:
lm(formula = grade ~ level, data = MacGrades)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5108 -0.3620  0.2088  0.4892  0.6380 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.3123669  0.0173797 190.588  < 2e-16 ***
level       0.0004962  0.0000776   6.394 1.75e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5925 on 5707 degrees of freedom
  (437 observations deleted due to missingness)
Multiple R-squared:  0.007112,  Adjusted R-squared:  0.006938 
F-statistic: 40.88 on 1 and 5707 DF,  p-value: 1.746e-10
  1. H_0: \beta_1 = 0
    H_a: \beta_1 \ne 0

  2. will vary

  3. The p-value is 0.001058, which is the same as the t-test on the enroll coefficient in summary(grades_model_A). The p-value is the same because they’re testing the same hypotheses!

anova(grades_model_B, grades_model_A)
Analysis of Variance Table

Model 1: grade ~ level
Model 2: grade ~ enroll + level
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1   5707 2003.2                                
2   5706 1999.7  1    3.5176 10.037 0.001542 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exercise 7: But is it a “good” model?

Solution
  1. The enroll predictor doesn’t seem practically meaningful (to me!). There’s a small decrease in average grade per enrollment (0.002) when controlling for level. Given the small class sizes at Mac, this per-enrollment decrease wouldn’t “add up” to a meaningful drop in average grade in most courses.

    But the level predictor does seem practically meaningful (to me!). For example, an increase of 0.6 in the average grade from a 100-level to a 600-level course, when controlling for enrollment, is a big jump on the GPA scale.

  2. Only 1.271%. (This is the R-squared value.)

  3. It seems that course grades are associated with course level and enrollments, but these factors are not “good” predictors of grades.

Exercise 8: Reference categories

Solution

H0: There is no relationship between course grade and department.

H1: There is some relationship between course grade and department.

mod <- lm(grade ~ dept, data = MacGrades)
summary(mod)

Call:
lm(formula = grade ~ dept, data = MacGrades)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5268 -0.2749  0.1475  0.4515  0.9039 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.5000000  0.4088676   8.560   <2e-16 ***
deptb       -0.2536066  0.4155163  -0.610    0.542    
deptB       -0.2614286  0.4278948  -0.611    0.541    
deptC        0.0268398  0.4106338   0.065    0.948    
deptd       -0.1825995  0.4098240  -0.446    0.656    
deptD        0.0617623  0.4105399   0.150    0.880    
depte       -0.0162500  0.4122607  -0.039    0.969    
deptE        0.1391667  0.4416275   0.315    0.753    
deptF       -0.2251373  0.4104679  -0.548    0.583    
deptg        0.0956522  0.4176615   0.229    0.819    
deptG       -0.3303306  0.4105537  -0.805    0.421    
deptH       -0.0307812  0.4120495  -0.075    0.940    
depti       -0.1474011  0.4111711  -0.358    0.720    
deptI        0.0011538  0.4243020   0.003    0.998    
deptj       -0.0625248  0.4108867  -0.152    0.879    
deptJ       -0.2815172  0.4116777  -0.684    0.494    
deptk        0.0124521  0.4104311   0.030    0.976    
deptK       -0.2431624  0.4123474  -0.590    0.555    
deptL        0.0468000  0.4169648   0.112    0.911    
deptm        0.0224798  0.4099682   0.055    0.956    
deptM       -0.4039440  0.4099067  -0.985    0.324    
deptn        0.1487363  0.4111080   0.362    0.718    
deptN        0.1670000  0.4139469   0.403    0.687    
depto       -0.3616667  0.4255629  -0.850    0.395    
deptO        0.0066856  0.4100242   0.016    0.987    
deptp       -0.0139370  0.4120744  -0.034    0.973    
deptP       -0.1140000  0.4249077  -0.268    0.788    
deptq        0.0001132  0.4104076   0.000    1.000    
deptQ       -0.0644094  0.4120744  -0.156    0.876    
deptR       -0.0381579  0.4110139  -0.093    0.926    
depts       -0.0154839  0.4218507  -0.037    0.971    
deptS       -0.0247500  0.4139469  -0.060    0.952    
deptt       -0.0126923  0.4166562  -0.030    0.976    
deptT       -0.2733962  0.4165106  -0.656    0.512    
deptU       -0.1136842  0.4298486  -0.264    0.791    
deptV       -0.0242500  0.4189646  -0.058    0.954    
deptW       -0.0507164  0.4100863  -0.124    0.902    
deptX       -0.1189865  0.4116209  -0.289    0.773    
deptY        0.0814894  0.4174763   0.195    0.845    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5782 on 5670 degrees of freedom
  (437 observations deleted due to missingness)
Multiple R-squared:  0.06037,   Adjusted R-squared:  0.05407 
F-statistic: 9.586 on 38 and 5670 DF,  p-value: < 2.2e-16

We have statistically significant evidence of a relationship between department and course grades at a significance level of 0.05 (p-value < 2.2 x 10^{-16}). We reject the null hypothesis that there is no relationship between course grade and department.

None of the individual p-values for department are significant! These p-values tell us about whether or not there is a statistically significant difference in course grades between each respective department and the reference department (Department “A”). This doesn’t contradict our answer to part (b) because there are different hypothesis tests that answer different questions!