11
$\begingroup$

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
Improve this question
$\endgroup$

1 Answer 1

Reset to default
8
$\begingroup$

It is a lot easier to do using the lsmeans package

library(lsmeans)
lsmeans(mod, pairwise ~ tension | wool)
lsmeans(mod, pairwise ~ wool | tension)
Improve this answer
$\endgroup$
7
  • $\begingroup$ 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). $\endgroup$
    – user122338
    Jul 6, 2016 at 8:27
  • $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Jul 6, 2016 at 17:07
  • $\begingroup$ How can I generalize to more predictors? if for example I have 3 predictors how does it work? $\endgroup$
    – have fun
    Aug 27, 2017 at 14:58
  • 1
    $\begingroup$ @havefun Please look at help("lsmeans", package = "lsmeans") and vignette("using-lsmeans"). There is a lot of documentation and many examples. $\endgroup$
    – Russ Lenth
    Aug 27, 2017 at 15:05
  • 1
    $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Aug 28, 2017 at 13:11

Your Answer

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

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