Search code examples
rglmmtmbglmm

Maximum estimate curve value to temporal GLMM

I'd like to find a better way to calculate the maximum value of a temporal (my tempora variable is months mes) mixed GLM model using the glmmTMB package. But I'm not sure if it's the best way to do it. I'm using the following code:

#Packages
library(glmmTMB)
library(tidyverse)
library(inflection)
library(ggeffects)
theme_set(theme_bw())

# Dataset
d <- read.csv("https://raw.githubusercontent.com/Leprechault/trash/main/p_lineares.csv",h=T)


# Better GLMM comparing with several others eg. poisson, nbinom2 and ziGamma
m.c <- glmmTMB(N ~ poly(mes,2) + poly(prec,2) + (1 | mes), data = d,
               family = "nbinom1",
               ziformula = ~ 1)

# Find the max inflexion using bede: Bisection Extremum Distance Estimator Method 
ds_F <- cbind(x=d$mes,y=exp(predict(m.c)))
ds_F <- as.data.frame(ds_F)
bb=bede(ds_F$x,ds_F$y,1);bb
bb$iplast

# Plota o gráfico 
ggpredict(m.c, terms = "mes [all]") %>% plot(add.data = TRUE) + geom_vline(xintercept = bb$iplast, colour="red", linetype = "longdash") + labs(
      title = " ",
      x     = "Meses",
      y     = "Insects counts")

bede

Please, any help will be appreciated.

Solution

  • I'm not sure how you want to account for the non-focal variables. I'll use emmeans.

    library(emmeans)
predfun <- function(mval) {
   ee <- emmeans(m.c, ~mes, at = list(mes=mval))
   as.data.frame(ee)$emmean  ## extract predicted value
}
optimize(predfun, interval = c(2,15), maximum = TRUE)

    If you want to be a little bit sloppier, you can predict at a large number of intermediate values and use which.max():

    mval_vec <- seq(2, 12, length = 1000)
ee <- emmeans(m.c, ~mes, at = list(mes=mval_vec)) |> as.data.frame()
with(ee, mes[which.max(emmean)])

    (By the way, using your provided data and code produces a completely different graph for me, one with the opposite curvature - it has its minimum near 7.5 and its maxima at either end, at 2 and 12)