Problems with the output of mixed-effect model in R using lmer

Question

sorry if this is kind of a basic question, I’m no expert in stats. I am a PhD student (biologist) working on my thesis and I’ve been dealing with some complicated data that I obtained from an experiment. My animals are fruit flies and I want to study whether their sleep changes after a treatment (sleep deprivation which means not letting them sleep for a whole night) and if that depends on the genotype of the animal.

My animal groups are three different genotypes: "experimental" are the animals which have a genetic manipulation in particular that might affect the effect of sleep deprivation on this animal. To achieve this genetic manipulation I have to use two different genetic lines of flies combined, so I have to define 2 control groups "control 1" and "control 2", each one accounts for one of the lines that I combined to achieve the "experimental" group. So my genotypes are "experimental" (might respond differently to the sleep deprivation treatment) and "control 1" and "control 2" (I expect these two controls to respond in the same way to the treatment, but I'm actually not sure that that will happen). I take half of the flies in each genotype and I apply the treatment (sleep deprivation) during one night and the other half I use it as a control (I don't apply sleep deprivation during the same night). After that, I measure the amount of time they sleep during the first six hours of the day after that night of sleep deprivation or control (depending on the group). I have two data tables, one (named "sleep_before") for the sleep before the moment of sleep deprivation, that has a column with the amount of sleep minutes during the first six hours of the day before the night of sleep deprivation. The other table (named "sleep_after") has a column with the amount of sleep minutes during the first six hours of the day after sleep deprivation.

So for my results I need, for each genotype, to assess whether the sleep deprivation treatment had an effect. So for that, for each genotype, I want to compare the sleep minutes before and after the treatment for sleep deprived and control animals. And then I would like to see statistically if this difference (between sleep after and sleep before) I found between sleep deprived and control for each genotype, is different between genotypes. In other words, I want to see if the genotype affects the way the animals respond to sleep deprivation.

The residuals of my data don’t follow a normal distribution and they don’t meet the homoscedasticity assumption either. I inspected this visually by doing density plots, histograms and qq-plots of the residuals for the distribution and residuals vs. fitted values for the variances. I also performed the Shapiro-wilks test on the residuals for distribution and Levene test for the variance, in both cases p<0.005. I also tried to transform the data by adding a constant (to make all values positive, as I have some negative values in the data) and performing the sqrt (didn’t work), the log10 (didn’t work) and using the boxcox function in R (didn’t work). For all cases data still wasn’t normally distributed.

With some help I ended up using a mixed model that looks like this:

model <- lmer(sleep ~ genotype*treatment*time_point + (1|id), data = sleep_data)

sleep is the number of minutes slept

genotype is a fixed effect of gene variant

treatment is a fixed effect of whether or not that fly was kept awake (sleep deprived or control animals)

time_point is a fixed effect of time, regardless of whether sleep deprivation happened or not (defined as "before" and "after" depending on which table -"sleep_before" of "sleep_after"- the quantified sleep minutes come from, I combined the results in the "sleep_data" table)

After fitting the model I used the summary function from the lmerTest library and I obtained this:

In the table "G4" and "UAS" are genotypes "Control1" and "Control2". "SD" is the group of animals that was treated with sleep deprivation.

I understand that the summary function is taking the data and comparing to a reference value per each factor, but I need to understand whether the treatment had different effects across genotypes. For each genotype, I would need a p-value that shows if there's a difference in "sleep after" compared to "sleep before" in sleep deprived animals compared to animals that were not sleep deprived (control). And then I would like to see statistically if this difference (between sleep after and sleep before) I found between sleep deprived and control for each genotype, is different between genotypes.

I still don't understand how to extract this information from my model. Can someone help me please? I would appreciate it very much.

Answer 1

I think you might succeed with a simpler model, as you seem to have exactly 2 observations per individual. That might even allow you to avoid setting this up as a mixed model, as you then will have something that's a generalization of a paired t-test.

Such experiments sometimes are better handled by modeling just the second observation value, while including the corresponding first observation value as a predictor along with treatment and genotype. Something like:

lm(sleepAfter ~ sleepBefore + genotype*treatment)

or maybe even

lm(I(sleepAfter-sleepBefore) ~ genotype*treatment)

The latter is particularly easy to interpret: you would be modeling the change in sleep between the two observations as a function of genotype/treatment combinations. Those changes should be close to an average of 0 for the flies that weren't sleep deprived.

You would reshape your data into a wide format with one row per individual. With only 1 row per individual, you don't need the extra overhead of random effects. You won't have the troublesome 3-way interaction between timepoint and genotype and treatment. In your current model, it looks like the 3-way interaction will complicate the simple genotype/treatment comparisons that you seek.

Even if you stick with a mixed model, don't depend on a summary() output like you show when factors have more than 2 levels. The displayed coefficients are only for comparisons against the (arbitrary) choice of reference level. They don't directly indicate the overall significance of each factor when all levels are included.

You can get combined estimates of whether there are any differences associated with genotype, treatment, and their interaction with an appropriate anova() function. The standard R anova() function can lead to problems with an lm model if the experimental design isn't perfectly balanced. It's generally safer to use the Anova() function in the R car package if you use lm() as I recommend. If you do anova() on your mixed model of class lmerModLmerTest, however, a specially designed and appropriate function will be used.

For further comparisons between specific combinations of treatments and genotypes, consider the linearHypothesis() function in the car package, or the tools for post-modeling analysis provided by the emmeans package.