如何在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)