##============================================================================= ## 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/script17") options(width=90) ##----------------------------------------------------------------------------- ## Definições da sessão. ## Lista de pacotes a serem usados na sessão. pkg <- c("lattice", "latticeExtra", "gridExtra", "doBy", "multcomp", "reshape", "plyr", "nlme", "wzRfun") sapply(pkg, require, character.only=TRUE) ## source("http://dl.dropboxusercontent.com/u/48140237/apc.R") ## source("http://dl.dropboxusercontent.com/u/48140237/segplotby.R") ## source("http://dl.dropboxusercontent.com/u/48140237/equallevels.R") sessionInfo() trellis.device(color=FALSE) ##----------------------------------------------------------------------------- ## Leitura dos dados. da <- read.table("http://www.leg.ufpr.br/~walmes/data/soja_solo_tcc.txt", header=TRUE, sep="\t") str(da) ##----------------------------------------------------------------------------- ## Visualizar. xyplot(ctc~adubpk|sistema, groups=prof, data=da, type=c("p","smooth"), auto.key=TRUE) ##----------------------------------------------------------------------------- ## Análise considerando o modelo para experimentos em parcela ## subdividida sendo o sistema a parcela e a adubação a subparcela. ## Efeito categórico para adub para começar. da$adub <- factor(da$adubpk) ## Níveis das parcelas e subparcelas. da$parc <- with(da, interaction(bloco, sistema, drop=TRUE)) da$subp <- with(da, interaction(bloco, sistema, adub, drop=TRUE)) str(da) ## Adub categórico. m0 <- lme(ctc~bloco+sistema*adub*prof, random=~1|parc/subp, data=da, method="ML") ##----------------------------------------------------------------------------- ## Diagnóstico. r <- residuals(m0, type="pearson") f <- fitted(m0) bp <- unlist(ranef(m0)[1]) bs <- unlist(ranef(m0)[2]) grid.arrange(nrow=3, xyplot(r~f, type=c("p", "smooth")), xyplot(sqrt(abs(r))~f, type=c("p", "smooth")), qqmath(r), qqmath(bp), qqmath(bs)) ##----------------------------------------------------------------------------- ## Quadro para o teste dos efeitos fixos (Wald). anova(m0) ## Abandono das interações não significativas. m1 <- update(m0, .~bloco+sistema*prof+adubpk) ## m2 <- update(m0, .~bloco+sistema*prof) anova(m1, m0) anova(m1) ## Sem efeito da adubação. O que existe é uma diferença entre camadas ## para a que depende do sistema. ##----------------------------------------------------------------------------- ## Obter as médias de ctc em cada prof para cada sistema. Pode-se ficar ## a adubação em qualquer nível, bem como bloco. X <- LSmatrix(lm(formula(m1), data=da), effect=c("sistema","prof"), at=list(adubpk=120)) head(X) grid <- attr(X, "grid") rownames(X) <- apply(grid, 1, paste, collapse="|") L <- by(X, grid$sistema, apc, lev=c(20,40)) L <- lapply(L, as.matrix); L str(L) Xc <- do.call(rbind, L) rownames(Xc) <- paste0(names(L), "/", rownames(Xc)) Xc ## Teste para os contrastes. summary(glht(m1, linfct=Xc)) ## Médias com IC. ci <- confint(glht(m1, linfct=X), calpha=univariate_calpha()) grid <- cbind(grid, ci$confint) grid <- equallevels(grid, da) grid ##----------------------------------------------------------------------------- ## Representação. ## names(trellis.par.get()) l <- levels(grid$prof) key <- list(corner=c(x=0.9, y=0.05), type="o", divide=1, title="Camada do solo (cm)", cex.title=1.1, lines=list(pch=1:length(l), col=trellis.par.get("plot.symbol")$col), text=list(l)) segplot(sistema~lwr+upr, centers=Estimate, groups=grid$prof, f=0.05, data=grid, draw=FALSE, panel=panel.segplotBy, key=key, scales=list(y=list(rot=90)), pch=as.integer(grid$prof), xlab="CTC do solo", ylab="Sistema de cultivo") ##----------------------------------------------------------------------------- ## Os erros padrões de comparações de prof dentro de sistema e de ## sistema dentro de prof são diferentes pois são fatores com níveis em ## estratos diferentes no modelos misto de parcela subsubdividida. ## Desdobrar prof dentro (fixando) de sistema. L <- by(X, grid$sistema, apc, lev=c(20,40)) L <- lapply(L, as.matrix) Xc <- do.call(rbind, L) rownames(Xc) <- paste0(names(L), "/", rownames(Xc)) summary(glht(m1, linfct=Xc)) ## Desdobrar sistema dentro (fixando) de prof. L <- by(X, grid$prof, apc, lev=levels(grid$sistema)) L <- lapply(L, as.matrix) Xc <- do.call(rbind, L) rownames(Xc) <- paste0(names(L), "/", rownames(Xc)) summary(glht(m1, linfct=Xc)) summary(glht(m1, linfct=X))