Exemplo de análise
Experimento em delineamento de blocos casualizados
Definições da sessão
##-----------------------------------------------------------------------------
## Pacotes necessários.
require(latticeExtra)
require(plyr)
require(doBy)
require(tables)
require(multcomp)
## Remove objetos da mémoria para começar vazia.
rm(list=ls())
##-----------------------------------------------------------------------------
## Definições gráficas.
## Tábua de cores.
mycol <- brewer.pal(6, "Set1")
## Definições de preenchimentos, linhas, etc.
ps <- list(box.rectangle=list(col=1, fill=c("gray70")),
box.umbrella=list(col=1, lty=1),
dot.symbol=list(col=1, pch=19),
dot.line=list(col="gray50", lty=3),
plot.symbol=list(col=1, cex=1.2),
plot.line=list(col=1),
plot.polygon=list(col="gray95"),
superpose.line=list(col=mycol),
superpose.symbol=list(col=mycol, pch=c(1:5)),
superpose.polygon=list(col=mycol),
strip.background=list(col=c("gray80","gray50"))
)
trellis.par.set(ps)
## show.settings()Carregar os dados
Download dos dados: http://www.leg.ufpr.br/~walmes/data/DBCprod.txt
Download do script: DBCprod.R
##-----------------------------------------------------------------------------
## Se o arquivo não existir, descomente aqui para fazer o download.
## url <- "http://www.leg.ufpr.br/~walmes/data/DBCprod.txt"
## download.file(url=url, destfile=basename(url))
## Diretório de trabalho.
## getwd()
## Arquivos presentes no diretório.
list.files(pattern="*.txt")## [1] "DBCprod.txt"## Encoding do arquivo. Função do sistema LINUX.
system("file DBCprod.txt")
## Lendo o arquivo.
da <- read.table(file="DBCprod.txt", header=TRUE, sep="\t",
colClasses=c("factor","factor","integer"),
encoding="utf-8")
str(da)## 'data.frame': 80 obs. of 3 variables:
## $ bloc: Factor w/ 5 levels "3","5","6","7",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ trat: Factor w/ 16 levels "ADR1","ADR2",..: 7 7 7 7 7 9 9 9 9 9 ...
## $ prod: int 2819 2498 2812 3299 2474 2395 2479 2466 2535 2912 ...## É possível ler direto pelo link sem download.
## da <- read.table(file=url, header=TRUE, sep="\t",
## colClasses=c("factor","factor","integer"),
## encoding="utf-8")
## str(da)
## Primeiras linhas.
head(da)## bloc trat prod
## 1 3 Controle 2819
## 2 5 Controle 2498
## 3 6 Controle 2812
## 4 7 Controle 3299
## 5 8 Controle 2474
## 6 3 Kimberlit 2395Análise exploratória preliminar
##-----------------------------------------------------------------------------
## Níveis de tratamento.
levels(da$trat)## [1] "ADR1" "ADR2" "ADR3"
## [4] "ADR4" "ADR5" "ADR6"
## [7] "Controle" "Helena" "Kimberlit"
## [10] "Rigrantec" "SJ1-BJ-SATS" "SJ2-OMEX"
## [13] "SJ3-Quimifol Fenix" "SJ4-Blend" "Tradecorp"
## [16] "UPL"## Para controle ficar como primeiro nível.
da$trat <- relevel(da$trat, ref="Controle")
## Tabela resumo com médias, frequências, variâncias, etc.
tabular(trat+1~prod*(length+mean+median+sd+min+max), data=da)##
## prod
## trat length mean median sd min max
## Controle 5 2780 2812 333.6 2474 3299
## ADR1 5 2634 2589 344.2 2216 3039
## ADR2 5 2621 2498 448.5 2066 3275
## ADR3 5 2792 2902 303.8 2264 3013
## ADR4 5 2668 2688 432.7 2083 3291
## ADR5 5 2751 2719 291.0 2444 3157
## ADR6 5 2584 2533 508.7 2094 3271
## Helena 5 2817 2772 304.9 2404 3235
## Kimberlit 5 2557 2479 204.4 2395 2912
## Rigrantec 5 2887 2914 520.9 2304 3645
## SJ1-BJ-SATS 5 2811 2767 137.5 2728 3054
## SJ2-OMEX 5 2810 2978 385.1 2229 3134
## SJ3-Quimifol Fenix 5 2540 2530 328.7 2054 2906
## SJ4-Blend 5 2879 2989 553.2 1987 3489
## Tradecorp 5 2889 3111 476.6 2154 3347
## UPL 5 2758 2816 224.0 2377 2951
## All 80 2736 2768 361.7 1987 3645## Gráfico.
dotplot(prod~trat, data=da, groups=bloc,
scales=list(x=list(rot=90)),
xlab="Tratamentos",
ylab=expression("Produtividade"~(kg~ha^{-1})))
Análise dos dados
##-----------------------------------------------------------------------------
## Especificação do modelo.
## Dadas as condições em que o experimento foi realizado, ou seja, sob
## um delineamento de blocos casualizados, o modelo estatístico
## correspondente é
##
## y_{ij} = \mu+\alpha_i+\tau_j+\epsilon_{ij}
##
## em que y são os valores observados, \alpha são os efeitos de blocos,
## \tau são os efeitos de tratamentos, \epsilon é um termo aleatório com
## média 0 e variância constante, i e j representam níveis de bloco e
## tratamento.
## Ajuste do modelo.
m0 <- lm(prod~bloc+trat, data=da)
## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Não se tem fuga dos pressupostos.
## Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 2991980 747995 7.1436 8.99e-05 ***
## trat 15 1059542 70636 0.6746 0.7985
## Residuals 60 6282464 104708
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Sem feito dos níveis de tratamentos. Não se rejeita a hipótese nula
## de que as médias são iguais.
##-----------------------------------------------------------------------------
## Teste da aditividade dos blocos.
## Ver:
## http://www.unh.edu/halelab/BIOL933/labs/lab5.pdf
## http://www.stat.purdue.edu/~bacraig/notes1/topic11.pdf
## Valores ajustados ao quadrado.
haty <- fitted(m0)^2
xyplot(residuals(m0)~haty, type=c("p", "smooth"))
## Teste da aditividade de Tukey com 1 parâmetro.
anova(update(m0, .~.+I(haty^2)))## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 2991980 747995 7.2327 8.272e-05 ***
## trat 15 1059542 70636 0.6830 0.7904
## I(haty^2) 1 180790 180790 1.7481 0.1912
## Residuals 59 6101673 103418
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Obtenção das médias dos tratamentos retornadas pelo modelo
## (correspondem as médias amostrais por o experimento ser regular).
## Valores das médias.
LSmeans(m0, effect="trat")## estimate se df t.stat p.value trat
## 1 2780.4 144.7119 60 19.21334 2.563387e-27 Controle
## 2 2633.8 144.7119 60 18.20030 4.138344e-26 ADR1
## 3 2620.8 144.7119 60 18.11046 5.323829e-26 ADR2
## 4 2791.6 144.7119 60 19.29074 2.081712e-27 ADR3
## 5 2667.6 144.7119 60 18.43386 2.158499e-26 ADR4
## 6 2751.0 144.7119 60 19.01018 4.439757e-27 ADR5
## 7 2583.8 144.7119 60 17.85478 1.095632e-25 ADR6
## 8 2816.6 144.7119 60 19.46349 1.310985e-27 Helena
## 9 2557.4 144.7119 60 17.67235 1.841649e-25 Kimberlit
## 10 2887.4 144.7119 60 19.95274 3.596487e-28 Rigrantec
## 11 2811.0 144.7119 60 19.42480 1.453679e-27 SJ1-BJ-SATS
## 12 2810.0 144.7119 60 19.41789 1.480772e-27 SJ2-OMEX
## 13 2540.2 144.7119 60 17.55349 2.588270e-25 SJ3-Quimifol Fenix
## 14 2878.6 144.7119 60 19.89193 4.218316e-28 SJ4-Blend
## 15 2889.2 144.7119 60 19.96518 3.481215e-28 Tradecorp
## 16 2758.4 144.7119 60 19.06132 3.864955e-27 UPLRepresentar os resultados
##-----------------------------------------------------------------------------
## Para representar em gráfico com IC.
## Matriz de pesos para multiplicar os coeficientes e ter as médias
## ajustadas.
L <- LSmatrix(m0, effect="trat")
rownames(L) <- levels(da$trat)
## names(attributes(L))
## str(L)
## Tabela com um tratamento por linha.
grid <- attr(L, "grid")
## Estimativas das funções lineares.
g0 <- glht(m0, linfct=L)
## Cobertura individual de 95%.
ic <- confint(g0, calpha=univariate_calpha())$confint
grid <- cbind(grid, ic)
grid$trat <- factor(grid$trat, levels=levels(da$trat))
rownames(grid) <- NULL
grid## trat Estimate lwr upr
## 1 Controle 2780.4 2490.933 3069.867
## 2 ADR1 2633.8 2344.333 2923.267
## 3 ADR2 2620.8 2331.333 2910.267
## 4 ADR3 2791.6 2502.133 3081.067
## 5 ADR4 2667.6 2378.133 2957.067
## 6 ADR5 2751.0 2461.533 3040.467
## 7 ADR6 2583.8 2294.333 2873.267
## 8 Helena 2816.6 2527.133 3106.067
## 9 Kimberlit 2557.4 2267.933 2846.867
## 10 Rigrantec 2887.4 2597.933 3176.867
## 11 SJ1-BJ-SATS 2811.0 2521.533 3100.467
## 12 SJ2-OMEX 2810.0 2520.533 3099.467
## 13 SJ3-Quimifol Fenix 2540.2 2250.733 2829.667
## 14 SJ4-Blend 2878.6 2589.133 3168.067
## 15 Tradecorp 2889.2 2599.733 3178.667
## 16 UPL 2758.4 2468.933 3047.867## Gráfico de segmentos.
segplot(trat~lwr+upr, data=grid, centers=Estimate, draw=FALSE,
xlab="Tratamentos",
ylab=expression("Produtividade"~(kg~ha^{-1})))
## As médias não são diferentes. Usar a mesma letra para todas. Colocar
## as médias com letras sobre os pontos médios no gráfico. Arredondar as
## médias para uma casa decimal.
grid$txt <- sprintf("%.1f a", grid$Estimate)
## Reordenar os níveis de trat da menor para maior média.
grid$trat <- with(grid, reorder(x=trat, X=Estimate))
levels(grid$trat)## [1] "SJ3-Quimifol Fenix" "Kimberlit" "ADR6"
## [4] "ADR2" "ADR1" "ADR4"
## [7] "ADR5" "UPL" "Controle"
## [10] "ADR3" "SJ2-OMEX" "SJ1-BJ-SATS"
## [13] "Helena" "SJ4-Blend" "Rigrantec"
## [16] "Tradecorp"## Ampliar um pouco a amplitude vertical do gráfico para caber o texto.
nl <- levels(grid$trat)
ylim <- c(0.5, length(nl)+1)
scales <- list(y=list(at=seq_along(nl), labels=nl))
segplot(trat~lwr+upr, data=grid, centers=Estimate, draw=FALSE,
ylab="Tratamentos",
xlab=expression("Produtividade"~(kg~ha^{-1})),
txt=grid$txt,
ylim=ylim, scales=scales,
panel=function(z, centers, txt, ...){
panel.segplot(z=z, centers=centers, ...)
panel.text(x=centers, y=as.integer(z),
label=txt, pos=3, cex=0.8)
})
##-----------------------------------------------------------------------------
## Versões dos pacotes e data do documento.
sessionInfo()## R version 3.2.0 (2015-04-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid methods stats graphics grDevices utils
## [7] datasets base
##
## other attached packages:
## [1] multcomp_1.3-8 TH.data_1.0-4 mvtnorm_1.0-1
## [4] tables_0.7.79 Hmisc_3.14-5 Formula_1.1-2
## [7] doBy_4.5-12 survival_2.38-1 plyr_1.8.1
## [10] latticeExtra_0.6-26 lattice_0.20-31 RColorBrewer_1.0-5
## [13] rmarkdown_0.3.10 knitr_1.8
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.11.3 cluster_2.0.1 splines_3.2.0 MASS_7.3-39
## [5] stringr_0.6.2 tools_3.2.0 nnet_7.3-9 htmltools_0.2.6
## [9] yaml_2.1.13 digest_0.6.4 Matrix_1.2-0 formatR_1.0
## [13] acepack_1.3-3.3 rpart_4.1-9 evaluate_0.5.5 sandwich_2.3-2
## [17] foreign_0.8-63 zoo_1.7-11Sys.time()## [1] "2015-05-23 22:43:36 BRT"