Search code examples
mgcvgratia

gratia::draw(): "'length.out' must be a non-negative number"

I'm trying to plot a factor-smooth term bs='fs' from my model, but I keep getting this error from gratia::draw(mod):

Error in seq.default(from = lower, to = upper, length.out = n) :
'length.out' must be a non-negative number

I can plot all other individual terms, but it seems to be the fs term that it's having trouble with. I've detached and re-installed/updated gratia from github, but the error persists. I've been able to plot models with this term before, so I'm not sure what I've done wrong here.

The model converged and I get the same error with the full dataset. Here's a small subset:

library(dplyr)
library(gratia)
library(mgcv)

toad3 <- subset(toad2_small, fSite %in% c(1,10,20,30,40),
                select = c(CYR, fSeason, fSite, area_sampled, num))

dput(toad3)

> dput(toad3)
structure(list(CYR = c(2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 
2008L, 2008L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2012L, 2012L, 2012L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2013L, 2013L, 2013L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2014L, 2014L, 2014L, 
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2016L, 2016L, 
2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2017L, 
2017L, 2017L, 2017L, 2017L, 2017L, 2017L, 2017L, 2017L, 2018L, 
2018L, 2018L, 2018L, 2018L, 2018L, 2018L, 2018L, 2018L, 2018L, 
2019L, 2019L, 2019L, 2019L, 2019L, 2019L, 2019L, 2019L, 2019L, 
2019L, 2020L, 2020L, 2020L, 2020L, 2020L, 2020L, 2020L, 2020L, 
2020L, 2020L), fSeason = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L), levels = c("DRY", "WET"), class = "factor"), 
    fSite = structure(c(10L, 1L, 20L, 40L, 20L, 30L, 10L, 1L, 
    40L, 30L, 20L, 10L, 1L, 10L, 1L, 20L, 30L, 40L, 40L, 30L, 
    20L, 1L, 1L, 10L, 40L, 30L, 20L, 30L, 20L, 10L, 1L, 40L, 
    40L, 30L, 20L, 1L, 10L, 40L, 30L, 20L, 1L, 10L, 40L, 20L, 
    30L, 10L, 1L, 40L, 30L, 20L, 10L, 1L, 40L, 30L, 20L, 10L, 
    1L, 30L, 20L, 10L, 1L, 40L, 30L, 20L, 1L, 10L, 40L, 30L, 
    20L, 10L, 1L, 40L, 20L, 30L, 10L, 1L, 40L, 30L, 20L, 10L, 
    1L, 40L, 30L, 20L, 10L, 1L, 1L, 10L, 40L, 30L, 20L, 40L, 
    20L, 10L, 1L, 40L, 30L, 20L, 10L, 1L, 10L, 1L, 40L, 30L, 
    20L, 10L, 1L, 40L, 30L, 20L, 40L, 30L, 1L, 20L, 10L, 40L, 
    30L, 20L, 10L, 1L, 20L, 1L, 10L, 40L, 30L), levels = c("1", 
    "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", 
    "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", 
    "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", 
    "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", 
    "43", "44", "45", "46", "47"), class = "factor"), area_sampled = c(3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L), num = c(0L, 0L, 0L, 0L, 1L, 0L, 5L, 0L, 
    0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 6L, 1L, 0L, 0L, 0L, 0L, 
    5L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 6L, 0L, 0L, 1L, 1L, 0L, 
    0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 4L, 0L, 2L, 0L, 1L, 
    0L, 1L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 2L, 0L, 3L, 0L, 3L, 0L, 0L, 
    0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 4L, 0L, 0L, 
    0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 5L)), row.names = c(1L, 
10L, 20L, 22L, 31L, 41L, 42L, 52L, 70L, 78L, 86L, 104L, 110L, 
115L, 126L, 134L, 144L, 154L, 166L, 172L, 185L, 190L, 196L, 197L, 
213L, 227L, 239L, 248L, 258L, 263L, 271L, 281L, 296L, 309L, 319L, 
323L, 324L, 342L, 352L, 362L, 369L, 377L, 384L, 400L, 410L, 418L, 
424L, 432L, 444L, 456L, 462L, 470L, 477L, 486L, 500L, 509L, 517L, 
526L, 539L, 546L, 554L, 563L, 574L, 585L, 592L, 593L, 618L, 628L, 
632L, 634L, 642L, 653L, 670L, 680L, 681L, 691L, 709L, 713L, 728L, 
732L, 739L, 745L, 759L, 769L, 773L, 780L, 798L, 799L, 816L, 819L, 
829L, 855L, 857L, 865L, 874L, 887L, 894L, 904L, 912L, 921L, 929L, 
938L, 949L, 961L, 970L, 976L, 982L, 989L, 999L, 1011L, 1024L, 
1046L, 1047L, 1057L, 1061L, 1072L, 1081L, 1090L, 1099L, 1109L, 
1125L, 1128L, 1136L, 1139L, 1160L), class = "data.frame")

I get the same error with gam() and bam() as well.

mod <- bam(num ~  
              s(CYR) +
              s(CYR, fSeason, bs="sz") +
              
              s(fSite, CYR, bs="fs") + 
              
              offset(log(area_sampled)), 
            
            data = toad3, 
            
            method = 'fREML',
            discrete = TRUE,
            
            select = FALSE,
            family = nb(link = "log"),
            control = list(trace = TRUE),
            drop.unused.levels=FALSE)

> gratia::draw(mod)
Error in seq.default(from = lower, to = upper, length.out = n) : 
  'length.out' must be a non-negative number

> gratia::draw(mod, select = 3)
Error in seq.default(from = lower, to = upper, length.out = n) : 
  'length.out' must be a non-negative number

# These work:
> gratia::draw(mod, select = c(1,2))
> gratia::draw(mod, select = c(1,2), unconditional = T)

Other functions like plot() work, but I think it's only displaying the linear component of the fs term:

> plot(mod3)
Hit <Return> to see next plot: 
Hit <Return> to see next plot:

plot command works

Solution

  • TLDR: use the ordering for the fs smooth s(CYR, fSite, bs="fs") instead.

    It's due to the ordering of the terms in the fs smooth: s(fSite, CYR, bs="fs").

    Normally one would write this as s(CYR, fSite, bs = "fs"), as in you want a smooth of CYR for the levels of fSite. I guess mgcv::smooth.construct.fs.smooth.spec is more tolerant than I expected. I should fix this as if it produces a valid smooth and plot.gam() handles it. then there's no problem with specifying things the other way round. Indeed, I have already had to handle this situation for the sz basis, which also is documented to be specified in one order but allows the reverse.