#========================================================================================== # ESTATÍSTICA EXPERIMENTAL COM O R # Prof. Walmes Zeviani # LEG - DEST - UFPR # # Análise de experimentos simples com afastamento de pressupostos. #========================================================================================== #------------------------------------------------------------------------------------------ # definições da sessão require(lattice) require(agricolae) require(contrast) require(multcomp) require(doBy) require(latticeExtra) require(ExpDes) require(nlme) require(plyr) #------------------------------------------------------------------------------------------ # lendo arquivos de dados da <- read.table("http://www.leg.ufpr.br/~walmes/data/volume.txt", header=TRUE, sep="\t", colClasses=c("factor","factor","NULL","numeric")) str(da) da$dose <- factor(da$dose, levels=c("0","5","25")) #------------------------------------------------------------------------------------------ # ver xyplot(volu~dose|gen, data=da, type=c("p","a")) #------------------------------------------------------------------------------------------ # análise m0 <- lm(volu~gen*dose, data=da) par(mfrow=c(2,2)); plot(m0); layout(1) # fuga! #------------------------------------------------------------------------------------------ # transformação MASS::boxcox(m0) # raíz cúbica m1 <- lm(volu^(1/3)~gen*dose, data=da) par(mfrow=c(2,2)); plot(m1); layout(1) # melhorou #------------------------------------------------------------------------------------------ # modelar a variância apropos("^var") m2 <- gls(volu~gen*dose, data=da, weights=varExp(), method="ML") logLik(m2) m2 <- gls(volu~gen*dose, data=da, weights=varPower(), method="ML") logLik(m2) logLik(m0) anova(m2, m0) plot(sqrt(abs(residuals(m2, type="pearson")))~fitted(m2)) # ok qqnorm(residuals(m2, type="pearson")) # ok anova(m2) summary(m2) #------------------------------------------------------------------------------------------ # não tem interação então? m3 <- gls(volu~gen+dose, data=da, weights=varPower(), method="ML") anova(m3, m2) #------------------------------------------------------------------------------------------ # desdobrando a interação X <- popMatrix(m0, effect=c("dose","gen")) str(X) #------------------------------------------------------------------------------------------ # médias ajustadas pm <- summary(glht(m2, linfct=X), test=adjusted(type="fdr")) pm pm <- confint(pm, calpha=univariate_calpha()) pm <- cbind(attr(X, "grid"), pm$confint) str(pm) #------------------------------------------------------------------------------------------ # gráfico pm$dose <- factor(pm$dose, levels=levels(da$dose)) segplot(dose~lwr+upr|gen, centers=Estimate, data=pm, ylab="Dose", xlab="Volume", draw.bands=FALSE, segments.fun=panel.arrows, ends="both", angle=90, length=1, unit="mm") #------------------------------------------------------------------------------------------ # contraste entre médias ajustadas gen <- split(1:nrow(X), f=attr(X, "grid")$gen) source("funcoes.R") pmc <- lapply(gen, function(i){ XX <- X[i,] XXX <- apc(XX, lev=levels(da$dose)) c0 <- summary(glht(m2, linfct=XXX), test=adjusted(type="fdr")) c0$focus <- "dose" c0 }) lapply(pmc, print) let <- lapply(pmc, function(l) cld(l)$mcletters$Letters) let #------------------------------------------------------------------------------------------ # colocando no gráfico str(pm) pm$cld <- do.call(c, let) segplot(dose~lwr+upr|gen, centers=Estimate, data=pm, ylab="Dose", xlab="Volume", draw.bands=FALSE, segments.fun=panel.arrows, ends="both", angle=90, length=1, unit="mm", panel=function(x, y, z, centers, subscripts, ...){ panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...) panel.text(centers[subscripts], z[subscripts], labels=paste(format(centers, digits=2), pm$cld[subscripts]), pos=3, cex=0.8) }) #------------------------------------------------------------------------------------------