I have run into a problem with the posthoc comparison for my linear mixed effects model. I'll try to explain it with a quickly constructed unperfect example:
Here my example data:
Variable<-as.factor(rep(c(1,2,3),5))
Random<-rep(c(1,2,2),5)
Result<-rnorm(15,mean=10,sd=2)
Data<-as.data.frame(cbind(Variable,Random,Result))
I actually have several fixed and random effects included in my model, but this is sufficient to illustrate my problem:
library(lme4)
LME=lmer(Result~Variable+(1|Random))
summary(LME)
Looking at the fixed effects output I only get the significances for the different levels of the variable compared to the Intercept
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 9.5104 1.3685 12.0000 6.949 1.54e-05 ***
Variable2 0.9155 1.9354 12.0000 0.473 0.645
Variable3 1.7386 1.9354 12.0000 0.898 0.387
However, I would now like to compare Variable level 1 with level 2 and Variable level 2 with level 3, so I tried the following:
library(multcomp)
summary(glht(LME, linfct=c("Variable2-Variable1=0","Variable3-Variable2=0")))
Leaving me with this error:
Error in h(simpleError(msg, call)) :
error in evaluating the argument 'object' in selecting a method for function 'summary': multcomp:::chrlinfct2matrix: variable(s) ‘Variable1’ not found
If I exclude variable level 1 and only look at the comparison of 2 to 3, the code works fine:
summary(glht(LME, linfct=c("Variable3-Variable2=0")))
Simultaneous Tests for General Linear Hypotheses
Fit: lmer(formula = Result ~ Variable + (1 | Random))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Variable3 - Variable2 == 0 0.8231 1.6694 0.493 0.622
(Adjusted p values reported -- single-step method)
I can also run the linfct function with Tukey contrasts:
summary(glht(LME, linfct= mcp(Variable="Tukey")),test=adjusted("none"))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lmer(formula = Result ~ Variable + (1 | Random))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0 0.9155 1.9354 0.473 0.636
3 - 1 == 0 1.7386 1.9354 0.898 0.369
3 - 2 == 0 0.8231 1.6694 0.493 0.622
(Adjusted p values reported -- none method)
Seeing that I am not interessted in the comparison of 3 to 1, I would then only use the other 2 p-values and adjust them in a sperate step, but this is not really the solution I am looking for. My data results in more than just the two comparisons shown here, so the option with the Tukey contrasts would leave me with a lot of comparisons I am not really interessted in.
Is there a way to get Variable1 from LME? In the fixed effects it is included as Intercept, replacing Variable1 with Intercept or any combinations I could think of did not do the trick. Or is there generally a better way to achieve the comparisons I am looking for?
Any help would be greatly appreciated!
You've really already got what you want. In this case, since Variable=1 is the reference group, it's coefficient is fixed at 0 with a variance of 0. So, the test of whether Variable1=Variable2=0 is really just a test of Variable2=0. Likewise with Variable3. You can see this from the fact that the two pieces of code below produce the same output:
summary(glht(LME, linfct=c("Variable2=0","Variable3=0", "Variable3-Variable2=0")))
# Simultaneous Tests for General Linear Hypotheses
# Fit: lmer(formula = Result ~ Variable + (1 | Random))
# Linear Hypotheses:
# Estimate Std. Error z value Pr(>|z|)
# Variable2 == 0 -0.6524 2.0145 -0.324 0.942
# Variable3 == 0 -2.0845 2.0145 -1.035 0.545
# Variable3 - Variable2 == 0 -1.4321 1.1199 -1.279 0.396
# (Adjusted p values reported -- single-step method)
summary(glht(LME, linfct=mcp(Variable="Tukey")))
# Simultaneous Tests for General Linear Hypotheses
# Multiple Comparisons of Means: Tukey Contrasts
# Fit: lmer(formula = Result ~ Variable + (1 | Random))
# Linear Hypotheses:
# Estimate Std. Error z value Pr(>|z|)
# 2 - 1 == 0 -0.6524 2.0145 -0.324 0.942
# 3 - 1 == 0 -2.0845 2.0145 -1.035 0.545
# 3 - 2 == 0 -1.4321 1.1199 -1.279 0.396
# (Adjusted p values reported -- single-step method)
So, if you only want the adjusted comparisons with the reference, you can do:
summary(glht(LME, linfct=c("Variable2=0","Variable3=0")))
# Simultaneous Tests for General Linear Hypotheses
# Fit: lmer(formula = Result ~ Variable + (1 | Random))
# Linear Hypotheses:
# Estimate Std. Error z value Pr(>|z|)
# Variable2 == 0 -0.6524 2.0145 -0.324 0.889
# Variable3 == 0 -2.0845 2.0145 -1.035 0.404
# (Adjusted p values reported -- single-step method)