A correct mixed model for longitudinal data

Question

This is a rather basic question and there are many similar questions here on CV. However, I could not find answers to my specific three questions below. My apologies, if I may have overseen something.

Quite often, I see the following experimental setup: "we have measured weight over time in a number of mice, and we want to know whether there exists a difference between the genotypes, i.e. between knockout (KO) and wild type (WT)."

Usually, I think of a simple linear mixed effects model with a random intercept for each individual. Or, if they suspect different time trends, a random slope as well, which often might be a bit too much due to the rather small data sets. Further, I point them towards nlme since they won't stop asking for p-values ...

However, I'm afraid this model is actually too simplified. My concerns are:

1. Since my weight at one day depends on my weight of the day before, do I need to take some kind of autocorrelation into account? Is maybe some kind of ARIMA more appropriate?
2. A colleague of mine told me to always use the measurement at baseline time as covariate. While I see that this makes sense in a standard regression model, I'm wondering whether this is required when I already have a random intercept which is supposed to cover differences at baseline?
3. What to do if lineary does not seem to hold at all? For example, due to a strong peak in the middle with equally strong descend afterwards. Again, a colleague modelled time with restricted cubic splines (rcs) to account for this, but I'm afraid it may be an overkill to make inference from a mixed effects model in which each individual has its own rcs ... however, I do not have a better suggestion.

---

Just as illustration, an example for the first two questions:

library(dplyr)
library(nlme)
library(ggplot2)

# simulate some data to get started
set.seed(16)
n=10
baseline <- c(rnorm(n=n,mean=20,sd=1),rnorm(n=n,mean=23,sd=1))
data <- data.frame("weight"=baseline, time=0, "Genotype"=c(rep("WT",n), rep("KO",n)), "ID"=1:(2*n))
for(i in 1:20){
  next.week <- baseline + rnorm(n=20,mean=rnorm(1),sd=3)
  new.data <- data.frame("weight"=next.week, time=i, "Genotype"=c(rep("WT",n), rep("KO",n)), "ID"=1:(2*n))
  data <- rbind(data, new.data)
  baseline <- next.week
}

# plotting the data
data %>% group_by(time, Genotype) %>% 
  summarise(mean.weight = mean(weight, na.rm=T)) %>%
  ggplot(., aes(x=time, y=mean.weight, colour=Genotype)) + 
  geom_line(aes(y = mean.weight, x=time)) + 
  theme_bw() + scale_y_continuous(limits = c(15,40)) + 
  xlab("Time") + ylab("Weight in g") 

The data often looks a bit messy and not perfectly linear, although the assumption of linearity does not seem completely off since one would expect a more or less linear increase of weight over time in general.

# a simple random intercept model
lme.model <- lme(weight ~ Genotype + time, random = ~ 1|ID, data=data)
summary(lme.model)

Random effects:
 Formula: ~1 | ID
        (Intercept) Residual
StdDev:    9.203191   6.1417

Fixed effects: weight ~ Genotype + time 
                Value Std.Error  DF   t-value p-value
(Intercept) 22.958956  2.982354 399  7.698268  0.0000
GenotypeWT  -4.040438  4.159205  18 -0.971445  0.3442
time         0.245306  0.049491 399  4.956556  0.0000

We see a significant increase of weight over time, and no significant difference between the genotypes (due to the small sample size in this simulated example). From the random effects, we could calculate the ICC and see that there is a lot of variance captured by the random intercept, so we should not use a model without a random intercept.

Answer 1

Since my weight at one day depends on my weight of the day before, do I need to take some kind of autocorrelation into account? Is maybe some kind of ARIMA more appropriate?

An AR(1) structure is often a very good idea in these types of models, of the data support such a model, although in practice there is often insufficient data to estimate the AR parameter. So I would suggest trying it.

A colleague of mine told me to always use the measurement at baseline time as covariate. While I see that this makes sense in a standard regression model, I'm wondering whether this is required when I already have a random intercept which is supposed to cover differences at baseline?

Yes, your intuition is correct here.

What to do if lineary does not seem to hold at all? For example, due to a strong peak in the middle with equally strong descend afterwards. Again, a colleague modelled time with restricted cubic splines (rcs) to account for this, but I'm afraid it may be an overkill to make inference from a mixed effects model in which each individual has its own rcs ... however, I do not have a better suggestion.

Looking at your plot, I don't see a big problem with linearity, but if you do have a clear non-linear pattern ("trong peak in the middle with equally strong descend afterwards"), then a good starting point is to simply add a quadratic term - you may also want to consider centering the data in this case.