# Regression: distribution of residuals at each point

I'm stydying In the Wikipedia page on Errors and Residuals, it is written that

However, because of the behavior of the process of regression, the distributions of residuals at different data points (of the input variable) may vary even if the errors themselves are identically distributed. Concretely, in a linear regression where the errors are identically distributed, the variability of residuals of inputs in the middle of the domain will be higher than the variability of residuals at the ends of the domain[citation needed]: linear regressions fit endpoints better than the middle. This is also reflected in the influence functions of various data points on the regression coefficients: endpoints have more influence.

I don't get what are the "distribution of residuals at different data points": if I have one residual per point, how would I have a different distribution at each point?

• If you throw a dice once, than you get one result. Still that result ist taken from a distribution. even if there is only one value sampled from from that distribution. Every single residual stems from some distribution of possible residuals that could have emerged. Does that answer the question or was the question broader? – Bernhard Mar 2 '17 at 10:10
• @Bernhard Right, then why is it that these distributions change at each data point? – mar tin Mar 2 '17 at 10:13
• I was not aware of that and as your citation states: A citation is needed. As I understand it, points on the far left or far right of the x-axis will influence the regression line more than those in the middle. See the concept of leverage ( en.wikipedia.org/wiki/Leverage_(statistics) ). If a point has more power to draw the regression line towards it, there will be smaller residuals. – Bernhard Mar 2 '17 at 10:26
• Ok, I wish there were a citation indeed. I do understand this you say intuitively, and know about the leverages/hat matrix, but I'd like a demonstration of influence varying with the x-axis. – mar tin Mar 2 '17 at 11:25
• I needed a lot of code an more then 600 characters, so I continued this not in the comments but as an answer. See simulation below. – Bernhard Mar 2 '17 at 12:03

Ok, I wish there were a citation indeed. I do understand this you say intuitively, [...], but I'd like a demonstration of influence varying with the x-axis

For demonstration I wrote the following short R simulation function:

n = 20
test <- function(){
xleft   <- 3*0 + rnorm(n)
xmiddle <- 3*1 + rnorm(n)
xright   <-3*2 + rnorm(n)
model <- lm(c(xleft, xmiddle, xright) ~ rep(0:2, each=n))
return(mean(model$residuals[1:n])<mean(model$residuals[n+1:2*n]))
}


It draws a random sample with three distinct x values (0, 1, 2) and n normally distributed random values following $y_i \tilde{} Norm(3\times x_i, 1)$. It then performs a linear regression (lm) and returns, whether the mean of the residuals at $x=0$ are smaller than the mean of the residuals at $x=1$. The errors are the same at each point on the x axis and according to the statement you cited from wikipedia, the mean residuals should be smaller at the end ($x=0$) then in the middle ($x=1$). Now we can test that by replicating 100,000 times and counting, how often this is the case:

sum(replicate(10000, test()))


With $n = 20$ I got the following results:

> set.seed(111)
> sum(replicate(10000, test()))
[1] 4931
> sum(replicate(10000, test()))
[1] 5024
> sum(replicate(10000, test()))
[1] 5106
> sum(replicate(10000, test()))
[1] 4947
> sum(replicate(10000, test()))
[1] 5013
> sum(replicate(10000, test()))
[1] 4963
> sum(replicate(10000, test()))
[1] 5031
> sum(replicate(10000, test()))
[1] 5073
> sum(replicate(10000, test()))
[1] 4970
> sum(replicate(10000, test()))
[1] 5015


So always around 5000 times out of 10000 this is the case. The effect is either small or it depends on further assumptions about the particular regression problem, which are not mentioned in the wikipedia article. I am not shure, whether further investigation is worthwhile without further references.

E D I T: In the above code, I made a stupid coding error. The last line in my test function needs two more pairs of parentheses to bei correct:

n = 20
test <- function(){
xleft   <- 3*0 + rnorm(n)
xmiddle <- 3*1 + rnorm(n)
xright   <-3*2 + rnorm(n)
model <- lm(c(xleft, xmiddle, xright) ~ rep(0:2, each=n))
return(mean(model$residuals[1:n])<mean(model$residuals[(n+1):(2*n)]))
}


However, the result remains basically the same:

> set.seed(111)
> sum(replicate(100000, test()))
[1] 50187
> binom.test(x=50187, n= 100000)

Exact binomial test

data:  50187 and 1e+05
number of successes = 50187, number of trials = 1e+05, p-value = 0.2382
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.4987660 0.5049739
sample estimates:
probability of success
0.50187

• Thanks, liked the simulation way. I think all my problems stems from the fact that that sentence is poorly written/justified on Wikipedia. I didn't manage to find any theoretical proofs to the statement so will accept the numerical proof you gave. – mar tin Mar 2 '17 at 12:15