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?