##----------------------------------------------------------------------------- ## Definções do knitr. Não rodar. opts_chunk$set( cache=FALSE, tidy=FALSE, fig.width=6, fig.height=4.5, fig.align="center", dpi=100, dev="png", dev.args=list(png=list(family="Ubuntu Light", type="cairo"))) options(width=90) ##============================================================================= ## Curso de Estatística Experimental com aplicações em R ## 12 à 14 de Novembro - Manaus/AM ## Embrapa Amazônia Ocidental ## ## Prof. Walmes Zeviani ## LEG - DEST - UFPR ##============================================================================= **** # Análise de VCU de soja ##----------------------------------------------------------------------------- ## Definições da sessão. ## Lista de pacotes a serem usados na sessão. pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "lme4", "plyr", "gridExtra", "aod", "wzRfun") sapply(pkg, require, character.only=TRUE) sessionInfo() trellis.device(color=FALSE) ##----------------------------------------------------------------------------- ## Leitura dos dados. da <- read.table("http://www.leg.ufpr.br/~walmes/data/grupoexperimentos.txt", header=TRUE, sep="\t", na.string=".") str(da) ##----------------------------------------------------------------------------- ## Editar. da <- transform(da, bl=factor(bl), blin=interaction(bl, local, drop=TRUE)) str(da) sum(is.na(da$rend)) da <- subset(da, !is.na(rend), select=c("local","gen","blin","rend")) str(da) ## Nos níveis dos fatores não pode ter nome com traço, dá problema na ## cld(). Trocar por espaço. levels(da$gen) <- gsub("-", " ", levels(da$gen)) levels(da$local) <- gsub("-", " ", levels(da$local)) levels(da$gen) levels(da$local) ##----------------------------------------------------------------------------- ## Layout. x <- xtabs(~local+gen, data=da) ## dimnames(x) <- NULL mosaicplot(x, off=10, las=4) t(x) ##----------------------------------------------------------------------------- ## Ver. xyplot(rend~gen|local, data=da, as.table=TRUE) ##----------------------------------------------------------------------------- ## Ajuste do modelo que considera apenas o efeito de genótipo como ## fixo, os demais e interações são aleatórios. ## * fixo: gen; ## * aleatório: local, bl dentro de local e local interação com gen. ## da$y <- sqrt(da$rend) ## da$y <- log(da$rend) da$y <- da$rend m0 <- lmer(y~gen+(1|local)+(1|blin)+(1|local:gen), data=da, REML=FALSE) summary(m0) ## Abandona a interação. m1 <- update(m0, .~.-(1|local:gen), data=da) anova(m0, m1) ## Existe interação local:gen. ## Predição dos efeitos aleatórios. lapply(ranef(m0), c) ##----------------------------------------------------------------------------- ## Teste para os termos de efeito fixo. anova(m0) ## Não tem p-valor! É de propósito. Douglas Bates não concorda que o ## procedimento adequado para ser avaliar a estatística F seja por meio ## de ajustes no número de graus de liberdade do denominador. Para mais ## detalhes leia: ## https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html ## Alternativas (em order de recomendação): ## 1. Teste de razão de verossimilhanças entre modelos aninhados (um com ## outro sem o termo fixo de interesse, usar REML=FALSE). ## 2. Teste de Wald (inferência aproximada, pacote aod). V <- vcov(m0) b <- fixef(m0) nobars(formula(m0)) Terms <- which(attr(m0@pp$X, "assign")==1) ## Chi-quadrado de Wald (baseado na aproximação quadrática da ## verossimilhança). wt <- wald.test(Sigma=V, b=b, Terms=Terms) wt ## Chi-quadrado da razão de verossimilhanças (não baseado em aproximação ## de função, melhor opção). m1 <- update(m0, .~.-gen) anova(m0, m1) ##----------------------------------------------------------------------------- ## Diagnóstico. r <- residuals(m0, type="pearson") f <- fitted(m0) bec <- unlist(ranef(m0)[1]) beb <- unlist(ranef(m0)[2]) be <- unlist(ranef(m0)[3]) xyplot(r~f, type=c("p", "smooth")) xyplot(sqrt(abs(r))~f, type=c("p", "smooth")) grid.arrange(qqmath(r), qqmath(bec), qqmath(beb), qqmath(be), nrow=2) qqmath(~r|da$local) ##----------------------------------------------------------------------------- ## Médias ajustadas. ## Formula só da parte de efeito fixo. f <- nobars(formula(m0)) X <- LSmatrix(lm(f, da), effect=c("gen")) rownames(X) <- levels(da$gen) ## Estimativas das médias. summary(glht(m0, linfct=X), test=adjusted(type="none")) ## Contraste. Xc <- apc(X) dim(Xc) g <- summary(glht(m0, linfct=Xc), test=adjusted(type="fdr")); g ## Estimativas das médias com comparações. ## grid <- apmc(X=X, model=m0, focus="gen", test="bonferroni") grid <- apmc(X=X, model=m0, focus="gen", test="fdr") grid ##----------------------------------------------------------------------------- ## Gráfico. segplot(reorder(gen, estimate)~lwr+upr, centers=estimate, ylab="Cultivar de soja", xlab="Produtividade", data=grid, draw=FALSE, cld=grid$cld, panel=function(x, y, z, centers, subscripts, cld, ...){ panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...) panel.text(x=centers[subscripts], y=as.numeric(z)[subscripts], labels=cld[subscripts], pos=4) }) ##----------------------------------------------------------------------------- ## Para desdobrar a interação local:gen, tratar como efeito fixo. formula(m0) m0 <- lmer(y~local*gen+(1|blin), data=da, REML=FALSE) ##----------------------------------------------------------------------------- ## Fazer o desdobramento dentro do primeiro local. Os demais seguem a ## mesma regra. X <- LSmatrix(lm(nobars(formula(m0)), data=da), effect=c("local","gen")) dim(X) L <- by(X, attr(X, "grid")$local, as.matrix) L <- lapply(L, "rownames<-", levels(da$gen)) str(L) grid <- lapply(L, apmc, model=m0, focus="gen", test="fdr") str(grid) grid <- ldply(grid); names(grid)[1] <- "local" grid <- equallevels(grid, da) str(grid) ##----------------------------------------------------------------------------- ## Gráfico. segplot(gen~lwr+upr|local, centers=estimate, ylab="Genivar de arroz", xlab="Produtividade", data=grid, draw=FALSE, strip=FALSE, strip.left=TRUE, scales=list( y=list(relation="free", rot=0)), layout=c(2,NA)) ##----------------------------------------------------------------------------- ## Truque para ordenar gen dentro dos locais no gráfico. ## Função que remove o texto antes do ponto separador. yscale.components.dropend <- function(...){ ans <- yscale.components.default(...) lab <- ans$left$labels$labels ans$left$labels$labels <- gsub("^.*\\.", "", lab) ans } ## Número de gen em cada local. ngen <- apply(xtabs(~local+gen, da), 1, function(i) sum(i>0)) ngen grid$locgen <- with(grid, interaction(local, gen)) ## Ordem original dos níveis do fator locgen. on <- levels(grid$locgen) neworder <- with(grid, order(local, estimate)) orderin <- by(grid, INDICE=grid$local, function(i){ as.character(i$locgen[order(i$estimate)]) }) ## Uma cópia do fator com ondem diferente dos níveis. grid$locgen2 <- factor(grid$locgen, levels=unlist(orderin)) segplot(locgen2~lwr+upr|local, centers=estimate, ylab="Cultivar de soja", xlab="Produtividade", data=grid, draw=FALSE, scales=list(y=list(relation="free", rot=0, tck=0.5, cex=0.7)), yscale.components=yscale.components.dropend, between=list(y=0.2), layout=c(2,NA), as.table=TRUE)