Euriann, Considere o código abaixo. library(doBy) # Conjunto de ensaios em DBC em 6 locais (gerado por simulação). da <- expand.grid(loc = gl(6, 1), rpt = gl(3, 1), blc = gl(4, 1), trt = gl(5, 1)) da$prod <- rnorm(nrow(da)) # Análise conjunta. m0 <- aov(terms(prod ~ loc/rpt/blc + trt, keep.order = TRUE), data = da) anova(m0) # Médias amostrais. aggregate(prod ~ trt, data = da, FUN = mean) # Médias ajustadas. LSmeans(m0, effect = "trt") # Análises por local. fits <- lapply(split(da, f = da$loc), FUN = aov, formula = terms(prod ~ rpt/blc + trt, keep.order = TRUE)) lapply(fits, anova) # Médias amostrais por local. aggregate(prod ~ trt + loc, data = da, FUN = mean) # Médias ajustadas. lapply(fits, FUN = LSmeans, effect = "trt") # Para o delineamento usado nesse exemplo, as médias ajustadas coincidem # com as médias amostrais. À disposição. Walmes.