Experimento em delineamento interiramente casualizado
##-----------------------------------------------------------------------------
## 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
## TRUEsource("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] methods splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.3 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-7
## [5] TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-11 MASS_7.3-34
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] rmarkdown_0.3.3 knitr_1.7
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 grid_3.1.1 htmltools_0.2.6
## [6] Matrix_1.1-4 Rcpp_0.11.3 sandwich_2.3-0 stringr_0.6.2 tools_3.1.1
## [11] yaml_2.1.13 zoo_1.7-10## 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 <- "/home/walmes/Dropbox/CursoR/cnpaf1/dados/volume.txt"
##url <- "http://www.leg.ufpr.br/~walmes/data/volume.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t", nrow=27)
## Mostra a estrutura do objeto.
str(da)## 'data.frame': 27 obs. of 4 variables:
## $ gen : Factor w/ 9 levels "ATF06B","ATF40B",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ dose: int 0 0 0 0 0 0 0 0 0 0 ...
## $ rep : int 1 2 3 1 2 3 1 2 3 1 ...
## $ volu: num 1.19 1.17 1.12 1.27 1.85 ...## Os dados no formato amplo (não empilhado).
unstack(da, form=volu~gen)## ATF06B ATF40B ATF54B BR001B BR005B BR007B BR008B P9401 SC283
## 1 1.190 1.275 0.939 1.421 0.520 1.431 1.398 1.141 0.872
## 2 1.173 1.849 1.056 1.240 0.567 1.425 1.552 1.032 1.093
## 3 1.118 1.306 0.712 1.614 0.416 1.152 1.127 1.226 1.136## Os dados são de um experimento para verificar o ganho de peso em
## função de tipos de rações. Foram usandos 4 níveis do fator ração e 5
## repetições.
##-----------------------------------------------------------------------------
## Análise exploratória.
## xyplot(volu~gen, data=da, type=c("p","a"))
xyplot(volu~reorder(gen, volu), data=da, type=c("p","a"))
Especificação e estimação do modelo
\[ \begin{aligned} Y|\text{fontes de variação} &\,\sim \text{Normal}(\mu_i, \sigma^2) \newline \mu_i &= \mu+\alpha_i \end{aligned} \]
em que \(\alpha_i\) é o efeito do genótipo \(i\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito dos genótipos, \(\mu_i\) é a média para as observações em cada nível de genótipo e \(\sigma^2\) é a variância das observações ao redor desse valor médio. Para fazer a estimação desses parâmetros é necessário que seja considerada uam restrição paramétrica haja visto que o sistema admite infinitas soluções pois existem 10 parâmetros (\(\mu\) e \(\alpha_i\) com \(i=1,\ldots,9\) sendo que existem apenas 9 grupos. Em outras palavras, busca-se a solução para 10 icognitas com apenas 9 equações. O R considera o efeito do primeiro nível como zero (\(\alpha_1 = 0\)). O aplicativo da concorrência considera o último nível com efeito igual à zero. Além dessas possibilidades, ainda tem-se a famosa retrição do tipo soma zero (\(\sum \alpha_i = 0\)) e para essa situação \(\mu\) representa a média do experimento.
##-----------------------------------------------------------------------------
## Estimação.
## Especificação do modelo.
m0 <- lm(volu~gen, data=da)
## Estimativas e erros padrões.
summary(m0)##
## Call:
## lm(formula = volu ~ gen, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2320 -0.1313 0.0190 0.0910 0.3723
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.16033 0.10154 11.427 1.1e-09 ***
## genATF40B 0.31633 0.14360 2.203 0.040877 *
## genATF54B -0.25800 0.14360 -1.797 0.089198 .
## genBR001B 0.26467 0.14360 1.843 0.081854 .
## genBR005B -0.65933 0.14360 -4.591 0.000227 ***
## genBR007B 0.17567 0.14360 1.223 0.236998
## genBR008B 0.19867 0.14360 1.383 0.183448
## genP9401 -0.02733 0.14360 -0.190 0.851173
## genSC283 -0.12667 0.14360 -0.882 0.389371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1759 on 18 degrees of freedom
## Multiple R-squared: 0.8031, Adjusted R-squared: 0.7156
## F-statistic: 9.177 on 8 and 18 DF, p-value: 5.458e-05## Restrição paramétrica de zerar o primeiro nível.
contrasts(da$gen)## ATF40B ATF54B BR001B BR005B BR007B BR008B P9401 SC283
## ATF06B 0 0 0 0 0 0 0 0
## ATF40B 1 0 0 0 0 0 0 0
## ATF54B 0 1 0 0 0 0 0 0
## BR001B 0 0 1 0 0 0 0 0
## BR005B 0 0 0 1 0 0 0 0
## BR007B 0 0 0 0 1 0 0 0
## BR008B 0 0 0 0 0 1 0 0
## P9401 0 0 0 0 0 0 1 0
## SC283 0 0 0 0 0 0 0 1## model.matrix(m0)
##-----------------------------------------------------------------------------
## Existe efeito de genótipos? Ou seja, todos os \alpha_i podem ser
## considerados iguais à zero?
## Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: volu
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 8 2.27107 0.283884 9.1775 5.458e-05 ***
## Residuals 18 0.55679 0.030933
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Mas essa análise é confiável uma vez que tem-se a suspeita de
## violação do pressuposto de homogeneidade de variâncias?Diagnóstico do modelo
##-----------------------------------------------------------------------------
## Análise de diagnóstico dos pressupostos via resíduos.
## As inferências acima só passam a ser valiadas após confirmação de não
## haver afastamento dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
A tranformação Box-Cox baseia-se na escolha do valor \(\lambda\) de tal maneira que se aumente a compatibilidade com a suposição de normalidade. A transformação aplicada aos é
\[ \begin{aligned} y_{\text{transformado}} &= \frac{y^{\lambda}-1}{\lambda}, \quad \text{ se } \lambda\neq 0 \newline y_{\text{transformado}} &= \log(y), \quad \text{ se } \lambda=0. \end{aligned} \]
A transformação só está definida para \(y>0\). Caso tenha-se zeros ou valores negativos pode-se somar uma constante aos dados. Como o que interessa da transformação é parte não linear, quando a transformação indicada é a potência, na prática apenas se aplica a função potência sem subtrair 1 e dividir por \(\lambda\). No entanto, essa é a transformação de fato aplicada para calcular o valor da log-verossimilhança. Se houver a curiosidade de fazer o cálculo da log-verrossimilhança passo a passo, lembre-se de considerar o jacobiano da transformação.
##-----------------------------------------------------------------------------
## Transformar os dados pode diminuir a fuga dos pressupostos?
## Chama a função box-cox do pacote MASS sem precisar carregá-lo.
MASS::boxcox(m0)
abline(v=c(0, 1/3, 1/2), col=c(2:4), lty=2)
text(x=c(0,1/3,1/2), y=0,
## labels=c("log","raíz cúbica","raiz quadrada"),
labels=c(
expression(y[t]==log(y)),
expression(y[t]==y^{1/3}),
expression(y[t]==y^{1/2})),
srt=90, adj=c(0.5,-0.25))
##-----------------------------------------------------------------------------
## Refazer a análise com a variável transformada.
m1 <- lm(log(volu)~gen, data=da)
## Diagnóstico.
par(mfrow=c(2,2)); plot(m1); layout(1)
## Quadro de anova.
anova(m1)## Analysis of Variance Table
##
## Response: log(volu)
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 8 2.70970 0.33871 15.388 1.366e-06 ***
## Residuals 18 0.39621 0.02201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Passa m1 para m0.
m0 <- m1Comparações múltiplas
As comparações múltiplas são uma classe particular de hipóteses chamadas de hipóteses lineares sobre os parâmetros. São hipóteses sobre os contrastes entre médias (\(\mu_i - \mu_j = 0\)). Nesse material dá se preferência para os métodos disponíveis para avaliação dessa classe de hipóteses que sejam gerais, ou seja, que possam ser aplicadados à experimentos que apresentam estrutura irregular, como os desbalanceados e com perda de celas experimentais ou na presença de covariáveis (ancova) além de experimentos com respostas que não sejam gaussianas, como contagens (Poisson) e de proporção (binomial). Testes como o de Tukey e similares se baseiam em uma dms (diferença mínima significativa) constante para indicar diferença de médias. Em experimentos desbalanceados as médias ajustadas tem erro-padrão que depende do número de repetições e com respostas não gaussianas existe sempre alguma relação média-variância. Sendo assim, não se tem como ter uma dms constante. Os métodos gerais de avaliação de testes de hipótese linear são gearais e o controle da multiplicidade (o fato de se avaliar mais de uma hipótese, uma série/família) que implica na elevação do erro tipo I global é feito por métodos de correção do p-valor de cada hipótese.
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## Médias ajustadas (LSmeans).
p0 <- LSmeans(m0, effect="gen", level=0.95)
## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, da)
p0## gen estimate se df t.stat p.value lwr upr
## 1 ATF06B 0.14835308 0.08565725 18 1.7319384 0.1003860664 -0.03160612 0.32831229
## 2 ATF40B 0.37485339 0.08565725 18 4.3762015 0.0003641404 0.19489418 0.55481259
## 3 ATF54B -0.11604299 0.08565725 18 -1.3547364 1.8077373441 -0.29600220 0.06391621
## 4 BR001B 0.34839593 0.08565725 18 4.0673257 0.0007229191 0.16843673 0.52835514
## 5 BR005B -0.69946415 0.08565725 18 -8.1658488 1.9999998173 -0.87942336 -0.51950495
## 6 BR007B 0.28468163 0.08565725 18 3.3234972 0.0037799898 0.10472242 0.46464083
## 7 BR008B 0.29804877 0.08565725 18 3.4795510 0.0026754616 0.11808956 0.47800797
## 8 P9401 0.12238686 0.08565725 18 1.4287974 0.1701868335 -0.05757235 0.30234606
## 9 SC283 0.02649122 0.08565725 18 0.3092701 0.7606674493 -0.15346798 0.20645043## Comparações múltiplas, contrastes de Tukey.
## Método de correção de p-valor: False Discovery Rate (fdr).
g0 <- summary(glht(m0, linfct=mcp(gen="Tukey")),
test=adjusted(type="fdr"))
g0##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = log(volu) ~ gen, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## ATF40B - ATF06B == 0 0.22650 0.12114 1.870 0.134509
## ATF54B - ATF06B == 0 -0.26440 0.12114 -2.183 0.095751 .
## BR001B - ATF06B == 0 0.20004 0.12114 1.651 0.189821
## BR005B - ATF06B == 0 -0.84782 0.12114 -6.999 1.12e-05 ***
## BR007B - ATF06B == 0 0.13633 0.12114 1.125 0.366935
## BR008B - ATF06B == 0 0.14970 0.12114 1.236 0.334710
## P9401 - ATF06B == 0 -0.02597 0.12114 -0.214 0.856472
## SC283 - ATF06B == 0 -0.12186 0.12114 -1.006 0.421404
## ATF54B - ATF40B == 0 -0.49090 0.12114 -4.052 0.002989 **
## BR001B - ATF40B == 0 -0.02646 0.12114 -0.218 0.856472
## BR005B - ATF40B == 0 -1.07432 0.12114 -8.869 1.42e-06 ***
## BR007B - ATF40B == 0 -0.09017 0.12114 -0.744 0.559509
## BR008B - ATF40B == 0 -0.07680 0.12114 -0.634 0.620173
## P9401 - ATF40B == 0 -0.25247 0.12114 -2.084 0.103340
## SC283 - ATF40B == 0 -0.34836 0.12114 -2.876 0.027853 *
## BR001B - ATF54B == 0 0.46444 0.12114 3.834 0.004377 **
## BR005B - ATF54B == 0 -0.58342 0.12114 -4.816 0.000623 ***
## BR007B - ATF54B == 0 0.40072 0.12114 3.308 0.011735 *
## BR008B - ATF54B == 0 0.41409 0.12114 3.418 0.010029 *
## P9401 - ATF54B == 0 0.23843 0.12114 1.968 0.122470
## SC283 - ATF54B == 0 0.14253 0.12114 1.177 0.352619
## BR005B - BR001B == 0 -1.04786 0.12114 -8.650 1.42e-06 ***
## BR007B - BR001B == 0 -0.06371 0.12114 -0.526 0.680999
## BR008B - BR001B == 0 -0.05035 0.12114 -0.416 0.744656
## P9401 - BR001B == 0 -0.22601 0.12114 -1.866 0.134509
## SC283 - BR001B == 0 -0.32190 0.12114 -2.657 0.041240 *
## BR007B - BR005B == 0 0.98415 0.12114 8.124 1.77e-06 ***
## BR008B - BR005B == 0 0.99751 0.12114 8.235 1.77e-06 ***
## P9401 - BR005B == 0 0.82185 0.12114 6.784 1.41e-05 ***
## SC283 - BR005B == 0 0.72596 0.12114 5.993 5.88e-05 ***
## BR008B - BR007B == 0 0.01337 0.12114 0.110 0.913355
## P9401 - BR007B == 0 -0.16229 0.12114 -1.340 0.295493
## SC283 - BR007B == 0 -0.25819 0.12114 -2.131 0.099732 .
## P9401 - BR008B == 0 -0.17566 0.12114 -1.450 0.257059
## SC283 - BR008B == 0 -0.27156 0.12114 -2.242 0.090763 .
## SC283 - P9401 == 0 -0.09590 0.12114 -0.792 0.544826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)## Resumo compacto com letras.
l0 <- cld(g0)
l0## ATF06B ATF40B ATF54B BR001B BR005B BR007B BR008B P9401 SC283
## "bd" "d" "b" "d" "a" "cd" "cd" "bd" "bc"## Formatando média seguida de letras.
p0$cld <- l0$mcletters$Letters
p0$let <- sprintf("%.2f %s", p0$estimate, p0$cld)
p0$let## [1] "0.15 bd" "0.37 d" "-0.12 b" "0.35 d" "-0.70 a" "0.28 cd" "0.30 cd" "0.12 bd"
## [9] "0.03 bc"##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- "Volume das raízes (log)"
ylab <- "Genótipo de sorgo"
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
p0 <- transform(p0, gen=reorder(gen, estimate))
levels(p0$gen)## [1] "BR005B" "ATF54B" "SC283" "P9401" "ATF06B" "BR007B" "BR008B" "BR001B" "ATF40B"p0 <- p0[order(p0$gen),]
segplot(gen~lwr+upr, centers=estimate,
data=p0, draw=FALSE,
xlab=xlab, ylab=ylab, sub=sub, letras=p0$let,
panel=function(x, y, z, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers, ...)
panel.text(y=as.numeric(z), x=centers,
label=letras, pos=3)
})
##-----------------------------------------------------------------------------
## Representando na escala original.
sel <- c("estimate","lwr","upr")
p0[,sel] <- exp(p0[,sel])
xlab <- "Volume das raízes"
p0$let <- sprintf("%.2f %s", p0$estimate, p0$cld)
segplot(gen~lwr+upr, centers=estimate,
data=p0, draw=FALSE,
xlab=xlab, ylab=ylab, sub=sub, letras=p0$let,
panel=function(x, y, z, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers, ...)
panel.text(y=as.numeric(z), x=centers,
label=letras, pos=3)
})
Por se tratar de um experimento regular (não desbalanceado), pode-se fazer essa análise com o pacote ExpDes. Deve ser mecionado que é aconselhável antes declarar o modelo na função lm() e fazer a análise dos resíduos. Apenas depois de confirmada a adequação dos pressupostos, feita as transformações, caso necessárias, fazer a análise com o pacote ExpDes.
##-----------------------------------------------------------------------------
## Análise com o pacote ExpDes.
require(ExpDes)## Loading required package: ExpDes
##
## Attaching package: 'ExpDes'
##
## The following object is masked from 'package:MASS':
##
## ginv## Funções do pacote ExpDes.
ls("package:ExpDes")## [1] "ccboot" "crd" "duncan" "fat2.ad.crd" "fat2.ad.rbd"
## [6] "fat2.crd" "fat2.rbd" "fat3.ad.crd" "fat3.ad.rbd" "fat3.crd"
## [11] "fat3.rbd" "ginv" "lastC" "latsd" "lsd"
## [16] "lsdb" "order.group" "order.stat.SNK" "rbd" "reg.poly"
## [21] "scottknott" "snk" "split2.crd" "split2.rbd" "tapply.stat"
## [26] "tukey"## Complete randomized design.
with(da, crd(treat=gen, resp=log(volu), quali=TRUE, mcomp="tukey"))## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Treatament 8 2.70970 0.33871 15.388 1.3657e-06
## Residuals 18 0.39621 0.02201
## Total 26 3.10591
## ------------------------------------------------------------------------
## CV = 169.51 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.300618
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a ATF40B 0.3748534
## a BR001B 0.3483959
## ab BR008B 0.2980488
## ab BR007B 0.2846816
## ab ATF06B 0.1483531
## ab P9401 0.1223869
## ab SC283 0.02649122
## b ATF54B -0.116043
## c BR005B -0.6994642
## ------------------------------------------------------------------------with(da, crd(treat=gen, resp=log(volu), quali=TRUE, mcomp="sk"))## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## Treatament 8 2.70970 0.33871 15.388 1.3657e-06
## Residuals 18 0.39621 0.02201
## Total 26 3.10591
## ------------------------------------------------------------------------
## CV = 169.51 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.300618
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## Scott-Knott test
## ------------------------------------------------------------------------
## Groups Treatments Means
## 1 a ATF40B 0.37485339
## 2 a BR001B 0.34839593
## 3 a BR008B 0.29804877
## 4 a BR007B 0.28468163
## 5 b ATF06B 0.14835308
## 6 b P9401 0.12238686
## 7 b SC283 0.02649122
## 8 b ATF54B -0.11604299
## 9 c BR005B -0.69946415
## ------------------------------------------------------------------------Embora seja comum aplicar testes de normalidade aos resíduos de um modelo, deve ser mencionado que os resíduos não são uma amostra aleatória independente como é exigido pelo teste de normalidade. Os resíduos a serem considerados para tal devem ser os resíduos studentizados. Ainda assim, o número de parâmetros estimados para chegar aos resíduos não é uma informação passada para o teste. Será que é valido aplicar o teste aos resíduos então? Qual teste aplicar não parece subjetivo? E se o teste rejeitar a normalidade, não serão necessários fazer gráficos para saber que providência tomar? Sendo assim, por que não conduzir apenas uma análise gráfica? Argumentos semenlhantes são considerados para os testes de homogenidade de variâncias.