# Different predictions for GAMMs with ordered vs. unordered factors

I've fitted 2 different binomial GAMM models, both with a smooth for a predictor Timesec that varies by Condition (a categorical factor). In the first, Condition was a nominal, unordered factor:

bam(cbind(FixTarget, NFixations - FixTarget)
~ 0 + Condition
+ s(Timesec, by=Condition, bs="tp")
+ s(Participant, by=Condition, bs="re")
+ s(Participant, Timesec, by=Condition, bs="re"),
data = by.subject,
family = "binomial",
discrete = TRUE)


In the second, Condition was an ordered factor, and as such, I have also included a 'simple' smooth for Timesec:

bam(cbind(FixTarget, NFixations - FixTarget)
~ 1 + Condition
+ s(Timesec, bs="tp") + s(Timesec, by=Condition, bs="tp")
+ s(Participant, bs="re") + s(Participant, by=Condition, bs="re")
+ s(Participant, Timesec, bs="re") + s(Participant, Timesec, by=Condition, bs="re"),
data = by.subject,
family = "binomial",
discrete = TRUE)


The first model estimates separate smooths for Timesec for the levels of Condition, whereas in the second we estimate a smooth for the reference level of Condition and additional difference smooths, relative to the reference level (see GAM factor smooth interaction--include main effect smooth?).

What I didn't expect was that the predictions from these two models differ (see image): in the model with an ordered factor (bottom left), the predictions for the non-reference level has a much wider CI than in the model with an unordered factor (top left).

I was expecting these two models to be equivalent but with different parametrizations. But it seems the uncertainty about the difference smooth in the ordered factor model is carrying to the non-reference level. Why so and what am I doing wrong?

Also note: In both of these models, there were by-Participant random effects. I've fitted separate models without random slopes, in which I removed s(Participant, Timesec, by=Condition, bs="re"). These models show consistent CIs for ordered vs. unordered models (see right column in the image). This seems to be a relevant piece of the puzzle.

Finally, for completeness, predictions were obtained excluding the random effects, as such (after creating a data frame with new data called predictions.df):

mterms <- sapply(m\$smooth, "[[",  "label")
excluded.terms <- mterms[grep("Participant", mterms)]

predictions <- predict(m, newdata = predictions.df,
se.fit = TRUE,
exclude = excluded.terms)


## EDIT

Following Gavin's suggestion below, I have run models with factor smooths, instead of with linear adjustments.

The unordered factor model:

bam(cbind(FixTarget, NFixations - FixTarget)
~ 0 + Condition
+ s(Timesec, by=Condition, bs="tp")
+ s(Timesec, Participant, by=Condition, bs="fs", m=1),
data = by.subject,
family = "binomial",
discrete = TRUE,
nthreads = parallel::detectCores() - 1)


The ordered model:

bam(cbind(FixTarget, NFixations - FixTarget)
~ 1 + Condition
+ s(Timesec, bs="tp") + s(Timesec, by=Condition, bs="tp")
+ s(Timesec, Participant, bs="fs", m=1)
+ s(Timesec, Participant, by=Condition, bs="fs", m=1),
data = by.subject,
family = "binomial",
discrete = TRUE,
nthreads = parallel::detectCores() - 1)


The problem remained. In the ordered model, the reference level still yields smooths with smaller CIs, the non-reference level yields larger CIs:

I'm not convinced you should have random slopes in a model where you also include smooth functions of the covariate that is part of the random slopes. Why not allow for random smooths (say via a bs = "fs" smooth)?