How to code contrasts in lme4 or nlme?

Question

I am looking for the most efficient way to run contrasts in R when using lme4. I have been working with a stats consultant that I really trust and she has given me the following code. I have contrasts between 6 treatments and I run these contrasts for 6 different years. So I end up writing 90 contrasts out. Now I am about to include another factor (sampling depth) in the model which will result in me writing 450 contrasts.

There must be a better way?

I have been reading ways to run contrasts in R, but haven't turned up much related to lme4. nlme would work for me too, but it is also not clear to me how it works with contrasts.

Here is my data:

https://www.dropbox.com/s/2ho6phfxhz6xlsy/Root%20biomass%2C%20whole%20core.csv

Here is the simplest form of the code, for just one year:

lm1 <- lmer(mass_sum ~ block.f+ trt + (1|block.f:trt), data = roots2)
coefs <- fixef(lm1)
varb <- vcov(lm1)

##CC vs CCW
c1 <- as.matrix(c(0,0,0,0,1,0,0,0,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccvccw <- (1-pt(abs(t1), df = 15))*2

##CC vs CS
c1 <- as.matrix(c(0,0,0,0,0,1,0,0,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccvcs <- (1-pt(abs(t1), df = 15))*2

##CC vs P
c1 <- as.matrix(c(0,0,0,0,0,0,1,0,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccvp <- (1-pt(abs(t1), df = 15))*2

##CC vs PF
c1 <- as.matrix(c(0,0,0,0,0,0,0,1,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccvpf <- (1-pt(abs(t1), df = 15))*2

##CC vs SC
c1 <- as.matrix(c(0,0,0,0,0,0,0,0,1))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccvsc <- (1-pt(abs(t1), df = 15))*2

##CCW vs CS
c1 <- as.matrix(c(0,0,0,0,1,-1,0,0,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccwvcs <- (1-pt(abs(t1), df = 15))*2

##CCW vs P
c1 <- as.matrix(c(0,0,0,0,1,0,-1,0,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccwvp <- (1-pt(abs(t1), df = 15))*2

##CCW vs PF
c1 <- as.matrix(c(0,0,0,0,1,0,0,-1,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccwvpf <- (1-pt(abs(t1), df = 15))*2

##CCW vs SC
c1 <- as.matrix(c(0,0,0,0,1,0,0,0,-1))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
ccwvsc <- (1-pt(abs(t1), df = 15))*2

##CS vs P
c1 <- as.matrix(c(0,0,0,0,0,1,-1,0,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
csvp <- (1-pt(abs(t1), df = 15))*2

##CS vs PF
c1 <- as.matrix(c(0,0,0,0,0,1,0,-1,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
csvpf <- (1-pt(abs(t1), df = 15))*2

##CS vs SC
c1 <- as.matrix(c(0,0,0,0,0,1,0,0,-1))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
csvsc <- (1-pt(abs(t1), df = 15))*2

##P vs PF
c1 <- as.matrix(c(0,0,0,0,0,0,1,-1,0))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
pvpf <- (1-pt(abs(t1), df = 15))*2

##P vs SC
c1 <- as.matrix(c(0,0,0,0,0,0,1,0,-1))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
pvsc <- (1-pt(abs(t1), df = 15))*2

##PF vs SC
c1 <- as.matrix(c(0,0,0,0,0,0,0,1,-1))
est1 <- t(c1)%*%coefs
varc1 <- t(c1)%*%varb%*%c1
t1 <- as.numeric(est1/sqrt(varc1))
pfvsc <- (1-pt(abs(t1), df = 15))*2

ccvccw
ccvcs
ccvp
ccvpf
ccvsc
ccwvcs
ccwvp
ccwvpf
ccwvsc
csvp
csvpf
csvsc
pvpf
pvsc
pfvsc

for what it's worth, `lme4` and `nlme` (and just about every other R package built on a linear modeling framework) pass the contrasts specification through to `?model.matrix`, so they all essentially work the same with respect to contrasts. — Ben Bolker, May 14 '14 at 21:07
if you're trying to compute all pairwise contrasts, you might consider the `lsmeans` package or the `multcomp` package ... — Ben Bolker, May 14 '14 at 21:08
http://mindingthebrain.blogspot.ca/2013/04/multiple-pairwise-comparisons-for.html — Ben Bolker, May 14 '14 at 21:12
FWIW: http://cran.r-project.org/web/packages/lmerTest/lmerTest.pdf; http://cran.r-project.org/web/packages/LMERConvenienceFunctions/LMERConvenienceFunctions.pdf. Both have some posthoc functions. — Henrik, May 14 '14 at 21:14
Great! I am going to try the method from the blog (by @Dan M.) first as the example is very easy to follow. — Nazer, May 14 '14 at 21:22

score 1 · Answer 1 · answered Jun 19 '15 at 13:47

For all pairwise contrasts, this site is helpful: http://www.ats.ucla.edu/stat/r/faq/testing_contrasts.htm

For example:

library(multcomp)
HUC04_model = lmer(Mean_Wr ~ HUC04 + (1|FIN), data, REML=F)
test = glht(HUC04_model,linfct=mcp(HUC04="Tukey"))
summary(test)

How to code contrasts in lme4 or nlme?

1 Answers1