Is it appropriate to compare gamm (mgcv) models to get an interaction effect and how to do it?

Question

I want to know if the development of volume over time (and I mean age, not waves) is different between groups. I also have some covariates. I made a simulation dataset:

library(simstudy)

def <- defData(varname = "volume", dist = "normal", formula = 1000, 
                      variance = 100)
def <- defData(def, varname = "group", formula = "1;2;3", dist = "categorical")
def <- defData(def, varname = "age", dist = "normal", formula = 14, variance = 2)
def <- defData(def, varname = "scanner", formula = "1;2", dist = "categorical")

set.seed(1)
sim.data <- genData(220, def)

I came up with the following model:

library(mgcv)

sim.data$group <- as.factor(sim.data$group)

m1 <- gamm(volume ~ group + s(age, by = group, k = 20, bs = "tp") + scanner , data = sim.data, random=list(id=~1))

summary(m1$gam)

which results in the following output:

Family: gaussian 
Link function: identity 

Formula:
volume ~ group + s(age, by = group, k = 20, bs = "tp") + 
    scanner

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 997.3096     2.7303 365.280   <2e-16 ***
group2        0.5983     1.9557   0.306    0.760    
group3        1.2819     1.7592   0.729    0.467    
scanner       1.3806     1.3627   1.013    0.312    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
              edf Ref.df     F p-value
s(age):group1   1      1 1.850   0.175
s(age):group2   1      1 0.001   0.969
s(age):group3   1      1 0.194   0.660

R-sq.(adj) =  -0.0133   
  Scale est. = 2.1375e-07  n = 220

I hoped that I could interpret the difference in development between groups with my s(age,by=group) term, however it gives me the smooth parameters of age, per group. How can I best evaluate the interaction effect (although I also read somewhere that interaction may not be the appropriate term in this additive model)?

I thought of using using anova(), comparing a model with by term to a model without a by term.

m2 <- gamm(volume ~ group + s(age, k = 20, bs = "tp") + scanner , data = sim.data, random=list(id=~1))

anova(m1$gam,m2$gam)

...but this does not give me the output I expected based on my experience comparing lm models with anova().

Family: gaussian 
Link function: identity 

Formula:
volume ~ group + s(age, by = group, k = 20, bs = "tp") + 
    scanner

Parametric Terms:
        df     F p-value
group    2 0.297   0.743
scanner  1 1.026   0.312

Approximate significance of smooth terms:
              edf Ref.df     F p-value
s(age):group1   1      1 1.850   0.175
s(age):group2   1      1 0.001   0.969
s(age):group3   1      1 0.194   0.660

Does anyone know whether comparing models is an appropriate way to study the interaction effect and how to get the actual effect? Preferably a method that I can also use for pair-wise analyses as well.

All the help would be greatly appreciated!

Answer 1

Someone offline offered me the solution. You need to make an ordered factor of your grouping variable and than add a contrast like this:

sim.data$group<-as.factor(sim.data$group)
sim.data$group<-as.ordered(sim.data$group)
contrasts(sim.data$group) <- 'contr.treatment'

you can than do your analysis as usual

m1 <- gamm(volume ~ group + s(age, by = group, k = 20, bs = "tp") + scanner , data = sim.data, random=list(id=~1))

summary(m1$gam)

Hope this helps anyone struggling with the same question that I had.