Experimento em arranjo fatorial duplo com repetições
Walmes Zeviani
Definições da sessão
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
## lattice latticeExtra doBy multcomp reshape plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUE
source("lattice_setup.R")
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()
## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.1 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-3
## [5] TH.data_1.0-3 mvtnorm_0.9-99992 doBy_4.5-10 MASS_7.3-33
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 lme4_1.1-6
## [5] Matrix_1.1-4 methods_3.1.1 minqa_1.2.3 nlme_3.1-117
## [9] Rcpp_0.11.1 RcppEigen_0.3.2.1.2 sandwich_2.3-0 stringr_0.6.2
## [13] tools_3.1.1 zoo_1.7-11
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
Importação dos dados
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/banana_culttecido.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")
## Tem perda de observações? Quantas?
cc <- !complete.cases(da)
sum(cc)
## [1] 4
da[cc,]
## cult meio rept plant altplan compraiz nfolha nraiz
## 297 t wpm 5 1 NA NA NA NA
## 298 t wpm 5 2 NA NA NA NA
## 299 t wpm 5 3 NA NA NA NA
## 300 t wpm 5 4 NA NA NA NA
## Perdeu-se uma unidade experimental inteira.
## Criando níveis para unidade experimental.
da$ue <- with(da, interaction(cult, meio, rept, drop=TRUE))
nlevels(da$ue)
## [1] 75
## Remove linhas que contenham NA.
da <- droplevels(na.omit(da))
## Mostra a estrutura do objeto.
str(da)
## 'data.frame': 296 obs. of 9 variables:
## $ cult : Factor w/ 5 levels "c","d","f","p",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ meio : Factor w/ 3 levels "bds","ms","wpm": 1 1 1 1 1 1 1 1 1 1 ...
## $ rept : int 1 1 1 1 2 2 2 2 3 3 ...
## $ plant : int 1 2 3 4 1 2 3 4 1 2 ...
## $ altplan : num 7 9 9.5 9 5.5 2 2 2 7 8.5 ...
## $ compraiz: num 11.5 20 18 14 9 0 0 0 19 14 ...
## $ nfolha : int 4 4 3 4 3 0 0 0 4 0 ...
## $ nraiz : int 5 8 11 6 6 0 0 0 7 10 ...
## $ ue : Factor w/ 74 levels "c.bds.1","d.bds.1",..: 1 1 1 1 16 16 16 16 31 31 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:4] 297 298 299 300
## .. ..- attr(*, "names")= chr [1:4] "297" "298" "299" "300"
## Os dados são de um experimento para verificar o efeito a adubação com
## potássio nos componentes de produção de soja sob diferentes níveis de
## umidade.
##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(altplan~meio|cult, data=da,
type=c("p","a"), jitter.x=TRUE)
Especificação e estimação do modelo
##-----------------------------------------------------------------------------
## Reduzindo os dados para média e número de elementos.
db <- aggregate(cbind(altplan, compraiz, nfolha, nraiz)~cult+meio+rept,
data=da, FUN=mean)
db <- merge(db, aggregate(cbind(n=altplan)~cult+meio+rept,
data=da, FUN=length))
str(db)
## 'data.frame': 74 obs. of 8 variables:
## $ cult : Factor w/ 5 levels "c","d","f","p",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ meio : Factor w/ 3 levels "bds","ms","wpm": 1 1 1 1 1 2 2 2 2 2 ...
## $ rept : int 1 2 3 4 5 1 2 3 4 5 ...
## $ altplan : num 8.62 2.88 5 6.75 8.38 ...
## $ compraiz: num 15.88 2.25 10.25 10.5 20 ...
## $ nfolha : num 3.75 0.75 1 2.75 4.75 0.75 3 3.75 2.5 3 ...
## $ nraiz : num 7.5 1.5 6.25 5 6.25 1 5.75 6 3.5 5 ...
## $ n : int 4 4 4 4 4 4 4 4 4 4 ...
## Modelo (os pesos são iguais e poderiam ser omitidos).
m0 <- lm(altplan~cult*meio, data=db, weights=n)
par(mfrow=c(2,2)); plot(m0); layout(1)
MASS::boxcox(m0)
m0 <- update(m0, log(altplan)~.)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: log(altplan)
## Df Sum Sq Mean Sq F value Pr(>F)
## cult 4 28.16 7.04 40.9 2.2e-16 ***
## meio 2 8.96 4.48 26.0 7.9e-09 ***
## cult:meio 8 2.89 0.36 2.1 0.05 *
## Residuals 59 10.16 0.17
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Comparações múltiplas
##-----------------------------------------------------------------------------
## Comparações múltiplas para o desdobramento.
LSmeans(m0, effect=c("meio","cult"))
## estimate se df t.stat p.value lwr upr meio cult
## 1 1.7710 0.0928 59 19.084 6.364e-27 1.5853 1.957 bds c
## 2 1.6699 0.0928 59 17.994 1.237e-25 1.4842 1.856 ms c
## 3 1.1219 0.0928 59 12.090 1.315e-17 0.9362 1.308 wpm c
## 4 1.0511 0.0928 59 11.326 1.972e-16 0.8654 1.237 bds d
## 5 1.0862 0.0928 59 11.705 5.104e-17 0.9005 1.272 ms d
## 6 0.9246 0.0928 59 9.963 2.953e-14 0.7389 1.110 wpm d
## 7 1.6213 0.0928 59 17.471 5.386e-25 1.4356 1.807 bds f
## 8 1.6493 0.0928 59 17.773 2.297e-25 1.4636 1.835 ms f
## 9 1.2911 0.0928 59 13.912 2.760e-20 1.1054 1.477 wpm f
## 10 2.1278 0.0928 59 22.929 4.487e-31 1.9421 2.314 bds p
## 11 1.9777 0.0928 59 21.311 2.126e-29 1.7920 2.163 ms p
## 12 1.8199 0.0928 59 19.611 1.580e-27 1.6342 2.006 wpm p
## 13 1.4318 0.0928 59 15.429 2.252e-22 1.2461 1.618 bds t
## 14 1.6472 0.0928 59 17.750 2.447e-25 1.4616 1.833 ms t
## 15 0.9715 0.1038 59 9.363 2.854e-13 0.7638 1.179 wpm t
X <- LSmatrix(m0, effect=c("meio","cult"))
grid <- equallevels(attr(X, "grid"), da)
## Divide a matriz de acordo com os níveis de A.
L <- by(X, INDICES=grid$cult, FUN=as.matrix)
L <- lapply(L, "rownames<-", levels(da$meio))
t(L[[1]])
## bds ms wpm
## (Intercept) 1 1 1
## cultd 0 0 0
## cultf 0 0 0
## cultp 0 0 0
## cultt 0 0 0
## meioms 0 1 0
## meiowpm 0 0 1
## cultd:meioms 0 0 0
## cultf:meioms 0 0 0
## cultp:meioms 0 0 0
## cultt:meioms 0 0 0
## cultd:meiowpm 0 0 0
## cultf:meiowpm 0 0 0
## cultp:meiowpm 0 0 0
## cultt:meiowpm 0 0 0
## Faz comparações múltiplas de médias entre os níveis de K dentro de
## cada nível de A.
grid <- lapply(L, apmc, model=m0, focus="meio", test="fdr")
grid <- ldply(grid); names(grid)[1] <- "cult"
grid <- equallevels(grid, da)
grid
## cult meio estimate lwr upr cld
## 1 c bds 1.7710 1.5853 1.957 a
## 2 c ms 1.6699 1.4842 1.856 a
## 3 c wpm 1.1219 0.9362 1.308 b
## 4 d bds 1.0511 0.8654 1.237 a
## 5 d ms 1.0862 0.9005 1.272 a
## 6 d wpm 0.9246 0.7389 1.110 a
## 7 f bds 1.6213 1.4356 1.807 a
## 8 f ms 1.6493 1.4636 1.835 a
## 9 f wpm 1.2911 1.1054 1.477 b
## 10 p bds 2.1278 1.9421 2.314 a
## 11 p ms 1.9777 1.7920 2.163 a
## 12 p wpm 1.8199 1.6342 2.006 a
## 13 t bds 1.4318 1.2461 1.618 a
## 14 t ms 1.6472 1.4616 1.833 a
## 15 t wpm 0.9715 0.7638 1.179 b
##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
ylab <- expression("(log) Altura de plantulas")
xlab <- expression("Meio de cultura")
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
## Gráfico final.
segplot(meio~lwr+upr|cult, centers=estimate, horizontal=FALSE,
data=grid, draw=FALSE, layout=c(NA,1),
xlab=xlab, ylab=ylab, sub=sub,
strip=strip.custom(
factor.levels=sprintf("Cultivar: %s", levels(da$cult))),
letras=sprintf("%0.2f %s", grid$estimate, grid$cld),
panel=function(x, y, z, subscripts, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers,
subscripts=subscripts, ...)
panel.text(x=as.numeric(z)[subscripts],
y=centers[subscripts],
label=letras[subscripts],
srt=90, adj=c(0.5,-0.5))
})
##-----------------------------------------------------------------------------
## A diferença em termos de distribuição para os dados originais e para
## as médias.
## Alguns zeros.
p1 <-
xyplot(compraiz~meio|cult, data=da, layout=c(NA,1),
type=c("p","a"), jitter.x=TRUE)
p2 <-
xyplot(compraiz~meio|cult, data=db, layout=c(NA,1),
type=c("p","a"), jitter.x=TRUE)
plot(p1, split=c(1,1,1,2), more=TRUE)
plot(p2, split=c(1,2,1,2), more=FALSE)
## Discreta.
p1 <-
xyplot(nfolha~meio|cult, data=da, layout=c(NA,1),
type=c("p","a"), jitter.x=TRUE)
p2 <-
xyplot(nfolha~meio|cult, data=db, layout=c(NA,1),
type=c("p","a"), jitter.x=TRUE)
plot(p1, split=c(1,1,1,2), more=TRUE)
plot(p2, split=c(1,2,1,2), more=FALSE)
## Discreta.
p1 <-
xyplot(nraiz~meio|cult, data=da, layout=c(NA,1),
type=c("p","a"), jitter.x=TRUE)
p2 <-
xyplot(nraiz~meio|cult, data=db, layout=c(NA,1),
type=c("p","a"), jitter.x=TRUE)
plot(p1, split=c(1,1,1,2), more=TRUE)
plot(p2, split=c(1,2,1,2), more=FALSE)
