Random effects in Generalised Additive Model - mgcv package

Question

I am confused about how random effects work in a Generalised Additive Model, specifically in mgcv. I want to specify a factor as a random effect. (Specifically, a binomial outcome for the type of job, with a factor for the region the job is in)

Do they work the same way as in a linear mixed effects model (such as lme4)?

Reading through the package guide I see they specify a random effect in the s() function as a smooth. This does not make sense for a factor, I thought it would only make sense to use s() on a continuous variable like time?

I have fit a model with the random effects below. Could some one interpret the effects for me so I can gain an understanding of what they do here?

My data structure if relevant, I am modelling Secondary_yes.

 $ Month        : num [1:1141] 1 1 1 1 1 1 1 1 1 1 ...
 $ Year         : num [1:1141] 2011 2011 2011 2011 2011 ...
 $ Primary      : Factor w/ 10 levels "Foot & Ankle",..: 10 10 7 8 8 9 9 7 3 1 ...
 $ LocSub       : Factor w/ 2 levels "Locum","Substantive": 2 2 1 2 2 2 2 2 2 2 ...
 $ Region       : Factor w/ 14 levels "East Midlands",..: 8 11 5 2 2 8 8 2 8 8 ...
 $ Secondary_yes: Factor w/ 2 levels "0","1": 2 2 1 2 2 1 1 1 1 1 ...

I have the following model which seems to fit quite well, does this make sense?:

mod = gam(Secondary_yes ~ s(Year) + s(Region, bs = "re") + 
    LocSub + s(Primary, bs = "re"), data = df , family = binomial, 
    method = "REML")

Summary:

Family: binomial 
Link function: logit 

Formula:
Secondary_yes ~ s(Year) + s(Region, bs = "re") + LocSub + 
    s(Primary, bs = "re")

Parametric coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -2.2682     0.5081  -4.464 8.05e-06 ***
LocSubSubstantive   0.8329     0.3031   2.748    0.006 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
             edf Ref.df Chi.sq  p-value    
s(Year)    2.229   2.78   3.16    0.314    
s(Region)  8.282  13.00  38.92 5.26e-05 ***
s(Primary) 8.105   9.00 166.46  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.208   Deviance explained =   20%
-REML = 518.17  Scale est. = 1         n = 1141

Answer 1

The random effects in mgcv are proper random effects; there is a way to view random effects, penalized smooths, Gaussian processes, & Gaussian Markov Random Fields all as a Gaussian random field. From that viewpoint, all these things are the same general type of thing and it is only the small details of each specific thing that gives rise to the diversity of terms in this general class.

In mgcv, GAMs are fitted via penalized likelihood where the penalty takes the form of

$$\boldsymbol{\lambda} \boldsymbol{\beta}^{\prime}\boldsymbol{S}\boldsymbol{\beta}$$

where $\boldsymbol{\lambda}$ is a vector of smoothness parameters, $\boldsymbol{\beta}$ are the estimated model coefficients, and $\boldsymbol{S}$ are penalty matrices.

In the normal smooth types in mgcv, with the default settings, the penalty matrix is defined such that $\boldsymbol{\beta}^{\prime}\boldsymbol{S}\boldsymbol{\beta}$ measures the wiggliness of the smooth, i.e. the integrated squared second derivative of the smooth.

But we can imagine different forms for $\boldsymbol{S}$ that measure other types of "complexity" and random effects use one such form for the penalty matrix, an identity matrix. The identity matrix means we have a ridge penalty on the model coefficients (or those associated with the random effect "smooth" in a model with multiple terms) associated with the random effect term, that penalty measures how far the coefficients deviate from 0, which implies shrinkage in the model estimates towards 0. Just like in random effects in a mixed model.

The smoothness parameter for random effects smooths is proportional to the variance parameter (IIRC it is actually related to the inverse of the variance parameter, the precision, rather than the variance itself).

And the Bayesian view of smoothing (with model = "REML" or "ML") means that everything above is really a fancy set of random effects, whether you're using spline smooths, random effects, MRFs, etc.

Why s() for random effects? I think this really just comes down to keeping the UI simple. Also, you wouldn't be able to combine, through tensor products, random effects if they had to be specified via their own function re() say. From the point of view of the model, these different types of effects are just due to variations in the form of the penalty matrix $\boldsymbol{S}$. As such, it makes sens I think to just view the s(), te(), etc as setting up these more general Gaussian fields and not to think of them as smooths in the sense of splines.

We can only set up a "smooth" using a factor if we declare that the basis type is the random effect basis (bs = "re"): if you don't declare this, you'll get an error for a factor.

Do they work the same way as in a linear mixed effects model (such as lme4)?

Yes, but there is a difference if you start adding terms for say random intercepts and slopes. The default syntax in lme results in correlated random effects and you have to use a more complex syntax/formulation to get uncorrelated effects. mgcv* doesn't do correlated random effects, as far as I can tell.

Assuming Region and Primary are factors, then these terms are adding random intercepts (deviations from the overall mean, the model constant term) for each level of the factors. If you plot the model (plot(mod)) you'll see QQ-plots for the two random effects terms, showing the estimated intercepts for the different levels, just as you would if you looked at a QQ plot of the random effects in a model fitted using lme4. The edf column in the output is showing you an estimate of the complexity, effective degrees of freedom, for each term. The Ref.df column shows the reference values for the tests; the difference between the two for the random effects terms is due to the shrinkage of the random intercepts towards 0 (i.e. towards the model intercept).

The test shown here for the random effects are specialised tests developed for random effects and where the reference distribution is corrected for testing on the boundary of the parameter space, unlike the GLRT tests in lme4 when comparing models that differ in terms of the random effects structure only.