##============================================================================= ## Ajuste de modelos lineares e mistos no ambiente R ## 08 e 09 de Novembro de 2014 - Esalq/USP ## Pós Graduação de Genética e Melhoramento de Plantas ## ## Prof. Walmes Zeviani ## LEG - DEST - UFPR ##============================================================================= ##----------------------------------------------------------------------------- ## Definções do knitr. Não rodar. opts_chunk$set( cache=FALSE, tidy=FALSE, fig.width=7, fig.height=6, fig.align="center", dpi=150, dev="png", dev.args=list(png=list(family="Ubuntu Light", type="cairo"))) opts_chunk$set(fig.path="fig/script20") options(width=90) ##----------------------------------------------------------------------------- ## 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) #source("http://dl.dropboxusercontent.com/u/48140237/apc.R") #source("http://dl.dropboxusercontent.com/u/48140237/equallevels.R") #source("http://dl.dropboxusercontent.com/u/48140237/desdobrglht.R") sessionInfo() trellis.device(color=FALSE) ##----------------------------------------------------------------------------- ## Leitura dos dados. da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf2/VCU0910.csv", header=TRUE, sep=";") str(da) ##----------------------------------------------------------------------------- ## Editar. da <- transform(da, bl=factor(bl), blin=interaction(bl, ensaio, drop=TRUE)) str(da) sum(is.na(da$prod)) da <- na.omit(da) 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$cult) <- gsub("-", " ", levels(da$cult)) levels(da$cult) levels(da$ensaio) ##----------------------------------------------------------------------------- ## Layout. x <- xtabs(~ensaio+cult, data=da) ## dimnames(x) <- NULL mosaicplot(x, off=10, las=4) x ##----------------------------------------------------------------------------- ## Ver. xyplot(prod~cult|ensaio, data=da, as.table=TRUE) ## xyplot(sqrt(prod)~cult|ensaio, data=da, as.table=TRUE) ##----------------------------------------------------------------------------- ## Ajuste do modelo que considera apenas o efeito de cultivar como ## fixo, os demais e interações são aleatórios. ## * fixo: cult; ## * aleatório: ensaio, bl dentro de ensaio e ensaio interação com cult. ## * Será análisado a raíz da produtividade porque nessa escala os ## pressupostos são melhor atendidos. da$y <- sqrt(da$prod) m0 <- lmer(y~cult+(1|ensaio)+(1|ensaio:bl)+(1|ensaio:cult), data=da, REML=FALSE) summary(m0) ## Abandona a interação. m1 <- update(m0, .~.-(1|ensaio:cult), data=da) anova(m0, m1) ## Existe interação ensaio:cult. ## 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, .~.-cult) 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$ensaio) ##----------------------------------------------------------------------------- ## Médias ajustadas. ## Formula só da parte de efeito fixo. f <- nobars(formula(m0)) X <- LSmatrix(lm(f, da), effect=c("cult")) rownames(X) <- levels(da$cult) ## 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="cult", test="bonferroni") grid <- apmc(X=X, model=m0, focus="cult", test="fdr") grid ##----------------------------------------------------------------------------- ## Gráfico. segplot(reorder(cult, estimate)~lwr+upr, centers=estimate, ylab="Cultivar de arroz", xlab="raíz da 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 ensaio:cult, tratar como efeito fixo. formula(m0) m0 <- lmer(y~ensaio*cult+(1|ensaio:bl), data=da, REML=FALSE) ## Aviso devido cultivares que não ocorrem em alguns ensaios. Com isso a ## estrutura não é completamente cruzada. length(fixef(m0)) ## Combinações previstas se fosse completamente cruzado. with(da, nlevels(interaction(ensaio, cult, drop=TRUE))) ## Combinações que de fato existem. with(da, nlevels(interaction(ensaio, cult, drop=FALSE))) ## Cria variável combinando níveis de ensaio com cult. da$enscult <- with(da, interaction(ensaio, cult, drop=TRUE)) ## O mesmo modelo escrito com uma variável combinada. m0 <- lmer(y~enscult+(1|ensaio:bl), data=da, REML=FALSE) head(fixef(m0)) ##----------------------------------------------------------------------------- ## Fazer o desdobramento dentro do primeiro ensaio. Os demais seguem a ## mesma regra. X <- LSmatrix(lm(nobars(formula(m0)), data=da), effect="enscult") rownames(X) <- levels(da$enscult) dim(X) ## O nome do primeiro ensaio. split <- levels(da$ensaio)[1] split ## Seleciona as linhas correspondentes à esse ensaio. Xs <- X[grep(split, rownames(X)),] dim(Xs) summary(glht(m0, linfct=Xs), test=adjusted(type="fdr")) confint(glht(m0, linfct=Xs), calpha=univariate_calpha()) ##----------------------------------------------------------------------------- ## Comparações multiplas dentro do primeiro ensaio. ## Ordenar as linhas da matriz Xs pela estimativa das médias. a <- confint(glht(m0, linfct=Xs), calpha=univariate_calpha()) Xs <- Xs[order(a$confint[,"Estimate"], decreasing=FALSE),] Xc <- apc(Xs) a <- confint(glht(m0, linfct=Xs), calpha=univariate_calpha()); a Xc <- apc(Xs) s <- summary(glht(m0, linfct=Xc), test=adjusted(type="fdr")) ret <- cld2(s) ci <- cbind(data.frame(a$confint), cld=ret$mcletters$Letters) arrange(ci, -Estimate) ##----------------------------------------------------------------------------- ## Obter todas as médias das cult em cada ensaio. g <- confint(glht(m0, linfct=X), calpha=univariate_calpha()) grid <- attr(X, "grid") ## Quebrando o nome no ponto. quebra <- do.call(rbind, strsplit(grid$enscult, split="\\.")) colnames(quebra) <- c("ensaio","cult") grid <- cbind(quebra, grid, g$confint) grid <- equallevels(grid, da) str(grid) ##----------------------------------------------------------------------------- ## Gráfico. segplot(cult~lwr+upr|ensaio, centers=Estimate, ylab="Cultivar 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 cultivares dentro dos ensaios 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 cult em cada ensaio. ncult <- apply(xtabs(~ensaio+cult, da), 1, function(i) sum(i>0)) ncult ## Ordem original dos níveis do fator enscult. on <- levels(grid$enscult) neworder <- with(grid, order(ensaio, Estimate)) orderin <- by(grid, INDICE=grid$ensaio, function(i){ as.character(i$enscult[order(i$Estimate)]) }) ## Uma cópia do fator com ondem diferente dos níveis. grid$enscult2 <- factor(grid$enscult, levels=unlist(orderin)) segplot(enscult2~lwr+upr|ensaio, centers=Estimate, ylab="Cultivar de arroz", xlab="raíz da 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)