Equivalence of random effects via likelihood and smoothed splines

Question

Some fake data:

X = runif(1000)
ff = rep(1:10,100)
E = rnorm(1000)
y = x+e+f
f = as.factor(ff)

When you fit a model like

library(mgcv)
m = gam(y~x+s(f,bs='re'),method='REML')

and

library(lme4)
mr = lmer(y~x+(1|f))

you are basically fitting the same model via different methods. The former specifies the model like a least-squares dummy variable fixed-effects regression, then picks a smoothing parameter to "shrink" the intercept estimates. The latter estimates the variance of the (assumed normally-distributed) variance of the effect of the factor.

Numerically, the two are roughly equivalent (but not exactly). Running the following and letting it spin for a while shows (very) small differences, which increase and decrease at weird intervals, but generally seem to be decreasing in N.

d = c()
X = runif(10000)
ff = rep(1:10,1000)
E = rnorm(10000)
for (i in 20:10000){
    N=i
    x = X[1:N]
    f = ff[1:N]
    e = E[1:N]
    y = x+e+f
    f = as.factor(f)
    m = gam(y~x+s(f,bs='re'),method='REML')
    summary(m)
    library(lme4)
    mr = lmer(y~x+(1|f))
    d[i] = mean(abs(predict(m)-predict(mr)))
    plot(d)
    n = (1:N)
    if (N>20){abline(lm(d~n))}
    }

So my questions:

1. Under what conditions are the two methods exactly identical? Is it only as $N$ or $n$ goes to infinity? Where can I find a (hopefully well-explained) proof of any equivalence?
2. I'm guessing that those weird jumps I get when running that little simulation are artifacts of the various optimizers (or something) used in the software. Is that more or less the case?

EDIT: It turns out that this has been shown by Harville, 1976 and 1977.

Answer 1

1. Generalized additve models (GAM) including smooth terms can be represented by generalized linear mixed models (GLMM) including random effects after appropriate reparameterization (see for proof, for example: Fahrmeir/Kneib/Lang 2004: 742f.)

2. According to the author of the mgcv package the two models should be equivalent. Maybe you are picking up some artifact due to numerical precision of the different routines.