# In R/mgcv, what's the difference between s(x) + s(y) + ti(x, y), and ti(x) + ti(y) + ti(x, y)

I have a continuous variable and a factor variable, and would like to build a tensor product interaction between the two with a cubic regression spline bs="cr" and a random effect smoother bs = "re".

I understand ti(x) re-parameterised marginal smooth, and but NOT to bs = "cr".

However, when I applied s(x) + s(y) + ti(x, y), and ti(x) + ti(y) + ti(x, y) separately to the data, the Approximate significance of smooth terms is bit different especially at the interaction term, no matter whether I turned on or off np option in ti() for bs = "re".

One gave

Approximate significance of smooth terms:
edf Ref.df  Chi.sq  p-value
s(age)           1.001  1.003 1669.45  < 2e-16 ***
s(region)        8.305 14.000   27.76 3.93e-06 ***
ti(age,region)  10.929 59.000   25.29   0.0878 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.921   Deviance explained = 93.5%
-REML = 523.72  Scale est. = 1         n = 285


and the other gave

Approximate significance of smooth terms:
edf Ref.df  Chi.sq p-value
ti(age)          1.002  1.003 1711.03 < 2e-16 ***
ti(region)       8.640 14.000   30.06 < 2e-16 ***
ti(age,region)  11.662 56.000   19.25 0.00583 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.922   Deviance explained = 93.6%
-REML = 523.83  Scale est. = 1         n = 285


The text book example in Section 7.3 (2nd) uses the first option, while the help manual refers to the second.

When including an interaction term, should we use s(x) or ti(x) as the main effect?, and what are the other differences between s(x) and ti(x).

A major difference is that the default basis in a tensor product, whether created via te(), t2(), or ti() is a cubic regression spline (bs = "cr"). In contrast, the default basis for s() is a low rank thin plate regression spline (bs = "tp"). The principle difference, however, is the default for k in s() and ti(). And it turns out this latter difference is the main reason for the discrepancies.

For s(), the default of k = -1 will result in k = 10 basis functions in the univariate case for both TPRS and CRS bases, before any identifiability constraints are applied. For ti(), the default k = NA sets k = 5 and results in k^d basis functions (d is the dimension of the tensor product, the number of marginal smooths, which in the univariate case if d = 1), and hence 5^1 = 5 basis functions generated.

Setting a common basis type and size should make the models comparable.

Here is a reproducible example:

library("gratia")
library("mgcv")
library("dplyr")
library("ggplot2")
library("tidyr")
library("GGally")

df <- data_sim("eg2", n = 1000, seed = 2, scale = 0.25)


If we plot the basis functions generated by s() and ti() in the univariate setting, we immediately see where the differences arise:

b_s  <- basis(s(x, bs = "tp"), data = df, constraints = TRUE)
b_ti <- basis(ti(x, bs = "tp"), data = df, constraints = TRUE)
draw(b_s) + draw(b_ti) + patchwork::plot_layout(ncol = 2)


Clearly, the basis functions formed using s() and ti() are very different here; there are also many fewer created by ti().

Likewise if we do the above using the default basis for ti(), the CRS basis:

b_s  <- basis(s(x, bs = "cr"), data = df, constraints = TRUE)
b_ti <- basis(ti(x, bs = "cr"), data = df, constraints = TRUE)
draw(b_s) + draw(b_ti) + patchwork::plot_layout(ncol = 2)


This suggests that setting k and bs appropriately should render the models equivalent. Let's check the basis functions, for the CRS basis as it is easiest to see the equivalency once we have the same number of basis functions

b_s  <- basis(s(x, bs = "cr", k = 10), data = df, constraints = TRUE)
b_ti <- basis(ti(x, bs = "cr", k = 10), data = df, constraints = TRUE)
draw(b_s) + draw(b_ti) + patchwork::plot_layout(ncol = 2)


Whatever ti() is doing to remove the "main effects" (which aren't there in the univariate ti() setting) By default, tensor product splines are reparameterised in such a way that the coefficients for the basis functions are in terms of the value(s) of the tensor product smooth at evenly spaced values of the covariate(s), which results in different basis functions for the TPRS basis. To turn off this reparameterisation, use ti(x, ..., np = FALSE).

b_s  <- basis(s(x, bs = "tp", k = 10), data = df, constraints = TRUE)
b_ti <- basis(ti(x, bs = "tp", k = 10), data = df, constraints = TRUE)
draw(b_s) + draw(b_ti) + patchwork::plot_layout(ncol = 2)


But the important thing is that there are now as many basis functions as with the s() version.

Now the model fits are essentially equal, first with the TPRS basis:

m1a <- gam(y ~ s(x, k = 10) + s(z, k = 10) +
ti(x, z, k = 5, bs = "tp"),
data = df, method = "REML")
m1b <- gam(y ~ ti(x, bs = "tp", k = 10) + ti(z, bs = "tp", k = 10) +
ti(x, z, k = 5, bs = "tp"),
data = df, method = "REML")
edf(m1a) |>
bind_rows(edf(m1b))

# A tibble: 6 × 2
smooth    edf
<chr>   <dbl>
1 s(x)     3.90
2 s(z)     4.24
3 ti(x,z) 11.1
4 ti(x)    3.90
5 ti(z)    4.24
6 ti(x,z) 11.1


and with the CRS basis:

m2a <- gam(y ~ s(x, k = 10, bs = "cr") + s(z, k = 10, bs = "cr") +
ti(x, z, k = 5, bs = "cr"),
data = df, method = "REML")
m2b <- gam(y ~ ti(x, bs = "cr", k = 10) + ti(z, bs = "cr", k = 10) +
ti(x, z, k = 5, bs = "cr"),
data = df, method = "REML")
edf(m2a) |>
bind_rows(edf(m2b))

# A tibble: 6 × 2
smooth    edf
<chr>   <dbl>
1 s(x)     3.93
2 s(z)     4.23
3 ti(x,z) 11.1
4 ti(x)    3.93
5 ti(z)    4.23
6 ti(x,z) 11.1


By definition of the way the TPRS basis functions are created by s(), the individual functions (i.e. the columns of the model matrix) are orthogonal. The ti() construction produces highly correlated basis functions, which might be undesirable in some settings.

bs_tidy <- tidyr::pivot_wider(b_ti, names_from = bf, values_from = value,
names_prefix = "bf")
ggpairs(bs_tidy, cols = 5:8) # e.g. for basis funs 1 through 4


I personally recommend using s(x) + s(z) + ti(x,z), but only because Simon has threatened to make asking for univariate tensor product smooths an error. I probably wouldn't use the ti() form with the TPRS basis because of the highly correlated basis functions. If you use the CRS basis then it doesn't really matter, but you do need to set bs = "cr" if you want exactly the same model in s(x) + s(z) + ti(x,z) and ti(x) + ti(z) + ti(x,z) forms.

After I posted I noticed this Answer from non other than Simon Wood, wherein he addresses the last point about the reprarameterised basis functions in the TPRS case with ti(). This is done by default so that the coefficients for the smooth can be interpreted as the values of the smooth at evenly spaced covariate values. If you turn this off with np = FALSE in the constructor, we recover the same basis function as with the s() constructor.

b_s  <- basis(s(x, bs = "tp", k = 10),
data = df, constraints = TRUE)
b_ti <- basis(ti(x, bs = "tp", k = 10, np = FALSE),
data = df, constraints = TRUE)
draw(b_s) + draw(b_ti) + patchwork::plot_layout(ncol = 2)


The reparameterisation is not needed for the CRS, hence we saw equivalent functions earlier with both constructors.

The main difference between the two constructors for univariate smooths, namely the default number of basis functions created, is also addressed in a comment to Simon's Answer.

In hindsight, this Q could have been closed as this Question (which Simon answered) essentially asks the same thing but in a very different context.

• This is the most comprehensive answer I have ever read. Thanks Gavin.
– Fred
Commented Jun 20, 2023 at 13:31
• Thanks @Fred I learned something while writing this so doing so helped me too Commented Jun 20, 2023 at 13:58