# load packages
library(tidyverse) # for data wrangling and visualization
library(tidymodels) # for modeling
library(openintro) # for the duke_forest dataset
library(scales) # for pretty axis labels
library(knitr) # for pretty tables
library(kableExtra) # also for pretty tables
library(patchwork) # arrange plots
# set default theme and larger font size for ggplot2
ggplot2::theme_set(ggplot2::theme_bw(base_size = 20))SLR: Mathematical models for inference
Announcements
HW 01 due TODAY at 11:59pm
Lab 01 due Thursday at 11:59pm
- Questions?
Statistics experience - due Mon, Nov 20 at 11:59pm
Course policy reminders
-
- HW and labs accepted up to 2 days late.
- 5% deduction for each 24-hour period the assignment is late.
-
- Can use on HW and individual labs
Lowest HW and lowest lab grade dropped at the end of the semester.
Topics
Define mathematical models to conduct inference for the slope
Use mathematical models to
calculate confidence interval for the slope
conduct a hypothesis test for the slope
construct intervals for predictions
Computational setup
The regression model, revisited
df_fit <- lm(price ~ area, data = duke_forest)
tidy(df_fit) |>
kable(digits = 3)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 116652.325 | 53302.463 | 2.188 | 0.031 |
| area | 159.483 | 18.171 | 8.777 | 0.000 |
Inference, revisited
- Earlier we computed a confidence interval and conducted a hypothesis test via simulation:
- CI: Bootstrap the observed sample to simulate the distribution of the slope
- HT: Permute the observed sample to simulate the distribution of the slope under the assumption that the null hypothesis is true
- Now we’ll do these based on theoretical results, i.e., by using the Central Limit Theorem to define the distribution of the slope and use features (shape, center, spread) of this distribution to compute bounds of the confidence interval and the p-value for the hypothesis test
Mathematical representation of the model
\[ \begin{aligned} Y &= \text{Model} + \text{Error} \\[8pt] &= f(X) + \epsilon \\[8pt] &= E(Y|X) + \epsilon \\[8pt] &= \mu_{Y|X} + \epsilon \\[8pt] &= \beta_0 + \beta_1 X + \epsilon \end{aligned} \]
where the errors are independent and normally distributed:
. . .
- independent: Knowing the error term for one observation doesn’t tell you anything about the error term for another observation
- normally distributed: \(\epsilon \sim N(0, \sigma_\epsilon^2)\)
Side note: Expected value and variance
Let \(X\) be a random variable and \(a\) and \(b\) be constants.
Expected value
\[ E(aX+b) = E(aX) + E(b) = aE(X) + b \]
Variance
\[ Var(aX + b) = a^2Var(X) \]
Mathematical representation, visualized
\[ Y|X \sim N(\beta_0 + \beta_1 X, \sigma_\epsilon^2) \]
- Mean: \(\beta_0 + \beta_1 X\), the predicted value based on the regression model
- Variance: \(\sigma_\epsilon^2\), constant across the range of \(X\)
- How do we estimate \(\sigma_\epsilon^2\)?
Regression standard error
Once we fit the model, we can use the residuals to estimate the regression standard error, the average distance between the observed values and the regression line
\[ \hat{\sigma}_\epsilon = \sqrt{\frac{\sum_\limits{i=1}^n(y_i - \hat{y}_i)^2}{n-2}} = \sqrt{\frac{\sum_\limits{i=1}^ne_i^2}{n-2}} \]
. . .
Why divide by \(n - 2\)?
Why do we care about the value of the regression standard error?
Standard error of \(\hat{\beta}_1\)
The standard error of \(\hat{\beta}_1\) quantifies the sampling variability in the estimated slopes
\[ SE_{\hat{\beta}_1} = \hat{\sigma}_\epsilon\sqrt{\frac{1}{(n-1)s_X^2}} \]
. . .
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 116652.33 | 53302.46 | 2.19 | 0.03 |
| area | 159.48 | 18.17 | 8.78 | 0.00 |
Mathematical models for inference for \(\beta_1\)
Hypothesis test for the slope
Hypotheses: \(H_0: \beta_1 = 0\) vs. \(H_A: \beta_1 \ne 0\)
. . .
Test statistic: Number of standard errors the estimate is away from the null
\[ T = \frac{\text{Estimate - Null}}{\text{Standard error}} \\ \]
. . .
p-value: Probability of observing a test statistic at least as extreme (in the direction of the alternative hypothesis) from the null value as the one observed
\[ p-value = P(|t| > |\text{test statistic}|), \]
calculated from a \(t\) distribution with \(n - 2\) degrees of freedom
Hypothesis test: Test statistic
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 116652.33 | 53302.46 | 2.19 | 0.03 |
| area | 159.48 | 18.17 | 8.78 | 0.00 |
\[ T = \frac{\hat{\beta}_1 - 0}{SE_{\hat{\beta}_1}} = \frac{159.48 - 0}{18.17} = 8.78 \]
Interpreting the test statistic
The test statistic is 8.78. What is the best interpretation?
- The estimated slope of 159.48 is 8.78 standard errors away from the mean.
- The estimated slope of 159.48 is 8.78 standard errors above the mean.
- The estimated slope of 159.48 is 8.78 standard errors above 0, the hypothesized mean.
- The estimated slope of 159.48 is 8.78 standard errors away from 0, the hypothesized mean.
Hypothesis test: p-value
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 116652.33 | 53302.46 | 2.19 | 0.03 |
| area | 159.48 | 18.17 | 8.78 | 0.00 |
Hypothesis test: p-value
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 116652.33 | 53302.46 | 2.19 | 0.03 |
| area | 159.48 | 18.17 | 8.78 | 0.00 |
A more exact p-value
2 * pt(q = 8.78, df = 96, lower.tail = FALSE)[1] 6.19602e-14Interpreting the p-value
What does the p-value mean in the context of the data?
- The probability there is no linear relationship between area and price is approximately 0.
- Given there is no linear relationship between area and price, the probability of observing a slope of 159.48 is approximately 0.
- Given there is a linear relationship between area and price, the probability of observing a slope of 159.48 is approximately 0.
- Given there is no linear relationship between area and price, the probability of observing a slope of 159.48 or more extreme is approximately 0.
Understanding the p-value
| Magnitude of p-value | Interpretation |
|---|---|
| p-value < 0.01 | strong evidence against \(H_0\) |
| 0.01 < p-value < 0.05 | moderate evidence against \(H_0\) |
| 0.05 < p-value < 0.1 | weak evidence against \(H_0\) |
| p-value > 0.1 | effectively no evidence against \(H_0\) |
Important
These are general guidelines. The strength of evidence depends on the context of the problem.
Hypothesis test: Conclusion, in context
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 116652.33 | 53302.46 | 2.19 | 0.03 |
| area | 159.48 | 18.17 | 8.78 | 0.00 |
- The data provide convincing evidence that the population slope \(\beta_1\) is different from 0.
- The data provide convincing evidence of a linear relationship between area and price of houses in Duke Forest.
Confidence interval for the slope
\[ \text{Estimate} \pm \text{ (critical value) } \times \text{SE} \]
. . .
\[ \hat{\beta}_1 \pm t^* \times SE_{\hat{\beta}_1} \]
where \(t^*\) is calculated from a \(t\) distribution with \(n-2\) degrees of freedom
Confidence interval: Critical value
# confidence level: 95%
qt(0.975, df = nrow(duke_forest) - 2)[1] 1.984984# confidence level: 90%
qt(0.95, df = nrow(duke_forest) - 2)[1] 1.660881# confidence level: 99%
qt(0.995, df = nrow(duke_forest) - 2)[1] 2.628016
95% CI for the slope: Calculation
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 116652.33 | 53302.46 | 2.19 | 0.03 |
| area | 159.48 | 18.17 | 8.78 | 0.00 |
\[\hat{\beta}_1 = 159.48 \hspace{15mm} t^* = 1.98 \hspace{15mm} SE_{\hat{\beta}_1} = 18.17\]
. . .
\[ 159.48 \pm 1.98 \times 18.17 = (123.50, 195.46) \]
95% CI for the slope: Computation
tidy(df_fit, conf.int = TRUE, conf.level = 0.95) |>
kable(digits = 2)| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 116652.33 | 53302.46 | 2.19 | 0.03 | 10847.77 | 222456.88 |
| area | 159.48 | 18.17 | 8.78 | 0.00 | 123.41 | 195.55 |
Intervals for predictions
Intervals for predictions
- Suppose we want to answer the question “What is the predicted sale price of a Duke Forest house that is 2,800 square feet?”
- We said reporting a single estimate for the slope is not wise, and we should report a plausible range instead
- Similarly, reporting a single prediction for a new value is not wise, and we should report a plausible range instead
Two types of predictions
Prediction for the mean: “What is the average predicted sale price of Duke Forest houses that are 2,800 square feet?”
Prediction for an individual observation: “What is the predicted sale price of a Duke Forest house that is 2,800 square feet?”
. . .
Which would you expect to be more variable? The average prediction or the prediction for an individual observation? Based on your answer, how would you expect the widths of plausible ranges for these two predictions to compare?
Uncertainty in predictions
Confidence interval for the mean outcome: \[\large{\hat{y} \pm t_{n-2}^* \times \color{purple}{\mathbf{SE}_{\hat{\boldsymbol{\mu}}}}}\]
. . .
Prediction interval for an individual observation: \[\large{\hat{y} \pm t_{n-2}^* \times \color{purple}{\mathbf{SE_{\hat{y}}}}}\]
Standard errors
Standard error of the mean outcome: \[SE_{\hat{\mu}} = \hat{\sigma}_\epsilon\sqrt{\frac{1}{n} + \frac{(x-\bar{x})^2}{\sum\limits_{i=1}^n(x_i - \bar{x})^2}}\]
. . .
Standard error of an individual outcome: \[SE_{\hat{y}} = \hat{\sigma}_\epsilon\sqrt{1 + \frac{1}{n} + \frac{(x-\bar{x})^2}{\sum\limits_{i=1}^n(x_i - \bar{x})^2}}\]
Standard errors
Standard error of the mean outcome: \[SE_{\hat{\mu}} = \hat{\sigma}_\epsilon\sqrt{\frac{1}{n} + \frac{(x-\bar{x})^2}{\sum\limits_{i=1}^n(x_i - \bar{x})^2}}\]
Standard error of an individual outcome: \[SE_{\hat{y}} = \hat{\sigma}_\epsilon\sqrt{\mathbf{\color{purple}{\Large{1}}} + \frac{1}{n} + \frac{(x-\bar{x})^2}{\sum\limits_{i=1}^n(x_i - \bar{x})^2}}\]
Confidence interval
The 95% confidence interval for the mean outcome:
new_house <- tibble(area = 2800)
predict(df_fit, new_house, interval = "confidence", level = 0.95) |>
kable()| fit | lwr | upr |
|---|---|---|
| 563205.5 | 529351 | 597060.1 |
. . .
We are 95% confident that mean sale price of Duke Forest houses that are 2,800 square feet is between $529,351 and $597,060.
Prediction interval
The 95% prediction interval for an individual outcome:
new_house <- tibble(area = 2800)
predict(df_fit, new_house, interval = "prediction", level = 0.95) |>
kable()| fit | lwr | upr |
|---|---|---|
| 563205.5 | 226438.3 | 899972.7 |
. . .
We are 95% confident that predicted sale price of a Duke Forest house that is 2,800 square feet is between $226,438 and $899,973.
Comparing intervals
Extrapolation
Using the model to predict for values outside the range of the original data is extrapolation.
. . .
Calculate the prediction interval for the sale price of a “tiny house” in Duke Forest that is 225 square feet.
. . .
No, thanks!
Next class
Multiple linear regression