Parametric term evaluation for ziplss models about gratia HOT 5 CLOSED

Hi Philipp, yeah, handling these extended formula models (distributional regression models) is not something gratia does at the moment; Things work for the smooths as Simon just stores them as if they were just a single model with many smooths. The parametric terms code was written without considering those more general models.

This is on the ToDo list that I hope to get to once I get a new course written/start of term, whichever happened first.

from gratia.

I worked out what was going on; a combination of parametric_terms() returning a list for models with more than one linear predictor, and then needing to rename terms according to mgcvs weird convention of add .1, .2 to all smooths and terms in the higher linear predictors.

This should now work for all LSS models, but I have thus far only check the ziplss examples you provided.

Thanks for reporting this.

from gratia.

Hi Gavin,
thanks so much for looking into this and expanding the functionality of this great package!

Unfortunately I discovered that there's a further complication in my use case that wasn't covered by my original reprex. It appears that categorical predictors cause additional problems.

Amended reprex:

#gratia error reprex
library(mgcv)
library(gratia)
## simulate some data...
f0 <- function(x) 2 * sin(pi * x); f1 <- function(x) exp(2 * x)
f2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 * 
  (10 * x)^3 * (1 - x)^10
n <- 500;set.seed(5)
x0 <- runif(n); x1 <- runif(n)
x2 <- runif(n); x3 <- runif(n)
x4 <- sample(factor(c('a','b','c')), size = n, replace = TRUE)

## Simulate probability of potential presence...
eta1 <- f0(x0) + f1(x1) - 3
p <- binomial()$linkinv(eta1) 
y <- as.numeric(runif(n)<p) ## 1 for presence, 0 for absence

## Simulate y given potentially present (not exactly model fitted!)...
ind <- y>0
eta2 <- f2(x2[ind])/3
y[ind] <- rpois(exp(eta2),exp(eta2))

## Fit ZIP model... 
b <- gam(list(y~s(x2)+x3,~s(x0)+x1),family=ziplss())
draw(b)

#ZIP model with a categorical predictor
b0 <- gam(list(y~s(x2)+x4,~s(x0)+x1),family=ziplss())
draw(b0)# works but warning issued
plot(b0,pages=1, all.terms = TRUE) #works

#ZIP model with linear and categorical predictor
b1 <- gam(list(y~s(x2)+x3+x4,~s(x0)+x1),family=ziplss())
draw(b1) # fails
draw(b1, parametric = FALSE) # works
plot(b1,pages=1, all.terms = TRUE) #works

Having a single categorical predictor yields a warning, but the partial effect plot appears.

Warning messages:
1: In model.matrix.default(Terms[[i]], mf, contrasts = object$contrasts) :
  variable 'x4' is absent, its contrast will be ignored

Having a categorical predictor and another term leads to a failure

Error: Columns 3, 4, 5, 6, 7, ... (and 3 more) must be named.
Use .name_repair to specify repair.
Call `rlang::last_error()` to see a backtrace
In addition: Warning messages:
1: In model.matrix.default(Terms[[i]], mf, contrasts = object$contrasts) :
  variable 'x4' is absent, its contrast will be ignored
2: In model.matrix.default(Terms[[i]], mf, contrasts = object$contrasts) :
  variable 'x4' is absent, its contrast will be ignored
3: In predict.gam(object, newdata = mf, type = "terms", terms = mgcv_names[ind],  :
  non-existent terms requested - ignoring

with the following traceback

<error>
message: Columns 3, 4, 5, 6, 7, ... (and 3 more) must be named.
Use .name_repair to specify repair.
class:   `rlang_error`
backtrace:
  1. gratia::draw(b1)
  2. gratia:::draw.gam(b1)
  4. gratia:::evaluate_parametric_term.gam(object, term = terms[i])
  6. tibble:::as_tibble.data.frame(evaluated)
  7. tibble:::as_tibble.list(unclass(x), ..., .rows = .rows, .name_repair = .name_repair)
  8. tibble:::lst_to_tibble(x, .rows, .name_repair, col_lengths(x))
  9. tibble:::set_repaired_names(x, .name_repair)
 14. tibble:::repaired_names(names(x), .name_repair = .name_repair)
 15. tibble:::check_unique(new_name)
Call `rlang::last_trace()` to see the full backtrace

which to me appears to indicate that names are being lost as they are being passed from the model frame to the prediction frames, but my knowledge of programming with tibbles is very limited. Any clues?

from gratia.

Thanks for the revised example; I was likely missing some subtleties when implementing this. mgcv can fit such a range of models and this was worked on for simpler models; guess I need to expand the test suites to include some more complex examples.

I’ll take a look at why this is failing.

from gratia.

Finally found some time to fix this after I'd completely changed the way I handle parametric effects (evaluate_parametric_term() has bee soft deprecated for a little while now). Sorry it took so long.

from gratia.