2
$\begingroup$

I'm sorry if this question has been clarified in another post. I've looked around for some time and have been unable to find an answer.

I'm conducting research to evaluate heart function before and after a surgical procedure. Subjects, divided into three conditions, have two measurements made: one pre- and another post-op. I'm trying to determine first whether heart function changes between the conditions with an ANOVA and then to quantify those changes with selected post-hoc pairwise comparisons.

A statistical consultant recommended I use a linear mixed-effects model from the lme4 package, but they haven't been able to clarify to me how it differs from a simple lm model. It does not seem to make any difference whether I run the ANOVA using lme4 or lm function.

How can I correctly estimate these group differences accounting for repeated measures?

Here are my data and R code:

library(tidyverse)
library(ggpubr)
library(lme4)
library(lmerTest)
library(emmeans)

# Reading in the data.
Echo_Long <- tibble(Subject = c("5801","5801","5802","5802","5803","5803","5804","5804","5805","5805","5806","5806","5807","5807","5808","5808","5809","5809","5810","5810","5817","5817","5818","5818","5819","5819","5820","5820","5821","5821","5822","5822","5823","5823","5824","5824","5871","5871","5872","5872","5873","5873","5874","5875","5875","5876","5876","5877","5878","5878","5879","5879","5880","5881","5881","5882","5882","5883","5883","5884","5884","5885","5885","5886","5886","5887","5887","5888","5888","5889","5889","5890","5890","5891","5891","5892","5892","5893","5893","5894","5894","5895","5895","5896","5896"),
                    Condition = factor(c("CLP","CLP","CLP","CLP","Sham","Sham","Control","Control","CLP","CLP","Sham","Sham","Control","Control","Control","Control","CLP","CLP","Sham","Sham","CLP","CLP","Sham","Sham","CLP","CLP","Control","Control","CLP","CLP","CLP","CLP","Sham","Sham","Control","Control","CLP","CLP","CLP","CLP","Sham","Sham","Control","CLP","CLP","Sham","Sham","Control","CLP","CLP","Sham","Sham","Control","CLP","CLP","CLP","CLP","Sham","Sham","Control","Control","CLP","CLP","Sham","Sham","Sham","Sham","Control","Control","CLP","CLP","Control","Control","CLP","CLP","Sham","Sham","CLP","CLP","Sham","Sham","Control","Control","Control","Control"),
                                       levels = c("Control","Sham","CLP")),
                    PrePost = factor(c("Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Pre","Post","Pre","Post","Pre","Pre","Post","Pre","Post","Pre","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post","Pre","Post"),
                                     levels = c("Pre","Post")),
                    FS = c(25.856713,30.169482,30.909077,30.546638,41.328412,36.42172,32.867139,34.471001,39.436617,39.626171,39.71118,30.272102,29.824546,42.21454,41.221393,32.624105,35.451506,47.49036,38.589217,30.612237,30.258312,33.224753,32.706771,31.249994,31.97492,30.344834,32.423221,41.218641,31.329113,35.438585,33.892605,42.295094,35.877867,31.561452,45.756452,25.751061,36.42385,39.676124,39.862529,32.208585,34.982328,31.76472,44.308943,36.781606,41.947568,32.989691,36.824324,24.303333,42.248064,32.713739,41.155224,47.767851,48.9712,31.598498,33.823523,28.861782,36.666664,34.96504,30.744348,47.15448,29.757777,29.823996,44.53442,40.30419,27.722751,32.996624,24.172197,31.929807,39.350181,27.984326,20.588239,28.040558,33.436527,37.818172,38.790033,28.807946,27.814569,27.586225,44.047609,37.282238,27.814569,39.393924,34.432243,37.596906,25.964909))

# Building the figure.
Echo_Long %>%
  ggplot(aes(x = PrePost, y = FS, fill = interaction(PrePost, Condition)), show.legend = FALSE) +
  stat_boxplot(geom = "errorbar", width=0.3, position = position_dodge(1), lwd=1, show.legend = FALSE) +
  geom_boxplot(position = position_dodge(1), outlier.shape = NA, color = "black", lwd=1, show.legend = FALSE) +
  geom_point(aes(x = PrePost, y = FS, color = Condition), shape=21, alpha = 0.6, color="black", size=3,
             show.legend = FALSE) +
  geom_path(aes(x = PrePost, y = FS, group = interaction(Condition, Subject)), color="grey", show.legend = FALSE) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, (Echo_Long$FS %>% max())*1.05)) +
  scale_fill_brewer(palette="Paired") +
  facet_wrap(~ Condition) +
  ggtitle("Fractional Shortening")

# A paired t-test gives a lower p value on a subset of the data.
TTestResults <- Echo_Long %>% filter(Condition == "Sham") %>% t.test(FS ~ PrePost, data = .)
PairedTTestResults <- Echo_Long %>% filter(Condition == "Sham") %>% t.test(FS ~ PrePost, data = ., paired = TRUE)

TTestResults
PairedTTestResults

# I find no difference between these models.
Mod1 <- Echo_Long %>%  
  lmer(FS ~ Condition*PrePost + (1|Subject), data = .)
Mod2 <- Echo_Long %>% 
  lm(FS ~ Condition*PrePost, data = .)

anova(Mod1)
anova(Mod2)
anova(Mod1, Mod2)

Results1 <- Mod1 %>% emmeans(list(pairwise ~ Condition + PrePost), adjust = "none")
Results2 <- Mod2 %>% emmeans(list(pairwise ~ Condition + PrePost), adjust = "none")  

Results1
Results2

Thank you.

$\endgroup$

1 Answer 1

2
$\begingroup$

If you check the warning posted by R after fitting your lmer model, Mod1, you'll notice the dreaded singular fit message. You can google that error message to learn more about it but suffice it to say that you should not trust a model which comes with this warning. It seems that there's virtually no variation among the random intercepts in your model so this throws R off. (This is likely the reason you are concluding that you can't see a difference between your lmer and lm models.)

I would like to suggest that you consider fitting your models using the gamlss package of R so that you can evade the singular fit warning produced by lmer and compare the various competing models.

Here is the R code you need for this:

library(gamlss)

Echo_Long$Subject <- factor(Echo_Long$Subject)


Model1  <- gamlss(FS ~ re(fixed = ~ Condition*PrePost, 
                          random = ~ 1|Subject), 
                  data = Echo_Long, 
                  family = NO) 

summary(getSmo(Model1))


Model2 <- gamlss(FS ~ Condition*PrePost, 
                 data = Echo_Long, 
                 family = NO)

summary(Model2) 

GAIC(Model1, Model2)

The model with the lowest (generalized) AIC value should be preferred.

Model2 is a linear regression model which assumes independence of FS values within the same subject as well as across different subjects. (NO stands for the Normal distribution.)

Model1 is a linear mixed effects model which assumes that the FS values within the same subject are possibly correlated (e.g., if the FS value for that subject is high at Pre, it will tend to also be high at Post if the within-subject correlation of FV values is positive; this correlation is assumed to be the same for all 3 conditions).

If you compare the confidence intervals for the fixed effects of interest in both models, you should expect the linear mixed effects model to produce wider confidence intervals because it adequately reflects that you have less information coming from two correlated values of FS coming from the same subject than from two independent values.

intervals(getSmo(Model1))

confint(Model2)

You can plot model diagnostics for the two models - it seems like the linear mixed effects model diagnostics look a bit better than the linear model ones.

plot(Model1) 
plot(Model2)

Some people would say that if there is no variation in the random intercepts, then the random intercept for subject should be excluded from the model. But I think there is no harm in keeping it in the model provided you can fit the model without getting any warning messages from R. The model should reflect the study design as close as possible.

I don't think the emmeans recognizes gamlss models so the post-hoc comparisons of interest would have to be coded manually if you decided to stay within the gamlss framework. (Recently, I posted several answers on gamlss on this forum - if possible, you can check them out for further clues.)

$\endgroup$
2
  • 2
    $\begingroup$ I think emmeans can handle some gamlss models, but not ones involving splines and such. $\endgroup$
    – Russ Lenth
    Oct 17, 2020 at 13:01
  • $\begingroup$ This is great. Thank you very much! $\endgroup$
    – SStandage
    Oct 22, 2020 at 16:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.