5
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ It sounds like you could turn your edit into a brief answer. $\endgroup$ – Glen_b Sep 20 '14 at 11:13
2
$\begingroup$
  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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.