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/vindoura2014-15.txt
Download do script: fazVindoura.R
##-----------------------------------------------------------------------------
## Se o arquivo não existir, descomente aqui para fazer o download.
## url <- "http://www.leg.ufpr.br/~walmes/data/vindoura2014-15.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" "vindoura2014-15.txt"## Encoding do arquivo. Função do sistema LINUX.
system("file vindoura2014-15.txt")
## Lendo o arquivo.
da <- read.table(file="vindoura2014-15.txt", header=TRUE, sep="\t",
colClasses=c("factor","factor","numeric"),
encoding="utf-8")
str(da)## 'data.frame': 100 obs. of 3 variables:
## $ bloc: Factor w/ 10 levels "1","10","2","3",..: 1 3 4 5 6 7 8 9 10 2 ...
## $ trat: Factor w/ 10 levels "Celere","Noble",..: 9 9 9 9 9 9 9 9 9 9 ...
## $ prod: num 4694 3956 4368 4662 4036 ...## É 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 1 Testemunha 4694
## 2 2 Testemunha 3956
## 3 3 Testemunha 4368
## 4 4 Testemunha 4662
## 5 5 Testemunha 4036
## 6 6 Testemunha 4231Análise exploratória preliminar
##-----------------------------------------------------------------------------
## Níveis de tratamento.
levels(da$trat)## [1] "Celere" "Noble" "Produquímica" "SJ1-BJ-SATS"
## [5] "SJ2-OMEX" "SJ3-Quimifol" "SJ4-Blend" "Stoller"
## [9] "Testemunha" "Ubyfol"## Para controle ficar como primeiro nível.
da$trat <- relevel(da$trat, ref="Testemunha")
da <- arrange(da, trat, bloc)
## 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
## Testemunha 10 4478 4515 363.7 3956 5146
## Celere 10 4521 4410 292.0 4197 5145
## Noble 10 4559 4561 271.1 4092 4983
## Produquímica 10 4606 4577 426.4 3877 5567
## SJ1-BJ-SATS 10 4643 4663 335.0 3887 5029
## SJ2-OMEX 10 4696 4644 230.9 4410 5135
## SJ3-Quimifol 10 4361 4310 517.8 3725 5257
## SJ4-Blend 10 4560 4421 299.3 4268 5126
## Stoller 10 4386 4361 147.6 4208 4696
## Ubyfol 10 4537 4575 254.6 4086 4850
## All 100 4535 4529 329.4 3725 5567## Gráfico.
dotplot(prod~trat, data=da, groups=bloc, type="b",
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 9 1889947 209994 2.17 0.033 *
## trat 9 995459 110607 1.14 0.345
## Residuals 81 7854552 96970
## ---
## 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 9 1889947 209994 2.15 0.035 *
## trat 9 995459 110607 1.13 0.351
## I(haty^2) 1 27951 27951 0.29 0.594
## Residuals 80 7826601 97833
## ---
## 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 lwr upr trat
## 1 4478 98.47 81 45.47 1.962e-59 4282 4674 Testemunha
## 2 4521 98.47 81 45.91 9.400e-60 4325 4716 Celere
## 3 4559 98.47 81 46.30 4.815e-60 4363 4755 Noble
## 4 4606 98.47 81 46.77 2.194e-60 4410 4801 Produquímica
## 5 4643 98.47 81 47.15 1.167e-60 4447 4839 SJ1-BJ-SATS
## 6 4696 98.47 81 47.69 4.804e-61 4500 4892 SJ2-OMEX
## 7 4361 98.47 81 44.29 1.532e-58 4166 4557 SJ3-Quimifol
## 8 4560 98.47 81 46.30 4.807e-60 4364 4755 SJ4-Blend
## 9 4386 98.47 81 44.54 9.962e-59 4190 4582 Stoller
## 10 4537 98.47 81 46.08 7.024e-60 4341 4733 UbyfolRepresentar 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 Testemunha 4478 4282 4674
## 2 Celere 4521 4325 4716
## 3 Noble 4559 4363 4755
## 4 Produquímica 4606 4410 4801
## 5 SJ1-BJ-SATS 4643 4447 4839
## 6 SJ2-OMEX 4696 4500 4892
## 7 SJ3-Quimifol 4361 4166 4557
## 8 SJ4-Blend 4560 4364 4755
## 9 Stoller 4386 4190 4582
## 10 Ubyfol 4537 4341 4733## 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" "Stoller" "Testemunha" "Celere"
## [5] "Ubyfol" "Noble" "SJ4-Blend" "Produquímica"
## [9] "SJ1-BJ-SATS" "SJ2-OMEX"## 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-6 TH.data_1.0-3 mvtnorm_0.9-9997
## [4] tables_0.7.79 Hmisc_3.14-4 Formula_1.1-1
## [7] doBy_4.5-10 MASS_7.3-39 survival_2.38-1
## [10] plyr_1.8.1 latticeExtra_0.6-26 lattice_0.20-31
## [13] RColorBrewer_1.0-5 rmarkdown_0.5.1 knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.11.0 cluster_2.0.1 splines_3.2.0 minqa_1.2.3
## [5] stringr_0.6.2 tools_3.2.0 nlme_3.1-120 htmltools_0.2.6
## [9] yaml_2.1.13 lme4_1.1-7 digest_0.6.4 Matrix_1.2-0
## [13] nloptr_1.0.4 formatR_1.0 evaluate_0.5.5 sandwich_2.3-0
## [17] zoo_1.7-10Sys.time()## [1] "2015-06-03 19:21:03 BRT"