11

I am performing post-hoc tests on a linear mixed-effects model in R (lme4 package). I am using multcomp package (glht() function) to perform the post-hoc tests.

My experimental design is repeated-measures, with a random block effect. The models are specified as:

mymod <- lmer(variable ~ treatment * time + (1|block), data = mydata, REML = TRUE)

Rather than attaching my data here, I am working off of the data called warpbreaks within the multcomp package.

data <- warpbreaks
warpbreaks$rand <- NA

I have added an extra random variable to mimic my "block" effect:

warpbreaks$rand <- rep(c("foo", "bar", "bee"), nrow(warpbreaks)/3)

This mimics my model:

mod <- lmer(breaks ~ tension * wool + (1|rand), data = warpbreaks)

I am aware of the the example in the "Additional Multcomp Examples- 2 Way Anova" This example leads you to a comparison of levels of tension within the levels of wool.

What if I want to do the opposite - compare the levels of wool within the levels of tension? (In my case, this would be comparing the levels of treatment (two - 0, 1) within the levels of time (three - June, July, August).

I have come up with the following code to do so, but it doesn't seem to work (see error message below).

First, from the example (with wool and tension swapped places):

tmp <- expand.grid(wool = unique(warpbreaks$wool), tension = unique(warpbreaks$tension))
X <- model.matrix(~ tension * wool, data = tmp)
glht(mod, linfct = X)

Tukey <- contrMat(table(warpbreaks$wool), "Tukey")

K1 <- cbind(Tukey, matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(warpbreaks$tension)[1], rownames(K1), sep = ":")

K2 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(warpbreaks$tension)[2], rownames(K2), sep = ":")

From here to bottom, my own code:

K3 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(warpbreaks$tension)[3], rownames(K3), sep = ":")

K <- rbind(K1, K2, K3)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))

> summary(glht(mod, linfct = K %*% X))
Error in summary(glht(mod, linfct = K %*% X)) : 
  error in evaluating the argument 'object' in selecting a method for function 'summary': Error in K %*% X : non-conformable arguments
luchonacho
  • 2,568
  • 3
  • 21
  • 38
Ashley Asmus
  • 113
  • 1
  • 1
  • 5

1 Answers1

7

It is a lot easier to do using the lsmeans package

library(lsmeans)
lsmeans(mod, pairwise ~ tension | wool)
lsmeans(mod, pairwise ~ wool | tension)
Russ Lenth
  • 15,161
  • 20
  • 53
  • Great, it works! Thanks. Note: this code only worked for my data after changing my repeat variable from numerical values (3 & 6) to alphabetical values (A & B). –  Jul 06 '16 at 08:27
  • Well, that matters a lot! Because it's a *different model* with `time` as a numeric predictor. I suspect you wanted it as a factor. – Russ Lenth Jul 06 '16 at 17:07
  • How can I generalize to more predictors? if for example I have 3 predictors how does it work? – have fun Aug 27 '17 at 14:58
  • 1
    @havefun Please look at `help("lsmeans", package = "lsmeans")` and `vignette("using-lsmeans")`. There is a lot of documentation and many examples. – Russ Lenth Aug 27 '17 at 15:05
  • One more question, why didn't you suggest `lsmeans(mod, pairwise ~wool*tension)` ? I tried with my data and the comparisons and estimated values are the same BUT p-values are very different! – have fun Aug 28 '17 at 12:02
  • 1
    Count the number of comparisons you obtain with each method, they are not the same. Also read up on multiple-testing adjustments. When you have a larger family of tests, the adjusted P values are different then for a smaller family. When you use a by variable, the adjustments are applied separately to each set. – Russ Lenth Aug 28 '17 at 13:11
  • I thought it was something in that direction...but which one is the correct way? – have fun Aug 28 '17 at 16:08