Conflicting p-values of regressors in simple regression vs multiple regression? [duplicate]

Question

Consider the results of the following code.

x_cm = 3*rnorm(100,50,3)
x_in = x_cm/2.5 + rnorm(100,0,0.1);
y = x_cm + rnorm(100,0,4)

mod <- lm(y ~ x_cm); summary(mod)
mod <- lm(y ~ x_in); summary(mod)
mod <- lm(y ~ x_in + x_cm); summary(mod)

The p-values given for $x_\text{cm}$ and $x_\text{in}$ are extremely small when we look at $y$ regressed on each of those predictors, individually, e.g. they are <2e-16. This indicates that we have a statistically significant result that the associated regression coefficients $\beta_\text{cm}$ and $\beta_\text{in}$ are non-zero.

However, when we perform a regression of $y$ on both of these variables simultaneously, and check their individual p-values in the summary, these p-values are now very large: 0.547 and 0.975 (although the p-value for the F-test is very small ?). This indicates that we have a statistically insignificant result that the regression coefficients $\beta_\text{cm}$ are $\beta_\text{in}$ are non-zero.

These two results regarding the hypotheses that the regression coefficients are non-zero are in direct conflict with each other.

Why is that when we look at the individual p-values after performing a multiple regression, that we have statistically insignificant results? Why are the individual p-values in multiple regression not the same as p-values obtained by performing two separate simple regressions?

While I am already aware of issues with p-values in statistical inference, e.g. p-hacking and data snooping, the phenomenon observed in this simple example makes it seems like we cannot trust p-values at all when it comes to multiple regression, is that in fact the case?

It seems like the results from the multiple regression case tell us that $x_\text{cm}$ and $x_\text{in}$ are not useful in our model, due to the insignificant p-value and thus our conclusion should be that these variables can't be used to predict $y$. But of course they can be used to predict $y$, since $y$ was directly generated by $x_\text{cm}$ and $x_\text{in}$ is highly correlated with $x_\text{cm}$ so it could also be used to predict $y$.

Finally, when I perform a multiple regression in general and I get large p-values for some coefficients, how should I interpret this situation and what should be my next steps?

Answer 1

It is a known issue in regression with correlated predictors that the standard errors on the coefficients get inflated, resulting in larger p-value and decreased power to reject a null hypothesis that a particular parameter is zero. You still can wind up with the Gauss-Markov theorem in effect to give you that unbiased estimator, but inferential ideas (p-values, standard errors, confidence intervals) are larger than they would be with uncorrelated predictors.

Note that this only happens when the predictor you're examining is correlated with the other predictors, meaning that if you have three predictors, two that are correlated with each other and a third that is uncorrelated with the others, inference on the uncorrelated predictor's parameter goes the same as it would if the three variables were uncorrelated. I address this idea in my self-answer here.

The F-test p-value is very small despite the individual parameter p-values being large because the F-test is comparing your model to the intercept-only model. Your conclusion would be that the included predictors do influence the response variable, even if you can't pin down which variable is doing the influencing. In other words, while you can't reject either of $H_0: \beta_{cm} = 0$ or $H_0: \beta_{in} = 0$, you can reject $H_0: \beta_{cm} = \beta_{in} = 0$, which is what R is testing in that F-test.

Answer 2

This is the problem of multicollinearity. When one variable lies close to the space spanned by the others then the p value will be split among them diluating its value on any one. The variance inflation factor (VIF) can be used to diagnose this. Typically a value above 5 is regarded as problematic and that is the case here.

library(faraway)

vif(mod1)
## x_cm 
##    1 

vif(mod2)
## x_in 
##    1 

vif(mod3)
##     x_in     x_cm 
## 1155.331 1155.331 

We can also see this by comparing mod2 and mod3. mod3 is not significantly different from mod2 (p = 0.7242) so adding the extra variable seems pointless.

anova(mod2, mod3)
## Analysis of Variance Table
## 
## Model 1: y ~ x_in
## Model 2: y ~ x_in + x_cm
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     98 1406.3                           
## 2     97 1404.5  1    1.8126 0.1252 0.7242

If you wish to pursue diagnostics further the mctest R package contains numerous collinearity diagnostics. See this article which discusses it.

Note

The input in reproducible form:

set.seed(123)

x_cm = 3*rnorm(100,50,3)
x_in = x_cm/2.5 + rnorm(100,0,0.1);
y = x_cm + rnorm(100,0,4)

mod1 <- lm(y ~ x_cm); summary(mod1)
mod2 <- lm(y ~ x_in); summary(mod2)
mod3 <- lm(y ~ x_in + x_cm); summary(mod3)