# Problems interpreting GAM output

I have been advised to run General Additive Models to be able to describe trends in my data, my data being animal harvest numbers by year. I have done so, but have a problem with interpreting the model output. Hopefully it will be sufficient to show you parts of my script, and then three of my graphs with corresponding model output.

The synthax below, fitting both a linear and a nonlinear component:

library(mgcv)
model1 <- gam(Tot~Year+s(Year),family=poisson,data=ds)
model2 <- gam(Tot~Year+s(Year),family=poisson,data=ls)
model3 <- gam(Tot~Year+s(Year),family=poisson,data=nw)


Then I get the summary output for the models:

summary(model1)
summary(model2)
summary(model3)


Then I group the data in 5-year periods, before generating the data for the plots:

YearP=seq(1975,2015,by=5)
model1.pred=predict(model1,newdata=data.frame(Year=YearP),type="response",se.fit=T)
model2.pred=predict(model2,newdata=data.frame(Year=YearP),type="response",se.fit=T)
model3.pred=predict(model3,newdata=data.frame(Year=YearP),type="response",se.fit=T)


I graph the three models in the similar manner, showing both the data and the model output:

library(gplots)
plotCI(x=YearP, y=model1.pred$$fit,uiw=2*model1.pred$$se.fit,
col="red",lwd=3,cex=1.2,las=1,
xlab="", ylab="Observed and fitted numbers")
points(ds$$Year,ds$$Tot,pch=19,cex=0.9)


Here's my three graphs xlab is years and ylab is "observed and fitted numbers", the data are the black dots and the model predictions are the red bars:

And here are the outputs for the three models:

I am trying to use the best method possible to explain my data, and I believe this is it. I have tried to google this, and I have read the relevant answers at CrossValidated, but have yet to find explanations that make sense to me. I need someone to tell me what the output means.

I specifically need help to understand the estimate of the parametric coefficient and the p-value, what do they tell me? E.g. in graph 1 there is obviously an increase, and in graph 3 a decrease, but coefficient estimates for both are positive. This is difficult to understand.

Then I also need to know what the "Approximate significance of the smooth terms" tell me, especially the edf and the p-value here?

Extremely grateful for help to understand.

• Welcome to CV, Dag. Your question would be very much improved if you edited it (using the "edit" link at lower left) to clarify which specific portions of the output you want help interpreting. You have posted literally dozens of lines of output and dozens of numbers along with your three graphs. Narrow it down. – Alexis Apr 19 '19 at 19:31

Without data or reproducible example (nor time right now to create one of my own) I'm going to speculate that the cause of the 0 estimate and NAs in the summary() output are due to model identifiability problems.

This is most likely due to your model not being the correct way (in mgcv at least) to fit a model like this. First off, we need to note that the basis expansion of Year contains a linear function. Hence the parametric Year effect plus the basis expansion of Year include two of the same terms. I thought mgcv had some way of identifying these kinds of things and handling it, but it may not be able to do something about it unless it essentially ends up messing with the parametric parts of the model (I'm speculating - I need to dig into this more to be sure).

This linear function in the basis is part of the null space of the basis. So is a flat function, which would be confounded with the model intercept, so it is removed from the initial basis expansion as we could add a scalar value to the intercept and remove that same value from the coefficient for the flat basis function and recover the same fit. Hence there are an infinity of models to choose from and that way madness lies.

What we need to do, if you want to fit a model of the form

$$\hat{y_i} = \beta_0 + \beta_1 \text{Year}_i + f(\text{Year}_i)$$

is to exclude the linear function from the basis, which we do by requesting that the basis expansion for s(Year) not have a null space. In mgcv-speak this is done as follows:

m <- gam(Tot ~ Year + s(Year, bs = 'tp', m = c(2,0)),
family = poisson, data = ds, method = 'ML')


Here we're being explicit about wanting to use a thin plate regression spline basis (bs = 'tp') — this is the default basis but I don't think the other bases all allow the same control over the null space, so I'm being explicit.

The m argument is how we pass information to the function that generates the basis on how we want to control the generation. Here, the form c(2, 0) says that we want a penalty on the second derivative (2, which is the default), whilst the 0 indicates that we do not want any null space in the resulting basis. We have to specify the 2; we can't use the single integer form as that is interpreted as the order of the basis. Hence we specify both aspects:

c(penalty_order, null_space_size)


Now you can interpret this model as having an intercept, plus a linear parametric term for a linear trend, plus some non-linear deviation from the linear trend.

You could compare that with a model that only included the parametric terms and compare them via a generalised likelihood ratio test via anova()(although do heed the warnings in ?anova.gam about the limits of this test). Or more simply, just look at the summary() output and the test for the smooth function is a test for the non-linearity. If the non-linear effect is small relative to the parametric trend, then it would suggest that the extra non-linearity is not needed. This is encoded in the test shown in the output from summary() for the smooth term.

Note: you might consider centring the Year variable for model fitting so that Year = 0 corresponds to mean(Year), which would put the intercept in the middle(-ish) of the time series, rather than extrapolating all the way out to Year 0

## Examples

First a Gaussian example, for a slightly non-linear true function F

set.seed(1)
F <- function(x) exp(2 * x)
scale <- 2
x0 <- runif(100)
f <- F(x0)
y <- f + rnorm(100, 0, scale)

m1 <- gam(y ~ s(x0), family = gaussian, method = 'REML')
m2 <- gam(y ~ x0 + s(x0), family = gaussian, method = 'REML')
m3 <- gam(y ~ x0 + s(x0, m = c(2,0)), family = gaussian, method = 'REML')


Here m1 does what we want but without the decomposition into linear and non-linear terms:

> summary(m1)

Family: gaussian

Formula:
y ~ s(x0)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.1928     0.1897   16.83   <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(x0) 2.044  2.563 34.4 1.39e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.47   Deviance explained = 48.1%
-REML = 207.16  Scale est. = 3.5973    n = 100


m2 doesn't work in the sense that we can't uniquely identify the parametric linear term, which gets aliased and as reported in the summary output, the model matrix is rank deficient

> summary(m2)

Family: gaussian

Formula:
y ~ x0 + s(x0)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.1928     0.1897   16.83   <2e-16 ***
x0            0.0000     0.0000      NA       NA
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df     F  p-value
s(x0) 2.031  2.547 34.62 1.37e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Rank: 10/11
R-sq.(adj) =   0.47   Deviance explained = 48.1%
-REML =  205.6  Scale est. = 3.5979    n = 100


m3 works as we get effectively the same fitted values as m1

> all.equal(fitted(m1), fitted(m3))
[1] "Mean relative difference: 6.15775e-07"


and we can identify a parametric linear term plus some non-linearity

> summary(m3)

Family: gaussian

Formula:
y ~ x0 + s(x0, m = c(2, 0))

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.5069     1.1486  -0.441  0.65996
x0            7.1443     2.1875   3.266  0.00151 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df     F p-value
s(x0) 1.044      8 0.304  0.0963 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.47   Deviance explained = 48.1%
-REML = 205.84  Scale est. = 3.5973    n = 100


Here's the same idea but for Poisson data

set.seed(1)
scale <- 0.5
ff <- F(x0)
g <- exp(ff * scale)
f <- log(g)
y <- rpois(100, g)

m4 <- gam(y ~ s(x0), family = poisson, method = 'REML')
m5 <- gam(y ~ x0 + s(x0), family = poisson, method = 'REML')
m6 <- gam(y ~ x0 + s(x0, m = c(2,0)), family = poisson, method = 'REML')


Wherein we see the same rank-deficient model matrix in the case of m5, but otherwise the fits are the same

> summary(m4)

Family: poisson

Formula:
y ~ s(x0)

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.63304    0.05125   31.86   <2e-16 ***
---
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(x0) 2.918   3.64  476.1  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.885   Deviance explained = 86.8%
-REML = 218.43  Scale est. = 1         n = 100
> summary(m5)

Family: poisson

Formula:
y ~ x0 + s(x0)

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.63304    0.05125   31.86   <2e-16 ***
x0           0.00000    0.00000      NA       NA
---
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(x0) 2.918   3.64  476.1  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Rank: 10/11
R-sq.(adj) =  0.885   Deviance explained = 86.8%
-REML = 217.51  Scale est. = 1         n = 100
> summary(m6)

Family: poisson

Formula:
y ~ x0 + s(x0, m = c(2, 0))

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.1134     0.4377  -0.259    0.796
x0            3.3724     0.8254   4.086 4.39e-05 ***
---
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(x0) 1.918      8  16.82 2.21e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.885   Deviance explained = 86.8%
-REML =  217.1  Scale est. = 1         n = 100

• can your answer explain why x + s(x) works without further futzing in the Gaussian family but not the Poisson case? – Andrew M Apr 28 '19 at 13:20
• Though I can verify that s(x, m = c(2, 0)) does at least formally identify the parametric terms from the non-parametric terms in the poisson family, though the resulting estimates are rather strange, which seems to relate to the TPRS splines. Cubic splines perform as expected. Will try to post some code. – Andrew M Apr 28 '19 at 13:28
• Hi Andrew, x + s(x) doesn't work "without further futzing" in the Gaussian case, not reliably in the sense that the model matrix can be rank deficient. I'll add an example where it fails. – Gavin Simpson Apr 29 '19 at 18:53

I'm not entirely clear on why your intercept is 0 for each of your models, as Andrew M pointed out. The default gam() function provides an intercept, which is probably why your intercept did not change when you added + 1 to your model formula.

Why do you feel the need to specify both a parametric term and a smooth term for Year? Including both terms in your model seems to complicate the interpretation unnecessarily. A smooth term for Year by itself would allow you to visualize both the general trend in your data as well as any anomalies. I may be missing something important regarding your field, but based off of what I know so far, I would specify my model as follows:

model <- gam(Tot ~ s(Year, k = 10), family = poisson, data = data)


Since you are trying to understand trends over the course of 40 years, I would also pay attention to the value of k for your smooth term. k specifies the maximum degrees of freedom your smooth can take (which is k-1), and the gam() function automatically constrains k to be reasonably flexible given the number of observations you have. The main thing is that you don't want k to be too small. This can be easily checked with k.check(model).

The easiest way to display the resulting smooth (and hopefully draw conclusions from it) is to use the plot.gam() function:

plot(model, shaded = T)


This will plot the smooth terms in your model with shaded standard errors.

In my opinion, smooth terms are generally difficult to make hard conclusions about. It's not like a parametric term where you can make a statement like: "animal harvests increased by 15% a year." However, if the association with a response is highly nonlinear, smooth terms can provide highly intuitive visualizations. And if you desire to make bold, concrete statements about the predictor, you can always use smooths as a preliminary analysis to a parametric model, where you could specify higher order terms as needed.

• It is perfectly reasonable to fit the model the OP wants as a decomposition into a linear trend plus potentially some non-linear component. That is, assuming that the linear basis function is removed from the basis expansion for the smoother. The theory for confidence/credible intervals and other forms of inference in the penalized spline / empirical Bayesian interpretation of splines starts to have issues where estimated functions are close to the null space; say when the s(Year) in your model is only slightly non-linear. The OP's model, suitably modified gets around this. – Gavin Simpson Apr 24 '19 at 17:46
• @Gavin Simpson That sounds good. I had never done that. I do wonder if that approach would be useful in the OP's situation, however, since it seems like decomposing into linear trends and nonlinear trends would be most useful if you expect there to be seasonal changes. I am completely ignorant about animal harvests, but it seems unlikely that there would be any seasonality on the level of years. Maybe I'm misunderstanding you though. Would it still be useful to decompose a time trend into linear and nonlinear components if you don't anticipate seasonality? – dante Apr 25 '19 at 18:06
• It's perfectly fine to test for non-linearity that is not seasonal - consider annual climate data in certain regions, they exhibit somewhat non-linear behaviour but it is reasonable to ask if this nonlinearity is sufficiently large so that the model with the non-linear trend fits significantly better than the linear-trend model given the extra complexity of the non-linear version. You can decompose data with seasonality too; you just add a cyclic spline for month or day-of-year of observation for example and you can still have the linear + non-linear decomposition. – Gavin Simpson Apr 25 '19 at 22:14

Your model is performing a regression through the origin -- the (Intercept) term is 0. It is unclear why mgcv is doing this, since as is shown below in an example, in the Gaussian case there is no difficulty to include both parametric and smoothed versions of a covariate.

The fact that the model, as fit, is regression through the origin explains why the Year coefficients as are all positive: at year 0, the poisson log-rate is 0, or 1 events occur per year. Given that constraint, by year 1980, where the observed rate is approximately 30 events per year, the parametric "year" coefficients are all positive, and approximately $$10^{-3}$$, since $$\exp(2000 \times 0.0013) \approx 35$$.

It's important to understand how the intercept and parametric year coefficient are identified. The family of smooth functions generated by s certainly includes linear functions. Identification of the parametric terms is achieved by constraining the smooth (s(Year)) term to sum to zero over the observed values of the covariates. This can lead to some counter-intuitive behavior. Consider two smooth functions with a linear component: \begin{align*} Y_1 & = 4 + 2x + \sin(4\pi x), \\ Y_2 & = 4 + 2x + \cos(4\pi x). \end{align*} In both cases, the non-linear part of the function integrates to zero for $$x \in [0,1]$$. But when we fit the gams (dashed lines), and compare the parametric portions to the estimate you'd get from just estimating the linear trend in y (with lm, solid lines), for the sin function, gam overestimates the linear trend, both in terms of the coefficients of the true function, and compared to the (mis-specified) least-squares linear trend. For the cos function, the parametric estimates are within a standard error of the true forms, and the linear part of the gam matches the linear trend estimated from the least-squares linear trend. I don't have a good explanation for why.

## Code

ex = data_frame(x = seq(0, 1, length.out = 100)) %>% mutate(
epsilon = rnorm(100, sd = .05),
y1 = 4 + 2 * x + sin(x * 4 * pi) + epsilon,
y2 = 4 + 2 * x + cos(x * 4 * pi) + epsilon
)

plot(ex$$x, ex$$y1, col = 'black')
points(ex$$x, ex$$y2, col = 'red')

gfit1 = gam(y1 ~ 1 + x + s(x), data = ex)
gfit2 = gam(y2 ~ 1 + x + s(x), data = ex)
lin1 = lm(y1 ~ x, data = ex)
lin2 = lm(y2 ~ x, data = ex)
gam_linear = function(fit) {
mm = model.matrix(fit)[, c('(Intercept)', 'x')]
beta = coef(fit)[c('(Intercept)', 'x')]
mm %*% beta
}

ex = ex %>% mutate(
p1 = predict(gfit1),
p2 = predict(gfit2),
glin1 = gam_linear(gfit1),
glin2 = gam_linear(gfit2),
lin1 = predict(lin1),
lin2 = predict(lin2)
)

lines(ex$$x, ex$$p1, col = 'black', lty = 2)
lines(ex$$x, ex$$p2, col = 'red', lty = 2)
lines(ex$$x, ex$$glin1, col = 'black', lty = 2)
lines(ex$$x, ex$$glin2, col = 'red', lty = 2)

lines(ex$$x, ex$$lin1, col = 'black', lty = 1)
lines(ex$$x, ex$$lin2, col = 'red', lty = 1)
legend(
'topleft',
legend = c(
'4 + 2 * x + sin(x * 4 * pi) (gam)',
'4 + 2 * x + cos(x * 4 * pi) (gam)',
'4 + 2 * x + sin(x * 4 * pi) (lm)',
'4 + 2 * x + cos(x * 4 * pi) (lm)'
),
col = rep(c('black', 'red'), 2),
lty = c(2, 2, 1, 1)
)


• Thank you Andrew M. I am not sure if I got this right, but I have now rerun my script after specifying Tot ~ 1 + Year + s(Year), but all coefficients are still exactly the same. – Dag Apr 20 '19 at 11:14
• Hmm, yeah. I just tried this with a poisson family and see similar behavior, though in my case the parametric x term is getting dropped. Setting s(x, bs = 'cr') seems to "fix" it, but this honestly seems like a bug to me. – Andrew M Apr 22 '19 at 4:21
• It's not a bug per se (unless the side constraint code is not working) as OP's model is trying to fit contains two versions of a linear trend: i) in the parameter Year term, and ii) one in the null space of the TPRS basis for Year. This is making the model potentially unidentifiable — I would expect the effect that you are seeing, that the parametric Year effect is aliased, rather than that the OP describes, which seems to be aliasing the intercept. The CRS basis in mgcv, IIRC, doesn't contain a linear function; if so that might explain the better handling of the identifiability issue. – Gavin Simpson Apr 24 '19 at 17:32
• As @Dag implies, this has nothing to do with regression through the origin; all model formulae in R include the intercept term implicitly. You would need to explicitly remove it with + 0 or - 1 in the formula if you wanted to fit through the origin. The sum-to-zero constraint is there to get rid of the flat function from the basis (confounded with the intercept). To fit the model the OP and you are trying to, you also need to get rid of the other function in the null space of the penalty, the linear function, which will be confounded with the linear parametric term. – Gavin Simpson Apr 24 '19 at 17:39
• @GavinSimpson I agree that the model formulae doesn't explain the effect that Dag sees and will modify my answer. But I don't follow why the sum-to-zero constraint would not identify $s(x)$ from functions of the form $bx$, since it identifies the intercept, and more generally additive smooths $s(x) + s(z)$ of correlated predictors $x, z$. Are you saying $bx + 0$ is in the null-space, $0x + c$ is not? This would make for a very strange spline. In any case, $x + s(x)$ works fine the Gaussian family... – Andrew M Apr 26 '19 at 20:55