Generalized additive mixed model (GAMM) reference level

Question

To investigate the non-linear age trajectories of various brain outcome measures (obtained using structural MRI), we want to run generalized additive mixed model (GAMM) analyses with R's gamm() from the mgcv package on our dataset. Our dataset consists of three groups that we want to statistically compare within one model, much like one would linearly using linear mixed modelling with lmer() or lme().

The formula looks like this:

gamm(Outcome ~ GroupOrdered + Sex + Scanner + Type + s(Age, k = 4) + 
    s(Age, by = GroupOrdered, k = 4),
  random = list(FamilyID = ~1, SubjectID = ~1), data = DF)

where:

• GroupOrdered includes three groups, ordered as follows, such that "Co" (i.e., Controls) is set as the reference level:

DF$GroupOrdered <- ordered(DF$GroupCat, levels = c("Co","Bdo","Szo"))

• Sex (F/M), Scanner (1/2) and Type (DICOM/PARREC) are categorical variables we need to control for
• s(Age, k = 4) is the smooth term for the Age variable
• s(Age, by = GroupOrdered, k = 4) is the interaction of Group with the smoothed Age term
• FamilyID (some individuals are from the same family) and SubjectID (a good number of individuals has been scanned twice) are added as random effects, with SubjectID nested within FamilyID

Now, this model seems to run nicely. For the (interactions with the) smooth terms, this gives the following output:

Approximate significance of smooth terms:
                           edf Ref.df      F p-value    
s(Age)                   1.856  1.856 10.543 0.00056
s(Age):GroupOrderedCoBdo 1.000  1.000  0.728 0.35127
s(Age):GroupOrderedCoSzo 1.000  1.000  0.184 0.58648

To reveal the statistics for the third group comparison, i.e., s(Age):GroupOrderedBdoSzo, you would think that, as with an e.g., lmer(), you simply run the model again, but this time with the reference level of the group variable changed to "Bdo". However, this changing of the reference level changes the statistical output of the Bdo vs Co comparison as well, which is not what I expected, and doesn't happen when running an lmer() or an lme().

Approximate significance of smooth terms:
                          edf   Ref.df      F p-value    
s(Age)                    1.654  1.654  2.954 0.02546
s(Age):GroupOrderedBboCo  1.763  1.763  2.898 0.05875
s(Age):GroupOrderedBDoSZo 1.000  1.000  0.001 0.96854

The same issue occurs when running the analysis on a dataset that only includes data of two of the groups. I've tried numerous things to find out what's happening here, but to no avail. I would like to be able to run the model on the entire dataset, producing smoothed age * group interactions for each of the three group comparisons. If anyone knows how I could achieve this, I would be happy to hear the solution.

I hope I've been clear enough in describing our problem. If any more information is required, please let me know.