Solved – How to report and interpret the output from linear mixed models with interaction terms

experiment-designgroup-differencesinterpretationlme4-nlmemixed model

So, we are doing a linear mixed effects model for analyzing some results of our study. In short, we have performed two different meal tests (i.e., two groups), and measured the response in various biomarkers at baseline as well as 1, 2, 3, and 4 hours after the meal.

I had a meeting with a statistician who explained that we should use linear mixed models for this and as such, using the nlme package in R the syntax looks like this:

model <-lme(biomarker~ as.factor(group)*visit, random = ~1|ID, data=data, method="ML")
summary(model)

The output (abbreviated for readability):

Linear mixed-effects model fit by maximum likelihood
 Data: data 
      AIC      BIC    logLik
  137.593 149.0651 -62.79649

Random effects:
 Formula: ~1 | ID
        (Intercept)  Residual
StdDev:    1.462879 0.6039689

Fixed effects: biomarker ~ as.factor(group) * visit 
                            Value Std.Error DF   t-value p-value
(Intercept)              7.869766 0.7157143 38 10.995681  0.0000
as.factor(group)3        1.295118 1.0121729  8  1.279542  0.2366
visit                   -0.096024 0.0679003 38 -1.414191  0.1654
as.factor(group)3:visit -0.358905 0.0960255 38 -3.737606  0.0006
 Correlation: 
                        (Intr) as.()3 visit 
as.factor(group)3       -0.707              
visit                   -0.247  0.174       
as.factor(group)3:visit  0.174 -0.247 -0.707

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.107751422 -0.303320567  0.004573801  0.377750437  1.967646127 

Number of Observations: 50
Number of Groups: 10

To summarize:

Exposure = one of two meal tests (group in the syntax)
Outcome = Biomarker
Time variable = Visit (5 in total for each participant, continuous)

My questions are:

Am I correct in interpreting this that there is an interaction between group and visit?
I am unclear as to how I should interpret the estimates here. Am I correct in saying that at time = 0, then the group difference (3 vs. 2) is 1.29? And further that this effect depends on the visit? What about the other timepoints?
Is it sufficient to report this model or should we also include a model without the interaction term that is just including group and visit?

Best Answer

1) Am I correct in interpreting this that there is an interaction between group and visit?

Yes.

2) I am unclear as to how I should interpret the estimates here. Am I correct in saying that at time = 0, then the group difference (3 vs. 2) is 1.29? And further that this effect depends on the visit? What about the other timepoints?

Yes, that is the group difference at time = 0 and yes, it depends on time. You can calculate the estimate at any combination of the variables by using the formula:

$7.87 + 1.30*(I(\text{group} = 3)) - 0.10*\text{visit} - 0.36*(I(\text{group} = 3))*\text{visit}$

3) Is it sufficient to report this model or should we also include a model without the interaction term that is just including group and visit?

That depends on what you are interested in, but when there is an interaction, the model with only main effects can be misleading.

Related Solutions

R Mixed Model – Interaction Term in a Linear Mixed Effect Model in R

Here's what I would do:

First, I would have a look here on how to specify the random term in your model1. I am not quite sure what you are trying to fit. There is also a lot of info on linear mixed effects models here on CV. Click on the lme4-nlme tag, which you also provided. It would also help if you could provide an example dataset, or at least the structure of your data.

Then, you most likely only need one model, which is presumably in the form of:

my_model <- lmer(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)

I specified the random effect to be + (1|site), because you said:

Study sites are included as the random effect in the model.

To get the ANOVA table you can either do:

library(car)
Anova(my_model)

or:

library(afex)
mixed(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)

or instead of running lmer() through the lme4 package, load the lmerTest package and run:

my_model <- lmer(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)
anova(my_model)

This will give you the ANOVA table you probably need eventually. Make sure to have a look at those functions and their arguments (?Anova, ?mixed, ?lmerTest::anova).

I don't quite understand why would want to exclude species if the interaction is significant and run separate models for all species?!

However, if your main effects are not significant you could consider tossing them out and re-running the model with the interaction only. However, if one or both main effects are significant, I would keep them both in the model and report this together with a potential significant interaction.

In any case, if you have a significant interaction you should focus on interpreting the interaction and not the main effects since their interpretation could now be misleading. The interpretation of the interaction should start by visualizing it. You could do this for example using the emmip() function in the emmeans package:

library(emmeans)
emmip(my_model, landuse ~ species)

Regarding the adjustment of p-values, you only need to do that if you are following up with post-hoc tests.

This could be done with the emmeans() function (also from the emmeans package):

emmeans(my_model, pairwise ~ species : landuse)

Model Selection in Longitudinal Data – Testing the Need for Random-Effects Terms in Longitudinal Data Analysis

The likelihood ratio test is slightly incorrect (in general, conservative) for testing the significance of a random effect, because the null value ($\sigma^2=0$) is at the boundary of the feasible space, but in this case there is overwhelmingly strong evidence against the null hypothesis. The model with random effects of individual is 15713-6772=8941 log-likelihood units better; twice the log-likelihood value is $\chi^2$ distributed, so the direct p-value calculation would give you ...

pchisq(2*8941,df=1,lower.tail=FALSE,log.p=TRUE)/log(10)
## -3885.251

... a p-value of approximately $10^{-3885}$.

You should really consider a random-slope model (random = ~time|id) as well.

Update: relative to the random-intercept model, the random-slopes model is again much better. The improvement is now 935 log-likelihood units, which doing the equivalent calculation as above corresponds to a rejection of the null hypothesis (among-individual variation in slope is equal to zero) with a p-value of "only" $10^{-408}$.

Best Answer

Related Solutions

R Mixed Model – Interaction Term in a Linear Mixed Effect Model in R

Model Selection in Longitudinal Data – Testing the Need for Random-Effects Terms in Longitudinal Data Analysis

Related Question