Mixed effects model with splines

Question

I am fitting a mixed effects model with a spline term in an application where the trend over time is known to be curvi-linear. However, what I would like to assess is whether the curvi-linear trend occurs due to individual deviation from linearity, or is it an effect at the group level that makes a group level fit appear curvi-linear. I give a reproducible example boring a dataset from the JM package.

library(nlme)
library(JM)
data(pbc2)

fitLME1 <- lme(log(serBilir) ~ ns(year, 2), random = ~ year | id, data = pbc2)
fitLME2 <- lme(log(serBilir) ~ year, random = ~ ns(year, 2) | id, data = pbc2)

Essentially I want to know which one of these better fits my data. However comparison by anova gives me an ominous warning:

        Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fitLME1     1  7 3063.364 3102.364 -1524.682                        
fitLME2     2  9 2882.324 2932.472 -1432.162 1 vs 2 185.0399  <.0001
Warning message:
In anova.lme(fitLME1, fitLME2) :
  fitted objects with different fixed effects. REML comparisons are not meaningful.

Now I am aware that there are difficulties making these kinds of comparisons via maximum likelihood methods - but what is the alternative?

Answer 1

As mdewey says then refit the model without the REML estimation method. As the warning says, comparisons are not meaningful when you have different fixed effects structures.

The next issues is that the models are not nested so the F-test presumably does not makes sense. You could look at the information criteria. Both favor fitLME2.