# mgcv::gam(). ANOVA decomposition of a factor-by-variable interaction

Suppose: y is a continuous response variable. z is a factor variable, and x is a continuous variable.

For example,

  #library(mgcv)
#data <- gamSim(eg=4, n =400, dist="normal")
n <- 400
scale <- 2
x <- runif(n, 0, 1)  # x is a continuous variable
z <- as.factor(sample(1:3, n, replace = TRUE))  # z is a factor varible
z.1 <- as.numeric(z == 1)  # dummy coding
z.2 <- as.numeric(z == 2)  # dummy coding
z.3 <- as.numeric(z == 3)  # dummy coding

f1 <- 2 * sin(pi * x)  # effect of x, when z = 1
f2 <- exp(2 * x) - 3.75887   # effect of x, when z = 2
f3 <- 0.2 * x^11 * (10 * (1 - x))^6 + 10 * (10 * x)^3 * (1 - x)^10   # effect of x, when z = 3

e <- rnorm(n, 0, scale)  # noise
y <- f1 * z.1 + f2 * z.2 + f3 * z.3 + e   # response variable; there is a z-by-x interaction on the mean of y.


I’m interested in performing an ANOVA decomposition of the following gam() model:

gam1 <- gam(y  ~ s(x, by =z) )


That is, I want to decompose the model gam1 into something of the form:

gam2  <- gam(y  ~ s(x) + z + ti(x, by =z) )


where “s(x)” corresponds to the main effect of x, the term “z” corresponds to the main effect of z, and the term ti(x, by =z) corresponds to the interaction between x and z.
Here, the problem is that the identifiability condition to be used in the "decomposition" model gam2 is not very clear to me.

I want the term ti(x, by =z) to model “pure” interactions (i.e., the component ti(x, by =z) does not include the main effect of x and z).

If z were a continuous variable, we could do an ANOVA decomposition of the smooth te(x,z) from gam(y~ te(x,z)) into something like: gam(y ~ ti(x) + ti(z) + ti(x, z)).

From Section 5.6.3 of the book “Generalized additive models, 2nd edition (Wood, 2017)” I understand that mgcv::ti() uses a particular identifiability condition that leaves ti(x) + ti(z) orthogonal to ti(x, z), which is great for interpretability and for testing the significance of the interaction term.

I understand that the identifiability constraint enabling such a decomposition is absorbed into the tensor product basis matrix associated with the smooth ti(x, z). In other words, the reparametrized tensor product basis associated with ti(x, z) will make the smooth ti(x, z) automatically exclude the main effects from the marginal terms of x and z (i.e., ti(x) and ti(z)). This is great because ti(x,z) will then specify “pure” interactions between x and z, and the main effect can be separated from the interaction effects.

But in my case, z is a factor variable.

I would like to do an analogous ANOVA decomposition of the smooth s(x, by =z) (this has a smooth per each level of z) obtained from gam(y ~ s(x, by =z)) into something like: gam(y ~ ti(x) + z + ti(x, by =z)).

And for this decomposition, I want to impose an identifiability condition that makes the components ti(x) + z orthogonal to the component ti(x, by =z).

Is there any method in mgcv, for a factor variable z (as in when we specify the smooth ti(x,z) for a continuous variable z), that imposes a particular identifiability condition on the basis matrix associated with the smooth ti(x, by =z), which makes the resulting component ti(x, by =z) orthogonal to ti(x) + z ?

Thank you very much for any help in advance!

## 1 Answer

Your first model y ~ s(x, by = z) need z as a factor also so it should be y ~ z + s(x, by = z).

Would you settle for the ANOVA decomposition afforded by a model similar to a proper ANOVA with a factor covariate? If so, look at the ordered factor variation for factor by smooths. Such a model would look like

oz <- ordered(z)
gam(y ~ oz + s(x) + s(x by = oz), method = 'REML')


The first s(x) now pertains to the reference level of z and the other s(x, by = oz) refers to a difference smooth between the main s(x) and and the other levels.

Another option could be to change the order of the penalty on the by smooth, so that you have

gam(y ~ z + s(x) + s(x by = z, m = 1), method = 'REML')


where the penalty is now on the first derivative of the second function, which penalises deviation from the a flat/horizontal line, which implies deviation from the "global" smooth of x, s(x).

• Thank you very much. gam.fit <- gam(y ~ z + s(x) + s(x, by = z, m = 1), method = 'REML') seems to work pretty well, although it doesn’t always give the desired decomposition. Anyways it’s very interesting to see that the identifiability condition (which enables an ANOVA decomposition) can be absorbed, to some extent, into the penalty matrix in estimation. Thanks! – syhyunpark Apr 1 at 19:02