Why do I see a pattern in the residuals in this well specified model? I made a model using simulated data.
I then fitted an OLS on it.
I know the assumptions of OLS are honored since this is simulated data.
Regardless there is a pattern in the residuals, they don't seem to be zero-mean.
Why is that?
Here is my code and the diagnostic plots:
df <- data.frame(id=seq(1, 12, 1))
df$age <- c(18, 19, 20, 40, 41, 42,
            60, 61, 62, 40, 41, 42)
df$treat <- c(rep(1,6), rep(0,6))
df$rec <- 2*df$age + rnorm(nrow(df), 0, 2)

mod2 <- lm(df$rec ~ df$treat+df$age)
print(summary(mod2))
par(mfrow=c(2,2))
plot(mod2)


 A: Thomas answer is great (+1), I just wanted to clear up a particular confusion in the wording of your question:

there is a pattern in the residuals, they don't seem to be zero-mean.

The mean is zero. You can easily check this:
set.seed(1)
df <- data.frame(id=seq(1, 12, 1))
df$age <- c(18, 19, 20, 40, 41, 42,
            60, 61, 62, 40, 41, 42)
df$treat <- c(rep(1,6), rep(0,6))
df$rec <- 2*df$age + rnorm(nrow(df), 0, 2)
mod2 <- lm(df$rec ~ df$treat+df$age)

# Mean value of the residuals
mean(residuals(mod2))

This is equal to -6.473289e-17, which is $0.0000000000000000647 \approx 0$. The only difference from zero here is due to (lack of) precision.
Note that you don't even have to set the second argument of rnorm() to zero:
set.seed(1)
df <- data.frame(id=seq(1, 12, 1))
df$age <- c(18, 19, 20, 40, 41, 42,
            60, 61, 62, 40, 41, 42)
df$treat <- c(rep(1,6), rep(0,6))
df$rec <- 2*df$age + rnorm(nrow(df), 1000, 2) # large mean
mod2 <- lm(df$rec ~ df$treat+df$age)

# Mean value of the residuals
mean(residuals(mod2))

Which returns -1.853384e-17... Still practically zero. So what happened? The 1000 just got added to the intercept.
A: How close the residuals at specific $x$ values are to zero depends on the sample size. Now, the sample size in real examples will be whatever it is, so there's not much use saying it should be bigger, but you do need to calibrate your expectations of what sort of departures from zero mean are detectable with small samples.
A useful step in this direction is to simulate multiple realisations, not just one. Here are four realisations of the top left plot from your diagnostics

Looking at just one of them, you might think there was a pattern. Looking at all four of them shows the sort of 'pattern' that arises just by chance.  You could do this for multiple sample sizes, and get more idea of the sort of patterns that arise by chance from well-specified models vs the sort that mean something.
