Bi-variate Regression in R
Bi-variate Regression
Introduction
Much of what we need to know regarding relationships boils down to correlations. Correlation tests, as we discussed, allow us to know what relationships are real. But, problematically, they don’t give us all the information we may want. For instance, they don’t tell us how much change in the y-variable is cause by the x-variable per unit increase. I know that sounds complicated, but its likely something you did in high school. In this document, we are going to cover bi-varite regression - which is regression with two variables. On the way, we’ll compare it to correlation tests, and we’ll see that while they’re largely similar, they do differ. We’ll also see that, in many ways, regression is similar to ANOVA - and in fact we can answer largely similar questions.
The reason these similarities exist is because of their reliance (aside from correlation tests) on the general linear model (GLM). I’ll do a separate document on that model, but know that what we do here will set us up to understand the GLM, and therefore the rest of regression for the semester.
Correlation Tests vs. Linear Regression
The good news is, there really isn’t a lot of new stuff when we talk about linear regression - just new ways of using some things you’re already familiar with. Often, when we’re looking at bi-variate regression and correlation, the question of their difference arises. The reality is, they’re not really that different. Mostly, it’s in what we’re attempting to find out. In the subsections below, I’ll run over a couple of easy examples of comparision between correlation and bi-variate regression. From there, we’ll move on to a larger example of bi-variate regression, which will set us up for a conversation on multi-variate regression - which we’ll address separately.
Correlation Tests
You’ll recall from the previous lesson that correlation tests are test to see if things “move together.” In other words, as one variable increases (x), what happens to the other (y). An example of a simple correlation is in the code below.
cor.test(mtcars$wt, mtcars$mpg)
Pearson's product-moment correlation
data: mtcars$wt and mtcars$mpg
t = -9.559, df = 30, p-value = 1.294e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.9338264 -0.7440872
sample estimates:
cor
-0.8676594 In this example, we can see that, based on a t-test (null hypothosis is that Pearson’s r = 0), we can reject the null and say that the variables vary together. In this case, because it’s a negative correlation, we can that the the sign on Pearson’s r (cor, in cor.test()) is negative. This means, as wight (wt) increases, miles-per-gallon (mpg) decreases. Because Pearson’s r is -.87, we would say this relationship is strong.
In addition, we can calculate the r-squared value for the correlation test by squaring the Pearson’s r. In this case, it would be roughly .75, meaning we explain about 75% of the variation in mpg with wt.
Bi-variate Regression
So, given that we can figure out how the variables are related, and how strong the relationship is, why would we bother with another analytical technique? Two reasons. First, in correlation tests you’re limited to understanding a single pair of variables at a time. While we’re going to focus on the same thing in regression to start, regression can be extended (theoretically) to an infinite number of variables. Additionally, regression can handle multiple levels of measurement, while (regular) correlation can really only handle interval/ratio data. Finally, while correlation gives you the size of the relationship, regression allows you to calculate the exact change in y based on the change in x. This may not seem that important to begin with, but knowing this allows us not only to know what relationships are important, but to predict what values would be - even beyond the sample we have available.
We’re going to keep the first example simple, by using exactly the same data that we use above.
model<-lm(mtcars$mpg~mtcars$wt)
summary(model)
Call:
lm(formula = mtcars$mpg ~ mtcars$wt)
Residuals:
Min 1Q Median 3Q Max
-4.5432 -2.3647 -0.1252 1.4096 6.8727
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
mtcars$wt -5.3445 0.5591 -9.559 1.29e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10You can see right away that the code we use is a bit different. First, we’re using the formula format to specify the “model” (which is just a way to say our regression). We then use the summary() function to pull the information that we need to interpret the model. We will cover residuals and checking the model later, so let’s start with some things that look familiar.
Coefficients
Examining the summary for the initial model, you can see that there is a section titled “coefficients.” This is really the primary are of interest in the model, for most purposes. It tells us what our individual regression values are (estimates), what their standard error is (how far the line misses by - more on this in a moment), t-values, and Pr(>|t|) - which is just a way to designate p-values. The asterics are shortcuts to know what level you have significance at (with the Signif. codes underneath).
One thing you’ve no doubt noticed is that, rather than a single row with the information, we have two rows. The first of these is designated (Intercept). This is the value of y (in our case, MPG), when the value of x is 0. Now, if you think about what that really means, it’s the MPG we would expect if a car weight 0 pounds. Of course that’s impossible, so we generally ignore it. It’s important for the calculation of the regression (and for plotting), but ultimately it’s not a parameter that makes a lot of sense for most models. In some (rare) instances, we may care what the Intercept value is, and whether it’s significant, but generally we do not.
So, since we can essentially dismiss the Intercept, we can move to our actual variable of interest - wt. As with most statistical tests, what we care about is the p-value. Like usual, in order to reject the null hypothesis (in the case of regression, that there is no relationship between the x and y variables), we need the p-value to be < .05. Helpfully, the three astericks let us know that the p-value is far below this (p = .00000000000000002). So, as with correlation tests, we know from this test that the relationship is significant. If it were not significant, we would ignore it.
Speaking of correlation tests, take a second to compare the t-value and the p-value to the correlation test we did in the section above. What do you notice? They’re exactly the same! That’s because we’re using the exact same test! However, if you look at the “Estimate” column you’ll see that the estimate for mtcars$wt is not the same as the Pearson’s r. What gives? Why the difference? And what is that number?
Regression coefficients and their interpretation
While correlation coefficients are useful, the way they’re calculated limits their ability to take into account more than one variable. Regression coefficients (insofar as we are concerned with here) are estimated using Ordinary Least Squares (OLS). We will get into this method in a moment, but suffice to say, the goal is to minimize error of a line representing the relationship between two variables. One way to think of it (though it’s not exact) is as the “average” relationship between two variables. By calculating that line - called the regression line - we are able to do a couple of things that are not possible with correlation which, in turn, affect how we can interpret the coefficients.
Specifically, the coefficient estimate for the variables (except the Intercept) is the specific amount y changes based on a one-unit increase in x. To put it into concrete terms based on the above, for every 1000lbs1 increase in a car’s weight (which is x) we see a -5.3445 mpg decrease (which is y).
Why might we care about this? Well, imagine we’re designing a car and we anticipate its weight at 1000lbs (we’re into tiny cars). We can use the regression above to predict its weight.
\(y = intercept + -5.3445(1)\) which is \(31.9406 = 37.2851 -5.3445\)
So, what this suggests is that a 1000lb car would get roughly 32mpg, based on our current dataset. We can also predict what a 10000lb car would get in terms of mpg.
\(y = intercept + -5.3445(10) = -16.1599\) which is \(-16.1599 = 37.2851 - 53.445\)
You’ll notice that this leads to an impossible value, -16.1599mpg, but the point is that we can predict exact values for any amount of weight. This is something we cannot do using correlation tests, and can be very useful if we want to, say predict crime based on population or other types of crime.
Examining the next column, we can see that there are standard errors for the coefficient. This is the average amount that the coefficient “misses” by. What this means is easiest to see in the visualization below.
library("ggpubr")
ggscatter(mtcars, x = "wt", y = "mpg",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Weight", ylab = "Miles per Gallon")
The gray area surrounding the regression (aka “best fit”) line is the confidence interval - which is the area we can be 95% sure that the real regression line falls within (based on our sample). The standard error is the average amount the line misses by across the full line. The larger the standard error, the worse the estimate provided by the coefficient. Thus, when you see large errors, there is a much lower liklihood that the variable will be significant.
The t-value is the same exact t-test as the cor.test() provides - it’s essentially how we know if the relationship between the variable represented by the coefficient and the y variable (in this case mpg) is real. This is the same way we tell if the Pearson’s r’s difference from 0 is real.
R-squared and Other Values
(I’m going to reproduce the information in the model from above here. It’s just for visual help, nothing has changed.)
summary(model)
Call:
lm(formula = mtcars$mpg ~ mtcars$wt)
Residuals:
Min 1Q Median 3Q Max
-4.5432 -2.3647 -0.1252 1.4096 6.8727
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
mtcars$wt -5.3445 0.5591 -9.559 1.29e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10Risidual standard error is something we will cover when we look at how to tell how “good” your model is. For now, think of it as a way to understand how well your model “fits” the data.
More important for our purposes is “Multiple R-squared,” which is the same as the r-squared we calculated manually above. First, notice that it includes “multiple.” While right now we only have one independent variable (wt), as we add variables we will continue to use this measure. So, while we have one independent variable, we just treat it exactly tthe same as we did r-squared. In this case, that means that 75% of the dependent variable is explained by the independent variable. Pretty good!
Lastly, we can take a look at the F-statistic. If you’re wondering if this is the same statistic as the F we used for ANOVA - the answer is yes! It’s exactly the same test (though caculated slightly differently) for exactly the same reason. When we have more than one independent variable (so not yet), this represents a single test of all the variables in the model. More on this when we do multiple regression.
Weight, in the mtcars dataset, is by 1000lb increments.↩
Comments
Post a Comment