14.1 The standard error of a regression coefficient
Of course, our primary interest in this course is hypothesis testing. We know that each and every time we generate an estimate from our data, we have to think about it as one estimateof an infinite number of estimates we could have generated if we repeatedly resampled. For all of our estimates, therefore, there is a sampling distribution which describes what sort of range of estimates we will get if we repeatedly sample.
Let’s go back to our old standby of playing god and knowing a population to look at what the sampling distribution for a regression will look like. The mvrnorm() function from the MASS package allows us to draw from the “multivariate normal” we don’t need to know what this is, but just know that we can control the covariance between the two variables, as well as the variance of each variable. In this case we are setting the variance of each of x and y to 1, and the covariance between the two variables as 0.2. Here’s one draw:
Call:
lm(formula = d$y ~ d$x)
Residuals:
Min 1Q Median 3Q Max
-2.2904 -0.7694 0.1109 0.6806 2.5906
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.03485 0.10287 0.339 0.736
d$x 0.08855 0.10909 0.812 0.419
Residual standard error: 1.026 on 98 degrees of freedom
Multiple R-squared: 0.006678, Adjusted R-squared: -0.003458
F-statistic: 0.6588 on 1 and 98 DF, p-value: 0.4189
Focusing just on the slope, the coefficient here is .08855. The hypothesis test that we should fail to reject the null hypothesis that there are no relationship between these two variables, even though I set them up with a positive covariance. So this is a false negative.
What happens when we repeatedly re-sample from this distribution, capturing the beta coefficient each time? Let’s do 6 times:
We get a pretty wide variety of outcomes. Now let’s look at doing this 1000 times. I’m going to record the beta each time, and stack on top of each other all of the regression lines.
The graph that has all of the slopes on it is kind of a “double fan” pattern. All the lines seem to converge at 0,0. Note that both of these variables have, in the population, an expected value of 0.
It is actually true that every single regression line passes through the mean of the x variable and the mean of the y variable. This is very simple to show mathematically.
The regression line is:
\[
\hat{y} = \hat{\alpha} + \hat{\beta}x_i
\]
And we want to know the value of \(\hat{y}\) for when \(x_i = \bar{x}\). Subbing in the equation we generated for \(\hat{\alpha}\):
So when x is equal to it’s mean, \(\hat{y}\) is equal to {y}. The inverse must be true because we are talking about a point on a line.
So every regression line passes through the sample means of the two variables. Looking at the graph to the left then, this makes it clear we will be more certain about the location of the regression line when we are close to the mean of x and y. This will have big ramifications when we start thinking about our ability to use regression for prediction!
Our main interest right now is the sampling distribution to the right. This appears to be a normal distribution centered around the true population slope coefficient.
To prove to ourselves that this is a normal distribution:
eval <-seq(-.2, .6, .001)plot(density(betas), main="Sampling Distribution of Slope coefficient")points(eval, dnorm(eval, mean=mean(betas), sd=sd(betas)), type="l", col="firebrick")
Yeah looks about right!
When we derived the standard error of the mean we learned that two components influenced the size of the standard error and thus the width of the sampling distribution: sample size and the randomness of the underlying variables. Is the same thing true here?
Yes for sure! We can see this in both graphs. The regression lines are much tighter around the “true” slope, and the sampling distribution is much smaller.
What about the randomness of the underlying variables? Again, we’ve seen before that if our underlying variables are more random it’s harder to make inferences about our samples. Does that matter here?
Which version of the above will lead to a lower variance in the slop of the regression line? It’s going to be the top. When x has less variation, small, random, changes in y are going to have large ramifications about where the regression line is. The best way to think about it is that more variation in x leads to more leverage in determining where the regression line is.
For the big reveal, here is the equation that generates the standard error of a regression coefficient:
The standard error is determined by the average squared residual, divided by the variance of x.
summary(m)
Call:
lm(formula = d$y ~ d$x)
Residuals:
Min 1Q Median 3Q Max
-2.59804 -0.62515 0.04112 0.67740 1.70046
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.10725 0.09155 1.172 0.2442
d$x 0.25168 0.12728 1.977 0.0508 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9079 on 98 degrees of freedom
Multiple R-squared: 0.03837, Adjusted R-squared: 0.02856
F-statistic: 3.91 on 1 and 98 DF, p-value: 0.05081
sqrt(mean(m$residuals^2)/sum((d$x-mean(d$x))^2))
[1] 0.1249118
So what makes for a smaller standard error? As we’ve seen, more variance in x will lead to a smaller standard error, all else equal.
What else? Because of the \(n\) in the average formula, all things being equal more \(n\) will make the standard error smaller, but that’s only true if the \(n\) you add have relatively small residuals. In total the numerator of the equation tells us that more and smaller residuals will lead to a low standard error. What leads to smaller residuals? Less variation in \(y\) around the regression line. In other words if the points show a consistent relationship then there is going to be a smaller standard error, which makes a lot of sense!
I want to re-emphasize that once you have an estimate and a standard error literally everything that we’ve done so far is exactly the same.
For the regression \(m\) the standard error is .127, which is the standard deviation of a t distributed sampling distribution with 98 degrees of freedom. Where does 98 come from? The degrees of freedom for a regression is \(n\) minus the number of coefficients you are estimating, including the intercept, which here is 2. So the t distribution looks like:
Same thing would go for a power calculation, or the confidence interval. Once you have the standard error, you are golden.
14.2 Regression with more than one independent variable
First, let’s install and load the “AER” package, which has some interesting built in data. We are just using this package for the data, it’s not a necessary package for regression or anything.
#install.packages("AER")library(AER)
The AER package has some built in data on California schools
Let’s set out a research question for these data. I want to know if the student to teacher ratio in a classroom is associated with better test scores. Specifically, as the student to teacher ratio increases, test scores should decrease.
The null hypothesis of this relationship is that student to teacher ratio has 0 association with test scores
Let’s first create two variables that capture these concepts. STR will be the student to teacher ratio:
Well, these look much messier than the fake data we used earlier, which is about right. even still, there does seem to be a generally negative relationship here, though it’s not a slam dunk, and could be due to random chance.
Let’s use our new tool of regression to describe the relationship we are seeing here.
m <-lm(score ~ STR, data=dat)summary(m)
Call:
lm(formula = score ~ STR, data = dat)
Residuals:
Min 1Q Median 3Q Max
-47.727 -14.251 0.483 12.822 48.540
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 698.9329 9.4675 73.825 < 2e-16 ***
STR -2.2798 0.4798 -4.751 2.78e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.58 on 418 degrees of freedom
Multiple R-squared: 0.05124, Adjusted R-squared: 0.04897
F-statistic: 22.58 on 1 and 418 DF, p-value: 2.783e-06
The coefficient on STR is negative. When STR increases by one unit, test scores decrease by -2.3 units, approximately.
Let’s put this in english. A one unit change in STR means one additional student for every teacher. So…A classroom with one additional student for every teacher, is associated with a 2.3 point decline in test scores on average.
Looking at the hypothesis test for this, the probability of obtaining a relationship this extreme if the truth was that there was no association between STR and test scores approaches zero. We reject the null hypothesis here.
What is the intercept telling us in this case ? When STR is zero, the average test score is 698.
But this is where we have to be smarter than the regression: can STR ever be zero??
14.3 Ommited variable bias
You’ve heard time and time again: correlation does not equal causation. Honestly I’m getting a bit tired of hearing it, as people often use it to mean “Can’t trust those scientists these days!”. Over the next few weeks we’ll formalize what this means and determine situations where correlation does equal causation. For right now, what we can think about doing is ruling out other plausible reasons for why these two variables might be correlated.
Theoretically, what we are interested in are variables that cause both a high student to teacher ratio and low test scores. In other words, might it be the case that some third, unseen, variable is driving this relationship?
There are some classic absurd omitted variable bias examples: There is a negative association between the number of firemen at a fire and property damage. It would be a mistake to conclude that firemen cause property damage. The truth is that both are driven by the severity of the fire.
Or above, we estimated the (fake) impact of years of schooling on income. What might be a third variable that causes someone to have more schooling and a higher income?
In the case we are dealing with here, a potential explanation for why a school might have a high student to teacher ratio and lower test scores might be because it is a school with a high number of children who speak english as a second language. Schools with a high number of children who speak english as a second language may be more likely to be in areas with a poor tax base and therefore less money to hire teachers. Standardized tests have been found to not capture the innate intelligence of those who do not speak english as their primary language.
In our original data, we have a variable that captures the percent of the school who has English as second language. To start let’s create a simple dummy variable that is equal to 1 if a school is above the median for ESL children, and 0 for schools that are below the median.
Just using the naked eye, there does seem to be a large number of schools which have both a high number of children who have English as a second language and also have lower test scores.
Here is what we want to do: Taking into account that there are two different types of schools, can we summarize what the relationship between STR and test scores are?
Mechanically, this is what regression is going to do: Choose a slope line for STR and TWO intercepts, one for each type of school, such that the sum of the squared residuals is minimized.
It’s much shallower. Why? The relationship between STR and test scores was driven, in part, by the fact that schools with high numbers of ESL children have both a high STR and low test scores. Once we take into account that fact (by focusing on the relationship within school types) the association between STR and test scores is much lower.
What does this look like in a regression equation?
m <-lm(score ~ STR + high.esl, data=dat)summary(m)
Call:
lm(formula = score ~ STR + high.esl, data = dat)
Residuals:
Min 1Q Median 3Q Max
-37.857 -10.853 -0.626 10.198 43.834
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 691.3267 8.1315 85.019 < 2e-16 ***
STR -1.3963 0.4171 -3.348 0.000889 ***
high.esl -19.4911 1.5763 -12.365 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.91 on 417 degrees of freedom
Multiple R-squared: 0.3058, Adjusted R-squared: 0.3025
F-statistic: 91.84 on 2 and 417 DF, p-value: < 2.2e-16
Helpfully, we are going to read these coefficients in the same way as we did with univariate regression. Each variable still represents how a one unit change in each independent variable is associated with a coefficient level change in the dependent variable. The major difference is that now these numbers are the effect of this variable holding the other variable constant. Another way to think of this is: accounting for the effect of the other variable, what’s the effect of this variable?
So for STR, holding constant whether a school has a high number of ESL students or not, the effect of an additional student for every teacher is to reduce test scores by 1.4 points.
The other variable high.esl is an indicator, which makes interpretation easier. Just as before a “one unit change” for an indicator variable is just changing categories, here from a low to a high ESL school.
Holding constant the student teacher ratio, schools with an above average number of children who have english as their second language perform nearly 20 points worse than schools with a below average number of childen who have english as their second language, holding constant the student to teacher ratio.
What about the intercept? For univariate regression we learned that the intercept is the average value of y when x is 0. Helpfully, the same thing applies here, with the caveat that now the intercept is the average value of y variable when all variables are equal to zero. In this case, the intercept is the average test scores when STR is zero (not real) and when high.esl is equal to zero.
What then, is the average test scores when STR is equal to zero in high esl schools? It would be 691-19.49! Indeed, when we look at the graph we’ve made, 19.49 is precisely the gap between the red and blue lines.
14.4 Multiple continuous variables.
What about when we have a regression with two continous variables?
The variable “engligh” is the percent of the students who have english as a second language. We think that “controlling” for this variable should have a similar effect on the relationship between STR and test scores as controlling for the the high.esl variable we created above. Schools with a high percentage of ESL students are likely to have higher STR and lower test scores.
But now, thinking about our original scatterplot, we can’t easily classify these dots into “high” and “low” as we did before
But I pretty much never do think about this way, and if you add just one additional variable there is no way to visualize!
At this point, just thinking about this through the regression output is much more helpful.
m <-lm(score ~ STR + english, data=dat)summary(m)
Call:
lm(formula = score ~ STR + english, data = dat)
Residuals:
Min 1Q Median 3Q Max
-48.845 -10.240 -0.308 9.815 43.461
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 686.03224 7.41131 92.566 < 2e-16 ***
STR -1.10130 0.38028 -2.896 0.00398 **
english -0.64978 0.03934 -16.516 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.46 on 417 degrees of freedom
Multiple R-squared: 0.4264, Adjusted R-squared: 0.4237
F-statistic: 155 on 2 and 417 DF, p-value: < 2.2e-16
Again, what we are seeing here is the association between each of these variables and the dependent variable while holding constant the other.
So in plain language: Holding constant the impact of the percent of students who have English as a second language, the impact of one additional student per teacher is to reduce test scores by 1.1 points.
Further, holding constant the impact of the student to teacher ratio, the impact of an additional 1% of students having English as a second language is to reduce test scores by .65 points.
14.5 Many variables
While visualization becomes harder/impossible with additional variables, the logic of what regression coefficients mean is infinitely scalable. The coefficients will always be the impact of a one-unit change in the variable on the dependent variable, holding the other variables constant.
So if we want to assess the impact of STR while controlling for % of students getting subsidized lunch as well as the income of (something, i’m not sure what this data is)
m <-lm(score ~ STR + english + lunch + income, data=dat)summary(m)
Call:
lm(formula = score ~ STR + english + lunch + income, data = dat)
Residuals:
Min 1Q Median 3Q Max
-30.7085 -4.9977 -0.2014 4.8411 28.7000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 675.60821 5.30886 127.261 < 2e-16 ***
STR -0.56039 0.22861 -2.451 0.0146 *
english -0.19433 0.03138 -6.193 1.42e-09 ***
lunch -0.39637 0.02741 -14.461 < 2e-16 ***
income 0.67498 0.08333 8.100 6.19e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.448 on 415 degrees of freedom
Multiple R-squared: 0.8053, Adjusted R-squared: 0.8034
F-statistic: 429.1 on 4 and 415 DF, p-value: < 2.2e-16
Note that regression does not care bout the order of variables on the “right hand” side of the equation:
m <-lm(score ~english + income + STR + lunch, data=dat)summary(m)
Call:
lm(formula = score ~ english + income + STR + lunch, data = dat)
Residuals:
Min 1Q Median 3Q Max
-30.7085 -4.9977 -0.2014 4.8411 28.7000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 675.60821 5.30886 127.261 < 2e-16 ***
english -0.19433 0.03138 -6.193 1.42e-09 ***
income 0.67498 0.08333 8.100 6.19e-15 ***
STR -0.56039 0.22861 -2.451 0.0146 *
lunch -0.39637 0.02741 -14.461 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.448 on 415 degrees of freedom
Multiple R-squared: 0.8053, Adjusted R-squared: 0.8034
F-statistic: 429.1 on 4 and 415 DF, p-value: < 2.2e-16
It’s important to note here that this is just enough information to get yourself in trouble. Specifying a regression is hard, theoretical, work. There are some big pitfalls to be aware of.
Primarily: more variables isn’t necessarily better! Indeed, adding certain types of variables can actually make your estimates worse. Regression should be tricky, and I wouldn’t go out applying this to all sorts of things without learning more about the theory.
Here’s an example: Imagine had a dataset that had whether people smoked or not, whether they got pneumonia or not and whether they died.
You think: better be careful, let’s do
lm(die ~ smoke + pneumonia)
But people who smoke die of pneumonia. The effects of smoking within those who get pneumonia is probably close to zero, or greatly reduced.
This is what we call “post-treatment bias”
Another example is looking at the effects of gender on wages controlling for length of parental leave. You’re controlling away (one of) the reasons why women get paid lower wages: they reasonably take more time off when they have a child.
