# Confusion about what is treated as the "reference" level in GAM using by= ordered factor in smooth

I have some confusion when trying to follow the example (https://fromthebottomoftheheap.net/2017/12/14/difference-splines-ii/) of using ordered factors in mgcv to get p-values for difference smooths. In that blog it states that when data is modeled using the formula Y~Factor+s(X)+s(X, by=Ordered_Factor), where "Factor" and "Ordered_Factor" define the same groups, one is just ordered and the other isn't, that the s(X) term should model the smooth on the reference level of the ordered factor. I'm curious about that "reference level" refers to here. For instance if my factor has ordered levels A,B,C,D,E,G is the reference level A? Or is the reference level the overall smooth for all the data without separating them by factor?

The reason I ask is that if I create some hypothetical data, it appears that the s(X) term fits the smooth for the data overall, not for factor level A:

library(mgcv)
library(ggplot2)
library(gratia)

set.seed(42)

x<-seq(0,20, by=2)

A<-data.frame(Treatment=rep("A",length(x)),Time=x,Expression=rep(rnorm(length(x),0,0.1)))

B<-data.frame(Treatment=rep("B", length(x)),Time=x,Expression=(100*(-exp(-x*0.5)+exp(-x*0.4))+rnorm(length(x), 0 , 0.1)))

C<-data.frame(Treatment=rep("C", length(x)),Time=x,Expression=(100*(-exp(-x*0.5)+exp(-x*0.44))+rnorm(length(x), 0 , 0.1)))

D<-data.frame(Treatment=rep("D", length(x)),Time=x,Expression=(100*(-exp(-x*0.20)+exp(-x*0.174))+rnorm(length(x), 0 , 0.1)))

E<-data.frame(Treatment=rep("E", length(x)),Time=x,Expression=(100*(-exp(-x*0.50)+exp(-x*0.4))+rnorm(length(x), 0 , 0.1)))

G<-data.frame(Treatment=rep("G", length(x)),Time=x,Expression=(100*(-exp(-x*0.50)+exp(-x*0.4))+rnorm(length(x), 0 , 0.1)))

df<-rbind(A,B,C,D,E,G)

df$$Treatment<-factor(df$$Treatment)

df$$Treatmentord<-factor(df$$Treatment, ordered=T)


The data look like this:

I then fit a GAM with mgcv using ordered factors in the model as described in the linked blog and use the gratia package to plot the smooth fits:

md_ordered<-gam(Expression ~ Treatment + s(Time) + s(Time, by=Treatmentord), data=df, method = "REML")

gratia::draw(md_ordered, residuals=T)


The output from gratia::draw() looks like this:

This output indicates to me that the s(Time) smooth is the smooth for Time when the factor levels aren't separated, rather than the smooth for time at the "reference" factor level.

If instead I fit the model without the s(Time) term and using unordered factors and plot the gratia::draw() output:

md_unordered<-gam(Expression ~ Treatment + s(Time, by=Treatment), data=df, method = "REML")
gratia::draw(md_unordered, residuals=T)


Then the output for factor level A of course makes sense, but now I don't get the benefit of being able to statistically test for differences in the smooths.

Can anyone help clarify if it is indeed possible to fit a model that has a smooth on the reference level alone and then difference smooths on the remaining factor levels, or if I'm simply misunderstanding how to use the mgcv and gratia packages. If it's not possible to do what I want, is it appropriate to test for the significance of a smooth but subsetting the data to only have two of the unordered factor levels and comparing the following two models (for example):

df_sub<-rbind(A,B)
df_sub$$Treatment<-factor(df_sub$$Treatment)

md1<-gam(Expression ~ Treatment + s(Time, by=Treatment), data=df_sub, method="REML")
md2<-gam(Expression ~ Treatment + s(Time), data=df_sub, method="REML")

anova(md1,md2)


This is an infelicity in gratia that I haven't thought enough about as yet. Your ordered factor code is correct; it's gratia that doesn't know enough to label the first smooth as belonging to A and the other smooths as difference smooths.

In the first set of plots, the figure in the top-left is the estimated smooth for the reference level A. The other panels show the estimated deviance/difference smooths, indicating how the stated level smoothly differs from the smooth for the reference level A.

These models are somewhat sensitive to the choice of reference level. Here I would argue that from the plots of the data, any level but A should be the reference, as all the other levels are much closer in shape/magnitude to one another compared to A.

FYI, something I hadn't appreciated sufficiently when I wrote that blog post is that you can reset the contrasts on the ordered factor such that it uses the default treatment contrasts and thus acts like a normal factor though because it is still ordered, gam does the fancy ANOVA-like decomposition for the smooths. If you were to add this line to the end of your example data code:

contrasts(df\$Treatmentord) <- "contr.treatment"


then use this factor for both the parametric term and the by factor:

gam(Expression ~ Treatmentord +
s(Time) +
s(Time, by = Treatmentord),
data = df,
method = "REML")


you will get a more useful set of parametric terms.

I'll see what can be done to make it clearer in the plots what each smooth refers to.

• Thank you for the answer! I actually tried changing the code to: x<-seq(0,20, by=1.7) so that each timecourse had 12 data points instead of 10 and then the smooth from the ordered factor model and the unordered factor model looked exactly the same. So in the ordered factor model is the reference level borrowing some information from the other levels? I guess I had imagined the model fitting that one independently of the others and then fitting the difference smooths to the reference. Commented Jul 22 at 20:32
• I think the basis for the reference smooth is formed using all the observations across all levels of the factor — likewise that for the difference smooths — because that seems most sensible and necessary and is how other by smooths work. The smooths are estimated independently (given the other terms in the model) but the basis used is using all the data. Commented Jul 23 at 6:08
• That makes sense. In the model with the unordered factor and no s(Time) term is there a difference in how the basis is formed as compared to the ordered factor model with the s(Time) term? I'm still a little bit confused about why the results for the smooths are different between the two models. Ideally I would want to be able to compare the smooths from each level B-G to level A using the ordered factor approach but maybe it's not possible here. Seems like it's time for me to read that GAM textbook by Simon Wood! Also, I really appreciate all the blog and video resources you've published Commented Jul 23 at 17:35
• The basis for the factor by version is also made from all the data across all levels. These models differ because in the factor by model, we are actually modelling the individual smooths, while for the ordered by we only model the actual smooth for the reference level; the smooths for other levels are a linear combination of the weighted basis fuctions for two smooths. There's another infelicity in gratia here; as it doesn't really know that the s(time) is for the reference level, when you add partial residuals it adds the data for all levels. Commented Jul 26 at 8:47