Equivalence of random effects via likelihood and smoothed splines

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.

• It sounds like you could turn your edit into a brief answer. – Glen_b Sep 20 '14 at 11:13