# When to choose GAM over GLMM and how to include random effect into GAM

This is a follow-up question to this question.

Here is a description of the dataset: the outcome variable is the number of contacts per participant for two periods before lockdown and under lockdown. My research topic is to compare the intra-individual variance of contact numbers before and during the lockdown. The interested variables include age, household size, sex, original country, education level, employment status, marital status, transportation method, whether this contact is recorded during the weekend or not the weekend, region, and province.

The outcome variable (number of contacts) is right skewed. After removing the outliers, the overdispersion is solved. So now I have three models that seem to be good to go. The first model is GLMM with negative binomial distribution. The 2nd model is the Zero-Inflated Negative Binomial model. The 3rd model is the GAM model.

For the first two models, random effects and interaction terms are included. For the GAM model, I could add the interaction term but when I added the random variables, it kept dropping errors. (update: solved)

My question is there seems to be no linear relationship between the outcome variable and two numeric variables, age and household size, so I categorized the age and household size variables according to plots of the relationships. So in this case, should I still consider the GAM?

Here are the plots of the relationship between Age with outcome, and Household size with outcome. Followed by the models, the DHARMa residual plots, the comparisons between the AIC and BIC of the four models (GLMM with Poisson distribution, GLMM with NB distribution, Zero-Inflated Negative Binomial model, GAMM model with NB distribution)and the codes for the plots.

# 1, GLMM with NB distribution
GLMM_nb <- glmer.nb(TotalNum ~ Wave * Age_3cat + Wave * Sex  + Wave * hh_3cat + Wave * (Country2) + Wave * (Employment2) + Wave * (Education2) + Wave * (Marital3) + Wave * (Transport2) + Wave * (region) +  Wave * (Weekend) + (1 | Token) + (1 | province), data = data, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))

# 2, GLMM with Poisson distribution
GLMM_ps <- glmer(TotalNum_nonhh ~ Wave * Age_3cat + Wave * (Sex)  + Wave * hh_3cat + Wave * (Country2) + Wave * (Employment2) + Wave * (Education2) + Wave * (Marital3) + Wave * (Transport2) + Wave * (region) +  Wave * (Weekend) + (1 | Token) + (1 | Bundesland),
data = data, family = poisson,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))

# 3, Zero-Inflated Negative Binomial model
zinb <- glmmTMB(TotalNum ~ Wave * Age_3cat + Wave * Sex + Wave * hh_3cat +
Wave * Country2 + Wave * Employment2 + Wave * Education2 + Wave * Marital2 + Wave * Transport2 + Wave * region + Wave * Weekend +  (1 | Token) + (1 | Bundesland),
ziformula = ~1, family = nbinom2, data = data)

# 4, GAMM
GAMM <- gam(TotalNum_nonhh ~ s(Age, k = 5) + s(hhsize, k = 5) +
Wave * Sex + Wave * region +  Wave* Employment2 + Wave* Marital3 +  Wave * (Education2) + Wave * (Country2)  + Wave * (Transport2) + Wave * (region) + Wave * (Weekend) +
s(Age, by = Wave, k = 5) + s(hhsize, by = Wave, k = 5) +
s(Token, bs = "re") + s(Bundesland, bs = "re"),
family = nb(),
data = data,
method = "REML" )


> AIC(GLMM_nb, GLMM_p, zinb, GAMM) #
df      AIC
GLMM_nb                               31.0000 2186.022
GLMM_p                                30.0000 2286.327
zinb                                  30.0000 2185.038
GAMM                                  93.7624 2184.484

> BIC(GLMM_nb, GLMM_p, zinb, GAMM) #
df      BIC
GLMM_nb                               31.0000 2311.344
GLMM_p                                30.0000 2407.606
zinb                                  30.0000 2306.317
GAMM                                  93.7624 2563.531

ggplot(data, aes(x = Age, y = TotalNum, fill = Wave)) +
geom_jitter(width = 0.3, height = 0.3, alpha = 0.5) + geom_smooth(method = "gam", formula = y ~ s(x, k = 5), col = "blue") +
labs(title = "Scatter Plot of Age vs TotalNum",
x = "Age",
y = "Total Number of Contacts") +
theme_minimal() +
facet_grid(Wave ~ ., scales = "free_y") +
theme(plot.title = element_text(hjust = 0.5))

ggplot(data, aes(x = hhsize, y = TotalNum, fill = Wave)) +
geom_jitter(width = 0.3, height = 0.3, alpha = 0.5) +
geom_smooth(method = "gam", formula = y ~ s(x, k = 5), col = "blue") +
labs(title = "Scatter Plot of Household size vs TotalNum",
x = "HH",
y = "Total Number of Contacts") +
theme_minimal() +
facet_grid(Wave ~ ., scales = "free_y") +
theme(plot.title = element_text(hjust = 0.5))

• There are a couple of issues here. First, Poisson regression makes assumptions about the conditional mean and variance, not the overall mean and variance. Second, what do you mean by "closer to reality"? It is not clear how any observations could tell you what the regression should be (or what that means). That's why you do regression. Commented Jul 14 at 11:40
• Btw @Chao it is usually good practice to retain the original title and content of the question on CV. There are two reasons for this. First, answers produced for the original question will look nonsensical (as mine now does) if the original content is missing. Second, if the query keeps changing, it makes it more and more difficult for answerers to provide an answer, as the goalpost will continuously shift with each change to the original content of the question. Edits to the original question are still okay, so long as the original context remains. Commented Jul 15 at 6:45
• @ShawnHemelstrand sorry, I will be aware next time to keep the original question!
– Chao
Commented Jul 15 at 7:33

#### Answer to OP's Original Question

This part of the question originally addressed OP's question about overdispersion and conventional residual plots from separate count models. See edit to see response to revised question.

Regarding this comment:

Also, when I compare the results of the Poisson and NB models and the residual plots between the two models, it seems that the NB models are a good decision.

Using these two plots do not provide overwhelming evidence of a fix, other than a new range of fitted and residual values, and with it some change in the distribution. It's usually not useful to analyze the typical residuals you get from GLMMs anyway. As an example from this vignette, finding the mis-specification of a Poisson regression using three different residual types leads to very difficult-to-interpret plots, with only five "banded" patterns to go off of:

Consequently, it is not clear from the plots here what is "wrong" about this regression. A better alternative for Poisson regression, if one wants to inspect residuals, is to use simulated ones. The link above shows how to get these from the DHARMa package, which simulates response data from fitted models, calculates their empirical CDF, then standardizes them. Overdispersion is typically much easier to visually identify with DHARMa residuals, as shown below, where a simple QQ plot of the standardized residuals shows an obviously poor pattern.

As to your more pressing concern, we cannot derive what the dispersion is simply by looking at the mean and variance, nor from a residual plot alone. It is first good to at least check if there is any apparent overdispersion problem (Hilbe, 2014, p.39, see also p.82), as there may be other underlying issues with the model. Apparent overdispersion can be determined by any one of the following causes:

• Non-random missingness.
• Outliers.
• Omitted variables / interactions.
• Problematic non-linearity.
• Too many zeroes than modeled by the regression (zero-inflation).
• Censored data (likely not the issue in your case).

Since there is no formal check here for dispersion issues, there may also be underdispersion present, though this is typically uncommon and is usually due to "clumping" of the data together (which doesn't seem to be the case for you). Hilbe's book on count data modeling (referenced below) has a lot of good content on this topic if you want to explore further.

#### Reference

Hilbe, J. M. (2014). Modeling count data. Cambridge Univ. Press.

#### Edit for Revised Question

Note that the original question has been substantially revised and the original plots/question have changed, so that my original answer does not make a ton of sense now. I am editing this answer to provide some help with the follow-up questions, but OP really should add back in all the original content to provide context for this answer.

Now that the jittered scatterplots, syntax, and DHARMa residual plots have been added, it is easier to see what the model is doing. Its not clear which smoothing is being applied in the scatterplot, but it seems apparent to me that most of the curving is attributed to the outlying points, which do not demonstrate clearly the trend they show. I would investigate those data points and see what may be contributing to their departure from the rest of the data. If one is worried about it’s influence on the rest of the fit, a more robust method could be applied, but that may not be necessary in your case. The linear fits with the GLMMs are likely more parsimonious and efficient compared to any spline-based fit in any case. The DHARMa residuals seem to at least convey this for the linear fits added to this question, as they seem much more reasonable now (compared to the original ones that were included in your question). That is my opinion, and others may chime in about the utililty of a spline over a linear fit here.

The GAMM that was posted was clearly mis-specified here. The GAM without random effects doesn’t appear to have any spline terms, so it does not achieve anything major over the other fits. The GAMM has both linear and nonlinear fits to every term, which greatly magnifies the degrees of freedom to be estimated. It also looks like you are using wave as both an interaction term as well as a by-factor smooth. I would just use one or the other, otherwise the model won’t make sense (I’m guessing that is what is going on here, but the model terms are fairly hard to piece apart, which I come back to later).

In any case, the GAM(M) eats up a lot more degrees of freedom, so the information criteria (particularly BIC, which has steeper penalties for complex fits) is consequently poor. The Poisson AIC/BIC predictably does poorly because it simply never achieves the unrealistic assumptions one needs for it to work. I would vote for the more parsimonious negative-binomial GLMM over the other fits, though others may have a different take on that.

As a stylistic note, it appears you include a very redundant interaction term in your model which makes it very difficult to read (and makes it look like you have more parameters than you actually have). The syntax Main Effect * (All Other Effects) allows you to fit all the interaction terms that are shared with one main effect so that it is easier to write. I would also organize your linear terms in the beginning and your nonlinear terms after rather than mixing them up. Otherwise anybody reading your code has to sift through all the terms to figure out what is going on.

• Thank you very much for your explanation, I further checked the data with the GAM model and found that there seems to be not a linear relationship between the outcome variable (number of contacts) and the interested variable (Age), s(Age) with a p-value 4.14e-07 ***, so I am trying the GAM models with interaction term and random effects.
– Chao
Commented Jul 14 at 19:50
• I use the DHARMa package and got three plots (1, Poisson with random effect and interaction term, 2, NB with random effect and interaction term, 3, GAM without random effect and interaction term. I updated the plots in the question. )
– Chao
Commented Jul 14 at 20:06
• Why did you remove the random effect for the GAM? What does it look like if you include it? It may help answerers if you provide some visualization of your data, as it's not immediately clear why some of the nonlinearity is present in the rank-transformed predictions vs residuals in the second pair of plots, which otherwise look okay. R syntax for your models may also be helpful. Commented Jul 15 at 4:47
• Because I tried since yesterday afternoon but still could not add the random effect to the GAM. It keeps dropping errors: Error: grouping factors must have > 1 sampled level. I will add some plots.
– Chao
Commented Jul 15 at 5:15
• Thanks for the edits. Just a few things. First, don't simply throw out the random effects for the GAM...they are necessary to deal with the non-independence of the errors in your model. I would troubleshoot what is the problem. If its working in glmer and not in gam, then it could simply be how the factors are coded (I know for example that factors in gam() require coding of the factors as non-character values). Since some of your plots are very discrete, I would jitter the points to see how the data is dispersed at each value. The dots overlap a lot, making it hard to see. Commented Jul 15 at 6:35