# GAM: Smooth and factor interaction

I am currently working on a GAM. It includes a smooth-factor interaction. I am trying to decide which function is the one I need.

Model_gamlog2 <- mgcv::gam(perpetration_all ~ dar_scale + s(dar_scale, by = age_binary) + age_binary, data = merged_all_complete, family=binomial, method = "REML")

> summary(Model_gamlog2)

Family: binomial

Formula:
perpetration_all ~ dar_scale + s(dar_scale, by = age_binary) +
age_binary

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   1.7862     0.3707   4.818 1.45e-06 ***
dar_scale    -0.2610     7.0429  -0.037    0.970
age_binary2   1.7874     1.2895   1.386    0.166
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(dar_scale):age_binary1 1.000  1.000  0.010   0.920
s(dar_scale):age_binary2 1.474  2.122  2.038   0.421

Rank: 20/21
R-sq.(adj) =  0.0173   Deviance explained = 9.94%
-REML = 39.827  Scale est. = 1         n = 134


Here I use the By= argument. I have some questions regarding this output:

1. Do I need s() as smoothing parameter or are ti() etc. also possible. I just want the smoothing parameters for each level of IV to vary if needed.
2. How do I see whether there is an interaction? As far as I understand smooth terms only tell me that per level of IV age there is no significant correlation to DV, not whether they differ.

Alternatively I have the following code:

Model_gamlog2 <- mgcv::gam(perpetration_all ~ dar_scale + s(dar_scale, age_binary, bs = "fs") + age_binary, data = merged_all_complete, family=binomial, method = "REML")

> summary(Model_gamlog2)

Family: binomial

Formula:
perpetration_all ~ dar_scale + s(dar_scale, age_binary, bs = "fs") +
age_binary

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   1.8761     0.3907   4.802 1.57e-06 ***
dar_scale     0.5246     0.4452   1.178   0.2386
age_binary2   1.2737     0.7561   1.685   0.0921 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(dar_scale,age_binary) 3.346     17  5.558   0.125

R-sq.(adj) =  0.0508   Deviance explained = 13.9%
-REML = 42.792  Scale est. = 1         n = 134


Does s(dar_scale,age_binary) show the interaction effect?

• By interaction effect, what specifically are you trying to capture with your model? Are you just trying to see how the association between dar_scale and perpetration_all looks between ages? Also, what is the nature of the age_binary variable? The naming seems to convey it has been dichotomized in some way, but there is not much to go off of in the question. Commented Jun 23 at 11:59
• Yes the interaction is between dar, which measures anger dysregulation, and the age groups. I dichotomised age with a median split. Commented Jun 23 at 13:52

#### What Not to Do

I will start by saying that it is generally not wise to dichotomize a numeric variable. Generally it results in a loss of information and produces bias in your standard errors. It also just doesn't make sense to do if you are trying to get a realistic approximation of how age works, particularly when age is likely to have nonlinear effects on the outcome anyway. As an example, I have seen in my own work that age has a nonlinear effect on word reading ability in Chinese, where students who get older tend to stop improving in reading ability and need to be tested on more complex reading skills (e.g. reading comprehension). Binning this data would cause some problems. You can see a simulated example below:

#### Simulate Data ####
set.seed(123)
n <- 300
age <- runif(n,0,6)
Age <- age + 6

#### Plot ####
plot(
Age,
pch=21,
bg="gray"
)


The plotted data shows that age seems associated with initial increases in reading, but then rapidly declines after around 7 or 8 years old.

Now suppose we bin this data and plot it again:

#### Visualize ####
bin.age <- ifelse(Age > mean(Age),1,0)

plot(
jitter(bin.age,factor = .1),
xlab="Age",
pch=21,
bg="gray",
)


Now we have completely lost the nonlinear trend, and its not clear now why one group varies a lot more than the other (without knowing how this data was generated).

#### What to Do

Now as far as what you can do instead, I would suggest running the same model with age as continuous, but fit it as a main effect and interaction like so:

gam(
perpetration_all # DV
~ s(dar_scale) # main effect smooth of scale
+ s(age) # main effect smooth of age
+ ti(dar_scale, age), # interaction between age and dar scale
data = merged_all_complete, # data
family=binomial, # using binomial family
method = "REML" # restricted maximum likelihood estimation
)


The s() function fits an independent smooth to each predictor, and ti() creates an interaction between age and your other predictor after conditioning for the main effects. When plotting, you should get back a regression line for each main effect, and a 3D plot of the interaction between them. Importantly, I have removed the linear effect you added to the model (where it was just listed as dar_scale), as this just fits a linear trend (when you already have a smooth in the model).

#### Edit

Per the comments, ti() is usually used for smooth interactions between multiple continuous variables and won't produce any plots if you don't have s() fitted main effects anyway (because there is nothing to build without them). The bs = "fs" argument in a typical s() smooth creates a smooth where each level of a factor shares the same smoothing parameter and is in a sense like the smooth version of random slopes. (See additional edit below).

The s() function fits an independent smooth to each predictor, and ti() creates an interaction between age and your other predictor after conditioning for the main effects. When plotting, you should get back a regression line for each main effect, and a 3D plot of the interaction between them. Importantly, I have removed the linear effect you added to the model (where it was just listed as dar_scale), as this just fits a linear trend (when you already have a smooth in the model).

#### Edit 2

To clarify, the ti() function can be used as both a main effect and interaction. It is however often used to construct interactions in the following way, where the s()/ti() function with a single variable estimates a main effect (by default thin plate regression splines for s() and cubic regression for ti()), and the combined variables in ti() create a separate interaction like the : operator in regular regressions in R.

s(x) + s(z) + ti(s,z)


Consequently, the model can also be constructed as so, where the ti() version has a smooth that is a tensor product but of a single marginal smooth (per Gavin's comments below):

ti(x) + ti(z) + ti(s,z)


Either one works and one can also just fit the interactions by themselves (if that is the goal). The by= argument simply fits a regression line to each level of a factor variable, so that if you have three levels in a factor, the spline produces three separate lines for each. The bs = "fs" argument instead estimates a smooth in the way I described earlier. I will note again that the by = and bs="fs" arguments in a smooth shouldn't apply in this case given you shouldn't be estimating different factors of dichotomized age. While one can compare differences between by-smooths, it shouldn't really apply here.

• Thank you very much. I did run this model before. Unfortunately, I had high concurvity. And as this is just a hypothesis generating analyses (I do not have many participants) I thought of adding another model to the appendix where age is dichotomised. That was done before with this variables. Commented Jun 23 at 16:01
• While trying to do so, I realised that I don´t quite understand how an interaction between a smooth and factor is done in gam. As shown before, by= results in smooths for every factor level but where is the interaction? And s() with 2 variables usually (as fas as I understood) does not discriminate between main effect and interaction. In a glm I would only get the interaction for one factor level with the continuous variable and would need to look further into the interaction effect. Here, I am a bit lost. Commented Jun 23 at 16:07
• The only thing that by= does is create a smooth regression line for each level of your factor variable. In that sense it is a bit like a nonlinear version of an interaction between a continuous and categorical variable, but I assume your issue is statistically testing a difference of smooths. One could perhaps fit a model with a simple age-only smooth and compare it to a by-factor smooth model and see if the LRT test or AIC values differ appreciably if that is the goal, though this only tests the term itself and not each smooth from each other. Commented Jun 23 at 16:28
• Thank you! I was just thinking, would this be applicable if I only want the interaction: ti(dar_scale, age_binary, bs = "fs"? This would use the tensor interaction (excludes the main effects) with a smoothing for the factor levels. Commented Jun 23 at 16:37
• @Laura while there's nothing stopping you modelling the main effects as linear, i) in this case you only knew they were linear because you fitted a model where you allowed them to be non-linear (so any further inference is biased in the wrong way), ii) if you think the interaction is non-linear, I would argue that there is no sensible scientific reason to model the main effects as linear; model them as smooths and if the estimated effect is linear, just live with it. Commented Jun 25 at 7:03