Multiple Regression

Multi-variate Regression AKA OLS Regression AKA Multiple Regression

Introduction

Up until this point (except briefly for Factorial ANOVA), we’ve been focused on examining the relationship of two variables at a time. Correlation, associated tests (cor.test(), bi-variate regression), and even t-tests all examine the relationship (or difference) between two variables - and only two - variables. While this can be helpful, obviously many relationships are much more complex. Crime, for instance, is not caused only by a single variable (e.g. poverty), and so simply analyzing crime relative to poverty is unlikely to yeild very complete answers (and would only inform a decision in a limited way).

Additionally, simply repeating analysis, like doing several bi-variate regressions, poses several problems. First, we have the repeated tests issue we discussed relative to t-tests. Second, when we do a bi-variate regression, there is an implicit assumption that “all else is equal” across cases outside of the model. If we have more than one model, we know that all else is not equal, as we are explicitly modeling the variation of another independent variable! This is called specification error (you model should be as complete as possible based on theoretical rationalle) and means that, while you might get an answer, your model is unlikely to be trustworthy. Additional problems include low model accuracy (there are not a lot of single-variable models that will explain a complex phenomena), and inefficiency in estimation.

Multiple Regression

Fortunately, multiple regression - which is a very simple expansion of bi-variate regression - can handle many independent variables. How many? Well, it depends on the size of your dataset, but generally it can handle n-1 predictors (we usually refer to the nuber of predictors, or independent variables, by k). So, for 10 cases, you could, in theory handle 9 predictors!¹

Additionally, how we add independent variables is also fairly straightforward. Before we had: \(\hat{y} = \alpha + \beta{x}\) as our model. To make it multiple regression, we simply add additional predictor variables: \(\hat{y} = \alpha + \beta_1{x_1} + \beta_2{x_2} ... + \beta_k{x_k}\).

The effects of this, on the back end, is to now be able to calculate the independent effects of each of the variables, while taking the other variables into account. Thus, when we look at our regression coefficients (which we will do in an example in a moment), we are seeing what the effect of that variable is while taking all of the other variables into account. This helps us understand a lot of information that we cannot get with regular bi-variate regressions, as will be seen below.

Example

As with most kinds of inferential statistics, an example can help. In this case, we’re going to examine the NYFS dataset we were using frequently during the beginning of the semester.

Introduction

In the current example we’re going to use marijuana use as the dependent variable of interest. To set it up further, let’s pretend that we’ve got a theory out there that suggests that marijuana usage can be explained by a combination of other drug use, minor offending, and a person’s age. We’ve measured all of those variables, and want to know both which variables are important - while accounting for the others - as well as how well our model performs as a whole at explaining marijuana use. This, in turn, will inform us regarding how well our theory holds up.

Data

The current analysis uses the National Youth and Family Study (NYFS) survey from 1974 (\(n = 578\)). It contains responses on a variety of crime-related variables, as well as control variables (like age, gender, etc.). A summary table with the relevant independent variables can be seen in Table 1, below.

	Overall (N=578)
marij
Mean (SD)	22.742 (85.397)
Range	0.000 - 900.000
age
Mean (SD)	20.024 (1.826)
Range	18.000 - 24.000
minor
Mean (SD)	7.163 (33.435)
Range	0.000 - 377.000
hdrugs
Mean (SD)	4.875 (29.500)
Range	0.000 - 438.000
alcohol
Mean (SD)	51.983 (73.321)
Range	0.000 - 500.000

We can see from Table 1, above, that several of our variables are widely dispersed, but aside from age, have relatively low means. This suggests that the variables have a strong negative skew. We can visualize this regarding them as well, as in Figure 1, below.²

In addition, we can also examine correlations among the variables, to get some sense of how well they move together - at least as pairs.

analysis<-subset(nyfs, select = c(hdrugs, marij, minor, alcohol, age))
cor(analysis)

##              hdrugs        marij      minor   alcohol          age
## hdrugs   1.00000000  0.327020741 0.10597577 0.1729969 -0.030382286
## marij    0.32702074  1.000000000 0.26204497 0.2405332 -0.006128634
## minor    0.10597577  0.262044975 1.00000000 0.2001076  0.044023770
## alcohol  0.17299693  0.240533216 0.20010758 1.0000000  0.141420490
## age     -0.03038229 -0.006128634 0.04402377 0.1414205  1.000000000

Here we can see that the correlation is relatively low among most variables. The exception seems to be hard drugs and marijuana, which have a Pearson’s r of .327 - which is generally considered moderate. A correlation plot gives us similar information in graphical form.

Hypotheses

Based on the review of the (pretend) literature, we have several hypotheses we can test. In the case of regression, we’re testing (for each independent variable, and then for the whole model) the hypotheses that the slopes of each variable are equal to 0. If this sounds familiar, it should! It’s basically the same test we use for correlation tests (cor.test()), only we can do it with a bunch of variables at once!³

H₀₍₁₎: \(\beta_{alcohol} = 0\)

H₀₍₂₎: \(\beta_{age} = 0\)

H₀₍₃₎: \(\beta_{hdrugs} = 0\)

H₀₍₃₎: \(\beta_{minor} = 0\)

And the associate research hypotheses are:

H₁₍₁₎: \(\beta_{alcohol} \neq 0\)

H₁₍₂₎: \(\beta_{age} \neq 0\)

H₁₍₃₎: \(\beta_{hdrugs} \neq 0\)

H₁₍₃₎: \(\beta_{minor} \neq 0\)

Model

Now that we have our hypotheses, we can test them in a single model. This is a very simple extension of the original bi-variate regression model. In formula form, for the bi-variate model, we had:

\[ \hat{y} = \alpha + \beta_{x} \]

For multiple regression, we have:

\[ \hat{y} = \alpha + \beta_{1}{x_1} + \beta_{2}{x_2} +...+\beta_{k}{x_k}\] This form is followed in how we define the model in R, which is seen below.

model<-lm(marij~age+minor+alcohol+hdrugs, data=analysis)
stargazer(model, type = 'html', title = "Table 2. Regression Results", dep.var.labels=c("Marijuana Use"), single.row = TRUE, align=TRUE)

Table 2. Regression Results

	Dependent variable:
	Marijuana Use
age	-1.341 (1.789)
minor	0.518*** (0.099)
alcohol	0.182*** (0.046)
hdrugs	0.804*** (0.112)
Constant	32.526 (35.713)
Observations	578
R2	0.182
Adjusted R2	0.176
Residual Std. Error	77.524 (df = 573)
F Statistic	31.790*** (df = 4; 573)
Note:	p<0.1; p<0.05; p<0.01

Brief Digression

Before, when we were talking about bi-variate regression, I skipped a line of the table - the \(F\)-statistic. The reason I did so is because it isn’t really as important for the bi-variate use-case (and honestly, is somewhat less important for regression as a whole). However, if you recall the conversation about multiple-test problems from the section on t-tests, you may wonder why we don’t have the same issue with regression, since we’re (basically) conducting a set of t-tests. The reasons is two-fold: First, we are testing the slopes independently (while holding the other variables constant), so it’s somewhat different than simply repeating t-tests. Second, we actually do use an omnibus test, the same \(F\)-test we use in ANOVA, in fact. And, indeed, that’s what’s reported in the \(F\)-Statistic line in the output.

What does that mean? Well, like ANOVA, it’s an omnibus test of all model effects at once. So, one way to treat it is exactly like the ANOVA \(F\). You can check it first and know that, if it’s significant, you’ll have at least one significant slope (coefficient) in the model. A second way to think of it is as the total-model test. If we’re less interested in specific variables’ effects, and more interested in straight prediction, the \(F\) is a useful statistic to see whether or not the model, as a whole, predicts the dependent variable significantly.

In all events, you’ll need to report the \(F\) statistic and associated degrees of freedom when you report regression model results. Chances are, for most of our use cases, that will be essentially your last interaction with it, but it’s helpful to know that it’s doing the same thing as the ANOVA \(F\) test.

Model Results

Examining the results in Table 2, above, we can see that the overall model is significant (\(F = 31.790(4,573), p < .05\)). Specifically, three variables have slopes significantly different from 0. Minor offending (\(\beta = .518, p < .05\)), alcohol use (\(\beta = .182, p < .05\)) and hard drug use (\(\beta = .804, p < .05\)) are all significant predictors of marijuana use. Age, on the other hand, is not a significant predictor of marijuana use. Thus, we reject three of our null hypotheses but must accept the null hypothesis regarding age.

Looking at the overall level of explanation, we can see that the model is generally weak with an Adj. \(R^2\) of just .176. This means that with the combination of the independent variables within the model, we have explained only 18% of the variation in the use of marijuana. Additionally, this only provides weak support for our theory, as while several of the variables were significant predictors, the overall model’s performace was generally poor (as judged by the \(R^2\) value).

[You’ll notice that I’m not visualizing the model. That’s because the model now has more dimensions (each predictor is a dimension) than can be realistically visualized. We will, on the other hand, use some visualization for model-checking in the next section.]

Checking the Model

So, off the bat, we know that several of these variables are not particularly good in terms of meeting the assumptions of Regression. The variables should be normally distributed (oops - look at the histograms), homoskedastic (maybe okay - we’ll see), and have a linear relationship (probably fine). We can check these, along with model accuracy, the same ways we checked it within bi-variate regression. We will use the function predict() to predict values for the model and then compare those with the actual values, to start.

analysis$predict<-predict(model) # predict values for y
qqplot(analysis$marij, analysis$predict)
abline(0,1)

We can see from the above plot that, well, the model’s not really fitting the data that well. That’s not surprising given the model’s overall \(R^2\). On average, we would expect a model that’s really good, to fit the data pretty well (though that’s not always the case). How to figure out what is wrong with the model is slightly more complicated. We can use the plot() function along with our model, to try to ascertain some answers. More information on interpreting these plots can be found here: https://data.library.virginia.edu/diagnostic-plots/.

par(mfrow=c(2,2))
plot(model)

There is a lot going on in the plots above. Generally, we wouldn’t include these in a write-up, but it’s important you figure out how your model is doing so you don’t give bad information. We’ll take the charts one-by-one.

Residuals vs. Fitted

The residuals vs. fitted plot is probably the easiest to understand. It take the distances between predicted points and the observed points (in the dependent variable) and plots them against the predicted points. It is a visual way of representing error. Think of it this way: a perfect model would have all of the points exactly on the regression line. So the distance from the line for every point would be 0. The prediction (because every point would be on the line) would also perfect. So, in the Residuals vs. Fitted plot, we would see a horizontal line at 0 with all of our points on it. Because we do have error (in every model), we know that not all of the points will fall on the line. But, if the model is good, and in particuar, if there is no problem with linearity in the variables, the error will be equally distributed around 0 across the full model.

Incidentally, in all of the plots, R tries to be helpful by showing outliers and cases that may be affecting the model. Regression is pretty robust to model assumptions, but outliers can definitely pull the model in one direction or another (more on that in a bit).

Normal QQ Plot

The Normal QQ Plot is slightly tougher to understand, but in some ways easier to read. Essentially, it’s a test to see if your residuals (one way to measure error) are normally distributed. We would expect them to be in a good model, and using the QQ plot, we should see the dots fall on a straight line. However, in our case, we definitely see the beginning and end of the distribution are well off the line. This can indicate that you are seeing some heteroskedasticity, which may in turn be caused by extreme values (of which we have some).

Scale-Location Plot

More specifically, we can check heteroskedasticity with the Scale-Locaiton plot. This examines whether the spread of your residuals is the same over the full range of your dependent variable. A perfect model would have an equal distribution of error across the full model - but we definitely don’t. This is unsurprising given the very high number of cases near 0 and the relatively small number of cases at the high end of the distribution of marijuana use (and the others, excepting age).

Residuals vs. Leverage

I’m going to replot this plot separately, because it needs to be larger to read:

plot(model, which = 5)

There is a lot going on in this plot, and you need to know a couple of things going in. First, Cook’s distance is a measure of how much an outlier a given datapoint is. But it’s not as easy as just saying “look at how high this value is!” You need to see it within the model output, which is what this is showing us. The dotted red lines are the areas we’re interseted in finding things outside of. To take an easy case, look at the point labeled ‘44’. That case has a cooks distance of over 1 (which is why it’s outside the red dotted lines), so we may want to think about removing it. Same with 99 and perhaps 15. Plotting it this way shows us which cases are likely most affecting the analysis.

To get an idea of how much of an impact the cases are having, let’s take out cases 44 and 99 and rerun the model.

analysis2<-analysis[-c(44,99),] # get rid of those cases
model2<-lm(marij~age+minor+alcohol+hdrugs, data=analysis2)
summary(model2)

## 
## Call:
## lm(formula = marij ~ age + minor + alcohol + hdrugs, data = analysis2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -388.18  -14.50   -9.47   -6.95  882.15 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.60863   33.31834   0.769 0.442446    
## age         -0.92056    1.66893  -0.552 0.581448    
## minor        0.22698    0.09886   2.296 0.022042 *  
## alcohol      0.14524    0.04325   3.358 0.000837 ***
## hdrugs       1.26116    0.13285   9.493  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 72.31 on 571 degrees of freedom
## Multiple R-squared:  0.1933, Adjusted R-squared:  0.1876 
## F-statistic:  34.2 on 4 and 571 DF,  p-value: < 2.2e-16

We can see that while the same variables are significant, that minor offending has a lower level of signifiance. Additionally, our \(R^2\) value jumped some, indicating a better model fit without those two cases. We can also see if it helped any of our other diagnostics.

While there are clearly still some issues in the model, we can see that some of these do look, albeit slightly, better (particularly the Residuals vs. Fitted plot). The Residuals vs. Leverage chart now has a couple of new potentially influential cases that we might want to deal with as well.

What else can we do?

We can see that our model here still is not great, in terms of how well its fitting the data. This (to me, and truly this is more art than science) looks like it’s because of how the data distributions are shaped - very long tails. There’s a couple of ways to deal with this (like getting rid of outliers, more on this in a bit), one of which is transformation.

Transformation

Transformation is a pretty straightforward idea. We take a variable (often the dependent, but can be others) and change the distribution with a known function. The most common version of this is what we call a “log-transform.” This is when we take the natural log of the depenedent variable and model against that. It basically squeezes the distribution so that the distance between 5-50 is the same as 50-500. This means that, while before we had a distribution that was really spread out, now it’s spread out, but less so.⁴

This is perhaps easist to see in an example:

par(mfrow=c(1,2))
hist(analysis$marij)
hist(log(analysis$marij+1)) # add 1 or you will get NA for any number that's a 0

You can see that the data not is much more condensed. Whereas before the spread was over 800, now its over about 7. Additionally, we can see that, though it’s still far from normal, it looks better than it did with the raw data. What happens when we model it?

analysis$logmarij<-log(analysis$marij +1)
model3<-lm(logmarij~age+minor+alcohol+hdrugs, data=analysis)
summary(model3)

## 
## Call:
## lm(formula = logmarij ~ age + minor + alcohol + hdrugs, data = analysis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2076 -0.7198 -0.5853  0.4915  5.7376 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.9040289  0.6667204   1.356 0.175654    
## age         -0.0177093  0.0333924  -0.530 0.596083    
## minor        0.0070207  0.0018447   3.806 0.000156 ***
## alcohol      0.0074383  0.0008579   8.671  < 2e-16 ***
## hdrugs       0.0119174  0.0020828   5.722  1.7e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.447 on 573 degrees of freedom
## Multiple R-squared:  0.226,  Adjusted R-squared:  0.2206 
## F-statistic: 41.82 on 4 and 573 DF,  p-value: < 2.2e-16

There’s a couple of important things to note here. First, look at the \(R^2\) value - it has gone up (and I’m using the original data, so the influential cases are still there). Second, look at the significant coefficients - they haven’t changed in terms of their significance. Why? Because, like correlation, regression is interested in movement and transformations do not change how variables move together. So, we can get better model performance without much of a trade-off in results.

What is the tradeoff? Interpretability. Now, rather than saying for a one-unit increase in minor offending, the level of marijuana usage goes up by .518 uses per year (which is what the original model had), we have to say for a one unit increase in minor ofending there is an associated increase of .007 in the log of marijuana use. Interpreting what a “log-use” of marijuana is is not possible, so we sacrifice that ability to interpret clearly from the model. We can get back to an interpretable coefficient if we exponentiate the coefficient, subtract one from this number, and multiply by 100. This gives the percent increase (or decrease) in the response for every one-unit increase in the independent variable. We can see how to do this for model3, below:

(exp(model3$coefficients)-1)*100

## (Intercept)         age       minor     alcohol      hdrugs 
## 146.9532705  -1.7553367   0.7045395   0.7466008   1.1988701

So, for every one unit increate in minor offending, we see a 70% increase in marijuana usage.

Whether or not you want to transform these back to an interpretable coefficient is largely a product of what it is you’re testing. If you’re interested in seeing which variables are best at explaining marijuana usage, and not in what effects a specific predictor has, you don’t need to move your coefficients back to an interpretable frame, but you do have to remember to speak about them in log as long as you haven’t.

Let’s look quickly and see if our diagnostic charts are any better after the log-transform.

par(mfrow=c(2,2))
plot(model3)

While still not great, they are definitely much improved over what they were before.

Outliers

We saw before that using the Residuals vs. Leverage plot how we can identify potentially influential cases. But we can use other plots to do that as well. One example is just a plain Cook’s Distance plot:

plot(model3, which = 4)

Here, we can see that there are three potential outliers. Let’s remove them and rerun the log transformed model we used before.

analysis3<-analysis[-c(44,289,420),]
model4<-lm(logmarij~age+minor+alcohol+hdrugs, data=analysis3)
summary(model4)

## 
## Call:
## lm(formula = logmarij ~ age + minor + alcohol + hdrugs, data = analysis3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7090 -0.7216 -0.5991  0.4773  5.7845 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.8534052  0.6551344   1.303    0.193    
## age         -0.0141271  0.0328162  -0.430    0.667    
## minor        0.0084380  0.0020465   4.123 4.30e-05 ***
## alcohol      0.0064784  0.0008677   7.466 3.11e-13 ***
## hdrugs       0.0214485  0.0031248   6.864 1.75e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.422 on 570 degrees of freedom
## Multiple R-squared:  0.2501, Adjusted R-squared:  0.2449 
## F-statistic: 47.54 on 4 and 570 DF,  p-value: < 2.2e-16

Once again, we can see an \(R^2\) increase and we’re now up to about 25% of the variation in (log of) marijuana use being explained. Additionally, we can see that, once again, it did not affect which variables were significant.

par(mfrow=c(2,2))
plot(model4)

Looking at our residual plots, we can see that, again they’re somewhat better (though still not perfect). If we wanted, we could examine the outliers presented within the plots above, and continue to refine the model and make it better.

The big outlier question

So the mechanics of detecting/removing outliers is pretty straightforward, but the bigger question is should you remove them? It’s really not always clear! One question, and I think perhaps the best question to start with, is: What is your goal? If it’s purely prediction, as we saw, you can improve your ability to explain variance (and thus predict better) by removing anythign designated as an outlier. On the other hand, if you’re interested in explaining a phenomena, it becomes a more challenging question. Should you remove a case which is simply extreme, even if it’s a possible value? If you do, you run the risk of (perhaps) removing information that actually is important for informing your investigation. For instance, let’s take case 44 (who is identified as an outlier in the original dataset).

head(analysis[44,])

When we examine case 44’s other variable values, it becomes clear why he’s an outlier. Look at the hdrugs variable! He reports using marijuana 100 times a year and using hard drugs 438 times in the past year. That’s a total of over 500 times a year, or taking drugs more than once a day, every day, for the whole year.

Now, is that possible? Maybe. The question is does it inform your analysis. Clearly this guy takes way more drugs than almost anyone else, so does he help explain marijuana usage for the “normal” person? Probably not. On the other hand, he might provide information relative to how much impact very heavy hard druge use is on marijuana use.

All that to say, it’s not very cut-and-dry when to eliminate outliers. You may chose to do so, but you always want to make sure you have a justification!

Findings Explained

So, what does all this mean? I would draw your attention to the fact that nothing really changed since the original analysis (outside of \(R^2\))! Why? Because even though it’s clear we violated some assumptions (even in our best model), generally speaking we know regression is pretty robust even when we do violate its assumptions. Additionally, though we transformed variables and elimitated influential cases, we didn’t really change the analysis, and the findings have been clear (in terms of statistical significance) from the beginning!

Conclusion

This. Was. A. Lot. We learned how to do a multiple regression, how to check the models and how to deal with certain types of problems through getting rid of outliers or transforming variables. Not every analysis includes all of these steps in a write-up, but every analyst should include them in building the models. As it stands, we didn’t get very close to a perfect model and…well it’s no really all that surprising. Real-life data is messy.

All that said, we got where we got, and as long as we make sure to show our work, and write up what we feel the limitations are, we are doing the work we are supposed to!

---

1. However, that’s generally too many predictors (indepdent variables) to try to include with such a small dataset, and your model would likely be pretty poor.↩

2. We will learn about transforming variables to help with some of this dispersion in a later lesson.↩

3. And as we’ll see shortly, we can do an omnibus test - \(F\) actually - to test them all at once.↩

4. We can actually log-transform any of our variables, or all of them. It just gets more complicated to interpret.↩