One X Variable

Module 4.4: Inference, Fit, and OLS in R

Author

Affiliation

Alex Cardazzi

Old Dominion University

All materials can be found at alexcardazzi.github.io.

Statistical Inference

Just like how we tested if our sample mean was statistically different from some hypothesized value in the last module, we can test if our regression coefficients are different from hypothesized values as well.

In the notes, we use zero, but the hypothesized value could have been anything.

Focusing on $\beta_1$, and skipping the math, we can write the standard error of $\beta_1$ as follows:

\[se_{\beta_1} = \sqrt{\frac{s^2}{\sum(x_i - \overline{x})^2}}\]

Here, $s^2 = \frac{SSR}{n-2}$ is the variance of the residuals, and $s$ is called the residual standard error. You will see it near the bottom of summary() output as “Residual standard error”. In words, it is the typical size of a miss: how far, in units of $Y$, a point usually sits from the regression line.

As a reminder, we need the standard error to calculate a confidence interval and/or perform a hypothesis test. Also, you can think of a standard error as a measure of how precisely we estimated $\beta$.

To test if your estimate of $\beta_i$ is different from some hypothesized number, we first need to generate a test statistic:

\[t = \frac{\beta_i - \#}{se_{\beta_i}}\]

Generally, people choose 0 for $\#$, since a coefficient of 0 would imply that the regressor (e.g., GPA) has no relationship with the outcome variable (e.g., SAT score). This way, the hypothesis test is set up such that it tests whether $X$ is related to $Y$.

Once we calculate this test statistic, we can calculate a $p$-value and decide whether to reject the null hypothesis that our regressor is not related to the outcome variable.

As a final note about inference (i.e., hypothesis tests), we need to make two important assumptions about our errors ($\hat{\epsilon}_i$) in order for the above t-statistic to be valid.

Normality: we assume that our errors are independent, random draws from a normal distribution.
Homoskedasticity: this big word means that the normal distribution we draw our errors from has a constant variance across all points. See below for examples of homoskedastic errors and heteroskedastic errors.

Plot

One graph of homoskedastic errors next to another graph of heteroskedastic errors.

Goodness of Fit

Once we estimate $\beta_0$ and $\beta_1$, we can also begin to talk about how well the model “fits” the data. In other words, what fraction of $Y$’s variance is explained by $X$.

Before getting too much further with this, I want to be very clear that goodness-of-fit measures do not determine whether your model is good. You may be able to explain a lot of $Y$’s variation, but still have a bad model. In addition, your model might have poor overall explanatory power, but still do a good job measuring the effect of $X$ on $Y$, which is what we care about. We will discuss this more as the course progresses.

Below are three equations for Total Sum of Squares (SST), Explained Sum of Squares (SSE), and Residual Sum of Squares (SSR). Hover over the equations for explanations.

\[\text{SST} = \sum (Y_i - \overline{Y})^2\]

\[\text{SSE} = \sum (\hat{Y}_i - \overline{Y})^2\]

\[\text{SSR} = \sum \hat{\epsilon_i}^2 = \sum (Y_i - \hat{Y}_i)^2\]

Once we have defined these three measures, we can connect all three by the following equation:

\[SST = SSE + SSR\]

If our model is very good at explaining the variation in $Y$, we would expect SSE to be very large relative to SSR. For a perfect model, we would see $SST = SSE$. Therefore, we can create a ratio between these two numbers, which will give us a percentage of the variation that is explained by the model. We call this ratio $R^2$, and define it as follows:

\[R^2 = \frac{SSE}{SST} = 1 - \frac{SSR}{SST}\]

Below, I demonstrate what I mean when I say $R^2$ should not be a be-all and end-all measure of how good a model is. I simulated a dataset of 100 $X$ values. Then, I created two $Y$ variables: $Y_1 = 5X + \epsilon$ and $Y_2 = 5X + 2.5\epsilon$. Both $Y$ variables have the same relationship with $X$, but the random errors are magnified in the second. The estimated coefficients are very similar, but the $R^2$ measure is quite different. Again, we care more about the coefficients being estimated accurately rather than the model itself having explanatory power.

Plot

OLS in R

So far, we have derived the OLS estimator, calculated statistics for hypothesis tests, and looked at measures for goodness of fit. Now, we are going to talk more about application.

First, to estimate a linear regression model, we will use the lm() function.

lm(), which stands for “linear model”, accepts the following arguments:

formula: Formulas in R take on the following form: y ~ x1 + x2. In the case for GPAs and SAT scores, we would write: sat_sum ~ hs_gpa.
data: This is the data you plan to use. This tells R where to look for the variables in your formula.
subset: An optional vector that specifies a subset of observations used to fit the line. For example, maybe we want to fit two lines: one for men and one for women. We would use subset = df$sex == 1 to estimate the model using only data from female students. Another reason to use this would be to eliminate outliers, etc.

We are going to keep examining the GPA data (data; documentation). Let’s read the data into R:

Code

df <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/openintro/satgpa.csv")

Let’s estimate the SAT and GPA model from the previous section of notes.

Code

reg_sat_gpa <- lm(sat_sum ~ hs_gpa, data = df)
reg_sat_gpa

Output


Call:
lm(formula = sat_sum ~ hs_gpa, data = df)

Coefficients:
(Intercept)       hs_gpa  
      66.99        11.36

The output of lm() returns the call you used to generate the output in addition to the estimated coefficients. These coefficients are exactly what we previously estimated with cov(), var(), cor(), etc. Not only does lm() output the coefficients, but it also returns both the fitted values (reg_sat_gpa$fitted.values) and the residuals/errors (reg_sat_gpa$residuals).

Code

head(data.frame(observed_y = df$sat_sum,
                fitted_y = reg_sat_gpa$fitted.values,
                residual_y = reg_sat_gpa$residuals))

Output

  observed_y fitted_y residual_y
1        127 105.6232  21.376776
2        122 112.4411   9.558875
3        116 109.6003   6.399667
4         95 109.6003 -14.600333
5        107 112.4411  -5.441125
6        111 112.4411  -1.441125

We can plot our regression line using the abline() function.

Code

library("scales")
plot(df$hs_gpa, df$sat_sum, las = 1, pch = 19,
     col = alpha("black", 0.2),
     xlab = "GPA", ylab = "SAT")
abline(reg_sat_gpa, lty = 1, col = "gold", lwd = 4)

Plot

Using the summary() function generates some important information about our regression and the individual coefficients.

Code

summary(reg_sat_gpa)

Output


Call:
lm(formula = sat_sum ~ hs_gpa, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-42.055  -8.669  -0.351   8.803  34.240 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  66.9885     2.4441   27.41   <2e-16 ***
hs_gpa       11.3632     0.7535   15.08   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.9 on 998 degrees of freedom
Multiple R-squared:  0.1856,    Adjusted R-squared:  0.1848 
F-statistic: 227.4 on 1 and 998 DF,  p-value: < 2.2e-16

In the table portion, each coefficient has an Estimate, Std. Errors, t value, and Pr(>|t|) In addition, you can find the Multiple R-Squared value, which is the measure for goodness of fit.

Regression for Male Students:

Code

reg_sat_gpa_m <- lm(sat_sum ~ hs_gpa,
                    data = df, subset = df$sex == 2)
summary(reg_sat_gpa_m)

Output


Call:
lm(formula = sat_sum ~ hs_gpa, data = df, subset = df$sex == 
    2)

Residuals:
    Min      1Q  Median      3Q     Max 
-31.997  -8.305  -0.455   8.574  33.287 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   55.851      3.579   15.61   <2e-16 ***
hs_gpa        13.715      1.079   12.71   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.33 on 482 degrees of freedom
Multiple R-squared:  0.2512,    Adjusted R-squared:  0.2496 
F-statistic: 161.7 on 1 and 482 DF,  p-value: < 2.2e-16

Regression for Female Students:

Code

reg_sat_gpa_f <- lm(sat_sum ~ hs_gpa,
                    data = df, subset = df$sex == 1)
summary(reg_sat_gpa_f)

Output


Call:
lm(formula = sat_sum ~ hs_gpa, data = df, subset = df$sex == 
    1)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.497  -8.323   0.153   8.877  31.153 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   71.280      3.183   22.39   <2e-16 ***
hs_gpa        11.019      1.003   10.98   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.55 on 514 degrees of freedom
Multiple R-squared:  0.1901,    Adjusted R-squared:  0.1885 
F-statistic: 120.6 on 1 and 514 DF,  p-value: < 2.2e-16

Take a look at the Pr(>|t|) column. This contains a p-value corresponding to the hypothesis test that $\beta_i \neq 0$. The value <2e-16 means that the p-value is smaller than 0.0000000000000002. R doesn’t report p-values that are smaller than this for brevity. In addition, R includes asterisks next to p-values that are statistically significant.

Nothing next to the p-value indicates $0.1 < p \leq 1$.
. next to the p-value indicates $0.05 < p \leq 0.1$
* next to the p-value indicates $0.01 < p \leq 0.05$
** next to the p-value indicates $0.001 < p \leq 0.01$
*** next to the p-value indicates $0 < p \leq 0.001$

This is helpful when looking to see if your coefficients are significant at different levels.

How can we interpret these results? Let’s start with $\beta_0$. For a male student with a GPA of 0, the expected SAT score is 55.85. For a female student with a GPA of 0, the expected SAT score is 71.28. In addition, for any GPA, we can compute the expected SAT score by plugging it into the equation.

Thankfully, R has a pre-defined function called predict(). This function accepts two arguments: a lm object and newdata. newdata must have the same column names as variables in the lm object.

Code

# Create a new data.frame of GPAs from 2 to 4 by 0.25 increments.
pred_df <- data.frame(hs_gpa = seq(2, 4, by = 0.25))
pred_df$sat_m <- predict(reg_sat_gpa_m, newdata = pred_df)
pred_df$sat_f <- predict(reg_sat_gpa_f, newdata = pred_df)
pred_df

Output

  hs_gpa     sat_m     sat_f
1   2.00  83.28204  93.31803
2   2.25  86.71088  96.07282
3   2.50  90.13971  98.82762
4   2.75  93.56855 101.58242
5   3.00  96.99738 104.33721
6   3.25 100.42622 107.09201
7   3.50 103.85505 109.84680
8   3.75 107.28389 112.60160
9   4.00 110.71272 115.35639

It seems that female students consistently outperform male students on the SAT, according to these data.

Another helpful way to look at this would be to plot both lines together with the underlying data.

In this plot, it’s very clear that the biggest difference between male and female students is when GPAs are relatively low. The gap then decreases as GPAs get larger.

While summary() provides helpful information about our regressions, the output is not very aesthetically pleasing. We are going to revisit the modelsummary package, and use the titular function modelsummary(). This function takes many arguments, but for simplicity, we will just discuss a few:

models: A named list of models we want to include in the table.
title: The title of our table.
estimate: This is a bit complicated, but we are going to set this equal to "{estimate}{stars}". This way, the table will display coefficients, standard errors, and stars to denote statistical significance. This is the norm. If you do not like the stars, you can remove them by not including estimate in the function call.
coef_map: A named vector variable names. This will rename the variable names from what they are labeled as in the data to something more human readable.
gof_map: A vector of desired goodness of fit measures for each model. We are going to stick with "nobs" and "r.squared" for the time being.

Code

library("modelsummary")
options("modelsummary_factory_default" = "kableExtra")

regz <- list(`All Students` = reg_sat_gpa,
             `Male` = reg_sat_gpa_m,
             `Female` = reg_sat_gpa_f)
coefz <- c("hs_gpa" = "High School GPA",
           "(Intercept)" = "Constant")
gofz <- c("nobs", "r.squared")
modelsummary(regz,
             title = "Effect of GPA on SAT Scores",
             estimate = "{estimate}{stars}",
             coef_map = coefz,
             gof_map = gofz)

Output

Effect of GPA on SAT Scores
	All Students	Male	Female
High School GPA	11.363***	13.715***	11.019***
	(0.754)	(1.079)	(1.003)
Constant	66.988***	55.851***	71.280***
	(2.444)	(3.579)	(3.183)
Num.Obs.	1000	484	516
R2	0.186	0.251	0.190