# A correct mixed model for longitudinal data

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.

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.

• Thank you very much. Regarding the AR(1) structure, would the following model be correct? lme(weight ~ Genotype + time, random = ~ 1|ID, correlation=corCAR1() data=data) I am actually impressed by how much the p-value of time is changed by adding the AR(1) structure. Oct 15, 2020 at 8:28
• Missing a comma of course after corCAR1(). And regarding the non-linearity, I must admit my example was formulated in an unfortunate way because it clearly indicated a quadratic relation. My apologies. If you see data in which there really isn't any simple transformation to achieve linearity, what would you recommend then? Or is it maybe a better option then to stick away from modeling the complete longitudinal data, and rather focus on specific (predefined) time points without consideration of the course over time? Oct 15, 2020 at 8:33
• Yes, that looks right, although questions about software are off topic here. As for linearity, I'm not suggesting a transformation (although that is another option), rather the addition of nonlinear terms. Splines would be a logical next step. Also, try not to be too concerned with p values. Oct 15, 2020 at 8:52
• Thanks! I always try to convince my clinical colleagues not to use p-values but it's an extremely difficult discussion. Often I achieve that they at least add confidence intervals to their papers :) Oct 15, 2020 at 10:28
• Ahh, I can feel your pain. I have been there. The struggle is real ! :) Oct 15, 2020 at 11:44