Global trend and random smooths do not match in HGAM

Hi,

I have fitted a beta-binomial hierarchical GAM in brms. Here are the questions I aim to answer with this model :

  • Is there an increase in success with experience ? (Global trend)
  • Do individuals differ in the relationship between success and experience? (Random smooths)

Here is the model structure in brms :

model_formula <- brmsformula(
    hunting_success | vint(4) ~
        s(Zcumul_xp) +
        s(Zcumul_xp, predator_id, bs = "fs") +
        Zgame_duration
)

Here is the model summary :

r$> summary(modA2)
 Family: beta_binomial2 
  Links: mu = logit; phi = identity 
Formula: hunting_success | vint(4) ~ s(Zcumul_xp) + s(Zcumul_xp, predator_id, bs = "fs") + Zgame_duration 
   Data: data (Number of observations: 100411) 
  Draws: 4 chains, each with iter = 1500; warmup = 500; thin = 4;
         total post-warmup draws = 1000

Smooth Terms:
                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sds(sZcumul_xp_1)                0.23      0.20     0.01     0.73 1.00      484      638
sds(sZcumul_xppredator_id_1)     1.29      0.09     1.13     1.46 1.00      603      672
sds(sZcumul_xppredator_id_2)     9.86      0.46     8.99    10.83 1.01      350      560
sds(sZcumul_xppredator_id_3)     2.59      0.17     2.27     2.93 1.01      775      951

Population-Level Effects:
               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          0.07      0.04    -0.01     0.14 1.01      163      348
Zgame_duration     0.60      0.01     0.59     0.61 1.00      900      830
sZcumul_xp_1       0.42      0.40    -0.22     1.38 1.00      810      936

Family Specific Parameters:
    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi     2.19      0.02     2.15     2.23 1.00      963      955

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Here is the code for my plots :

figA2 <- conditional_effects(
  modA2, method = "fitted",
  robust = FALSE, re_formula = NULL
 )

tabA2_a <- figA2[[4]] # group level smooths
tabA2_b <- figA2[[2]] # predator XP

# Transform as data.table
tabA2_a <- data.table(tabA2_a)
tabA2_b <- data.table(tabA2_b)

# Back transform x-axis values
sequence <- (seq(0, 500, 100) - mean(data$cumul_xp_pred))
standev <- sd(data$cumul_xp_pred)
scaled_breaks <- sequence / standev

# Compute non standardized cumulative XP
tabA2_a[, cumul_xp := (Zcumul_xp * standev) + mean(data$cumul_xp_pred)]
tabA2_b[, cumul_xp := (Zcumul_xp * standev) + mean(data$cumul_xp_pred)]

# Estimates and CIs back to probability scale
tabA2_a[, estimate_back := estimate__ / 4]
tabA2_b[, estimate_back := estimate__ / 4]

 tabA2_b[
   , c("lower_back", "upper_back") := lapply(.SD, function(x) {x / 4}),
   .SDcols = c("lower__", "upper__")
 ]

  tabA2_prey_b[
   , c("lower_back", "upper_back") := lapply(.SD, function(x) {x / 4}),
   .SDcols = c("lower__", "upper__")
 ]

# Plots
plotA2_a <- ggplot(tabA2_a,
                      aes(x = Zcumul_xp,
                          y = estimate_back,
                          color = predator_id)) +
    geom_line(linewidth = 1) +
    scale_color_viridis(discrete = TRUE, option = "D") + #B
    ylab("Hunting success\n") +
    scale_y_continuous(breaks = seq(0, 1, 0.25),
                       limits = c(0, 1)) +
    scale_x_continuous(breaks = scaled_breaks,
                       labels = seq(0, 500, 100)) +
    xlab("\nCumulative experience") +
    custom_theme

 plotA2_b <- ggplot(
   tabA2_b,
   aes(x = Zcumul_xp,
       y = estimate_back)
   ) +
   geom_ribbon(
     aes(
       ymin = lower_back,
       ymax = upper_back),
     fill = "gray") +
   geom_line(linewidth = 1) +
   ylab("Hunting success\n") +
   scale_y_continuous(breaks = seq(0, 1, 0.25),
                      limits = c(0, 1)) +
   scale_x_continuous(breaks = scaled_breaks,
                      labels = seq(0, 500, 100)) +
   xlab("\nCumulative experience") +
   custom_theme

When I use conditional_effects() to plot the global trend and the random smooths for every individual, it appears like the global trend is far lower along the y axis than the trends of every individual, which is surprising because the model should penalize any deviations from the global trend.

Here are the plots :

Moreover, when I look at the data generated by conditional_effects, I see that there is a “hunting_success” column where the value appears to be the mean hunting success, which is 2.05351 or 0.5133775 on the probability scale. If this is the general intercept, then why is my global trend so low? I see that there is also one predator_id, which are my individuals. Does this mean that conditional_effects is plotting the line for one individual and not the global trend?

     r$> head(tabA2_b)
   Zcumul_xp hunting_success Zgame_duration predator_id cond__ effect1__ estimate__      se__   lower__
1: -1.649585         2.05351   2.658875e-16  pred100494      1 -1.649585  0.5750221 0.1409850 0.3333904
2: -1.611829         2.05351   2.658875e-16  pred100494      1 -1.611829  0.5898576 0.1361832 0.3561160
3: -1.574074         2.05351   2.658875e-16  pred100494      1 -1.574074  0.6051662 0.1313861 0.3784752
4: -1.536318         2.05351   2.658875e-16  pred100494      1 -1.536318  0.6209535 0.1266540 0.3981343
5: -1.498563         2.05351   2.658875e-16  pred100494      1 -1.498563  0.6372188 0.1220532 0.4189676
6: -1.460807         2.05351   2.658875e-16  pred100494      1 -1.460807  0.6539533 0.1176481 0.4417914
     upper__     cumul_xp estimate_back lower_back upper_back
1: 0.8766665 -0.003868562     0.1437555 0.08334759  0.2191666
2: 0.8783981  5.046676101     0.1474644 0.08902900  0.2195995
3: 0.8834687 10.097220765     0.1512915 0.09461880  0.2208672
4: 0.8879595 15.147765429     0.1552384 0.09953358  0.2219899
5: 0.8934005 20.198310093     0.1593047 0.10474190  0.2233501
6: 0.9043724 25.248854756     0.1634883 0.11044785  0.2260931

I would really appreciate any help with this. Thank you very much!

Maxime

Ok so diving a bit more into the code, I see that indeed, conditional_effects is plotting the trend for individual “pred100494”.

plotA2_a <- ggplot(tabA2_a[predator_id == "pred100494"],
                      aes(x = Zcumul_xp,
                          y = estimate_back,
                          color = predator_id)) +
    geom_line(linewidth = 1) +
    scale_color_viridis(discrete = TRUE, option = "D") + #B
    ylab("Hunting success\n") +
    scale_y_continuous(breaks = seq(0, 1, 0.25),
                       limits = c(0, 1)) +
    scale_x_continuous(breaks = scaled_breaks,
                       labels = seq(0, 500, 100)) +
    xlab("\nCumulative experience") +
    custom_theme

Thus, my question would be : How can I use conditional_effects to plot the global trend and not the trend of an individual?

Thank you very much

Ultimately, I did some reading and found out that I need to use the conditions argument in the conditional_effects() function and set the grouping factor to NA. This post helped me get an answer.

I did this :

test <- conditional_effects(
  modA2, method = "fitted",
  robust = FALSE, re_formula = NULL,
  effects = "Zcumul_xp",
  conditions = data.frame(predator_id = NA)
 )

And now I think I have the proper trend :

Maxime