2

lme4 を使用するときに R でコントラストを実行する最も効率的な方法を探しています。私は本当に信頼している統計コンサルタントと協力しており、彼女は私に次のコードを提供してくれました。私は6回の治療の間にコントラストを持っており、これらのコントラストを6年間実行しています. だから私は90のコントラストを書き出すことになります. ここで、別の要素 (サンプリング深度) をモデルに含めようとしています。これにより、450 のコントラストが作成されます。

もっと良い方法があるはずですか？

私はRでコントラストを実行する方法を読んできましたが、lme4. nlme私にとってもうまくいくでしょうが、それがコントラストでどのように機能するかは私には明らかではありません.

ここに私のデータがあります：

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

以下は、1 年間だけの最も単純な形式のコードです。

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
4

1 に答える 1

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)
于 2015-06-19T13:47:22.060 に答える