Why does checking normality of residuals give a different result than checking bivariate normality of the two variables?

Question

I am checking the conditions for hypothesis testing a Pearson correlation as significant or not, and also checking the residuals normality conditions for OLS. Why are the following methods giving different results (why are the variables not bivariate normal but the residuals of an OLS model are normal?)? Also, why does changing the regression from y~x to x~y alter the p value so much for normality of the residuals?

I have included a reproducible example below.

library(tidyverse)
library(mvnormtest)

x <- c(21,21,20,16,15,20,16,16,23,21,9,14,18,6,35,39,29,26,32,7)
y <- c(226,223,233,209,207,251,268,255,357,401,323,302,422,245,389,409,352,374,389,418)

data <- tibble(x,y)

# Histogram of residuals for visual inspection
data %>% 
  mutate(resid = residuals(lm(x ~ y))) %>% 
  ggplot(aes(x = resid)) + 
  geom_histogram()

# Univariate Shapiro Wilks of residuals
data %>% 
  mutate(resid = residuals(lm(x ~ y))) %>% 
  select(resid) %>% 
  t() %>%
  shapiro.test()

# Univariate Shapiro Wilks of residuals (switch x and y order)
data %>% 
  mutate(resid = residuals(lm(y ~ x))) %>% 
  select(resid) %>% 
  t() %>%
  shapiro.test()

# Bivariate Shapiro Wilks
data %>%
  t() %>%
  mshapiro.test()

I was also confused because I thought that hypothesis testing a Pearson correlation has similar assumptions to fitting OLS model, in that the variables should be bivariate normal, but is this mistaken? Do I need to check normality individually for the x and y variables?

Sources:

1. Bivariate normality is a necessary condition for testing Pearson correlation (but alternatively, the univariate normality of the two variables can be separately checked (?)): https://statistics.laerd.com/spss-tutorials/pearsons-product-moment-correlation-using-spss-statistics.php

2. Both variables must be individually normal (seems to contradict the previous source (1)): https://toptipbio.com/what-is-pearson-correlation/

3. Holding constant x, then y must be normally distributed (this assumption seems to be the same as (1)): https://courses.lumenlearning.com/introstats1/chapter/testing-the-significance-of-the-correlation-coefficient/

If anyone has an authoritative source that explains this, I'd much appreciate it.

Answer 1

confused because I thought that hypothesis testing a Pearson correlation has similar assumptions to fitting OLS model,

As simple regression, sure, and so it, too, is fairly insensitive to normality of errors in large samples

Bivariate normality and marginal normality are not the same and neither is strictly required for testing a Pearson correlation. (Bivariate normality is sufficient but not necessary. Marginal normality on its own is neither sufficient nor necessary)

Also, why does changing the regression from y~x to x~y alter the p value so much for normality of the residuals?

Because you're looking at the errors from a different line and conditioning on a different variable. They may be entirely different.

Testing normality is not necessarily particularly helpful and doesn't answer the question you need answered.

---

To return to the title question:

Why does checking normality of residuals give a different result than checking bivariate normality of the two variables?

They're different things. Specifically, while bivariate normality will result in normally distributed residuals from a regression, you can get normally distributed residuals while neither variable is marginally normal (and neither are the variables jointly normal).

Consider, for example, a simple regression where x is either 0 or 1 and the model $Y_i = \beta_0 + \beta_1 x_i +\epsilon_i$, $i=1,2,...,n$, where the $\epsilon_i$ are i.i.d. $N(0,\sigma^2)$. In this case, under $\rho(x,Y)=0$, which is equivalent to $\beta_1=0$, the sample Pearson correlation $r$ has the usual null distribution, and so the test statistic $t=\frac{r\sqrt{n-2}}{1-r^2}$ has a t-distribution with $n-2$ d.f. when $H_0$ is true.

Here's an illustration (not proof, but proofs of the necessary results can be found in pretty much any decent regression text):

The plot on the left is under H1 (the variables are correlated with 1000 values of x and y). Clearly they won't be bivariate normal whatever the slope, because one variable is binary. Under H0, I simulated 1000 sets of 20 values of x and y - i.e. smaller samples and population slope 0 - and computed the correlation for each set, then calculated the t-test statistic, and then finally transformed it by its own theoretical cdf (the t cdf with 18 df). The result looks like a random sample from a uniform distribution, exactly as it should; suggesting that the properties of significance levels and p-values will be what they are claimed to be under H0 -- the test works exactly as advertized even though (X,Y) was not bivariate normal, it just had normal errors on the conditional distribution of one of the variables.

We can use a similar approach to investigate sensitivity to non-normality under whatever your sample size and pattern of x's is. That is what I'd suggest is a good starting point for deciding how sensitive the test is to that assumption, by exploring plausible possible assumptions (and then seeing how far you have to push them before the test's properties become too far from the nominal properties for your own requirements.

---

Regarding your sources:

bivariate normality is a necessary condition for testing Pearson correlation (but alternatively, the univariate normality of the two variables can be separately checked)

Univariate normality of the margins is not sufficient for bivariate normality. Bivariate normality establishes both linearity of conditional expectation, and both conditional and marginal normality. As mentioned before bivariate normality, while sufficient, is not required for the significance levels to be correct, or the test to be consistent, or to get good power.

Both variables must be individually normal

Neither variable need be individually normal; In the simple regression example I gave, the x's were Bernoulli$(\frac12)$ and the y's were conditionally normal with constant variance, but the marginal distribution of y will only be normal when the population correlation is 0; under H1 it can be bimodal. Under other patterns of x-values one can get left skewed, right skewed, heavy-tailed or light tailed marginal distributions of y (when the correlation is non-zero).

Holding constant x, then y must be normally distributed

Correct, though additional things are required. However, note that holding x constant, you might only have a single y value there -- consider replacing my earlier example with x having a beta(0.1,0.1) distribution, say.

this assumption seems to be the same as (1)

It is not. It is the conditional normality assumption from regression. To be clear, that's the assumption that is made in deriving the t-distribution of the test statistic when H0 is true. The test is not especially sensitive to that assumption in moderately large samples, though.

---