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

Question

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.

Answer 1

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.