# How to extract the penalty term from the GAMM function in R and what is used to estimate the penalty terms

I'm working on some log GDP data and I am testing out a Generalized Additive Mixed Model for it with a penalized cubic regression spline. I am using the GAMM function from the mgcv package in R and I have two questions that have popped up.

1. Once the model has been generated, how can I view what the penalty term(s), $$\lambda$$, is/are for the model?
2. What type of method is used to select the penalty term(s), is it cross-validation or is it some other method as I'm having some difficulty trying to find some proper documentation about it.

Here's a small worked example of the model that I am working with,

data_GDP= c(5.493884, 5.505332, 5.519860, 5.559527, 5.582368, 5.608006, 5.631928, 5.636217,
5.616771, 5.603594, 5.609105, 5.600642, 5.637643, 5.671259, 5.730749, 5.768008, 5.817111,
5.840932, 5.861925, 5.875492, 5.885548, 5.888878, 5.907267, 5.942274, 5.961005, 5.970496,
5.969219, 5.955837, 5.954022, 5.956096, 5.968708, 5.990714, 6.023690, 6.043820, 6.064250,
6.080162, 6.086093, 6.100319, 6.111911, 6.132313, 6.152307, 6.156979, 6.172744, 6.163104,
6.147399, 6.156979, 6.185797, 6.213808, 6.234999, 6.259008, 6.263398, 6.270232, 6.296372,
6.293604, 6.301886, 6.291939, 6.300786, 6.319869, 6.341593, 6.364062, 6.386879, 6.397596,
6.411818, 6.417222, 6.432458, 6.445402, 6.468320, 6.482954, 6.506979, 6.520179, 6.539586,
6.547216, 6.576191, 6.593318, 6.619139, 6.648855, 6.679222, 6.690842, 6.708816, 6.725394,
6.738389, 6.744059, 6.762961, 6.781512, 6.812785, 6.839798, 6.857304, 6.875232, 6.901033,
6.916715, 6.937314, 6.945147, 6.957688, 6.972981, 6.990349, 6.992648, 7.034564, 7.052981,
7.071319, 7.081961, 7.115257, 7.143934, 7.162863, 7.192107, 7.228026, 7.254107, 7.268084,
7.297294, 7.307336, 7.333088, 7.352441, 7.377571, 7.387771, 7.409681, 7.444132, 7.474091,
7.506866, 7.524183, 7.542532, 7.567501, 7.595186, 7.628469, 7.658464, 7.679852, 7.697485,
7.754310, 7.781180, 7.814763, 7.834630, 7.859876, 7.888934, 7.909820, 7.933725, 7.936446,
7.957352, 8.001556, 8.046934, 8.059118, 8.089667, 8.095843, 8.093859, 8.111328, 8.121569,
8.132295, 8.152889, 8.182783, 8.213165, 8.241361, 8.270807, 8.296447, 8.314906, 8.330526,
8.350005, 8.365184, 8.386355, 8.399333, 8.413587, 8.421849, 8.435484, 8.446256, 8.460030,
8.477662, 8.493843, 8.518792, 8.531766, 8.554489, 8.572212, 8.594062, 8.614556, 8.632752,
8.647414, 8.656468, 8.678070, 8.692826, 8.702028, 8.700298, 8.705364, 8.720444, 8.733256,
8.742654, 8.758271, 8.775055, 8.789751, 8.806993, 8.814256, 8.825986, 8.836679, 8.855621,
8.870059, 8.888329, 8.899881, 8.916680, 8.925627, 8.933400, 8.946778, 8.958360, 8.970623,
8.991288, 9.003488, 9.019156)

n = length(data_GDP)
t =1:n
k = floor(n/3)
data_gdp_gamm = gamm(data_GDP ~ s(t, k=k,
bs="cs"), correlation=corAR1())
data_gdp_gamm

• Questions solely about how software works are off-topic here, but you may have a real statistical question buried here. You may want to edit your question to clarify the underlying statistical issue. You may find that when you understand the statistical concepts involved, the software-specific elements are self-evident or at least easy to get from the documentation. May 21, 2022 at 22:09

### Q1

See ?gamObject and in particular note the $sp and $full.sp components of the fitted GAM model object and note their description. As your data_gdp_gamm object is of class "gamm", which is a list with two components

1. $lme the equivalent mixed model object 2. $gam the equivalent GAM object (which won't have all of the components described in ?gamObject)

you'll want to find the smoothing parameters ($$\boldsymbol{\lambda}_j$$) within the $gam component; data_gdp_gamm$gam\$sp for example.

### Q2

This is pretty well covered in ?gamm and its Details section, although it is unreasonable to expect some quite hefty theory to be presented in the documentation of a function, so you'll need to do some further reading of the references cited in that section.

What is happening here is that the theory showing that penalized splines and mixed effects models are two sides of the same coin is being used to fit the GAM as a linear mixed effects model. This Bayesian view of smoothing is that the penalty matrix(es) $$\mathbf{S}_{\lambda}$$ can be viewed as imposing a Gaussian prior on the model coefficients, and the resulting Bayesian (smoothing model) marginal likelihood has exactly the same form as the REML criterion of a mixed model. The smoothing parameters $$\boldsymbol{\lambda}_j$$, in this framework, are proportional to the precision (inverse of the variance) of the random effect variances in the mixed model form.

As such, selection of the smoothing parameters becomes a problem of jointly estimating the model coefficients for part of the smooth (coefs for the basis functions that are in the penalty null space, which are treated as fixed effects) and the posterior modes of the the random effects (whose model matrix contains the wiggly basis functions of the smooth) using either REML or ML. So, it is the Bayesian marginal log likelihood (=== REML criterion)

$$\mathcal{V}_{\text{r}}(\boldsymbol{\lambda}) = \log \int f(\mathbf{y} | \beta) f(\beta) d \beta$$

is being maximised (from Wood, 2017, p263). As this is not a fully Bayesian model (no prior on the smoothing parameters), this kind of approach is also known as empirical Bayes.

Wood, S.N., 2017. Generalized Additive Models: An Introduction with R, Second Edition. CRC Press.

• Thank you for the response Gavin and for answering both questions. Q1 makes more sense as I want to be able to try different values of sp to see how the model compares with different smoothing parameter values. For Q2 I have the old version of Wood's book and wasn't aware of his second version, so I'll check that one out for further reference. Thank you. May 23, 2022 at 15:57