fit GAMM to test for unimodality

unimodal_GAMM.R

library(tidyverse)
library(lme4)
library(gamm4)

# simulate data
set.seed(238472)
simdat <- tibble(env = c(seq(0, 1, length.out = 20),
                       seq(2, 3, length.out = 20),
                       seq(4, 5, length.out = 20)),
                 site = rep(letters[1:3], each = 20),
                 ybar = case_when(
                   between(env, 0, 1) ~ 0 + 1 * env,
                   between(env, 2, 3) ~ 2 + 0 * env,
                   between(env, 4, 5) ~ 5 - 1 * env
                 ),
                 y = map_dbl(.x = ybar, .f = ~rnorm(1, mean = .x, sd = 0.2)))


ggplot(simdat, aes(x = env, y = y, color = site))+
  geom_point()+
  geom_smooth(method = "lm", se = FALSE)  # not the "right" way to do it, but a good shortcut

# fit GLMM
mod_glmm <- lmer(y ~ 1 + env + (1 | site), data = simdat)
summary(mod_glmm)  # no effect of environment

# fit GAMM
mod_gamm <- gamm4(y ~ 1 + s(env), random = ~(1 | site), data = simdat)
summary(mod_gamm$mer)
summary(mod_gamm$gam)  # significant effect of environment

# make predictions across the full range of env
pred_gamm <- as_tibble(predict(mod_gamm$gam,
                               newdata = data.frame(env = seq(0, 5, length.out = 100)),
                               se.fit = TRUE)) %>%
  mutate(env = seq(0, 5, length.out = 100),
         fit_uci = fit + qnorm(0.95) * se.fit,
         fit_lci = fit + qnorm(0.05) * se.fit)

# plot the env smooth to visualize the effect
ggplot()+
  geom_ribbon(data = pred_gamm, aes(x = env, ymin = fit_lci, ymax = fit_uci), fill = "grey")+
  geom_line(data = pred_gamm, aes(x = env, y = fit))+
  geom_point(data = simdat, aes(x = env, y = y, color = site))+
  theme_bw(base_size = 20)