How to break regression

Regression models are a cornerstone of modern social science. They’re at the heart of efforts to estimate causal relationships between variables in a multivariate environment and are the basic building blocks of many machine learning models. Yet social scientists can run into a lot of situations where regression models break.

Famed social psychologist Richard Nisbett recently argued that regression analysis is so misused and misunderstood that analyses based on multiple regression “are often somewhere between meaningless and quite damaging.” (He was mainly talking about cases in which researchers publish correlational results that are covered in the media as causal statements about the world.)

Below, I’ll walk through some of the potential pitfalls you might encounter when you fire up your favorite statistical software package and run regressions. Specifically, I’ll be using simulation in R as an educational tool to help you better understand the ways in which regressions can break.

Using simulations to unpack regression

The idea of using R simulations to help understand regression models was inspired by Ben Ogorek’s post on regression confounders and collider bias.

The great thing about using simulation in this way is that you control the world that generates your data. The code I’ll introduce below represents the true data-generating process, since I’m using R’s random number generators to simulate the data. In real life, of course, we only have the data we observe, and we don’t really know how the data-generating process works unless we have a solid theory (like Newtonian physics or evolution) where the system of relevant variables and causal relationships is well understood and to which there is really no analogous phenomenon in social science.

What I’ll do here is create a dataset based on two random standard normal variables by simulating them using the rnorm() function, which draws random values from a normal distribution with mean 0 and standard deviation 1, unless you specify otherwise. I’ll create a functional relationship between y and x such that a 1 unit increase in x will be associated with a .4 unit increase in y.

1 # make the code reproducible by setting a random number seed
2 set.seed(100)
3
4 # When everything works:
5 N <- 1000
6 x <- rnorm(N)
7 y <- .4 * x + rnorm(N)
8 hist(x)
9 hist(y)
10
11 # Now estimate our model:
12 summary(lm(y ~ x))
13
14 Call:
15 lm(formula = y ~ x)
16 Residuals:
17    Min      1Q  Median      3Q     Max 
18 -3.0348 -0.7013  0.0085  0.6212  3.1688 
19 Coefficients:
20            		Estimate Std. Error t value Pr(>|t|)    
21 (Intercept) 	0.003921   0.031039   0.126    0.899    
22 x           	0.413415   0.030129  13.722   <2e-16 ***
23 ---
24 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
25 Residual standard error: 0.9814 on 998 degrees of freedom
26 Multiple R-squared:  0.1587,	Adjusted R-squared:  0.1579 
27 F-statistic: 188.3 on 1 and 998 DF,  p-value: < 2.2e-16
28
29 # Plot it
30 library(ggplot2)
31 qplot(x, y) +
32  geom_smooth(method='lm') +
33  theme_bw() +
34  ggtitle("The Perfect Regression")

Notice that the model estimates the functional relationship between x and y that I simulated quite well. The plot looks like this:

What about omitted variables? Our machinery actually still works if there is another factor causing y, as long as it is uncorrelated with x.

The dreaded omitted variable bias

Omitted variable bias (OVB) is much feared, and judging by the top internet search results, not well understood. Some top sources say it occurs when “an important” variable is missing or when a variable that “is correlated” with both x and y is missing. I even found a university econometrics course that defined OVB this way.

But neither of those definitions are quite right. OVB occurs when a variable that causes y is missing from the model (and is correlated with x). Let’s call that variable w. Because w is in play when we consider the causal relationship between x and y, it’s often referred to as “endogenous” or a “confounding variable.”

The example below first demonstrates that w, our confounding variable, will bias our results if we fail to include it in our model. The next two examples are essentially a re-telling of the post I mentioned above on collider bias, but emphasizing slightly different points.

1 w <- rnorm(N)
2 x <- .5 * w + rnorm(N)
3 y <- .4 * x + .3 * w + rnorm(N)
4
5 m1 <- lm(y ~ x)
6 summary (m1) # Omitted variable bias
7
8 Call:
9 lm(formula = y ~ x)
10
11 Residuals:
12    Min      1Q  Median      3Q     Max 
13 -3.2190 -0.7025  0.0314  0.7120  3.1158 
14
15 Coefficients:
16            Estimate Std. Error t value Pr(>|t|)    
17 (Intercept)  0.01126    0.03310    0.34    0.734    
18 x            0.50179    0.03049   16.46   <2e-16 ***
19 ---
20 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
21
22 Residual standard error: 1.046 on 998 degrees of freedom
23 Multiple R-squared:  0.2135,	Adjusted R-squared:  0.2127 
24 F-statistic: 270.9 on 1 and 998 DF,  p-value: < 2.2e-16

There it is: classic omitted variable bias. We only observed x, and the influence of the omitted variable w was attributed to x in our model. If you re-rerun the regression with w in the model, you no longer get biased estimates.

1 m2 <- lm(y ~ x + w)
2 summary (m2) # No omitted variable bias after conditioning on w
3
4 Call:
5 lm(formula = y ~ x + w)
6
7 Residuals:
8    Min      1Q  Median      3Q     Max 
9 -3.2748 -0.6632 -0.0001  0.6933  2.9664 
10
11 Coefficients:
12            Estimate Std. Error t value Pr(>|t|)    
13 (Intercept)  0.02841    0.03141   0.905    0.366    
14 x            0.40627    0.03132  12.973   <2e-16 ***
15 w            0.32344    0.03439   9.405   <2e-16 ***
16 ---
17 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
18
19 Residual standard error: 0.9927 on 997 degrees of freedom
20 Multiple R-squared:  0.3024,	Adjusted R-squared:  0.301 
21 F-statistic: 216.1 on 2 and 997 DF,  p-value: < 2.2e-16
22

Note that the regression errors, also known as residuals, are correlated with w:

1 cor(w,m1$residuals)
2 [1] 0.2597859

Now, recall above that I wrote that it’s wrong to say that OVB occurs when our omitted variable is correlated with both x and y. And yet w, x and w and y are all correlated in this first example:

1 cormatrix <- cor(as.matrix(data.frame(x,y,w)))
2 round(matrix, 2)
3
4     x    y    w
5 x 1.00 0.49 0.41
6 y 0.49 1.00 0.43
7 w 0.41 0.43 1.00

So why can’t we just say that OVB occurs when our omitted variable is correlated with both x and y? As the next example will show, correlation isn’t enough — w needs to cause both x and y. We can easily imagine a case in which we don’t have causality but we still see this kind of correlation — when x and y both cause w.

Let’s make this a little more concrete. Suppose we care about the effect of news media consumption (x) on voter turnout (y). One factor that some researchers think may cause both news media consumption and turnout is political interest (w). If we only measure media consumption and voter turnout, political interest is likely to confound our estimates.

But another school of thought from social psychology — along the lines of self-perception theory and cognitive dissonance — suggests that the causality could be reversed: Voting behavior might be mostly determined by other factors, and casting a ballot might prompt us to be more interested in political developments in the future. Similarly, watching the news might prompt us to become moreinterested in politics. Let’s suppose that second school of thought is right. If so, our simulated data will look like this:

1 media_consumption_x <- rnorm(N)
2 voter_turnout_y <- .1 * media_consumption_x + rnorm(N)
3
4 # Political interest increases after consuming media and participating, and, 
5 # in this hypothetical world, does *not* increase media consuption or participation
6 political_interest_w <- 1.2 * media_consumption_x + .6 * voter_turnout_y + rnorm(N)
7
8 cormat <- cor(as.matrix(data.frame(media_consumption_x, voter_turnout_y, political_interest_w)))
9 round(cormat, 2)
10
11                     media_consumption_x voter_turnout_y 
12 political_interest_w
media_consumption_x                 1.00            0.11                 0.70
13 voter_turnout_y                     0.11            1.00                 0.46
14 political_interest_w                0.70            0.46                 1.00

As you can see, all factors are again correlated with each other. But this time, if we only include x (media consumption) and y (turnout) in the equation, we get the correct estimate:

1 summary(lm(voter_turnout_y ~ media_consumption_x))
2
3 Call:
4 lm(formula = voter_turnout_y ~ media_consumption_x)
5
6 Residuals:
7    Min      1Q  Median      3Q     Max 
8 -2.8460 -0.6972 -0.0076  0.6702  3.3925 
9
10 Coefficients:
11                    Estimate Std. Error t value Pr(>|t|)    
12 (Intercept)         -0.01202    0.03217  -0.374 0.708839    
13 media_consumption_x  0.11719    0.03321   3.529 0.000436 ***
14 ---
15 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
16
17 Residual standard error: 1.014 on 998 degrees of freedom
18 Multiple R-squared:  0.01233,	Adjusted R-squared:  0.01134 
19 F-statistic: 12.46 on 1 and 998 DF,  p-value: 0.0004359

What makes defining omitted variable bias based on correlation so dangerous is that if we now include w (political interest), we will get a different kind of bias — what’s called collider bias or endogenous selection bias.

1 summary(lm(voter_turnout_y ~ media_consumption_x +  political_interest_w))
2
3 Call:
4 lm(formula = voter_turnout_y ~ media_consumption_x + political_interest_w)
5
6 Residuals:
7    Min      1Q  Median      3Q     Max 
8 -2.1569 -0.5981 -0.0129  0.5701  2.8356 
9
10 Coefficients:
11                      Estimate Std. Error t value Pr(>|t|)    
12 (Intercept)           0.003155   0.027098   0.116    0.907    
13 media_consumption_x  -0.437084   0.039102 -11.178   <2e-16 ***
14 political_interest_w  0.444571   0.021928  20.274   <2e-16 ***
15 ---
16 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
17 
18 Residual standard error: 0.854 on 997 degrees of freedom
19 Multiple R-squared:  0.3007,	Adjusted R-squared:  0.2993 
20 F-statistic: 214.3 on 2 and 997 DF,  p-value: < 2.2e-16

Simpson’s paradox

Simpson’s paradox often occurs in social science (and medicine, too) when you pool data instead of conditioning it on group membership (i.e., adding it as a factor in your regression model).

Suppose that, all other things being equal, consuming media causes a slight shift in policy preferences toward the left. But, on average, Republicans consume more news than non-Republicans. And we know that generally Republicans have much more right-leaning preferences.

If we just measure media consumption and policy preferences without including Republicans in the model, we’ll actually estimate that the effect goes in the direction opposite of the true causal effect.

1 N <- 1000
2
3 # Let's say that 40% of people in this population are Republicans
4 republican <- rbinom(N, 1, .4)
5
6 # And they consume more media
7 media_consumption <- .75 * republican + rnorm(N)
8
9 # Consuming more media causes a slight leftward shift in policy
10 # preferences, and Republicans have more right-leaning preferences
11 policy_prefs <- -.2 * media_consumption + 2 * republican + rnorm(N)
12
13 # for easier plotting later
14 df <- data.frame(media_consumption, policy_prefs, republican)
15 df$republican = factor(c("non-republican", "republican")[df$republican + 1])
16
17 # If we don't condition on being Republican, we'll actually estimate
18 # that the effect goes in the *opposite* direction
19 summary(lm(policy_prefs ~ media_consumption))
20
21
22 Call:
23 lm(formula = policy_prefs ~ media_consumption)
24
25 Residuals:
26    Min      1Q  Median      3Q     Max 
27 -3.6108 -0.9559 -0.0198  0.9257  3.9537 
28
29 Coefficients:
30                  Estimate Std. Error t value Pr(>|t|)    
31 (Intercept)        0.68923    0.04323   15.94  < 2e-16 ***
32 media_consumption  0.15269    0.03966    3.85 0.000126 ***
33 ---
34 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
35
36 Residual standard error: 1.317 on 998 degrees of freedom
37 Multiple R-squared:  0.01463,	Adjusted R-squared:  0.01365 
38 F-statistic: 14.82 on 1 and 998 DF,  p-value: 0.0001257
39
40 # Naive plot
41 qplot(media_consumption, policy_prefs) +
42   geom_smooth(method='lm') +
43   theme_bw() +
44   ggtitle("Naive estimate (Simpson's Paradox)")

The estimate goes in the opposite direction of the true effect! Here’s what the plot looks like:

To resolve this paradox, we need to add a factor in the model that indicates whether or not a respondent is a Republican. Adding that factor lets us estimate separate slopes for Republicans and non-Republicans. Note that this is not like estimating an interaction term, where two explanatory variables are multiplied together. It’s not that the slopes are different, we just need to estimate separate ones for Republicans and non-Republicans.

1 # Condition on being a Republican to get the right estimates
2 summary(lm(policy_prefs ~ media_consumption + republican))
3
4 Call:
5 lm(formula = policy_prefs ~ media_consumption + republican)
6
7 Residuals:
8    Min      1Q  Median      3Q     Max 
9 -3.5518 -0.6678 -0.0186  0.6562  3.3009 
10
11 Coefficients:
12                  Estimate Std. Error t value Pr(>|t|)    
13 (Intercept)        0.05335    0.03904   1.366    0.172    
14 media_consumption -0.13615    0.03111  -4.376 1.34e-05 ***
15 republican         1.93049    0.06758  28.565  < 2e-16 ***
16 ---
17 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 
18
19 Residual standard error: 0.9774 on 997 degrees of freedom
20 Multiple R-squared:  0.4581,	Adjusted R-squared:  0.457 
21 F-statistic: 421.4 on 2 and 997 DF,  p-value: < 2.2e-16
22 
23 # Conditioning on being Republican
24 qplot(media_consumption, policy_prefs, data=df, colour = republican) +
25   scale_color_manual(values = c("blue","red")) +
26   geom_smooth(method='lm') +
27  theme_bw() +
28  ggtitle("Conditioning on being a Republican (Simpson's Paradox)")

Here’s what the plot looks like:

Correlated errors

Another cardinal sin — and one that we should worry a lot about because it often arises from social desirability bias in survey responses — is the phenomenon of correlated errors. This example is inspired by Vavreck (2007).

Here, self-reported turnout and media consumption are caused by a combination of social desirability bias and true turnout and true consumption, respectively:

1 N <- 1000
2
3 # The "Truth"
4 true_media_consumption <- rnorm(N)
5 true_vote <- .1 * media_consumption + rnorm(N)
6
7 # social desirability bias
8 social_desirability <- rnorm(N)
9 #what we actually observe from self reports:
10 self_report_media_consumption <- true_media_consumption + social_desirability
11 self_report_vote <- true_vote + social_desirability

Let’s compare the estimated effect sizes of the self-reported data and the “true” data:

1 # Self reports
2 summary(lm(self_report_vote ~ self_report_media_consumption))
3
4 Call:
5 lm(formula = self_report_vote ~ self_report_media_consumption)
6
7 Residuals:
8    Min      1Q  Median      3Q     Max 
9 -3.9604 -0.7766  0.0142  0.8465  4.1811 
10
11 Coefficients:
12                              Estimate Std. Error t value 
   Pr(>|t|)    
13 (Intercept)                    0.02020    0.03951   0.511    0.609    
14 self_report_media_consumption  0.54605    0.02716  20.102   <2e-16 ***
15 ---
16 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
17
18 Residual standard error: 1.248 on 998 degrees of freedom
19 Multiple R-squared:  0.2882,	Adjusted R-squared:  0.2875 
20 F-statistic: 404.1 on 1 and 998 DF,  p-value: < 2.2e-16
21
22 # "Truth"
23 summary(lm(true_vote ~ true_media_consumption))
24
25 Call:
26 lm(formula = true_vote ~ true_media_consumption)
27
28 Residuals:
29     Min      1Q  Median      3Q     Max 
30 -3.5814 -0.6677 -0.0077  0.6829  3.4799 
31
32 Coefficients:
33                       Estimate Std. Error t value Pr(>|t|)
34 (Intercept)             0.01372    0.03217   0.426    0.670
35 true_media_consumption  0.01313    0.03245   0.404    0.686
36
37 Residual standard error: 1.017 on 998 degrees of freedom
38 Multiple R-squared:  0.0001639,	Adjusted R-squared:  -0.000838 
39 F-statistic: 0.1636 on 1 and 998 DF,  p-value: 0.686

The self-reported data is biased toward over-estimating the effect size, a very dangerous problem. How could we fix this? Well, one way is to actually measure social desirability and include it in the model:

1 summary(lm(self_report_vote ~ self_report_media_consumption + social_desirability))
2
3 Call:
4 lm(formula = self_report_vote ~ self_report_media_consumption + 
5    social_desirability)
6
7 Residuals:
8     Min      1Q  Median      3Q     Max 
9 -3.6042 -0.6774 -0.0127  0.6899  3.4470 
10
11 Coefficients:
12                               Estimate Std. Error t value Pr(>|t|)    
13 (Intercept)                    0.01208    0.03220   0.375    0.708    
14 self_report_media_consumption  0.01220    0.03246   0.376    0.707    
15 social_desirability            1.02245    0.04547  22.487   <2e-16 ***
16 ---
17 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
18
19 Residual standard error: 1.017 on 997 degrees of freedom
20 Multiple R-squared:  0.5277,	Adjusted R-squared:  0.5268 
21 F-statistic:   557 on 2 and 997 DF,  p-value: < 2.2e-16

Note that this while most people think about social desirability as being a problem related to measurement error, it is essentially the same problem as omitted variable bias, as described above.

It’s important to remember that omitted variable bias and correlated errors are just two potential problems with regression analysis. Regression models are also not immune to issues associated with low levels of statistical power, the failure to account for the influence of extreme values, and heteroskedasticity, among others. But by simulating the data-generating process, researchers can get a good sense of some of the more common ways in which statistical models might depart from reality.