smooth_estimates() returns random smooths for levels on in the by smooth

[StefanoMezzini]

Hi Gavin,

I have a dataset with control and treatment animals (C_* and T_*, respectively). Plotting the partial effect of s(doy, animal_year, by = group, bs = 'fs') returns the random smooths for all levels for both groups rather than only the levels within each group. This results in unnecessary, flat smooth estimates with the sample-level CIs. (alpha is 0.1 for all geom_ribbons.)

> unique(filter(partials, term == 's(doy,animal_year):groupControl')$animal_year)
 [1] C_100 2023 C_100 2024 C_102 2023 C_102 2024 C_106 2023 C_107 2023 C_108 2023
 [8] C_161 2023 C_161 2024 C_200 2024 T_151 2023 T_151 2024 T_152 2023 T_153 2023
[15] T_153 2024 T_155 2023 T_155 2024 T_158 2023 T_158 2024 T_159 2023 T_160 2023
[22] T_163 2023 T_163 2024 T_166 2023 T_166 2024 T_169 2023
26 Levels: C_100 2023 C_100 2024 C_102 2023 C_102 2024 C_106 2023 ... T_169 2023

Cheers,
Stefano

[gavinsimpson]

Thanks @StefanoMezzini; I'll take a look at what's going wrong here, most likely is the special internal method use for fs smooths isn't doing the required subsetting of the levels for the by smooth.

[gavinsimpson]

Hmm @StefanoMezzini I can't reproduce this with a similar example.

# Analyse the Rat hormone data from Fahrmeir et al (2013) Regression: Models,
#  Methods and Applications. Springer
pkgs <- c("mgcv", "lme4", "ggplot2", "readr", "dplyr", "forcats", "tidyr",
          "gratia")

# load data
vapply(pkgs, library, logical(1), character.only = TRUE, logical.return = TRUE,
       quietly = TRUE)

rats_url <- "https://bit.ly/rat-hormone"
rats <- read_table(rats_url, col_types = "dddddddddddd-")
# ignore the warning - it"s due to trailing white space at the ends of each
#   row in the file

rats <- rats |>
    mutate(treatment = fct_recode(factor(group, levels = c(1, 2, 3)),
                                  Low = "1",
                                  High = "2",
                                  Control = "3"),
           treatment = fct_relevel(treatment, c("Control", "Low", "High")),
           subject = factor(subject))

rats |>
    na.omit() |>
    count(subject) |>
    count(n, name = "n_rats")

plt_labs <- labs(y = "Head height (distance in pixels)",
                 x = "Age in days",
                 colour = "Treatment")

ggplot(rats, aes(x = time, y = response,
   group = subject, colour = treatment)) +
   geom_point(size = 1) +
   geom_line() +
   facet_wrap(~ treatment, ncol = 3) +
   plt_labs

K <- 7

# fit model
m <- bam(
    response ~ treatment + s(time, subject, by = treatment, bs = "fs", k = 7),
  discrete = TRUE, nthreads = 6
)

# plot
draw(m)

Produces

It looks like you are using a pretty old version of the package if you need to index the term variable (in, I presume, output from smooth_estimates()?) as it hasn't been called that for a long time.

If I've misunderstood, can you show me the code you're using and the specific model you are fitting?

[StefanoMezzini]

Yes, I think I was using an old version because I don't have the issue after updating gratia. Sorry!

[gavinsimpson]

Thanks for confirming that the issue is fixed in the latest version @StefanoMezzini.