multcomp Tukey-Kramer

Question

Tôi có một thử nghiệm không cân bằng ở ba vị trí (L, M, H), chúng tôi đo một tham số (met) trong bốn loại thảm thực vật khác nhau (a, b, c, d). Tất cả các loại thực vật đều có mặt ở cả ba địa điểm. Các kiểu thảm thực vật được nhân đôi 4 lần tại L và M và 8 lần tại H.

Do đó, một anova đơn giản và TukeyHSD sẽ không hoạt động. Các gói Agricolae (HSD.test) và DTK (DTK.test) chỉ hoạt động cho các thiết kế một chiều, và sau đó có nhiều ... Kiểm tra Tukey trong hàm mcp có tính tương phản Tukey-Kramer hay không, tương đương với Tukey thông thường? Tôi cho rằng trường hợp đầu tiên là trường hợp vì gói được hướng tới thử nghiệm nhiều so sánh cho các thiết kế không cân bằng, nhưng tôi không chắc chắn vì các giá trị p được tạo ra với cả hai cách tiếp cận hầu như giống nhau. Thử nghiệm nào sau đó sẽ phù hợp?

Ngoài ra, có cách tiếp cận phù hợp hơn để thực hiện một cách anova hai chiều cho các tập dữ liệu không cân bằng không?

library(multcomp) 

(met  <- c(rnorm(16,6,2),rnorm(16,5,2),rnorm(32,4,2))) 
(site <- c(rep("L", 16), rep("M", 16), rep("H", 32))) 
(vtype <- c(rep(letters[1:4], 16), rep(letters[1:4], 16), rep(letters[1:4], 32))) 

dat <- data.frame(site, vtype, met) 

# using aov and TukeyHSD 
aov.000 <- aov(met ~ site * vtype, data=dat) 
summary(aov.000) 
TukeyHSD(aov.000) 

# using Anova, and multcomp 
lm.000  <- lm(met ~ site * vtype, data=dat) 
summary(lm.000) 
library(car) 
Anova.000 <- Anova(lm.000, data=dat) 

dat$int <- with(dat, interaction(site, vtype, sep = "x")) 
lm.000 <- lm(met ~ int, data = dat) 
summary(lm.000) 
summary(glht.000 <- glht(lm.000, linfct = mcp(int = "Tukey"))) 

Answer 1

Đối với dữ liệu không cân bằng, anova có loại III SS có thể được sử dụng thay cho loại I SS [1]. Tính toán loại III anova trong R [2]:

model <- (met ~ site * vtype) 
defopt <- options() 
options(contrasts=c("contr.sum", "contr.poly")) 
print(drop1(aov(model),~.,test="F")) 
options <- defopt 

Đối với dữ liệu không cân bằng, có thể sử dụng so sánh hai phương pháp điều chỉnh. Tính toán trong R [4]: ​​

library(lsmeans) 
print(lsmeans(model, list(pairwise ~ site)), adjust = c("tukey")) 
print(lsmeans(model, list(pairwise ~ vtype)), adjust = c("tukey")) 
print(lsmeans(model, list(pairwise ~ site | vtype)), adjust = c("tukey")) 
print(lsmeans(model, list(pairwise ~ vtype | site)), adjust = c("tukey")) 

Dòng 2 và 3 so sánh các mức hiệu ứng chính "trang web" và "vytpe". Các dòng 4 và 5 so sánh các mức của một yếu tố ở mỗi cấp của một yếu tố khác một cách riêng biệt.

Tôi hy vọng điều này sẽ hữu ích.

Tham chiếu

[1] Miliken và Johnsen. 2009. Phân tích dữ liệu lộn xộn. Khối lượng 1.

[2] http://www.statmethods.net/stats/anova.html

[3] http://cran.r-project.org/web/packages/lsmeans/vignettes/using-lsmeans.pdf