multcomp Tukey-Kramer

Question

Ho un esperimento sbilanciato dove in tre siti (L, M, H) misuriamo un parametro (met) in quattro diversi tipi di vegetazione (a, b, c, d). Tutti i tipi di vegetazione sono presenti in tutti e tre i siti. I tipi di vegetazione vengono replicati 4 volte a L e M e 8 volte a H.

Pertanto un anova e TukeyHSD semplici non funzioneranno. I pacchetti Agricolae (HSD.test) e DTK (DTK.test) funzionano solo per i progetti one-way, e quindi c'è multcomp ... Il test Tukey nella funzione mcp calcola i contrasti Tukey-Kramer o fornisce i normali contrasti Tukey? Presumo che il primo sia il caso perché il pacchetto è orientato a testare confronti multipli per progetti sbilanciati, ma non sono sicuro perché i valori di p prodotti con entrambi gli approcci sono praticamente gli stessi. Quale test sarebbe quindi appropriato?

Inoltre, esistono approcci più idonei per eseguire un'anova bidirezionale per insiemi di dati non bilanciati?

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

Per i dati sbilanciati, anova di tipo III SS possono essere usati al posto di tipo I SS [1]. Calcolo di tipo III ANOVA a [2] R:

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

Per i dati sbilanciati, può essere utilizzato confronti a coppie di mezzi adeguati. Calcolo in [4] R:

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")) 

linee 2 e 3 confrontare i livelli di effetti principali "del sito" e "vytpe". Le righe 4 e 5 confrontano i livelli di un fattore per ciascun livello di un altro fattore separatamente.

Spero che questo aiuti.

Riferimenti

[1] Miliken e Johnsen. 2009. Analisi di dati disordinati. Volume 1.

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

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