How to compare two levels of one factor

Question

I´ve run a GLMM model with the following variables: Response variable: continuous (with positive, negative and zero values; Explanatory variable: 1 factor with 6 levels (DB, DF, NDB, NDF, BB, BF); individual as a random effect variable

After running a lmer function, and checking the assumptions, I want to perform a posteriori specific comparisons. How do I do that? I want to test for example if there is any difference between DF vs. DB. How do I write this in R??

CODE IS AS FOLLOWS

#Model:
m1 <- lmer(Vueltasmin ~ Condicion + (1 | Bicho), Datos)
summary(m1)

#Checking asusmptions: OK 

#Comparisons:
model.matrix.gls <- function(object, ...) {
  model.matrix(terms(object), data = getData(object), ...)
}
model.frame.gls <- function(object, ...) {
model.frame(formula(object), data = getData(object), ...)
}
terms.gls <- function(object, ...) {
  terms(model.frame(object), ...)
}


#Comparisons desired:
#DB-DF
#NDB-NDF
#BB-BF
#DB-BB

Here I show how the data plot looks like, and the comparisons desired.

Thanks everyone for the help provided!

Answer 1

Please look at the documentation for emmeans::contrast; it provides for such comparisons:

library("emmeans")
emm <- emmeans(m1, "Condicion")
contrast(emm, list(
    `DB-DF` = c(-1,1,0,0,0,0),
    `NDB-NDF` = c(0,0,-1,0,0,0),
    `BB-BF` = c(0,0,0,0,-1,1),
    `DB-BB` = c(-1,0,0,0,1,0)
))

Note: This is based on the assumption that the levels of Condicion are in the order given in the question. If they are not, then they are probably in alphabetical order (the default) and you need to either permute the 1s and -1s accordingly, or create

Datos$Condicion <- factor(as.character(Datos$Condicion),
    levels = c("DB", "DF", "NDB", "NDF", "BB", "BF"))

and then re-fit the model before using the above code to obtain the comparisons.

Answer 2

Here is a small vignette. Note the lsmeans package has been deprecated in favor of the emmeans package, but they have very similar syntax for this simple example.

require(lme4)
require(emmeans)

N <- 1000
M <- 200
Condicion <- as.factor(sample(1:6, 1000, replace=T))
X <- model.matrix(~Condicion)
Bicho <- as.factor(sample(rep(1:200, each=5), 1000, replace=F))
Z <- model.matrix(~Bicho+0)
Vueltasmin <- X %*% matrix(rnorm(6, 0, 5), nrow=6, ncol=1) + Z %*% matrix(rnorm(M, 0, 10), nrow=M, ncol=1) + rnorm(1000, 0, 1)

mod <- lmer(Vueltasmin ~ Condicion + (1 | Bicho))
summary(mod)

mod.emm <- emmeans(mod, "Condicion")
mod.pair <- pairs(mod.emm)
mod.comp <- mod.pair[c(1, 4, 10, 15)]
mod.comp

They key is that after you call pairs on the emmobj, you index for the contrasts of interest (in your case, 1, 4, 10 and 15). The p-values are then readjusted for the smaller number of comparisons.