# GAMs and random effects: Significant differences between GAMM and GAMM4 outputs

When introducing a random effect in gamm and gamm4, I receive different p-values for the $gam parametric coefficients. Which seems to be due to a difference in SE calculation(?). Differences in the smooth term values as well. Computing a random effect structure in gam gives a comparable result to that of gamm$gam.

EDIT

Reproducible example with modified models:

library(tidyr)
library(dplyr)
library(mgcv)
library(gamm4)
library(performance)
library(ggplot2)

# Generate the data
leaf_rand <- data.frame(
Treatment = rep(factor(c("B", "C", "A", "C", "B", "A", "C", "A", "B")),
times = length(2015:2022)),
Year = rep(2015:2022, each = 9),
Plot = factor(rep(1:9, times = 8)),
N.P = round(c(
15.2, 13.6, 15.8, 14.5, 16.1, 17.5, 19.4, 13.8, 14.5,
16.4, 15.7, 14.1, 14.1, 16.1, 18.4, 19.3, 14.8, 16.9,
17.0, 15.1, 18.9, 17.0, 17.3, 17.0, 19.5, 17.0, 16.0,
18.2, 13.1, 19.7, 17.3, 17.3, 17.7, 22.0, 18.3, 14.9,
17.9, 14.4, 16.3, 15.7, 15.8, 20.0, 14.6, 15.4, 12.8,
17.2, 12.4, 13.9, 11.2, 16.3, 22.0, 18.5, 12.4, 12.1,
16.6, 13.6, 17.0, 10.5, 16.2, 18.3, 17.9, 14.5, 14.0,
22.0, 15.4, 20.0, 17.5, 20.8, 21.8, 21.4, 19.4, 17.1
))
)

### GAMM Models ---------------------------------------------------------------

gam_model <- gamm4(N.P ~ Treatment +
s(Year, Treatment, k = 5, bs = "fs"),
random = ~(1|Plot),
data = leaf_rand, REML = TRUE)

gam_model1 <- gamm(N.P ~ Treatment +
s(Year, Treatment, k = 5, bs = "fs"),
random = list(Plot=~1),
data = leaf_rand, REML = TRUE)

gam_model2 <- gam(N.P ~ Treatment +
s(Year, Treatment, k = 5, bs = "fs") +
s(Plot, k = 9, bs = 're'), data = leaf_rand, method="REML")

# Model summaries
summary(gam_model$$gam) summary(gam_model1$$gam)
summary(gam_model2)

# Compare the models
compare <- compare_performance(gam_model$$gam, gam_model1$$gam, gam_model2)
print(compare)

### ggplot by treatment
ggplot(leaf_rand, aes(x = Year, y = N.P, color = Treatment)) +
geom_point() +
geom_smooth(aes(group = Treatment), method = "loess", se = FALSE) +
theme_minimal() +
theme(legend.position = "bottom") +
labs(
x = "Year",
y = "N.P",
color = "Treatment"
)


What is wrong here?

Which one should I proceed with?

Is it correct to assume that the differences between the treatments are significant according to the $gam parametric coefficient p-values? For more context about the data in question - there are three treatments and three plots per treatment as the repetitions. The treatments have been applied annually during several years and each plot was sampled once per year. It is a first attempt at additive models trying to catch information from an inter annual fluctuating curve, where a linear model perhaps not completely satisfy. The overall goal is still to test for differences between the treatments. Along the entire treatment period, as well as within the years. ## 1 Answer The differences are likely due to the different approaches functions gamm4 and gamm use to approximate the likelihood. nlme (and thereby gam and gamm) uses PQL to approximate the integrands. lme4 (and thereby gamm4) uses Gauss-Hermite quadrature. (RE)ML estimation of GLMMs requires integrating the random effects out of the model likelihood. There is no closed-form solution or ways to solve this analytically, so numerical methods must be used to approximate the integrals. From the package documentation of function gamm4::gamm4: "gamm4 is based on gamm from package mgcv, but uses lme4 rather than nlme as the underlying fitting engine via a trick due to Fabian Scheipl. gamm4 is more robust numerically than gamm, and by avoiding PQL gives better performance for binary and low mean count data." Dimitris Rizopoulos gives a great explanation of PQL and the different ways to numerically approximate the integrals: https://stats.stackexchange.com/a/436711/173546 ### Edit As pointed out by Ben Bolker in a comment, PQL versus GHQ for Gaussian responses should not make a difference. However, as shown above, results of mgcv::gam, mgcv::gamm and gamm4::gamm4 differ. Differences can be due to different optimizers used, but such differences would generally be small. The standard error differences for the parametric terms seem substantial, I don't know why, could be due to instability due to relatively small sample size. If you want to test for differences between the treatments, a by smooth might be more appriate (although it takes up more df than a factor smooth, but I honestly don't know how to use those for testing differences). Assuming you want to take level A for treatment as the reference category, and check whether each of the other two levels differ, I would take the 'ordered factor' approach. This directly allows you to test whether the parametric and smooth terms differ between the different levels of Treatment: leaf_rand$$oTreatment <- ordered(leaf_rand$$Treatment) gam_model <- gamm4(N.P ~ Treatment + s(Year, k = 5) + s(Year, by = oTreatment, k = 5), random = ~(1|Plot), data = leaf_rand, REML = TRUE) gam_model1 <- gamm(N.P ~ Treatment + s(Year, k = 5) + s(Year, by = oTreatment, k = 5), random = list(Plot=~1), data = leaf_rand, REML = TRUE) gam_model2 <- gam(N.P ~ Treatment + s(Year, k = 5) + s(Year, by = oTreatment, k = 5) + s(Plot, bs = 're'), data = leaf_rand, method="REML") # Model summaries and plots summary(gam_model$$gam) summary(gam_model1$$gam) summary(gam_model2) plot(gam_model$$gam) plot(gam_model1$$gam) plot(gam_model2)  The summaries indicate a difference between the smooth effects of Year between Treatment levels A and C. At least, the p-values are $$< .05$$. Only the gamm4 model indicates a difference in the parametric effect of Treatment levels A and C. But point estimates are all in the same direction. Inspecting the plots (and edf values) shows that the difference between the two differing levels can be described by a negative linear effect over time. I do not see the use of comparing fit between the fitted gam, gamm and gamm4 models. They are essentially equivalent models with identical degrees of freedom (hence the warnings), they just used different estimation and/or optimization approaches. • PQL vs GHQ shouldn't be relevant when fitting a LMM (i.e., identity-link Gaussian response). I would love to see a reproducible example ... (this could also be tried with recent versions of glmmTMB, for further comparison ... Commented Aug 10 at 20:17 • @MarjoleinFokkema thank you. I am still not sure how to interpret this into practical means. On which should one rely, if at all. The difference, at least in my case, are quite large it seems coming down to one implying significance, and the other not. – mink Commented Aug 11 at 13:27 • @mink I would have difficulty choosing the single model to trust, and would try to interpret results in tandem. AFAIK, computing (valid) p-values in (G)LMMs and thereby GAMs is not straightforward and I would therefore interpret them with care. Otherwise, if a single model must be picked, I would follow "gamm4 is more robust numerically than gamm", because the package author knows better than I do. As Ben suggests above, a reproducible example would be nice. Commented Aug 11 at 22:32 • Here the differences in the estimated SDs associated with the smooth components are very small, attributable to unimportant computational differences, and likely to be practically unimportant as well. The big differences seem to be in estimated uncertainty of the scale parameters/SDs, which are considerably harder to get at because uncertainty is estimated in different ways (in particular the $outer.info$hessian component that gam.vcomp uses to get confidence intervals on the SDs is missing, and not trivial to reconstruct ... Commented Aug 11 at 23:15 • @MarjoleinFokkema , @BenBolker Added reproducible data with a modified model structure which appears better overall. I would like to add a question here - what would be the appropriate post-hoc for pairwise comparisons and validating the significance of the differences between the treatments indicated by the gamm4$gam?
– mink
Commented Aug 13 at 10:18