Post-Hoc Test on lmer Model – How to Perform and Interpret

lme4-nlmepost-hocr

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?

Best Answer

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")

The next option is difflsmeans. 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 default by emmeans, so it might be best to explicitly specify the method you prefer.

library(lmerTest)
difflsmeans(model, test.effs = "Group", ddf="Kenward-Roger")

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)
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."

Related Solutions

Solved – What does the anova() command do with a lmer model object

Use the Source, Luke. We can peek inside the ANOVA function by doing getAnywhere(anova.Mermod). The first part of that function is for comparing two different models. The anova on the fixed effects comes in the big else block in the second half:

 dc <- getME(object, "devcomp")
        X <- getME(object, "X")
        asgn <- attr(X, "assign")
        stopifnot(length(asgn) == (p <- dc$dims[["p"]]))
            ss <- as.vector(object@pp$RX() %*% object@beta)^2
        names(ss) <- colnames(X)
        terms <- terms(object)
        nmeffects <- attr(terms, "term.labels")[unique(asgn)]
        if ("(Intercept)" %in% names(ss)) 
            nmeffects <- c("(Intercept)", nmeffects)
        ss <- unlist(lapply(split(ss, asgn), sum))
        stopifnot(length(ss) == length(nmeffects))
        df <- vapply(split(asgn, asgn), length, 1L)
        ms <- ss/df
        f <- ms/(sigma(object)^2)
        table <- data.frame(df, ss, ms, f)
        dimnames(table) <- list(nmeffects, c("Df", "Sum Sq", 
            "Mean Sq", "F value"))
        if ("(Intercept)" %in% nmeffects) 
            table <- table[-match("(Intercept)", nmeffects), 
                ]
        attr(table, "heading") <- "Analysis of Variance Table"
        class(table) <- c("anova", "data.frame")
        table

object is the lmer output. We start to calculate the sum of squares in line 5: ss <- as.vector ... The code multiplies the fixed parameters (in beta) by an upper triangular matrix; then squares each term. Here is that upper triangular matrix for the irrigation example. Each row corresponds to one of the five fixed effects parameters (intercept, 3 degrees of freedom for irrigation, 1 df for variety).

zapsmall(irrigation.lmer@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]  [,5]
[1,] 0.813 0.203  0.203  0.203 0.407
[2,] 0.000 0.352 -0.117 -0.117 0.000
[3,] 0.000 0.000  0.332 -0.166 0.000
[4,] 0.000 0.000  0.000  0.287 0.000
[5,] 0.000 0.000  0.000  0.000 2.000

The first row gives you the sum of squares for the intercept and the last gives you the SS for the within-fields variety effect. Rows 2-4 only involve the 3 parameters for the irrigation levels, so pre-multiplication gives you three pieces of the SS for irrigation.

These pieces are not interesting in themselves since they come from the default treatment contrast in R, but in line ss <- unlist(lapply(split .... Bates scoops up bits of sums of squares according to the number of levels and which factors they refer to. There's a lot of book-keeping going on here. We get the degrees of freedom as well (which are 3 for irrigation). Then, he gets the mean squares which show up on the printout of anova. Finally, he divides all his means squares by the within-groups residual variance, sigma(object)^2.

So what's going on? The philosophy of lmer has nothing to do with the method of moments approach used by aov. The idea in lmer is to maximize a marginal likelihood obtained by integrating out the unseen random effects. In this case, the random fertility level of each field. Chapter 2 of Pinheiro and Bates describes the ugliness of this process. The RX matrix used to get the sums of squares is their matrix $R_{00}$ from equation 2.17, page 70 of the text. This matrix is obtained from a bunch of QR decompositions upon (amongst other things) the design matrix of the random effects and $\sqrt{\sigma^2/\sigma_f^2}$, where $\sigma_f^2$ is the variance of the field effect. This is that missing factor you were asking about, but it does not enter into the solution in a transparent or simple way.

Asymptotically, the estimates of the fixed effects have distribution:

$$ \hat{\beta} \sim \mathcal{N}(\beta, \sigma^2[R_{00}^{-1}R_{00}^{-T}]) $$

which means that the components of $R_{00}\hat{\beta}$ are independently distributed. If $\beta=0$, the squares of those terms (divided by $\sigma^2$) follow a chi-squared distribution. Dividing through by the estimate of $\sigma^2$ (another chi-squared when divided by $\sigma^2$) gives an F statistic. We don't divide by the mean square error for the between-groups analysis because that has nothing to do with what's going on here. We need to compare the sums of squares obtained through $R_{00}$ by comparing them to an estimate of $\sigma^2$.

Note that you would not have got the same F statistics had the data been unbalanced. Nor would you have obtained the same F statistics if you had used ML instead of REML.

The idea behind aov is that the expected mean squares for irrigation is a function of $\sigma^2$, $\sigma_f^2$ and the irrigation effects. The expected mean square for field residuals is a function of $\sigma^2$ and $\sigma_f^2$. When the irrigation effects are 0, both of these quantities estimate the same thing and their ratio follows an F distribution.

Interestingly, Bates and Pinheiro recommend using the ANOVA over fitting two models and doing a likelihood ratio test. The latter tends to be anti-conservative.

Of course, if the data are unbalanced, it is no longer clear exactly what hypothesis the ANOVA is testing. I removed the first observation from the irrigation data and refit. Here is the $R_{00}$ matrix again:

zapsmall(fit2@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]   [,5]
[1,] 0.816 0.205  0.205  0.205  0.457
[2,] 0.000 0.354 -0.119 -0.119 -0.029
[3,] 0.000 0.000  0.334 -0.168 -0.040
[4,] 0.000 0.000  0.000  0.288 -0.071
[5,] 0.000 0.000  0.000  0.000  1.874

As you can see, the sums of squares for the irrigation parameters now contain some of the variety effect as well.

Solved – Addressing “NOTE: Results may be misleading due to involvement in interactions” warning with Tukey post-hoc comparisons in lsmeans R package

My view is that the $F$ test of statistical significance of the interaction effect is less important than the subjective nature of the interaction, as evidenced by the plot. The plot tells me that it is reasonably sensible to compare the overall averages of Depression and Top, but it'd be silly to compare those averages with the overall average of Slope -- whether or not these comparisons are statistically significant. Basically, I'd say to avoid doing comparisons that don't make sense -- so my advice is do not ignore the warning note in this case. If the curve for Top were fairly parallel with the other two, that's when you could ignore it.

In general, I suggest looking at enough plots that you can tell what's going on, and then restrict your post-hoc testing to things that are sensible.

Since P is continuous, you're really fitting straight lines (they look curved because you chose unequally spaced points). You can compare the slopes of these lines:

R> lstrends(Dens.LMER, pairwise ~ Contour, var = "P")

$lstrends
 Contour        P.trend          SE    df    lower.CL     upper.CL
 Depression -0.00681143 0.004901195 39.68 -0.01671957  0.003096714
 Slope      -0.03376293 0.010533875 41.88 -0.05502295 -0.012502911
 Top        -0.01306992 0.010499548 41.97 -0.03425936  0.008119525

Confidence level used: 0.95 

$contrasts
 contrast               estimate         SE    df t.ratio p.value
 Depression - Slope  0.026951501 0.01161827 42.00   2.320  0.0639
 Depression - Top    0.006258486 0.01158716 41.81   0.540  0.8520
 Slope - Top        -0.020693015 0.01487290 41.99  -1.391  0.3545

P value adjustment: tukey method for a family of 3 tests

The comparison between the shallowest and largest slopes has an adjusted $P$ value of about $.06$.

Best Answer

Related Solutions

Solved – What does the anova() command do with a lmer model object

Solved – Addressing “NOTE: Results may be misleading due to involvement in interactions” warning with Tukey post-hoc comparisons in lsmeans R package

Related Question