如何在 lme4 或 nlme 中编码对比？

Question

我正在寻找使用 lme4 时在 R 中运行对比的最有效方法。我一直在与一位我非常信任的统计顾问合作，她给了我以下代码。我有 6 种治疗方法之间的对比，我运行这些对比 6 不同的年份。所以我最终写出了90个对比。现在我将在模型中加入另一个因素（采样深度），这将导致我写出 450 个对比。

一定有更好的办法吗？

我一直在阅读在 R 中运行对比的方法，但没有发现与 lme4 有太大关系。 nlme 也适用于我，但我也不清楚它是如何与对比一起工作的。

这是我的数据：

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

这是最简单的代码形式，仅一年：

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

Answer 1

对于所有成对的对比，这个网站很有帮助： http://www.ats.ucla.edu/stat/r/faq/testing_contrasts.htm

例如：

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