Multiple Regression with Categorical Variables
Regression with Non-Continuous Variables
Introduction
So, regression turns out to be a very useful tool for both understanding the impact of a set of independent variables and predicting values of \(\hat{y}\) for some dependent variable. But, thus far, we’ve only looked at continuous variables, so you might be wondering if we can use categorical variables, and the answer is “yes”!
However, using them requires a bit more care in interpretation! While unit increases are easy to understand with continuous variables, we have to think a bit differently about categorical variables. Additionally, we have to be very careful about how they’re coded! Make sure you’ve got your factors ordered in the appropriate way (if ordinal). Once you’ve done that, then the challenge (and it’s not a big one) is just in the interpretation.
As always, it’s easiest to see in an example!
Breif Note
By this point, I’m assuming you understand what we’re testing and how to generate a hypothesis for regression. If you don’t, please go back and rewatch the videos on multiple regression and look at that Rmd file. I’m not going to go over that again, though I may refer to “hypotheses” and will assume you know what I’m talking about.
Example
For this example, I’m going to use some capital punishment data I found on Kaggle. It containes information on individuals executed since the reinstatement of the Death Penalty in 1976 (https://www.kaggle.com/usdpic/execution-database). I’ve included s table of some of our variables of interest, below.
| Overall (N=1442) | |
|---|---|
| Age | |
| Mean (SD) | 41.480 (9.472) |
| Range | 22.000 - 77.000 |
| Federal | |
| No | 1439 (99.8%) |
| Yes | 3 (0.2%) |
| Victim.Count | |
| Mean (SD) | 1.537 (4.496) |
| Range | 1.000 - 168.000 |
| Race | |
| Asian | 6 (0.4%) |
| Black | 495 (34.3%) |
| Latino | 120 (8.3%) |
| Native American | 16 (1.1%) |
| Other | 2 (0.1%) |
| White | 803 (55.7%) |
| Method | |
| Electrocution | 158 (11.0%) |
| Firing Squad | 3 (0.2%) |
| Gas Chamber | 11 (0.8%) |
| Hanging | 3 (0.2%) |
| Lethal Injection | 1267 (87.9%) |
We can see from the above, that there are some categories with very small percentages of cases in them. I’m going to go ahead and collapse those categories - particularly in race and execution method.
| Overall (N=1442) | |
|---|---|
| race3 | |
| Black | 495 (34.3%) |
| Other | 144 (10.0%) |
| White | 803 (55.7%) |
| method3 | |
| Electrocution | 158 (11.0%) |
| Lethal Injection | 1267 (87.9%) |
| Other | 17 (1.2%) |
Now that the variables have been recoded to where there are at least 1% of cases in each, we can enter them into the regression. Importantly, I’ve imported them as factors from the original csv (see my code with the option stringsasfactors=T). Additionally, these do not have a particular order. If they were ordinal, I would need to make sure that they were ordered correctly. As these variables are nominal, I don’t have to worry about that (more on this later).
In this example, the dependent variable is age. That is, we’re asking the question about what variables are associated with increases or decreases of the age of the person executed. Let’s go ahead and run the model to see how the results look different from when we have only continuous variables.
model<-lm(Age~Federal+Victim.Count+race3+method3, data=capital)
summary(model)##
## Call:
## lm(formula = Age ~ Federal + Victim.Count + race3 + method3,
## data = capital)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.564 -6.879 -1.082 5.928 37.065
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.95228 0.79568 46.441 < 2e-16 ***
## FederalYes 6.57952 6.49731 1.013 0.311
## Victim.Count -0.07371 0.06587 -1.119 0.263
## race3Other -0.98941 0.88158 -1.122 0.262
## race3White 3.49214 0.53007 6.588 6.24e-11 ***
## method3Lethal Injection 3.19320 0.78511 4.067 5.02e-05 ***
## method3Other -2.05970 2.36383 -0.871 0.384
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.244 on 1435 degrees of freedom
## Multiple R-squared: 0.05154, Adjusted R-squared: 0.04758
## F-statistic: 13 on 6 and 1435 DF, p-value: 2.339e-14Here we have R’s basic output. We can see that, largely, it looks the same. The intercept is still at the top (and, per usual, we’re going to ignore it), and we can see that the omnibus test of model significance, \(F\), is significant, indicating we will have significant variables within the model.
When we check the specific variables, we can see that there are some differences in what the output looks like! We now have “race3Other” and “race3White” as race categories, and “method3Lethal Injection” and “method3Other” as our method variables. So what’s going on here?
Well, R has split the variable, because it’s a factor, into its component categories and tested them one at a time. You’ll notice that while we have 3 categories in each of our categorical variables, only two appear in the model. There’s a reason for that I’ll explain in a bit (called multicoliniarity), but for now, know that the fact that there are only two categories in the model has an important effect on interpretation. Specifically, we must interpret those categories in reference to the one missing (which we call the reference category). So, in the example above, we can say the White people have a higher age of execution in reference to African Americans. Specifically, compared to African Americans, whites have an higher age by about 3.5 years. Why compared to African Americans? Because that’s the category not in the model - the reference category.
Similarly. we can see that method3Lethal injection is significant. We must interpret this relative to electrocution, which is the reference category. Specifically, those killed by lethal injection are, on average, 3.2 years older than those killed by electrocution.
Breif ANOVA Digression
One way to think of this is like an ANOVA inside a regression. We’re looking at group differences on the dependent variable (age of execution) and the factors are the same kind of factors we would see in an ANOVA. In fact, if we were to run an ANOVA we would see the same answers!
## Df Sum Sq Mean Sq F value Pr(>F)
## race3 2 4719 2359.4 27.62 1.69e-12 ***
## method3 2 1820 909.8 10.65 2.56e-05 ***
## Residuals 1437 122737 85.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Age ~ race3 + method3, data = capital)
##
## $race3
## diff lwr upr p adj
## Other-Black -0.520202 -2.573160 1.532756 0.8230989
## White-Black 3.512661 2.273606 4.751715 0.0000000
## White-Other 4.032863 2.070635 5.995091 0.0000047
##
## $method3
## diff lwr upr p adj
## Lethal Injection-Electrocution 3.154198 1.324818 4.98357870 0.0001627
## Other-Electrocution -2.074528 -7.609042 3.45998583 0.6534016
## Other-Lethal Injection -5.228726 -10.522717 0.06526481 0.0537619In the above, we can see I calculated a quick Factorial ANOVA with only main effects and (because I knew it would be significant) just rant the post-hoc tests. You can see that White and Black, and Black and Other are different, as are Lethal Injection and Electrocution. This is similar to the findings in the regression above.
So what’s the difference? Honestly, not that much! But the reality is that regression is a lot more flexible in many instances, and there are some (slight) differences under the hood (specifically in how they assign explained variance). In any event, which you use can be dictated by the question, as well as the types of variables you’re dealing with.
Back to Regression
Interestingly, there’s not much else different about interpretation here. We see that there are no other significant variables in the model, so after reporting a very poor \(R^2\), we wouldn’t have much else to say!
Multicoliniarity
So, in the explanation above, I used the word multicoliniarity. While regression is robust to violation of a lot of its assumptions, one of the ones that will trip it up is multicoliniarity.
What is it? Well, it’s really when two of your variables are very closely related. You often run into it when you measure similar things in different ways. For instance, if I measured anger and frustration, chances are they would vary quite closely with one another. If I then included them in a regression model, one of two things might happen. First, R may just kick on the variables out and I might get NAs where I would expect coefficients. Second, I might get an error from R telling me that the model can’t run.
In practice, we run into multicoliniarity quite commonly, but mostly because we don’t fully understand our measurements. It happens often with secondary datasets, and generally is most easily solved by simply dropping one of the variables.
Here’s an example of when we run into it, and what it looks like, with some today data.
thing1<-c(1,2,3,4,5,6,7,8,9,10)
thing2<-c(2,3,4,5,6,7,8,9,10,11)
thing3<-c(12,13,44,22,5,4,55,66,4,3)
toy<-data.frame(thing1, thing3, thing3)
toy_model<-lm(thing3~thing1+thing2)
summary(toy_model)##
## Call:
## lm(formula = thing3 ~ thing1 + thing2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.545 -18.647 -8.748 16.573 42.230
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.6667 17.0255 1.214 0.259
## thing1 0.3879 2.7439 0.141 0.891
## thing2 NA NA NA NA
##
## Residual standard error: 24.92 on 8 degrees of freedom
## Multiple R-squared: 0.002492, Adjusted R-squared: -0.1222
## F-statistic: 0.01998 on 1 and 8 DF, p-value: 0.8911As you can see, it kicked out the thing2 variable. Why? Because it can’t separate out the variance explained by things 1 and 2 - they’re too closely related. Think of trying to separate the circles of a perfectly overlapping Venn diagram - you couldn’t do it! Rather than try to figure out which circle is which, R simply kicks one out.
Diagnostics
Diagnostics work exactly the same when you include factors - so all the same charts still apply!
par(mfrow=c(2,2))
plot(model)
You can see some differences because we’re dealing with age in years (with a not very large range, it looks like), so we do see some different grouping than before, but generally, the charts can be interpreted largely in the same way as before.
Ordinal Data
So what happens with ordinal data? Well, in terms of what you enter into the model, not much. However, R does something really different with ordinal factors in terms of interpretation. Let’s take a look at an example using NYFS data.
nyfs<-read.csv("/Users/joshu/Dropbox/Fall 2020/CJ 625/R Stuff/NYFS_1974_clean.csv")
nyfs$relatt<-ordered(nyfs$relatt, c("never", "several times a year", "once or twice a month", "once a week", "several times/week"))
model2<-lm(alcohol~relatt+minor, data=nyfs)
summary(model2)##
## Call:
## lm(formula = alcohol ~ relatt + minor, data = nyfs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -190.40 -39.78 -19.31 13.13 437.13
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.21327 3.30140 13.695 < 2e-16 ***
## relatt.L -33.83058 8.24285 -4.104 4.65e-05 ***
## relatt.Q -7.34384 7.81167 -0.940 0.348
## relatt.C -0.58876 6.48771 -0.091 0.928
## relatt^4 -0.03819 6.40038 -0.006 0.995
## minor 0.39039 0.08942 4.366 1.50e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 71.02 on 572 degrees of freedom
## Multiple R-squared: 0.07003, Adjusted R-squared: 0.0619
## F-statistic: 8.615 on 5 and 572 DF, p-value: 6.967e-08We can see that while we still get the same number of contrasts (that’s the “categories” in the model - always k-1), they now have weird attached letters. These stand for Linear, Quadratic, Cubic, and (because we have 5 categories, and therefore 4 show in the model ^4th power). This is R trying to be helpful, but it’s a bit confusing.
The one, usually, to pay attention to is the .L - or linear. This gives you the average linear trend across the categories of the variable. So when we set the order and then used the variable in the model, R knows that it’s ordinal and tries to take advantage of that information by estimating the average effects of each level of increase across all the categories. In other words, it tries to estimate the linear change between categories. That estimate is the .L. So, in the current case, for every increase in relatt, we see on average, a -33.83 decrease in alcohol consumption.
The quadratic is a bit more challenging to interpret. It essentially is asking whether there is a quadratic function within the independent variable (relatt in this case) helps predict differences between the categories. In other words, does a quadratic function of relatt increase cause a difference in alcohol consumption. In this case, no, because it’s not significant. If it was, the interpretation would be that the effects across categories of relatt decreased as they went higher. Cubic and fourth order polynomials (what these are called) are a lot harder to interpret, and honestly, don’t matter (usually) for our purposes.
Essentially, the first and second order polynomials (.L and .Q, respectively) are the ones we need to pay attention to, and mostly to the linear polynomial (the .L), because that’s generally what we’re interested in and is most interpretable.
Let’s compare it to what we we would get with them unordered:
nyfs$relatt<-factor(nyfs$relatt, ordered = F)
model3<-lm(alcohol~relatt+minor, data=nyfs)
summary(model3)##
## Call:
## lm(formula = alcohol ~ relatt + minor, data = nyfs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -190.40 -39.78 -19.31 13.13 437.13
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.86577 6.73673 9.332 < 2e-16 ***
## relattseveral times a year -5.34572 8.55389 -0.625 0.532256
## relattonce or twice a month -13.75444 9.47352 -1.452 0.147082
## relattonce a week -25.99733 9.12602 -2.849 0.004548 **
## relattseveral times/week -43.16504 12.33347 -3.500 0.000502 ***
## minor 0.39039 0.08942 4.366 1.5e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 71.02 on 572 degrees of freedom
## Multiple R-squared: 0.07003, Adjusted R-squared: 0.0619
## F-statistic: 8.615 on 5 and 572 DF, p-value: 6.967e-08Now we see that we have similar output to when we used a nominal categorial variable (e.g. race3) in the original example. The interpretation here would be similar, drinking by people who relattonce a week is lower than those who never attend (the reference category) by about 26 drinks a year. Relattseveral times/week is about 43 drinks lower than those who never attend.
In general, it’s better to use all the information at your disposal, which means if you have an ordered factor (ordinal variable), you should generally use it as an ordered factor rather than as a nominal variable (unordered factor). However, there’s nothing technically wrong with bringing it down a level to make it somewhat easier to interpret, particularly if you’re interested in specific comparisons - called contrasts (e.g. those who go every week vs. never) - rather than the effect of the variable across all categories.
Conclusion
Using regression with different level variables is, again, largely just an extension of the earlier model. It does involve some thought as to what kind of variable you have (e.g. ordered vs. unordered), but entering it into the model is pretty straightforward, and even the more complex interpretation is largely simple when you know what you’re looking for. Additionally, elements like diagnostics do not change relative to the variables you’ve entered into the models.
With this, you’ve now reached the point where you can really engage in some interesting modeling/research. You can, for example, find out how well religious attendance and religious belief, coupled with a person’s age, sex, and alcohol use, predict felony offending. These are real questions, that you now have the tools not only to answer, but to answer well.s
Comments
Post a Comment