# Model Comparison of GAMs using gamm4

I am using a GAM from package gamm4 in R to fit a varying coefficient model for longitudinal data. I have data consisting of roughly 400,000 total observations of 100,000 subjects (averaging roughly 4 observations times per subject), including a continuous response variable that changes over time, and two binary categorical/grouping variables which are time-invariant. I am interested in the time-varying effects of the two grouping variables and would like to include random intercepts for each subject.

I have two main questions:

1. When I first fit the model treating the binary covariates as numerical (0 or 1 valued) I get a somewhat sensible answer. Once I converted them to factors the model fit was totally different though, and I got a message that the fixed effect model matrix was rank deficient. I tried throwing a -1 into the model formula but that didn't solve the problem (I realize this addresses the constant intercept and the issue lies in the smooth terms). Is there a way to drop a smooth intercept, or another way around this issue? And is the model using the numeric-coded variables correct?
2. What is the best way to compare gamm4 models? I have seen some discussion on here surrounding comparison of various mgcv models but not much on gamm4. For instance, it seems that AIC does not work, and I am unsure if anova is reliable.

I can't share the actual data, but here is a small example:

library(mgcv)
library(gamm4)
set.seed(050318)

N <- 1000
ni <- 7
amax <- 20
fac1 <- factor(rep(sample(c("f1","f2"),N,T),each = ni))
fac2 <- factor(rep(sample(c("g1","g2"),N,T),each = ni))
id <- factor(rep(1:N,each = ni))
t <- rep(NA,N*ni)
for(i in 1:N){
t[(ni*(i-1)+1):(ni*i)] <- sort(sample(1:amax,ni,F))
}
y <- 1 + 0.4*(fac1=="f2") + 0.5*(1/log(t + 1)) +
0.3*(fac2=="g2")*(1/(log(t + 1))^2) + rnorm(ni*N,0.1,0.2)

# data frame with grouping variables coded as factors
df <- data.frame(y = y, id = id, t = t, fac1 = fac1, fac2 = fac2)

# data frame with grouping variables coded as numeric, 0/1
df2 <- df
df2$$fac1 <- as.numeric(df2$$fac1) - 1
df2$$fac2 <- as.numeric(df2$$fac2) - 1

# Compare model using factor to model using numeric grouping variables
m0 <- gamm4(y ~ s(t,by = fac1) + s(t,by = fac2),random = ~ (1|id),data = df)
summary(m0$$mer) summary(m0$$gam) # Error in if(pval > 0.5) {: missing value where TRUE/FALSE
needed}
plot(m0$gam,pages = 1) m0.1 <- gamm4(y ~ s(t) + s(t,by = fac1) + s(t,by = fac2), random = ~ (1|id),data = df2) summary(m0.1$$mer) summary(m0.1$$gam) plot(m0.1$gam,pages = 1)

# Compare these models
m1 <- gamm4(y ~ fac1 + fac2 + s(t,by = fac1) + s(t,by = fac2),
random = ~ (1|id),data = df)
summary(m1$$mer) summary(m1$$gam) # Error in if(pval > 0.5)
plot(m1$gam,pages = 1) m2 <- gamm4(y ~ fac1 + fac2 + s(t,by = fac2), random = ~ (1|id),data = df) summary(m2$$mer) summary(m2$$gam) plot(m2$gam,pages = 1)

anova(m1$$gam,m2$$gam) # Error in if(pval > 0.5)


EDIT: One work-around that seems to answer question 1 is to group fac1 and fac2 together into one 4-level variable, say newfac. For example:

df$$newfac <- "a" df$$newfac[which(df$$fac1=="f2" & df$$fac2=="g1")] <- "b"
df$$newfac[which(df$$fac1=="f1" & df$$fac2=="g2")] <- "c" df$$newfac[which(df$$fac1=="f2" & df$$fac2=="g2")] <- "d"
df$$newfac <- as.factor(df3$$newfac)

m3 <- gamm4(y ~ newfac + s(t,by = newfac),random = ~ (1|id),data = df)
summary(m3$$mer) summary(m3$$gam)

However I find the interpretation of the smooths estimated by model m0.1 to be more natural. In that model, the smooth for fac1 gives the contrast between f2 and f1, and the smooth for fac2 gives the contrast between g2 and g1, which is really what I'm looking for. So is there anything lost but using a numeric coding of the variables as in m0.1? And is there a simpler way to do this with factors than by combining multiple factors into a new factor variable?