Unintuitive result: Multiple Linear Regression coefficients

Question

I am working on the candy-power-ranking dataset and noted something very strange. When I calculate the pairwise correlation between the variable fruity (0=without fruity taste, 1=with fruity taste) and the target variable winpercent (from 0 to 100) I get a negative correlation. But when I use a multiple linear regression (winpercent ~ all other variables) the coefficient of the fruity term ends up beeing positive and significant (p < 0.01). But that would imply that adding fruity taste would improve the winpercent values.

I know the multiple linear regression is a little bit more complicated than a standard regression, but still this result seems to be very unintuitive. Maybe nonlinear effects also play a role her.

I hope someone can enlight me on this problem.

Answer 1

This is not uncommon. It indicates that fruity correlates with another variable.

Essentially, having taken all other variables into account, fruity will correlate positively with winpercent.

For example, maybe hard correlates with fruity. Say, we have the following table:

hard | fruity | # disliked | # liked 
------------------------------------
 no  |   no   |     10     |    90
 no  |   yes  |      0     |    20
 yes |   no   |     30     |     0
 yes |   yes  |     70     |    20

This table shows that fruity correlates negatively with liked -- but if we only look at hard candies, we see that fruity is actually a good thing; and the same goes for non-hard candies!

(The above is an example of Simpson's paradox.)

This kind of thing happens whenever the explanatory variables are strongly correlated. In my example, hard correlates strongly with fruity.

You can test this in your own data. Find the OLS model excluding fruity, and then see how strongly fruity correlates with the residual. It should have a positive correlation.

Answer 2

I found a very easy example that show the effect. The overall regression (including hard and soft candy) indicates a decreasing trend. But if you only look at the soft candy you will realize that the regression shows a positive trend.

Answer 3

I must add another answer. If you like proofs, let me agree with you that a single example, and a single picture, doesn't prove anything :) Your example depends on the number of soft/hard candies (there are much more hard candies than soft candies), on their position etc. Building another example where the overall regression matches the regression for soft candies would be easy. Moreover, looking at your example one could wrongly think that 'strange' coefficients may arise only when the regressors are negatively correlated.

Let's consider a general case, the model $y=\beta_0+\beta_1x^*+\beta_2z^*+\epsilon$, where $x^*,z^*$ are $x,z$ centered and scaled. Then: $$\hat\beta=(X^TX)^{-1}X^Ty$$ where $$X^TX=\begin{bmatrix}n & 0 \\ 0 & R\end{bmatrix},\qquad R=\begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}$$ In general, $\text{var}(\hat\beta)=(X^TX)^{-1}\sigma^2$, and $\sigma^2$ is estimated by the residual mean square (e.g. see here). As to $(X^TX)^{-1}$, if the regressors are centered and scaled it is: $$(X^TX)^{-1}=\begin{bmatrix} \frac{1}{n} & 0 & 0 \\ 0 & \frac{1}{1-\rho^2} & -\frac{\rho}{1-\rho^2} \\ 0 & -\frac{\rho}{1-\rho^2} & \frac{1}{1-\rho^2} \end{bmatrix}$$ Therefore: $$\text{var}(\hat\beta_0)=\frac{\sigma^2}{n},\quad \text{var}(\hat\beta_1)=\frac{\sigma^2}{1-\rho^2},\quad \text{var}(\hat\beta_2)=\frac{\sigma^2}{1-\rho^2}$$

As you can see, the accuracy of the scaled regression coefficients depends only on the error variance $\sigma^2$ (estimated by the residual mean square) and the correlation between the two regressors. The scaled coefficients cannot be estimated accurately if the correlation is close to $\pm 1$. Standard errors much larger than the residual standard error are a sign that you can't trust in your estimates, even if they look significant.

If two regressors are strongly correlated, since one of them can explain a large share of total variation, the other one can't add anything and its estimate is fickle, unstable. As Seber & Lee (Linear Regression Analysis, John Wiley & Sons, 2003, my source) say:

"Intuitively, when the data are well spread over the $(x,z)$ plane, the fitted regression plane is well supported by the data. When the correlation is high, and $x$ and $z$ are almost linearly dependent, the regression plane is supported by a narrow ridge of points, and is consequently unstable, with a small change in the data resulting in a big change in the fitted plane" (p. 252)

Therefore, the simplest example is whatever linear regression with some "almost linearly dependent" regressors.

Answer 4

I agree with Robby the Belgian and I only wish to add an example.[1]

y = number of deaths in London from Dec 1 to Dec 15, 1952
x = smog (mg/m3)
z = sulfur dioxide (parts per million)

> airpoll
     y    x    z
1  112 0.30 0.09
2  140 0.49 0.16
3  143 0.61 0.22
4  120 0.49 0.14
5  196 2.64 0.75
6  294 3.45 0.86
7  513 4.46 1.34
8  518 4.46 1.34
9  430 1.22 0.47
10 274 1.22 0.47
11 255 0.32 0.22
12 236 0.29 0.23
13 256 0.50 0.26
14 222 0.32 0.16
15 213 0.32 0.16

x and z correlate positively with y:

> attach(airpoll)
> cor(y,x)
[1] 0.7559842
> cor(y,z)
[1] 0.8308853

but when you fit the linear model $y=\beta_0+\beta_1x+\beta_2z+\epsilon$, the smog coefficient is negative:

> fit <- lm(y ~ x + z)
> coef(fit)
(Intercept)           x           z 
    89.5108   -220.3244   1051.8165 

Why? Because smog and sulfur dioxide are strongly correlated:

> cor(x,z)
[1] 0.9875641

QED :)

---

[1] This is an old dataset that I found on my hard disk, and I can't remember the source. Sorry.