Error or invalid behaciour in `draw.gam` when using factors in tensor product interaction

[asgersvenning]

Using factors in ti smooth terms can result in errors or invalid behavior with draw.gam.

1D: Error

Here is just a small reprex that produces an error:

library(tidyverse)
library(mgcv)
#> Loading required package: nlme
#> 
#> Attaching package: 'nlme'
#> The following object is masked from 'package:dplyr':
#> 
#>     collapse
#> This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
library(gratia)
#> 
#> Attaching package: 'gratia'
#> The following object is masked from 'package:stringr':
#> 
#>     boundary

# 1D-interaction: Error
iris_gam_1d <- gam(
  Sepal.Length ~ s(Sepal.Width) + ti(Sepal.Width, Species, bs="fs"),
  data = iris,
  method = "REML"
)
plt <- draw(iris_gam_1d)
#> Error in Summary.factor(structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, : 'min' not meaningful for factors

^{Created on 2025-08-27 with reprex v2.1.1}

2D: Invalid behaviour

And here is the invalid behaviour:

library(tidyverse)
library(mgcv)
#> Loading required package: nlme
#> 
#> Attaching package: 'nlme'
#> The following object is masked from 'package:dplyr':
#> 
#>     collapse
#> This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
library(gratia)
#> 
#> Attaching package: 'gratia'
#> The following object is masked from 'package:stringr':
#> 
#>     boundary

# 2D-interaction: Invalid behaviour
iris_gam_2d <- gam(
  Sepal.Length ~ s(Sepal.Width) + s(Petal.Width) + ti(Sepal.Width, Petal.Width, Species, bs="fs"),
  data = iris,
  method = "REML"
)
plt <- draw(iris_gam_2d)

^{Created on 2025-08-27 with reprex v2.1.1}

Fix

The behaviour is quite deeply integrated in gratia and the intended behaviour not entirely clear, but I took a stab at fixing it anyways.

The first issue is the error Error in Summary.factor(structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, : 'min' not meaningful for factors, which is not super helpful, but I tracked it down to this function call:

gratia/R/too_far.R

Lines 165 to 171 in d79f691

	ind <- too_far(
	x = input[[sm_vars[1L]]],
	y = input[[sm_vars[2L]]],
	ref_1 = reference[[sm_vars[1L]]],
	ref_2 = reference[[sm_vars[2L]]],
	dist = dist
	)

Which doesn't check if input[[sm_vars[1L]]] and input[[sm_vars[2L]]] are numeric, as assumed by too_far. Once, I fixed that the regex for the 1D case doesn't give an error, but instead produces this figure:

and these warnings:

Warning messages:
1: `stat_contour()`: Zero contours were generated 
2: In min(x) : no non-missing arguments to min; returning Inf
3: In max(x) : no non-missing arguments to max; returning -Inf
4: `stat_contour()`: Zero contours were generated 
5: In min(x) : no non-missing arguments to min; returning Inf
6: In max(x) : no non-missing arguments to max; returning -Inf
7: `stat_contour()`: Zero contours were generated 
8: In min(x) : no non-missing arguments to min; returning Inf
9: In max(x) : no non-missing arguments to max; returning -Inf

Which is an improvement I think, but is the same unintended behaviour as in the 2D case, only now it gives errors because they both call plot_smooth.bivariate_smooth, which accidentally works okay for the 2D case, but gives these not-so-easy-to-understand warnings for the 1D case.

In any case, I decided to try to fix the issue properly. I identified that when draw_smooth_estimates checks which type of plot to use, in particular in the first part (1D):

gratia/R/smooth-estimates.R

Lines 1198 to 1236 in d79f691

	if (sm_dim == 1L &&
	sm_type %in% c(
	"TPRS", "TPRS (shrink)", "CRS", "CRS (shrink)",
	"Cyclic CRS", "Cyclic P spline", "P spline", "B spline", "Duchon spline",
	"GP",
	"Mono inc P spline",
	"Mono dec P spline",
	"Convex P spline",
	"Concave P spline",
	"Mono dec conv P spline",
	"Mono dec conc P spline",
	"Mono inc conv P spline",
	"Mono inc conc P spline",
	"Mono inc 0 start P spline",
	"Mono inc 0 start P spline"
	)) {
	class(object) <- append(class(object), "mgcv_smooth", after = 0L)
	} else if (grepl("1d Tensor product", sm_type, fixed = TRUE)) {
	class(object) <- append(class(object), "mgcv_smooth", after = 0L)
	} else if (sm_type == "Random effect") {
	class(object) <- append(class(object),
	c("random_effect", "mgcv_smooth"),
	after = 0
	)
	} else if (sm_type == "Factor smooth") {
	class(object) <- append(class(object),
	c("factor_smooth", "mgcv_smooth"),
	after = 0
	)
	} else if (sm_type == "Constr. factor smooth") {
	class(object) <- append(class(object),
	c("sz_factor_smooth", "mgcv_smooth"),
	after = 0
	)
	} else if (sm_type == "SOS") {
	class(object) <- append(class(object),
	c("sos", "mgcv_smooth"),
	after = 0
	)

It has a if-else branch for "Factor smooth", but the ti term is not caught by it, as it has the type label "Tensor product int." instead. Unfortunately, this label is also shared with tensor product interactions between numeric-continuous variables, so it isn't as easy as just checking the label. My suggestion is to change the check to something like:

  } else if (
    sm_type == "Factor smooth" || (
      sm_type %in% c("Tensor product int.", "Tensor product") && 
      any(map_lgl(object[sm_vars], is.factor))
    )) {

Basically, manually checking if any of the smooth variables are a factor if it is either a "Tensor product int." (ti) or just "Tensor product" (te), at least that catches the most common use-cases.

The new plots for the two cases are then:

1D fixed

2D fixed

Can't currently plot multivariate 'fs' smooths.
Skipping: ti(Sepal.Width,Petal.Width,Species)

which at least seems more consistent, although it is a shame that the unintended behaviour which sort of worked, gives that error instead, although it seems (since I didn't add the error message) that it what is intended.

Code

I have forked the repository and created a feature branch here:

https://github.com/asgersvenning/gratia-pr/tree/fix-draw-factor-smooth-interaction

Which also includes a new unit test:

https://github.com/asgersvenning/gratia-pr/blob/fix-draw-factor-smooth-interaction/tests/testthat/test-draw-methods.R#L207-L214

that fails without the code changes, but passes after (at least on my machine).

I hope that makes sense, and if you think my suggested fix looks alright I will open a PR as well!

PS: I also ran all the other unit tests after the fix, and except for some small changes to some snapshots which I think are an artefact of some system differences, not code changes, they all passed.

[gavinsimpson]

Thanks for the report; the 1D example is a little odd usage and likely something I didn't anticipate, a tensor interaction with an fs marginal. I honestly didn't think you could put fs smooths in a tensor product.

I also don't quite understand the 2D behaviour of mgcv; that call is odd as bs is supposed to be specified as a vector for each marginal and fs is invalid for each marginal so mgcv must be doing something odd. I would expect ti(x,z,f, bs = c("cr", "cr", "re")) to be what you really are trying to do. But what mgcv turns that call into I don't know.

I'll have a look at what mgcv actually creates in these two examples and then see what you are proposing as fixes.

[asgersvenning]

I'm 99% sure that mgcv falls back to the default smooth term ("tp") when it isn't passed, so ti(x, fac, bs="fs") is identical to ti(x, fac, bs=c("tp", "fs")).

Regarding what it is attempting to do, I just want to be able to decompose the effect sizes of the "main" effect and the "individual smooth" effects as mentioned in mgcv::factor.smooth under "Sum to zero smooth interactions", so specifying s(x) + ti(x, fac, ...) is preferred over s(x, fac, ...). Gives much more interpretable summary table with summary or overview.

I am also not such a big fan of the way the basis for the continuous covariate is specified with s(x, fac, bs="fs", xt=list(bs="tp")), at least I find it is easier to understand s(x, bs="tp") + ti(x, fac, bs=c("tp", "fs")).

This example might be a bit more understandable:

library(tidyverse)
library(mgcv)
#> Loading required package: nlme
#> 
#> Attaching package: 'nlme'
#> The following object is masked from 'package:dplyr':
#> 
#>     collapse
#> This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
library(gratia)
#> 
#> Attaching package: 'gratia'
#> The following object is masked from 'package:stringr':
#> 
#>     boundary

fmt_formula <- function(fml) fml %>% 
  deparse() %>% 
  paste0(collapse="") %>% 
  str_replace_all("\\s+", " ") %>% 
  str_replace_all("(?<=~) ", "\n     ") %>% 
  str_replace_all(" (?=[\\+])", "\n  ")

tidy_birds <- gratia::bird_move %>% 
  as_tibble %>% 
  filter(species != "sp6") %>% 
  mutate(species = droplevels(species)) %>% 
  uncount(count)

bird_1d <- gam(
  latitude ~ s(week, bs="cc") + ti(week, species, bs = "fs", k=c(8, 5)),
  data = tidy_birds,
  method = "REML"
)

bird_1d_identical <- gam(
  latitude ~ s(week, bs="cc") + ti(week, species, bs = c("tp", "fs"), k=c(8, 5)),
  data = tidy_birds,
  method = "REML"
)

bird_1d_verbose <- gam(
  latitude ~ s(week, bs="cc") + ti(week, species, bs = c("cc", "fs"), k=c(8, 5)),
  data = tidy_birds,
  method = "REML"
)

bird_1d_independent <- gam(
  latitude ~ s(week, bs="cc", by=species, k=8),
  data = tidy_birds,
  method = "REML"
)

bird_s_fs <- gam(
  latitude ~ s(week, species, bs="fs", k=8, xt=list(bs="cc")),
  data = tidy_birds,
  method = "REML"
)

bird_te_fs <- gam(
  latitude ~ te(week, species, bs=c("cc", "fs"), k=c(8, 5)),
  data = tidy_birds,
  method = "REML"
)

map(
  c(
    "bird_1d",
    "bird_1d_identical",
    "bird_1d_verbose",
    "bird_1d_independent",
    "bird_s_fs",
    "bird_te_fs"
  ),
  function(mod_name) {
    mod <- get(mod_name, envir=.GlobalEnv) # I didn't know if I was gonna use the model names ;) 
    
    mod %>% 
      fitted_values(data_slice(., species=species, week=seq(0, 52, 1))) %>% 
      ggplot(aes(week, .fitted, color=species)) +
      geom_point(
        aes(y=latitude),
        alpha = 0.1,
        shape = 16,
        data=tidy_birds
      ) +
      geom_line() +
      geom_ribbon(aes(ymin=.lower_ci, ymax=.upper_ci, fill=after_scale(color)), alpha=0.1) +
      labs(
        title=fmt_formula(mod$formula),
        subtitle=paste0("AIC:", round(AIC(mod), 2))
      ) +
      scale_x_continuous(limits = c(0, 52), expand = expansion()) 
  }) %>% 
  patchwork::wrap_plots(guides = "collect", axes="collect", nrow=2) &
  theme(
    plot.title = element_text(size = 14, hjust = 0),
    legend.position = "none"
  )

^{Created on 2025-08-28 with reprex v2.1.1}

You can see my first point by the top-left and top-middle plots being identical in appearance and AIC (also when looking at model$smooth).

The bottom-left figure just shows that independent smooths is worse in this case.

The last three show the difference between splitting the terms with [s + ti], [s] and [te], as is to be expected using one term of splitting yields a bit more power, but the difference between the two are not super substantial.

[gavinsimpson]

What mgcv is doing with ti(x, f, bs = "fs") means you are relying on several fallbacks as fs is not valid here, but Simon accommodates this silently. As you say, "fs" gets extended to c("tp", "fs") in this case, but "fs" is still not valid for a single factor, which is what mgcv sees here. Instead, mgcv:::smooth.construct.fs.smooth.spec defaults to s(f, bs = "re") in the case of a single factor: https://github.com/cran/mgcv/blob/master/R/smooth.r#L2023-L2027

  ## deal with no factor case, just base smooth constructor
  if (is.null(fterm)) {
    class(object) <- object$base.bs
    return(smooth.construct(object,data,knots))
  }

So, you might as well use ti(x, f, bs = c("tp", "re")), which is both clearer of intent, and also explains why you get invalid behaviour or errors, because gratia doesn't yet handle tensor products involving random effect terms.

bird_1d_re <- gam(
      latitude ~ s(week, bs="cc") + ti(week, species, bs = c("tp", "re"), k=c(8, 5)),
      data = tidy_birds,
      method = "REML"
    )

which gives

> AIC(bird_1d, bird_1d_re)
               df      AIC
bird_1d  38.90097 8585.937
bird_1d3 38.90097 8585.937

But in this case you probably should set the basis to be "cc" for the marginal smooth of week to match that of the "main effect" smooth, assuming we really want to enforce the wrapping at the ends ofweek:

bird_1d_fs2 <- gam(
      latitude ~ s(week, bs="cc") + ti(week, species, bs = c("cc", "re"), k=c(8, 5)),
      data = tidy_birds,
      method = "REML"
    )
bird_1d_re2 <- gam(
      latitude ~ s(week, bs="cc") + ti(week, species, bs = c("cc", "re"), k=c(8, 5)),
      data = tidy_birds,
      method = "REML"
    )

yielding

> AIC(bird_1d_fs2, bird_1d_re2)
                  df      AIC
bird_1d_fs2 34.06218 8615.587
bird_1d_re2 34.06218 8615.587

Your models are really comparing apples with apples because the models where you just use a single s() term for example have fewer basis functions in the model for the week term.

Take the example above (your bird_1d); there will be 10-2 (IIRC) basis functions for the s(week) term and then another 8-1 per species (in your usage of bs = "fs" as that will have a "tp" marginal) or 8-2 per species if we have a "cc" marginal. But for the single s() term version, you allow just 8-1 basis functions per species. Giving more basis functions to the latter to make them more comparable

bird_1_re3 <- gam(
  latitude ~ s(week, k = 10, bs="cc") + ti(week, species, bs = c("cc", "re"), k=c(14, 5)),
    data = tidy_birds,
    method = "REML"
)
bird_s_fs2 <- gam(
  latitude ~ s(week, species, bs="fs", k=14, xt=list(bs="cc")),
    data = tidy_birds,
    method = "REML"
)

In this case, the partitioning into two smooths is actually slightly worse on the AIC metric

> AIC(bird_1_re3, bird_s_fs2)
                 df      AIC
bird_1_re3 57.26408 8448.175
bird_s_fs2 59.90174 8425.965

The model with completely separate smooths per species with the factor by term is also somewhat restricted because of the lower k per species

> bird_1d_independent2 <- gam(
      latitude ~ s(week, bs="cc", by=species, k=14),
      data = tidy_birds,
      method = "REML"
    )
r$> AIC(bird_1d_independent2)
[1] 8443.95

shows better performance when we give the individual smooths somewhat more flexibility. Though even in this case the bird_s_fs2 model is still "better" on AIC.

As described in #169, handling random effects in tensor products is a little more involved when we consider that we might have 3d or 4d smooths (so random 2d or 3d surfaces), and I don't think I'll be able to do much for random 3d surfaces or higher dimensional smooths. smooth_estimates() works for such models, it's just what to do when auto plotting.

Also note that because the "fs" basis is fully penalised, you can do almost as well in AIC terms as the best model in my post if you just include the average smooth and the random smooths

bird_s_fs3 <- gam(
  latitude ~ s(week, bs = "cc") + s(week, species, bs="fs", k=14, xt=list(bs="cc")),
  data = tidy_birds,
  method = "REML"
)

These models aren't quite the same because this latter one has an extra smoothing parameter

> AIC(bird_s_fs2, bird_s_fs3)
                 df      AIC
bird_s_fs2 59.90174 8425.965
bird_s_fs3 59.50597 8426.683

> bird_s_fs2$sp
s(week,species)1 s(week,species)2
      0.00968096       0.68114772
> bird_s_fs3$sp
         s(week) s(week,species)1 s(week,species)2
      5.66840675       0.01136233       0.68184089

which I think is what partially explains the higher AIC.

[asgersvenning]

Ahh, that makes sense, of course! :) I should have known that it might also have fall back on the basis for the factor, when I already knew it did for the continuous variable.

I must have assumed that because:

# These raise the error:
#> Error in quantile.default(xu, seq(0, 1, length = nk)) : 
#>   (unordered) factors are not allowed
gam(latitude ~ s(week, species), data=tidy_birds)
gam(latitude ~ te(week, species), data=tidy_birds)
gam(latitude ~ s(week) + ti(week, species), data=tidy_birds)

# No error
gam(latitude ~ s(week, species, bs="fs"), data=tidy_birds)
gam(latitude ~ te(week, species, bs="fs"), data=tidy_birds)
gam(latitude ~ s(week) + ti(week, species, bs="fs"), data=tidy_birds)

That all three would have similar behaviour, even thought that is indeed not the case:

suppressPackageStartupMessages({
  library(tidyverse)
  library(mgcv)
})

tidy_birds <- gratia::bird_move %>% 
  as_tibble %>% 
  filter(species != "sp6") %>% 
  mutate(species = droplevels(species)) %>% 
  uncount(count)

get_smooth_type <- function(mod) {
  sm <- mod$smooth
  .get_terms_type <- function(terms) {
    out <- lapply(terms, function(x) {
      type <- class(x)[1]
      names(type) <- "type"
      if ("margin" %in% names(x)) {
        margin_type <- .get_terms_type(x$margin)
        type <- c(type, margin_type)
      }
      type
    })
    names(out) <- sapply(terms, function(x) {
      term <- x$term
      if (is.character(term) && length(term) == 1) return(term)
      paste0("(", paste0(term, collapse=", "), ")")
    })
    out
  }
  .get_terms_type(sm)
}

gam(latitude ~ s(week, species, bs="fs"), data=tidy_birds) %>% 
  get_smooth_type() %>% 
  str
#> List of 1
#>  $ (week, species): Named chr "fs.interaction"
#>   ..- attr(*, "names")= chr "type"

gam(latitude ~ te(week, species, bs="fs"), data=tidy_birds) %>% 
  get_smooth_type() %>% 
  str
#> List of 1
#>  $ (week, species):List of 3
#>   ..$ type   : chr "tensor.smooth"
#>   ..$ week   : Named chr "tprs.smooth"
#>   .. ..- attr(*, "names")= chr "type"
#>   ..$ species: Named chr "random.effect"
#>   .. ..- attr(*, "names")= chr "type"

gam(latitude ~ s(week) + ti(week, species, bs="fs"), data=tidy_birds) %>% 
  get_smooth_type() %>% 
  str
#> List of 2
#>  $ week           : Named chr "tprs.smooth"
#>   ..- attr(*, "names")= chr "type"
#>  $ (week, species):List of 3
#>   ..$ type   : chr "tensor.smooth"
#>   ..$ week   : Named chr "tprs.smooth"
#>   .. ..- attr(*, "names")= chr "type"
#>   ..$ species: Named chr "random.effect"
#>   .. ..- attr(*, "names")= chr "type"

^{Created on 2025-08-29 with reprex v2.1.0}

And disregarded the form you suggest:

gam(latitude ~ s(week, bs = "cc") + s(week, species, bs="fs", k=14, xt=list(bs="cc")), data=tidy_birds, method="REML")
#> Warning message:
#> In gam.side(sm, X, tol = .Machine$double.eps^0.5) :
#>   model has repeated 1-d smooths of same variable.

because of the warning.

In any case, would you be open to a PR, considering that the behavior I described is the same for the correct form:

suppressPackageStartupMessages({
  library(tidyverse)
  library(mgcv)
  library(gratia)
})

tidy_birds <- gratia::bird_move %>% 
  as_tibble %>% 
  filter(species != "sp6") %>% 
  mutate(species = droplevels(species)) %>% 
  uncount(count)

gam(latitude ~ s(week, bs = "cc") + ti(week, I(week^2), species, bs=c("cc", "cc", "re")), data=tidy_birds, method="REML") %>%
  draw

gam(latitude ~ s(week, bs = "cc") + ti(week, species, bs=c("cc", "re")), data=tidy_birds, method="REML") %>%
  draw
#> Error in `map()`:
#> ℹ In index: 2.
#> Caused by error in `Summary.factor()`:
#> ! 'min' er ikke meningsfuld for faktorer

^{Created on 2025-08-29 with reprex v2.1.1}

[gavinsimpson]

In the failing case:

gam(latitude ~ s(week, bs = "cc") + ti(week, species, bs=c("cc", "re")), data=tidy_birds, method="REML") %>%
  draw

A PR that identified this specific case and redirected to call the plot_smooth.factor_smooth method would be very welcome, as there shouldn't be any plotting differences between the ti(week, species, bs = c("cc", "re")) with one continuous and one random effect margin and a s(week, species, bs = "fs", xt = list(bs = "cc")) smooth.