# Mixed effects model with splines

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?

• Try re-fitting with method = "ML" or whatever the correct phrase is and see if the warning goes away. – mdewey Dec 4 '16 at 14:17

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.