Normality of Residuals

Question

Say I have a multilinear regression $Y= a1X1+a2X2+a3X3$ ( having more than one variable; not necessarily having 3 independent variables)

We then want the residuals :{$ Y_i-(a1X_{1i}+a2X_{2i}+a3X_{3i} ) $}

to be normally-distributed (with mean =0, as pointed in the answer by Enzo). Question: If the $Y_i$ are normally-distributed, does it follow that each of the residuals {$Y-a_1'$}, {$Y-a_2'X_2$}, {$Y-a_3'X_3$} (i.e., we regress Y against each variable separately ) is also normally -distributed?

Also, does this depend on whether :{$ Y_i-(a1X_{1i}+a2X_{2i}+a3X_{3i} ) $} is jointly normal? Thanks in Advance.

Answer 1

(One step back first) Typically, the assumptions underlying a linear regression model $$y_i = x_i^T\beta + e_i,\,\,\, i=1,\dots,n$$ are:

1. The errors $e_i$ are i.i.d. with Normal distribution with mean zero and variance $\sigma^2$.
2. The covariates are either a sequence of deterministic vectors or they come from a joint distribution such that for large enough $n$ the matrix $X^TX$ is positive definite, where $X$ is the design matrix.
3. $x_i \bot e_i$, the covariates and the errors are independent.

Of course, there are all sorts of generalizations of these assumptions (e.g. heteroscedasticity).

Suppose that you remove some covariates and keep $z_i$ covariates, then $y_i-z_i^T\beta_z$ are not necessarily normal since $e_i = y_i-x_i^T\beta \neq y_i-z_i^T\beta_z$, and consequently nothing guarantees the normality of the residuals under the smaller model.

In practice, if you fit a model, and the residuals look normal, this does not imply that under a smaller model the residuals will also look normal. Have a look at the following example in R for instance:

# Simulated data
ns = 1000 # sample size
X = cbind(1,rgamma(ns,5,5),rgamma(ns,5,5,)) # design matrix
e = rnorm(ns,0,0.5) # errors
beta = c(1,2,3) # true regression parameters 
y = X%*%beta + e  # simulating the responses 

# fitting the model
lmr = lm(y~-1+X) 

# residuals
res = lmr$residuals

# histogram and normality test: nicely normal looking
hist(res)
shapiro.test(res)

# Using only two covariates (one of them is the intercept)
Z = X[,1:2]

# Fitting the smaller model
lmrz = lm(y~-1+Z) 

# residuals
resz = lmrz$residuals

# histogram and normality test: not normal looking and failing the test
hist(resz)
shapiro.test(resz)

Answer 2

Consider the following case (where $N(\mu,\sigma)$ is a normal distribution, and $I(s)$ is 1 if $s$ is true and 0 if $s$ is false):

$W=N(0,100)$

$X_1=W*I(W>0)$

$X_2=W*I(W<0)$

And suppose we have $Y=X_1+X_2+N(0,1)$. It should be trivial to establish that $Y$ is normally distributed.

If we regress $Y$ on $X_1$ and $X_2$ jointly, what will the resulting models and residuals look like? What about individually?

Answer 3

No, because every time you build a model, you must test the normal distribution and the zero mean of the residuals. So, in your case, if you drop a variable of the first model you must re-test the normal distribution of the residuals (and the zero mean of course).