18
$\begingroup$

This is my data frame:

Group   <- c("G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3")
Subject <- c("S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15","S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15","S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15")
Value   <- c(9.832217741,13.62390117,13.19671612,14.68552076,9.26683366,11.67886655,14.65083473,12.20969772,11.58494621,13.58474896,12.49053635,10.28208078,12.21945867,12.58276212,15.42648969,9.466436017,11.46582655,10.78725485,10.66159358,10.86701127,12.97863424,12.85276916,8.672953949,10.44587257,13.62135205,13.64038394,12.45778874,8.655142642,10.65925259,13.18336949,11.96595556,13.5552118,11.8337142,14.01763101,11.37502161,14.14801305,13.21640866,9.141392359,11.65848845,14.20350364,14.1829714,11.26202565,11.98431285,13.77216009,11.57303893)

data <- data.frame(Group, Subject, Value)

Then I run a linear-mixed effects model to compare the 3 Groups' difference on "Value", where "Subject" is the random factor:

library(lme4)
library(lmerTest)
model <- lmer (Value~Group + (1|Subject), data = data)
summary(model)

The results are:

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept) 12.48771    0.42892 31.54000  29.114   <2e-16 ***
GroupG2     -1.12666    0.46702 28.00000  -2.412   0.0226 *  
GroupG3      0.03828    0.46702 28.00000   0.082   0.9353    

However, how to compare Group2 with Group3? What is the convention in academic article?

$\endgroup$
16
$\begingroup$

You could use emmeans::emmeans() or lmerTest::difflsmeans(), or multcomp::glht().

I prefer emmeans (previously lsmeans).

library(emmeans)
emmeans(model, list(pairwise ~ Group), adjust = "tukey")

Note difflsmeans cannot correct for multiple comparisons, and uses the Satterthwaite method for calculating degrees of freedom as default instead of the Kenward-Roger method used by emmeans.

library(lmerTest)
difflsmeans(model, test.effs = "Group")

The multcomp::glht() method is described in the other answer to this question, by Hack-R.

Also, you can get the ANOVA p-values by loading lmerTest and then using anova.

library(lmerTest)
anova(model)

Just to be clear, you intended for the Value to be assessed three times for each subject, right? It looks like Group is within-subjects, not between-subjects.

$\endgroup$
  • 1
    $\begingroup$ I just want to add to the response of Kayle Sawyer that the package lsmeans is being deprecated in favor of emmeans. $\endgroup$ – Downhiller Jun 14 '18 at 19:52
  • $\begingroup$ Note if you specify the library, you must use lmerTest::lmer(), not lme4::lmer() for anova() to show the p-values. $\endgroup$ – Kayle Sawyer Jun 28 '18 at 19:20
11
$\begingroup$

After you've fit your lmer model you can do ANOVA, MANOVA, and multiple comparison procedures on the model object, like this:

library(multcomp)
summary(glht(model, linfct = mcp(Group = "Tukey")), test = adjusted("holm"))
   Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = Value ~ Group + (1 | Subject), data = data)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)  
G2 - G1 == 0 -1.12666    0.46702  -2.412   0.0378 *
G3 - G1 == 0  0.03828    0.46702   0.082   0.9347  
G3 - G2 == 0  1.16495    0.46702   2.494   0.0378 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- holm method)

As for the convention in academic papers, that's going to vary a lot by field, journal, and specific subject matter. So for that case just review related articles and see what they do.

$\endgroup$
  • $\begingroup$ Thank you. But which adjustment was actually used? Tukey or holm? Why both appears in the post-hoc test? $\endgroup$ – Ping Tang Sep 29 '16 at 2:41
  • $\begingroup$ @PingTang You're welcome. It's Bonferroni-Holm correction of all-pair multiple comparison. That's just one option, of course. You could also do summary(glht(model, linfct = mcp(Group = "Tukey"))). If you want to see the full academic / statistical descriptions of the various tests that can be performed check out the references in ?glht and multicomp more generally. I think Hsu 1996 would be the main one. $\endgroup$ – Hack-R Sep 29 '16 at 2:54
  • 3
    $\begingroup$ @PingTang , the mcp function, the Group = Tukey just means to compare all pairwise groups in the variable "Group". It doesn't mean a Tukey adjustment. $\endgroup$ – Sal Mangiafico Aug 6 '17 at 17:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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