Two methods of adding random effects to a GAM give very different results. Why is this and which one should be used?

Question

A particular section of the mgcv documentation gives multiple methods of incorporating random effects into a generalized additive model. Two methods are 1) to add a smooth term in the class labels using bs="re" in gam; 2) Use the function gamm, which includes similar facilities to lme, combined with the existing functions for gam. However, on simulated data, the two give pretty different model fits. Why is that and which one should be used?

x <- rnorm(1000) 
ID <- rep(1:200,each=5)
y <- x 
for(i in 1:200) y[which(ID==i)] <- y[which(ID==i)] + rnorm(1)
y <- y + rnorm(1000)
ID <- as.factor(ID)
m1 <- gam(y ~ x + s(ID,bs="re"))
m2 <- gamm(y ~ x, random=list(ID=~1) )
mean( (fitted(m1)-fitted(m2$gam))^2 ) 

Answer 1

I suspect the difference is in terms of what fitted values you are getting. If you look at what I would call the model fit, the coefficient estimates, variance terms, the models are identical. Compare summary(m2$lme) with summary(m1) and gam.vcomp(m1).

> summary(m1)

Family: gaussian 
Link function: identity 

Formula:
y ~ x + s(ID, bs = "re")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.05234    0.07932    0.66     0.51    
x            1.01375    0.03535   28.68   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
        edf Ref.df     F p-value    
s(ID) 167.1    199 5.243  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.674   Deviance explained = 72.9%
GCV = 1.2133  Scale est. = 1.0082    n = 1000
> summary(m2$lme)
Linear mixed-effects model fit by maximum likelihood
 Data: strip.offset(mf) 
       AIC     BIC    logLik
  3218.329 3237.96 -1605.165

Random effects:
 Formula: ~1 | ID
        (Intercept) Residual
StdDev:    1.025306 1.003452

Fixed effects: y.0 ~ X - 1 
                 Value  Std.Error  DF   t-value p-value
X(Intercept) 0.0523358 0.07922717 799  0.660578  0.5091
Xx           1.0137531 0.03535887 799 28.670404  0.0000
 Correlation: 
   X(Int)
Xx 0.014 

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.80375873 -0.67702485  0.04245145  0.64026891  2.59257295 

Number of Observations: 1000
Number of Groups: 200 

We see that the estimates for $\hat{\beta}_x$ are 1.01375 and $\hat{\beta}_0$ ~0.05 in both models. Note also the estimate for the between subject/group variance (as a standard deviation), $\hat{\sigma}_{\mathrm{ID}}$ is given as 1.025 in the output from summary(m2$lme). The same information can be computed from the gam model using gam.vcomp(), which gives for m1

> gam.vcomp(m1)
   s(ID) 
1.027795

which is near enough a match for us to not worry about it.

Hence the fitted methods must be returning different fitted values; if we generate fitted values from m2$lme, then we get the same values as that produced by fitted(m1):

> mean((fitted(m1)-fitted(m2$lme))^2) 
[1] 2.966927e-07

which is for all intents and purposes 0.

fitted.lme is documented to return contributions from the the population level (averaging over ID) and for the subject-specific components. This is what fitted.gam will be doing for m1 because it represents the random effect as a spline "fixed" effect. In the case of the gamm model, fitted.gam is returning fitted values for the "fixed" effect part of the model, which would explain the difference. (I'm writing "fixed" because with splines the terms "fixed" and "random" effects get a little blurred.)

In his book, Simon Wood mentions this issue in the example fit to the Sole data using gamm(). He talks about using resid() and fitted() on the $lme component of the gamm model as excluding and including the random effects, respectively. He says this is appropriate for model diagnostics in the specific example here is using.

Which you need will depend on the context of your specific usage/research question.

If all you need are simple random effects like this and you are familiar with GAMs and mgcv then it might be simpler all round to just use the random effect spline basis with gam() rather than having to deal with the weird output of the hybrid that is a GAMM model fitted via gamm(). As I've shown above, the two models are effectively equivalent and the difference you report is just down to whether the fitted values include or exclude the subject (or ID) specific effects.