How can we get different mean and median values for residuals when working with the same dataset?

Question

Following this R instruction

> fit <- lm(spending ~ sex + status + income + verbal, data=spending)

I would like to calculate the mean and median of the residuals. Both my friend and I get different answers for the mean for the same data.

> mean(resid(fit))
[1] -3.065293e-17

while with the other model, results are:

> fit1<- lm(spending ~ status + income + sex + verbal, data=spending) 
> mean(resid(fit1))
[1] 4.064605e-16

Why did we get the same median but different mean if we are using the same data set?

Answer 1

By definition, the mean of the residuals of both models should be zero. Deviations are due to rounding issues.

Here's an example:

> x <- rnorm(10)
> y <- rnorm(10)
> z <- rnorm(10)
> res1 <- residuals(lm(y ~ x+z))
> res2 <- residuals(lm(y ~ z+x))

Note that both models include the same variables, but the order is different.

Look at the mean:

> mean(res1)
[1] 3.8849e-17
> mean(res2)
[1] 1.525337e-17
> median(res1)
[1] 0.06405684
> median(res2)
[1] 0.06405684

The mean is not identical, but the median is.

The residuals are nearly identical (as expected):

> all.equal(sort(res1), sort(res2))
 [1] TRUE

Check ?all.equal.

But they are not completely identical:

> identical(as.numeric(sort(res1)), as.numeric(sort(res2)))
[1] FALSE

The deviance from zero is due to rounding issues in R. Hence, the difference in the order of model predictors may cause this kind of differences.