1
$\begingroup$

I am interested in fitting an additive model with tensor interaction terms. This construction is described in Wood (2006) and used in the mgcv R package.

Assume we want to fit a model like this:

$y = f_{x}(x) + f_{z}(z) + f_{xz}(x, z)$,

where $f_{xz}(\cdot, \cdot)$ is a pure interaction term, constructed by removing $f_{x}(\cdot)$ and $f_{z}(\cdot)$ from its basis. In mgcv, an example such a model is set up with

m0 <- gam(y ~ ti(x, k = 4) + ti(z, k = 4) + ti(x,z, k = c(4, 4)))

As far as I understand, it is critical that the bases for ti(x, ...) and ti(z, ...) are exactly the same as used for ti(x, z, ...). Otherwise the main effects will not be exactly cancelled out from the interaction term, which hence will be confounded. This means that if we instead use any of the following formulations, the ti(x, z, ...) term is strictly speaking not a "pure" interaction term, but will also have some elements of the main effects in it:

m1 <- gam(y ~ ti(x, k = 10) + ti(z, k = 4) + ti(x,z, k = c(4, 4)))
m2 <- gam(y ~ ti(x, k = 4) + ti(z, k = 4, bs = "tp") + ti(x,z, k = c(4, 4)))
m3 <- gam(y ~ s(x, k = 4) + s(z, k = 4) + ti(x,z, k = c(4, 4)))

However, mgcv accepts models like m1 and m2 without giving any messages or warnings. And sometimes these types of models seem necessary. For example, with a highly nonlinear x term and a close to linear z term, it might be necessary to increase k for the main effect. However, having the same k in the interaction term can be computationally very demanding, and lead to nonconvergence. To the best of my knowledge, in the case of model m1, a model which would yield a "pure" interaction term would be the following, but this one is much more computationally demanding due to k = c(10, 4).

m1b <- gam(y ~ ti(x, k = 10) + ti(z, k = 4) + ti(x,z, k = c(10, 4)))

My question is thus:

  • Is my reasoning above correct? If so: Does anyone have a suggested way of measuring the amount of confounding of the interaction term, resulting from using main effects with bases that are not exactly identical to the bases of the interaction term? I assume this should be some measure of colinearity between the terms, but I am not able to come up with an exact formulation.

For reference, here is a reproducible example showing the various setups.

# Data generation code excerpt from help("ti", "mgcv")
library(mgcv)
#> Loading required package: nlme
#> This is mgcv 1.8-27. For overview type 'help("mgcv-package")'.
test1 <- function(x,z,sx=0.3,sz=0.4) { 
  x <- x*20
  (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
                 0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
}
n <- 500
x <- runif(n)/20
z <- runif(n)
f <- test1(x,z)
y <- f + rnorm(n)*0.2

# In this model, ti(x, z) has the main effects completely removed
m0 <- gam(y ~ ti(x, k = 4) + ti(z, k = 4) + ti(x,z, k = c(4, 4)))

# In this model, ti(x, z) has some confounding
m1 <- gam(y ~ ti(x, k = 10) + ti(z, k = 4) + ti(x,z, k = c(4, 4)))

# Also in this mode, there is some confounding
m2 <- gam(y ~ ti(x, k = 4) + ti(z, k = 4, bs = "tp") + ti(x,z, k = c(4, 4)))

Created on 2019-04-08 by the reprex package (v0.2.1)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.