
Pessoal, Estes dados apresentam um padrão peculiar onde a solução transformação Box-Cox não resolve. O Wagner mostrou um ajuste que modela a variância com uma função. Aqui estou repetindo o que ele fez, usando funções do pacote nlme, para fazer teste de razão de verossimilhança. Como a relação média variância nos dados não é exatamente uma função monótona, o ajuste com função Exp para a variância não resultou em melhoria. Então eu modelei a variância em função de grupos. Pro caso de comparar os tratamentos com a testemunha, reordene os níveis do fator para a testemunha ser o primeiro nível. Assim, o summary vai exibir estes contrastes. Fenois <- c(337.311, 344.874, 342.353, 325.546, 333.950, 330.588, 328.067, 328.067, 318.824, 331.429, 333.950, 334.790, 336.471, 338.151, 342.353, 259.160, 252.437, 268.403, 265.882, 266.723, 287.731, 88.571, 88.571, 90.252, 41.513, 52.437, 49.076, 88.571, 88.571, 90.252, 64.202, 60.000, 61.681) Cor <- factor(c(rep("ambar",6), rep("ambar_claro",3), rep("ambar",6), rep("ambar_claro",6),rep("branco",6), rep("extra_ambar_claro",3),rep("branco",3))) dados <- data.frame(Fenois, Cor) xtabs(~Cor, data = dados) library(lattice) library(latticeExtra) xyplot(Fenois ~ Cor, data = dados) # ATENÇÃO. O padrão a relação média variância não é uma função # monótona, ou seja, pelo gráfico, a variância não aumenta com a média # pois o trat 3 tem mesma média do 4 mas com variância bem maior o trat # 1 que é o de maior média tem variância baixa. mv <- aggregate(Fenois ~ Cor, data = dados, FUN = function(y) { cbind(m = mean(y), v = var(y)) }) str(mv) xyplot(Fenois[, 2] ~ Fenois[, 1], data = mv, xlab = "Média", ylab = "Variância", type = c("p", "r")) + layer(panel.text(x = x, y = y, labels = mv$Cor, pos = 3)) m0 <- lm(Fenois ~ Cor, data = dados) par(mfrow = c(2, 2)) plot(m0) layout(1) MASS::boxcox(m0) # A Box-Cox sugere transformação quando existe relação média variância # por função monótona ou problema de simetria. Esse padrão ela corrige. library(nlme) m1 <- gls(Fenois ~ Cor, data = dados, weights = varExp(), method = "ML") logLik(m0) logLik(m1) anova(m1, m0) # Razão entre modelos aponta que a função de variância exponencial (que # é monótona) não difere do modelo nulo para a variância. anova(m1) # Criar um fator que represente tratamento com a mesma variância. Nesse # caso serão formados dois grupos. dput(levels(dados$Cor)) dados$grupo <- factor(ifelse(dados$Cor %in% c("ambar_claro", "branco"), "alta", "baixa")) xyplot(Fenois ~ Cor, data = dados, groups = grupo, auto.key = TRUE) m2 <- gls(Fenois ~ Cor, data = dados, weights = varIdent(form = ~1 | grupo), method = "ML") anova(m2, m0) anova(m2) summary(m2) À disposição Walmes.