10
$\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
$\endgroup$
6
$\begingroup$

It is a lot easier to do using the lsmeans package

library(lsmeans)
lsmeans(mod, pairwise ~ tension | wool)
lsmeans(mod, pairwise ~ wool | tension)
$\endgroup$
  • $\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 '16 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 '16 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 '17 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 '17 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 '17 at 13:11

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.