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