1
$\begingroup$

In the GAM package in R created by Simon Wood there is a selection of the smooth function basis. I sort of understand the options such as bs='tp', bs='cr', etc.

But bs='re' seems odd... that does not seem to be specifying a smooth function. Instead, it is specifying a type of effect. So what smooth basis functions are used for the bs='re'?

I was hoping this webpage would help but it doesn't seem to. https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/smooth.construct.re.smooth.spec.html

$\endgroup$
1
$\begingroup$

There aren't really basis functions in the same sense of the thin plate regression spline or cubic regression spline bases you mention. The basis functions are binary indicator functions that return a 1 if the observation belongs to the a specific group or 0 otherwise.

For these models, all we need to do is to evaluate the basis functions at the observed values of the covariate(s) that form the smooth function: these values create the model matrix for which coefficients will be estimated. We also need a penalty matrix $\mathbf{S}$ that when written in quadratic form, $\boldsymbol{\beta}^{\mathsf{T}}\mathbf{S}\boldsymbol{\beta}$, gives the desired wiggliness penalty; by default the penalty matrix is constructed such that it penalizes the integrated squared second derivative of the spline. This form is convenient as the entire penalty can be written as a function of the penalty matrix and the model coefficients only. Estimated values for $\boldsymbol{\beta}$ are sought which minimise the penalized log-likelihood criterion.

For a random effect smoother, the columns of the model matrix that pertain to the random effect smoother are just a set of binary indicator vectors, indexing whether an observation belongs to the $j$th level of the factor/category. The penalty matrix $\mathbf{S}$ is an identity matrix of size $J$ (where $J$ is the number of levels of the factor forming the random effect). This has the effect of penalizing each group-level coefficient for the random effect in proportion to its squared deviation from 0. In other words, the group level coefficients for the random effect are shrunk towards zero, just as they would be if you fitted a random effect in a mixed effects model.

To see this, here's what mgcv sets up for simple model with 1 random effect smooth and 1 regular smooth:

library('mgcv')
set.seed(1)
dat <- gamSim(1, n = 400, scale = 2, verbose = FALSE) ## 4 term additive truth

## add a factor for the random effect and include it's effect in the simulated data
fac <- as.factor(sample(1:10, 400, replace = TRUE))
dat$X <- model.matrix(~ fac - 1)
b <- rnorm(10) * 0.5
dat <- transform(dat, y = y + X %*% b)

## fit the model
m <- gam(y ~ s(fac, bs = "re") + s(x0), data = dat, method = "ML")

In this example we have a 10-level factor so we have 10 group-level coefficients (for s(fac)) plus 9 basis functions for the regular TPRS term (9 because one of the basis functions in confounded with the model intercept term):

> coef(m)
(Intercept)    s(fac).1    s(fac).2    s(fac).3    s(fac).4    s(fac).5 
 7.71917475  0.24762873  0.48122511 -0.17534482  0.03670067 -0.03800310 
   s(fac).6    s(fac).7    s(fac).8    s(fac).9   s(fac).10     s(x0).1 
-0.52197990  0.53661911 -0.23833765 -0.26088376 -0.06762439  0.02740419 
    s(x0).2     s(x0).3     s(x0).4     s(x0).5     s(x0).6     s(x0).7 
 0.15183166 -0.06347520 -0.15041864 -0.04918018  0.15609539 -0.03564953 
    s(x0).8     s(x0).9 
 0.78256472 -0.26510925

each of which corresponds to a basis function that is a binary indicator variable in the model matrix

> head(model.matrix(m)[, 2:11])
  s(fac).1 s(fac).2 s(fac).3 s(fac).4 s(fac).5 s(fac).6 s(fac).7 s(fac).8
1        1        0        0        0        0        0        0        0
2        1        0        0        0        0        0        0        0
3        0        0        0        0        0        0        0        0
4        0        0        0        0        0        0        0        1
5        1        0        0        0        0        0        0        0
6        0        0        0        0        0        1        0        0
  s(fac).9 s(fac).10
1        0         0
2        0         0
3        0         1
4        0         0
5        0         0
6        0         0

The penalty matrix for the s(fac, bs = 're') random effect smooth is an identity matrix

> m[['smooth']][[1L]][['S']]
[[1]]
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    0    0    0    0    0    0    0    0     0
 [2,]    0    1    0    0    0    0    0    0    0     0
 [3,]    0    0    1    0    0    0    0    0    0     0
 [4,]    0    0    0    1    0    0    0    0    0     0
 [5,]    0    0    0    0    1    0    0    0    0     0
 [6,]    0    0    0    0    0    1    0    0    0     0
 [7,]    0    0    0    0    0    0    1    0    0     0
 [8,]    0    0    0    0    0    0    0    1    0     0
 [9,]    0    0    0    0    0    0    0    0    1     0
[10,]    0    0    0    0    0    0    0    0    0     1
$\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.