Gavin, these are wonderful tools. Two questions: (1) can the difference_smooths function accommodate complex factor smooths, such as s(x1, x2, by="fac)? I was having trouble getting that to work with the current version of the package. (2) Ideally, this sort of marginalization / conditional differentiation would result in some equivalent of the lme4 prediction routine of being able marginalize out predictors, and factors from associated smooths, through some equivalent of re.form in the prediction code in lme4. Do you think something like that is currently possible or on the horizon? Thanks again for this great work.

Answering my first question: yes it does work for complex smooths, but then all components of said smooth, e.g., "s(x1,x2)" have to be entered without spaces.

Yeah; that's how mgcv refers to the smooths so I figured it was best to stick to that. You can see what mgcv calls your smooths using the smooths() function.

Regarding (2); I'm not sure what you mean. Could you give an example of what you want? I don't think you can marginalize out a covariate from a 2d TPRS smooth as it is needed for the basis to exist. But I could see you being able to specify the values at which say x2 is evaluated at while keeping a fine grid of evaluation points for x1 rather than having a fine grid created for both. You should be able to do that with the newdata argument.

These differences are always done at the smooth level so it you meant to excluded other terms from the model when computing differences, that always happens. In fact, I don't even include the by factors in the difference as this is computing actual differences in estimated smooth functions, not that plus the means of the pair of groups.

Gavin, I think with using the exclude option passed through from predict.gam, that works just fine for what I need.

The changes in .4.2 also look very helpful, particularly the ability to predict on the link function. I am having trouble understanding the functional difference, though, between fitted_samples and predicted_samples. Is the latter adding a random residual or variance component into the posterior estimate? My project goal is to understand both the mean and the variance variance around those mean predictions for various groups in my data, which would be summarized from the simulated predictions. I'd appreciate any clarity you can offer.

Yes, that's the intended differentiation between the two functions is that predicted_samples draws actual observations from the implied posterior distribution (assuming Gaussian approximation to the posterior of the parameters) while fitted_samples() returns expectations.

If you want to understand the mean and variance of your response, you might be best served modelling both (assuming a Gaussian response) or modelling the location, scale, (and possibly shape) of the distribution of the response using GAMLSS models. mgcv has several families that allow this kind of model.

Does the difference-smooth function also provide valid standard errors and confidence intervals for GAMs using non-gaussian distributions?

Yes, at least in the same sense that a credible interval on a single smooth is "valid" for a GLM/GAM. They are approximate intervals, of course, and I'm not aware of results showing that the across-the-function coverage frequentist interpretation of the Bayesian credible interval applies to these differences, but they are just a linear contrast of model terms that we would compute for contrasts with GLMs etc using packages such as emmeans or marginal_effects.

Tx for the quick response and the nice blogs. Suppose you first use the difference_smooth estimate function (or any other similar function in the gratia package; e.g. smooth_estimate,...) and then apply the add_confint function to calculate simultaneous CIs, how and where do you best apply the tranformations using the link function in these functions? Specifying fun = ilink does not seem to do the do the trick, nor does specifying tranform = T...

I would definitely benefit from seeing some more examples using non-gaussian distributions or location-scale type distributions.

add_confint() doesn't add a simultaneous interval, it just adds the usual credible interval, and there isn't a method for objects of class "difference_smooth" as that function already returns a confidence interval. And I have little idea where you are referencing fun and transform arguments from.

There's no point applying the inverse of the link function to a difference; I don't even know what that would mean.

What you seem to want is to have a difference on the response scale; is that correct? If so, this is on the todo list, and shouldn't be too far off, but it's not something you can just do out of the box with {gratia} currently.

If you are doing this to the output of smooth_estimates(), you need to apply the constant term and then transform:

df <- data_sim("eg1", dist = "poisson", scale = 0.5, seed = 1)
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, family = poisson, method = "REML")

smooth_estimates(m, "s(x2)") |> 
    add_constant(coef(m)[1L]) |> 
    transform_fun(inv_link(m))

More examples of what? You have referenced a large part of {gratia}'s functionality but I need specific questions or examples to address.

Many thanks for the answer. I thought add_confint had a similar functionality as confint in gratia. The latter did have the functionality to provide simultaneous CIs I thought.

A code for extracting smooth values with confidence interval for a non-gaussian GAM (e.g. beta distribution) should then look something like this? It seems there is something odd with this code as it doesnt contain the CI values between 0 and 1...

# GAM with beta distribution and factor-smooth interaction smoother
sm1 <- smooth_estimates(m1, "Doy", partial_match = TRUE) |>
add_constant(coef(m1)[1]) |>
transform_fun(inv_link(m1))
names(sm1)
View(sm1)

# Adding CIs
sm2 <- add_confint(sm1, level = 0.95,
parm = "Doy",
partial_match = TRUE)
View(sm2)

# Subset only P splines & smoothers for the P splines
m3 <- subset(sm2, type == "P spline")
View(m3)

I´m indeed referring to differences on the response scale. Im looking forward to any updates in the package concerning this topic. Any idea on where I might look to do this using a not-so-out-of-the box approach in the meantime?

You can't create the interval after you transform to the response scale; it's the link function and its inverse that forces things to stay between any lower / upper bounds, so you have to apply the inverse of the link function after forming the interval on the linear predictor scale:

sm1 <- smooth_estimates(m, "s(x2)") |>
  add_confint() |>
  add_constant(coef(m)[1]) |>
  transform_fun(inv_link(m))

You want `fitted_samples()` to draw from the posterior of the two smooths you want to difference, and then difference the fitted values for each smooth within each draw to get a posterior for th difference.

Thanks for the response but this code does not provide sensible results for me. It seems the actual estimate is not transformed. Perhaps you could pinpoint mistakes? If I just give in inv_link(model), the function also just returns NULL.

Consider the following model:

m1 <- gam(Ratio ~ s(Doy, bs = 'ps', by = Genotype) +
s(Doy, Individual, bs = "re") +
Genotype,
family = betar(),
method = "REML",
data = data2)

To plot the smooths ;

# Grad the link function
ilink <- family(m1)$linkinv
iilink <- inv_link(betar())

# Get the smooth estimates, confidence intervals and transform to the response scale
sm1 <- smooth_estimates(m1, "Doy", partial_match = TRUE) |>
add_confint() |>
add_constant(coef(m1)[1]) |>
transform_fun(iilink)
names(sm1)

# Subset only P splines & smoothers for the P splines
sm2 <- subset(sm1, type == "P spline")

Adding the constant before calculating the Confidence intrérvals and performing the transformations seems more reasonable:

sm1 <- smooth_estimates(m1, "Doy", partial_match = TRUE) |>
add_constant(coef(m1)[1]) |>
add_confint() |>
transform_fun(iilink)
names(sm1)

However, how could one explain the converging\diverging behavior of the confidence intervals?

It seems like Im missing something crucial here...

You should also upgrade to the development version of gratia too (from github) as I've had to update a few of these functions recently as some bugs got introduced that meant the transformation wasn't being applied to all the columns it needed to be.

I'll post a fully reproducible example tomorrow; your example has way too much extraneous information that is irrelevant from my point of view as I don't have your data.

Upgrading to the development version solved the discrepancy between the model fit and the confidence intervals. It also made the inv_link function work properly. Im still stuck with the converging\diverging behavior of the confidence intervals... Not sure whether this is due to some model misspecification or something else...

Are those actual model fits? (i.e. estimates of the smooth that you back-transformed to the response scale?) That's easy to explain - those effects are essentially linear in the link scale and because of the centring (sum to zero) constraint on the smooths, they have to pass through 0 on the y-axis (this is the overall model intercept) and thence have a zero-width confidence interval at that point. When plotting these things with `plot.gam()` and `draw.gam()` you can add on the uncertainty in the model constant term to avoid the bow-tie credible intervals. I don't believe you can do this via `smooth_estimates()` and the code I posted earlier.

Why are you including a random slope for each individual but not a random intercept? I think your model should have `+ s(Individual, bs = "re")` as well as the random slope ranef you already include.

If you want to plot the estimated smooths on the response scale, use `fitted_values()` to generate the data you want/need for the plots, as the intervals computed for the fitted / predicted values will include the constant term and all the other model terms too. You'll likely want to exclude the random effect, so add `exclude = "s(Doy,Individual)"` (`exclude = c("s(Doy,Individual)", "s(Individual)")` if you do want a random intercept) to the `fitted_values()` call so that you only get back the main fixed effects.

You'll like want to prepare data for use by `fitted_values()` to achieve a good plot. I would do:

ds <- data_slice(m3, oGeneotype = evenly(oGenotype), DOY = evenly(Doy, n = 100))

To get evenly spaced `Doy` values for each level of `oGeneotype`, while the `Individual` gets set to some representative value - but don't worry, we ignore the random effects when generating the fitted values if you use the `exclude` argument as mentioned above. Then you can do

fv <- fitted_values(m3, data = ds, scale = "response", exclude = "s(Doy,Individual)")

or
fv <- fitted_values(m3, data = ds, scale = "response", exclude = c("s(Doy,Individual)", "s(Individual)"))

depending on whether you do include the random intercept or not (I think you should...)

The trick of applying the constant and back-transforming partial effects works in some instances but it is better to get into the habit of producing model counterfactuals, where you create the data you want to evaluate the model at and then evaluate it at those values.

Any idea on why it works is some instances but not in other instances? I´ve had some success with the back-transforming partial effects of gams (using a gaulss or beta distribution) but now that I want to plot the confidence interval of a gam with poisson distribution (and factor-smooth interaction smoothers), the transformation doesn´t seem to be accurate anymore (as the confidence intervals surpass the upper limit).

sm1 <- smooth_estimates(m1, "Doy", partial_match = TRUE) |>
add_confint() |>
add_constant(coef(m1)[1]) |>
transform_fun(inv_link(m1))

I´m using the development version of gratia 0.8.1.34. Any thoughts would be appreciated. Thanks in advance!

I'm sorry, but {mgcv} isn't my software so if you have an issue with it (and you know it is an issue with mgcv and not a mistake you made) then you need to contact the maintain of that package. Flooding my comments with code that is almost useless to me (I don't have the data so I can't run it!) isn't very helpful.

Anyway, you didn't follow the correct procedure to fit difference smooths via the ordered-factor approach; you need to include the smooth for the reference level:


m3 <- gam(Ratio ~ s(Doy. bs = "ps") + # reference level smooth
s(Doy, bs = 'ps', by = oGenotype) + # smooth diffs of other levels to the reference level
s(Individual, bs = "re") + # random interpcept for individuals
s(Doy, Individual, bs = "re") + # random linear effect of Doy per individual
oGenotype, # group level effects
family = betar(),
method = "REML",
data = data3)


Because there wasn't a reference level smooth, `gam()` ignored you and fitted what it could, which is basically just a normal factor-by model.

Thanks for the efforts, I appreciate the advice.

Is adding the extra random effect a necessity? I dont have enough data for including that random effect. Unfortunately, there is only one measurement per individual per day of the year. I cant include the extra random effect due to the fact that the model has more coefficients than data, if I do so. For the same reason I cant use GAMMs with the random effect specified as:
random = list(Individual=~1, Doy=~Individual),

Hence, I also cannot apply weights to account for heteroscedasticity. Hence, using GAMLSS would be preferable.

The code you suggested concerning the fitted_values function results in horizontal lines in my case. I notice it also always returns the same Individual and a lot of repeated values.

I indeed overlooked adding the extra reference smoother. However, I still end up getting unreliable results. This might be due to my inability to plot the results decently. If I understand you correctly, the bow-tie credible intervals are not incorrect and could somewhat be corrected for using the draw function?

plot3 <- draw(m1,
unconditional = TRUE,
overall_uncertainty = TRUE, # by setting the overall_uncertainty to true
ci_col = 'blue',
smooth_col = 'red',
constant = coef(m1)[1],
fun = inv_link(m1))
plot3

I think I would benefit from a fully reproducible example of an original model with non-gaussian distribution, factor-interaction smoother and nested random effect, how to plot the term plots on the response scale. Likewise, how to do this for a similar model using the ordered-factor approach.

I am wondering if I run:

comp <- compare_smooths(m1, m2, m3)
draw(comp)

...if there is any chance to get this comparing graphic with an y-axis that corresponds to the response variable?

draw(comp,trans=exp) won't work.

As an alternative, I also could not find out how to compare plots with an y-axis corresponding to the response variable with the command plot.

Like this works perfectly for one plot showing the y-axis corresponding to the response variable

m1 <- gam()
plot(m1, seWithMean = T, trans = exp)

....but would not work with plot (m1, m2, m3, trans=exp) for the aim of comparing within the same graphic.

If you want to do things like this, you are better off generating predictions from the model. Create the data you want to evaluate a smooth at (use `data_slice()`), then do `fitted_values()` on each model you want to compare (but use argument `exclude` to remove the effects you're not interested in, combine the objects and add a `model` variable to distinguish them, and then plot yourself with `ggplot()`. If you're not sure about how to do this, perhaps ask here https://github.com/gavinsim... with a reproducible example. It's easier and more findable to have this stuff answered there.