## -----------------------------------------------------------------------------
## 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
## TRUE TRUE TRUE TRUE TRUE
## plyr wzRfun
## TRUE 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
## [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] splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.1 plyr_1.8.1 reshape_0.8.5
## [4] multcomp_1.3-3 TH.data_1.0-3 mvtnorm_0.9-99992
## [7] doBy_4.5-10 MASS_7.3-33 survival_2.37-7
## [10] 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
## [4] lme4_1.1-6 Matrix_1.1-4 methods_3.1.1
## [7] minqa_1.2.3 nlme_3.1-117 Rcpp_0.11.1
## [10] 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)
## -----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/zimmermann_ql11x11_prop.txt"
url <- "/home/walmes/Dropbox/CursoR/DADOS/zimmermann_ql11x11_prop.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file = url, header = TRUE, sep = "\t", dec = ",")
## Para facilitar, truncar os nomes.
names(da) <- substr(names(da), 1, 4)
## Converte de inteiro para fator e cria nível de unidade experimental.
da <- transform(da, linh = factor(linh), colu = factor(colu), trat = factor(trat),
ue = interaction(linh, colu, trat, drop = TRUE))
## Mostra a estrutura do objeto.
str(da)
## 'data.frame': 484 obs. of 6 variables:
## $ linh: Factor w/ 11 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ colu: Factor w/ 11 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat: Factor w/ 11 levels "1","2","3","4",..: 2 2 2 2 6 6 6 6 3 3 ...
## $ suba: int 1 2 3 4 1 2 3 4 1 2 ...
## $ prop: num 0.9 0.88 0.903 0.875 0.737 ...
## $ ue : Factor w/ 121 levels "9.1.1","5.2.1",..: 12 12 12 12 56 56 56 56 23 23 ...
## Vendo o quadrado no plano.
cast(da, linh ~ colu, value = "trat", fun = length)
## linh 1 2 3 4 5 6 7 8 9 10 11
## 1 1 4 4 4 4 4 4 4 4 4 4 4
## 2 2 4 4 4 4 4 4 4 4 4 4 4
## 3 3 4 4 4 4 4 4 4 4 4 4 4
## 4 4 4 4 4 4 4 4 4 4 4 4 4
## 5 5 4 4 4 4 4 4 4 4 4 4 4
## 6 6 4 4 4 4 4 4 4 4 4 4 4
## 7 7 4 4 4 4 4 4 4 4 4 4 4
## 8 8 4 4 4 4 4 4 4 4 4 4 4
## 9 9 4 4 4 4 4 4 4 4 4 4 4
## 10 10 4 4 4 4 4 4 4 4 4 4 4
## 11 11 4 4 4 4 4 4 4 4 4 4 4
## cast(da, linh~colu, value='prop', fun=mean)
## Tem registros perdidos?
sum(is.na(da$prop))
## [1] 0
## -----------------------------------------------------------------------------
## Análise exploratória.
xyplot(prop ~ trat, data = da, type = c("p", "a"), jitter.x = TRUE)
## A proporção é uma variável resposta limitada. Quando os valores tendem aos
## extremos do intervalo, 0 ou 1, é notável o surgimento de assimetria na
## distribuição. Em muitos casos a análise deve ser feita com a variável
## transformada.
xyplot(prop^5 ~ trat, data = da, type = c("p", "a"), jitter.x = TRUE)
## Ao elevar à uma potência pode-se perceber correção da assimetria. Não
## devemos escolher a transformação por esse gráfico pois está se esquecendo
## do efeito de linha e coluna e focando atenção apenas no efeito de
## tratamentos.
\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_{ijk}+ue_{ij}, \sigma^2) \newline \mu_{ijk} &= \mu+\alpha_i+\beta_j+\tau_k \end{aligned} \]
em que \(\alpha_i\) é o efeito da linha \(i\), \(\beta_j\) é o efeito da coluna \(j\) e \(\tau_j\) o efeito da variedade \(j\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito dos termos mencionados, \(\mu_{ijk}\) é a média para as observações de procedencia \(ijk\) e \(\sigma^2\) é a variância das observações ao redor desse valor médio. O termo \(\text{ue}_{ij}\) representa o efeito aleatório da unidade experimental representada pela combinação de índices de linha com de coluna.
## -----------------------------------------------------------------------------
## Estimação.
## O modelo abaixo declara o efeito de ue como sendo fixo. Na realiadade é
## mais adequado que ele seja aleatório. Está sendo feito para dar acesso aos
## resíduos do modelo e investigar a transformação. Depois disso, um modelo
## mais apropriado pode ser declarado.
m0 <- lm(prop ~ linh + colu + trat + ue, data = da)
deviance(m0)
## [1] 2.128
df.residual(m0)
## [1] 363
## É o modelo pois o espaço de ue contém o espaço de linh+colu+trat.
m0 <- lm(prop ~ ue, data = da)
deviance(m0)
## [1] 2.128
df.residual(m0)
## [1] 363
par(mfrow = c(2, 2))
plot(m0)
layout(1)
## Forte relação variância-esperança, algo que de fato já é esperado para
## proporções se considerarmos os modelos teóricos binomial e beta.
bc <- MASS::boxcox(m0, lambda = seq(3, 6, l = 20))
abline(v = 5, col = 2)
bc$x[which.max(bc$y)]
## [1] 5.091
## Conforme sugerido pela transformação Box-Cox, será analisada a resposta
## tranformada de forma que y_trans = y_orig^lambda, em que lambra = 5.
## -----------------------------------------------------------------------------
## Ajuste na escala transformada.
m1 <- lm((prop^5) ~ ue, data = da)
par(mfrow = c(2, 2))
plot(m1)
layout(1)
## -----------------------------------------------------------------------------
## Declarar o modelo com a aov() para que ue seja entendida como um termo de
## erro do estrato parcela.
m2 <- aov(prop^5 ~ linh + colu + trat + Error(ue), da)
summary(m2)
##
## Error: ue
## Df Sum Sq Mean Sq F value Pr(>F)
## linh 10 3.41 0.341 5.10 7.2e-06 ***
## colu 10 0.86 0.086 1.29 0.25
## trat 10 4.15 0.415 6.19 3.9e-07 ***
## Residuals 90 6.03 0.067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 363 14.3 0.0393
## Como visto anteriormente, objetos de classe 'aovlist' possuem poucos
## métodos dispoíveis. Essa análise pode ser conduzida com a média das
## parcelas e o número de observações como peso.
## -----------------------------------------------------------------------------
## Agregando os dados para médias e pesos. Não esquecer que existem NA.
db <- ddply(da, .(linh, colu, trat), summarise, n = sum(!is.na(prop)), propm = mean(prop,
na.rm = TRUE))
str(db)
## 'data.frame': 121 obs. of 5 variables:
## $ linh : Factor w/ 11 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ colu : Factor w/ 11 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ trat : Factor w/ 11 levels "1","2","3","4",..: 2 11 10 9 5 4 6 8 7 3 ...
## $ n : int 4 4 4 4 4 4 4 4 4 4 ...
## $ propm: num 0.89 0.857 0.869 0.934 0.81 ...
## -----------------------------------------------------------------------------
## Ajuste do modelo.
m2 <- lm(propm ~ linh + colu + trat, db, weight = n)
## Diagnóstico.
par(mfrow = c(2, 2))
plot(m2, which = 1:3)
MASS::boxcox(m2, seq(3, 8, l = 20))
abline(v = 6, col = 2)
layout(1)
m3 <- lm(propm^6 ~ linh + colu + trat, db, weight = n)
par(mfrow = c(2, 2))
plot(m3, which = 1:3)
layout(1)
## Ok!
## -----------------------------------------------------------------------------
## Quadro de análise de variância.
anova(m3)
## Analysis of Variance Table
##
## Response: propm^6
## Df Sum Sq Mean Sq F value Pr(>F)
## linh 10 4.64 0.464 5.54 2.2e-06 ***
## colu 10 1.10 0.110 1.31 0.24
## trat 10 5.48 0.548 6.53 1.6e-07 ***
## Residuals 90 7.55 0.084
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## -----------------------------------------------------------------------------
## Comparações múltiplas.
## Médias ajustadas (LSmeans).
p0 <- LSmeans(m3, effect = "trat", level = 0.95)
## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, db)
p0
## trat estimate se df t.stat p.value lwr upr
## 1 1 0.4695 0.04367 90 10.752 7.927e-18 0.3828 0.5563
## 2 2 0.5876 0.04367 90 13.456 2.890e-23 0.5009 0.6744
## 3 3 0.6328 0.04367 90 14.490 2.960e-25 0.5460 0.7195
## 4 4 0.6421 0.04367 90 14.704 1.168e-25 0.5553 0.7288
## 5 5 0.4577 0.04367 90 10.483 2.855e-17 0.3710 0.5445
## 6 6 0.3254 0.04367 90 7.453 5.389e-11 0.2387 0.4122
## 7 7 0.6073 0.04367 90 13.906 3.875e-24 0.5205 0.6940
## 8 8 0.5845 0.04367 90 13.385 3.986e-23 0.4977 0.6712
## 9 9 0.7050 0.04367 90 16.144 2.562e-28 0.6182 0.7917
## 10 10 0.4683 0.04367 90 10.724 9.042e-18 0.3816 0.5551
## 11 11 0.6374 0.04367 90 14.597 1.860e-25 0.5507 0.7242
## Comparações múltiplas, contrastes de Tukey. Método de correção de
## p-valor: single-step.
g0 <- summary(glht(m3, linfct = mcp(trat = "Tukey")), test = adjusted(type = "fdr"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = propm^6 ~ linh + colu + trat, data = db, weights = n)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 0.11809 0.06176 1.91 0.10743
## 3 - 1 == 0 0.16324 0.06176 2.64 0.02803 *
## 4 - 1 == 0 0.17256 0.06176 2.79 0.02330 *
## 5 - 1 == 0 -0.01177 0.06176 -0.19 0.94942
## 6 - 1 == 0 -0.14407 0.06176 -2.33 0.05730 .
## 7 - 1 == 0 0.13774 0.06176 2.23 0.06465 .
## 8 - 1 == 0 0.11496 0.06176 1.86 0.10987
## 9 - 1 == 0 0.23546 0.06176 3.81 0.00138 **
## 10 - 1 == 0 -0.00121 0.06176 -0.02 0.98439
## 11 - 1 == 0 0.16789 0.06176 2.72 0.02544 *
## 3 - 2 == 0 0.04516 0.06176 0.73 0.59676
## 4 - 2 == 0 0.05448 0.06176 0.88 0.53598
## 5 - 2 == 0 -0.12986 0.06176 -2.10 0.08096 .
## 6 - 2 == 0 -0.26216 0.06176 -4.25 0.00049 ***
## 7 - 2 == 0 0.01965 0.06176 0.32 0.86062
## 8 - 2 == 0 -0.00312 0.06176 -0.05 0.97754
## 9 - 2 == 0 0.11737 0.06176 1.90 0.10743
## 10 - 2 == 0 -0.11930 0.06176 -1.93 0.10721
## 11 - 2 == 0 0.04980 0.06176 0.81 0.56622
## 4 - 3 == 0 0.00932 0.06176 0.15 0.94942
## 5 - 3 == 0 -0.17501 0.06176 -2.83 0.02330 *
## 6 - 3 == 0 -0.30731 0.06176 -4.98 4.3e-05 ***
## 7 - 3 == 0 -0.02551 0.06176 -0.41 0.81372
## 8 - 3 == 0 -0.04828 0.06176 -0.78 0.57146
## 9 - 3 == 0 0.07222 0.06176 1.17 0.38550
## 10 - 3 == 0 -0.16445 0.06176 -2.66 0.02803 *
## 11 - 3 == 0 0.00465 0.06176 0.08 0.97564
## 5 - 4 == 0 -0.18434 0.06176 -2.98 0.01826 *
## 6 - 4 == 0 -0.31664 0.06176 -5.13 4.2e-05 ***
## 7 - 4 == 0 -0.03483 0.06176 -0.56 0.71774
## 8 - 4 == 0 -0.05760 0.06176 -0.93 0.51159
## 9 - 4 == 0 0.06290 0.06176 1.02 0.46254
## 10 - 4 == 0 -0.17377 0.06176 -2.81 0.02330 *
## 11 - 4 == 0 -0.00467 0.06176 -0.08 0.97564
## 6 - 5 == 0 -0.13230 0.06176 -2.14 0.07671 .
## 7 - 5 == 0 0.14951 0.06176 2.42 0.04809 *
## 8 - 5 == 0 0.12674 0.06176 2.05 0.08770 .
## 9 - 5 == 0 0.24723 0.06176 4.00 0.00088 ***
## 10 - 5 == 0 0.01056 0.06176 0.17 0.94942
## 11 - 5 == 0 0.17966 0.06176 2.91 0.02092 *
## 7 - 6 == 0 0.28181 0.06176 4.56 0.00017 ***
## 8 - 6 == 0 0.25904 0.06176 4.19 0.00050 ***
## 9 - 6 == 0 0.37953 0.06176 6.15 1.2e-06 ***
## 10 - 6 == 0 0.14286 0.06176 2.31 0.05746 .
## 11 - 6 == 0 0.31196 0.06176 5.05 4.2e-05 ***
## 8 - 7 == 0 -0.02277 0.06176 -0.37 0.83456
## 9 - 7 == 0 0.09773 0.06176 1.58 0.18935
## 10 - 7 == 0 -0.13895 0.06176 -2.25 0.06430 .
## 11 - 7 == 0 0.03016 0.06176 0.49 0.76575
## 9 - 8 == 0 0.12050 0.06176 1.95 0.10635
## 10 - 8 == 0 -0.11617 0.06176 -1.88 0.10858
## 11 - 8 == 0 0.05293 0.06176 0.86 0.54132
## 10 - 9 == 0 -0.23667 0.06176 -3.83 0.00138 **
## 11 - 9 == 0 -0.06757 0.06176 -1.09 0.42289
## 11 - 10 == 0 0.16910 0.06176 2.74 0.02544 *
## ---
## 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, decreasing = TRUE)
l0
## 1 2 3 4 5 6 7 8 9 10 11
## "bd" "abc" "a" "a" "cd" "d" "ab" "abc" "a" "bd" "a"
## -----------------------------------------------------------------------------
## Representação gráfica dos resultados na escala da proporção.
i <- c("estimate", "lwr", "upr")
p0[, i] <- p0[, i]^(1/6)
## Formatando média seguida de letras.
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
xlab <- "Proporção de hastes viáveis"
ylab <- "Tratamento das hastes"
sub <- list("Médias seguidas de mesma letra indicam diferença nula à 5%.",
font = 1, cex = 0.8)
p0 <- transform(p0, trat = reorder(trat, estimate))
p0 <- p0[order(p0$trat), ]
xlim <- c(0.75, 1)
segplot(trat ~ lwr + upr, centers = estimate, data = p0, draw = FALSE, xlim = xlim,
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.integer(z), x = xlim[1], label = letras, pos = 4,
cex = 0.8)
})
## -----------------------------------------------------------------------------
## Esse seria o resultado se não fosse feita a transformação da resposta.
p0 <- LSmeans(m2, effect = "trat", level = 0.95)
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, db)
g0 <- summary(glht(m3, linfct = mcp(trat = "Tukey")), test = adjusted(type = "fdr"))
l0 <- cld(g0, decreasing = TRUE)
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
p0 <- transform(p0, trat = reorder(trat, estimate))
p0 <- p0[order(p0$trat), ]
xlim <- c(0.75, 1)
segplot(trat ~ lwr + upr, centers = estimate, data = p0, draw = FALSE, xlim = xlim,
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.integer(z), x = xlim[1], label = letras, pos = 4,
cex = 0.8)
})
