Getting response scale plots for a binomial GAMM

[yksakana]

Hi, I'm currently stuck on how to plot the results of a binomial GAMM on a response scale using gratia. I specified the model with the following code:
mod <- mgcv::gam(Move ~ te(Temp, Length, by=Condition)+s(Time)+s(ID, bs="re"), family=binomial(link="logit"), data=data)

The goal of this model is to examine the interaction effects of temperature (Temp) and length (Length) on the movement (Move; 0 or 1) of animals under two different conditions (Condition; A and B ). I want to plot the effects of each term on a response scale, ranging from 0 to 1.

I've found that the smooth_estimates function can provide estimates of the partial effects, but I'm unsure how to transform these partial effects to the response scale in the case of a binomial GAMM. Does the following code achieve plotting the effects of the interaction term under condition A on a response scale? I initially added the value of intercept to the partial effect of that term to "est" and then converted it to a response scale using the inverse link function of logit.
Also, I'd like to ask if the "est" value provided by smooth_estimates assume that other terms are fixed at zero?

linkinv <- family(mod)$linkinv
ggplot(data=smooth_estimates(mod) %>%
                         add_confint() %>%
                         filter(smooth=="te(Temp,Length):ConditionA") %>%
                         mutate(response=(est+coef(mod)[1]) %>% linkinv())
            )+
  geom_raster(mapping=aes(x=Temp, y=Length, fill=response))+
  geom_contour(mapping=aes(x=Temp, y=Length, z=response), color="black")

Thank you for your help!

[gavinsimpson]

Firstly, you are using an out of date version of gratia; I suggest you update before proceeding.

Secondly, your model is incorrectly specified; you are not including the group means for Condition. As you have just two levels for Condition it is perhaps best to do this through a fixed parametric effect such that your model will become:

mod <- mgcv::gam(Move ~ Condition +    # Condition means
    te(Temp, Length, by = Condition) + # tensor product for each level of Condition
    s(Time) +                          # smooth of time
    s(ID, bs = "re"),                  # random subject-level intercepts
  family = binomial(link = "logit"),
  data = data)

I'm not sure it makes a lot of sense to plot these on the response scale and it gets complicated because you will now have the parametric effect for the Condition means to deal with (and thus the model constant term (intercept) is for Condition == A).

Because you have to bring in these other effects, the plots aren't really partial effects now, but really are just predictions from the model at some combination of predictors. As such, it might be easier to just generate predictions at the required values and ignore the terms you are not interested in. Assuming you update to the current CRAN version of gratia, this can be done as

ds <- data_slice(mod,
  Temp = evenly(Temp),
  Length = evenly(Temp),
  Condition = evenly(Condition))

# note that something will be in `ds` for `Time` and `ID`
# but we will exclude those effects when predicting

fv <- fitted_values(mod, data = ds, exclude = c("s(Time)", "s(ID)"))

Now you can plot these using ggplot() however you want to as the data in .fitted are on the response scale (as are the credible intervals) as inside fitted_values() I do everything on the link scale and then transform to the response scale if scale = "response", which it is by default.

[yksakana]

Thank you very much for the prompt answer! I really appreciate your help.

The code worked perfectly, but as you pointed out, the plot of the predictions on a response scale does not seem to make much more sense compared to the usual plot provided by gratia::draw()... Thanks anyway!

[gavinsimpson]

Not sure what you mean by "does not seem to make much more sense"? These are the fitted values from the model for both Condition == A and Condition == B, ignoring the effects of $f(\text{Time}$ and $\gamma_{\text{ID}}$ (the random effect for ID). They should be different to what you see for link scale predictions and the output from draw() on your model.

[yksakana]

Thank you for your help. I meant that I found that the plot produced by draw() was sufficient for interpreting the results. I had made some incorrect assumptions about the values provided by these functions. I appreciate your assistance.