Ajuste de modelos lineares e mistos no ambiente R |
|
| 09 e 10 de Outubro de 2014 - Piracicaba - SP |
Prof. Dr. Walmes M. Zeviani |
| Escola Superior de Agricultura “Luiz de Queiroz” - USP |
Lab. de Estatística e Geoinformação - LEG |
| Pós Graduação em Genética e Melhoramento de Plantas | Departamento de Estatística - UFPR |
##-----------------------------------------------------------------------------
## 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-6
## [5] TH.data_1.0-3 mvtnorm_1.0-0 doBy_4.5-10 MASS_7.3-34
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 lme4_1.1-7 Matrix_1.1-4
## [6] methods_3.1.1 minqa_1.2.3 nlme_3.1-117 nloptr_1.0.4 Rcpp_0.11.2
## [11] sandwich_2.3-1 stringr_0.6.2 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/ramalho_feijao1.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t",
colClasses=c("factor","factor",NA,NA))
## Para facilitar, abreviar os nomes.
names(da) <- abbreviate(names(da))
## Mostra a estrutura do objeto.
str(da)
## 'data.frame': 240 obs. of 4 variables:
## $ rptc: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ faml: Factor w/ 24 levels "1","10","11",..: 1 12 18 19 20 21 22 23 24 2 ...
## $ plnt: int 1 1 1 1 1 1 1 1 1 1 ...
## $ prdc: num 4.5 5 6.4 3.3 0.9 1.4 4.3 7.4 4.9 1.6 ...
da$faml <- with(da, {
neword <- levels(faml)[order(as.integer(levels(faml)))]
factor(faml, levels=neword)
})
## Tabela de frequência.
xtabs(~rptc+faml, data=da)
## faml
## rptc 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## 2 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## Não são esses blocos grandes demais para serem considerados
## homogêneos?
## Tem registros perdidos?
sum(is.na(da$prdc))
## [1] 3
##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(prdc~faml|rptc, data=da, type=c("p","a"))
\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_{ij}+\text{ue}_{ij}, \sigma^2) \newline \mu_{ij} &= \mu+\alpha_i+\tau_j \newline \text{ue}_{ij} &\sim \text{Normal}(0, \sigma_{ue}^2) \end{aligned} \]
em que \(\alpha_i\) é o efeito do bloco \(i\) e \(\tau_j\) o efeito da família \(j\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito de bloco e variedade, \(\mu_{ij}\) é a média para as observações em cada combinação de níveis de bloco e família e \(\sigma^2\) é a variância das observações ao redor desse valor médio. O termo \(\text{ue}_{ij(k)}\) representa o efeito aleatório da unidade experimental.
##-----------------------------------------------------------------------------
## Estimação.
## Esse modelo assume que Var(ue) == 0. Pode ser uma suposição
## forte. Na prática esse modelo não é útil pois não representa
## adequadamente às fontes de variação presentes. Será considerado aqui
## por razões didáticas.
m0 <- lm(prdc~rptc+faml, data=da)
## Declarando o modelo especificando termo de efeito aleatório que
## corresponde à um estrato de variação.
da$ue <- with(da, interaction(rptc, faml))
m0 <- aov(prdc~rptc+faml+Error(ue), data=da)
##-----------------------------------------------------------------------------
## Quadro de análise de variância estratificado.
summary(m0)
##
## Error: ue
## Df Sum Sq Mean Sq F value Pr(>F)
## rptc 1 3 2.77 0.13 0.72
## faml 23 631 27.43 1.30 0.27
## Residuals 23 487 21.18
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 189 2377 12.6
##-----------------------------------------------------------------------------
## A classe 'aovlist' possuí poucos metodos disponíveis, o que a torna menos
## atraente.
class(m0)
## [1] "aovlist" "listof"
methods(class="aovlist")
## [1] dummy.coef.aovlist* model.frame.aovlist* model.tables.aovlist*
## [4] print.aovlist* proj.aovlist* se.contrast.aovlist*
## [7] summary.aovlist*
##
## Non-visible functions are asterisked
## Análise dos resíduos e comparações múltiplas são coisas complicadas
## de fazer com objetos de classe 'aovlist'. Nesse caso, pode-se reduzir
## os dados à totais ou médias e usar o número de observações como
## ponderador na análise.
## Por curiosidade, pode-se declarar o modelo com os termos todos de
## efeito fixo. O quadro de anova será igual exceto que não será
## estratificado em com valores de F obtidos usando o QMR com
## denominador (o que está errado). Isso permite acessar os resíduos
## porque o modelo será de classe 'lm' e é só para isso que se deve
## usá-lo.
m1 <- lm(prdc~rptc+faml+ue, data=da)
## Diganóstico dos resíduos.
par(mfrow=c(2,2)); plot(m1, which=1:3)
MASS::boxcox(m1)
abline(v=1/3, col=2); layout(1)
## Aponta clara relação esperança-variância.
##-----------------------------------------------------------------------------
## Agregando os dados para médias e pesos. Não esquecer que existem NA.
db <- ddply(da, .(rptc, faml), summarise,
n=sum(!is.na(prdc)),
prdcm=mean(prdc, na.rm=TRUE))
str(db)
## 'data.frame': 48 obs. of 4 variables:
## $ rptc : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ faml : Factor w/ 24 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ n : int 5 5 5 5 5 5 5 5 5 5 ...
## $ prdcm: num 5.38 6.46 8.16 3.88 3.82 3.84 4.28 6.48 3.94 5.4 ...
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m2 <- lm(prdcm~rptc+faml, db, weight=n)
## Diagnóstico.
par(mfrow=c(2,2)); plot(m2, which=1:3)
MASS::boxcox(m2)
abline(v=1/3, col=2); layout(1)
## Cadê aquele afastamento todo?
## Lembre-se de uma propriedade envolvendo ordem de funções lineares e
## não lineares. sqrt(mean(x)) é diferente de mean(sqrt(x)). A ordem
## importa.
## Qual a correspondência entre os quadros de anova?
summary(m0)
##
## Error: ue
## Df Sum Sq Mean Sq F value Pr(>F)
## rptc 1 3 2.77 0.13 0.72
## faml 23 631 27.43 1.30 0.27
## Residuals 23 487 21.18
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 189 2377 12.6
anova(m2)
## Analysis of Variance Table
##
## Response: prdcm
## Df Sum Sq Mean Sq F value Pr(>F)
## rptc 1 3 2.77 0.13 0.72
## faml 23 631 27.43 1.30 0.27
## Residuals 23 487 21.18
## O primeiro estrato é igual ao quadro de anova usando as médias e
## pesos. A vantagem é ter um objeto de classe 'lm' que tem muito mais
## métodos disponíveis.
methods(class="lm")
## [1] add1.lm* addterm.lm* alias.lm* anova.lm*
## [5] attrassign.lm* boxcox.lm* case.names.lm* confint.lm
## [9] cooks.distance.lm* deviance.lm* dfbeta.lm* dfbetas.lm*
## [13] drop1.lm* dropterm.lm* dummy.coef.lm effects.lm*
## [17] esticon.lm extractAIC.lm* family.lm* formula.lm*
## [21] hatvalues.lm* influence.lm* kappa.lm labels.lm*
## [25] linest.lm logLik.lm* logtrans.lm* model.frame.lm*
## [29] model.matrix.lm nobs.lm* plot.lm* predict.lm
## [33] print.lm* proj.lm* qqnorm.lm* qr.lm*
## [37] residuals.lm rstandard.lm* rstudent.lm* simulate.lm*
## [41] summary.lm variable.names.lm* vcov.lm*
##
## Non-visible functions are asterisked
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## Haja visto que o F da análise de variância não rejeitou a hipótese do
## efeito de família ser nulo, não há necessidade de conduzir
## comparações múltiplas. Por razões didáticas, o código para isso será
## mantido.
## Médias ajustadas (LSmeans).
p0 <- LSmeans(m2, effect="faml", level=0.95)
## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, da)
p0
## faml estimate se df t.stat p.value lwr upr
## 1 1 5.600 1.455 23 3.848 8.197e-04 2.5896 8.610
## 2 2 7.760 1.455 23 5.332 2.053e-05 4.7496 10.770
## 3 3 7.303 1.534 23 4.760 8.473e-05 4.1291 10.477
## 4 4 3.970 1.455 23 2.728 1.199e-02 0.9596 6.980
## 5 5 4.780 1.455 23 3.285 3.248e-03 1.7696 7.790
## 6 6 4.540 1.455 23 3.120 4.817e-03 1.5296 7.550
## 7 7 6.859 1.534 23 4.470 1.744e-04 3.6846 10.033
## 8 8 7.010 1.455 23 4.817 7.348e-05 3.9996 10.020
## 9 9 4.930 1.455 23 3.388 2.533e-03 1.9196 7.940
## 10 10 5.150 1.455 23 3.539 1.754e-03 2.1396 8.160
## 11 11 3.675 1.534 23 2.395 2.516e-02 0.5007 6.849
## 12 12 5.020 1.455 23 3.450 2.180e-03 2.0096 8.030
## 13 13 10.420 1.455 23 7.160 2.722e-07 7.4096 13.430
## 14 14 5.820 1.455 23 3.999 5.636e-04 2.8096 8.830
## 15 15 4.030 1.455 23 2.769 1.091e-02 1.0196 7.040
## 16 16 4.990 1.455 23 3.429 2.292e-03 1.9796 8.000
## 17 17 6.270 1.455 23 4.308 2.610e-04 3.2596 9.280
## 18 18 6.280 1.455 23 4.315 2.566e-04 3.2696 9.290
## 19 19 9.170 1.455 23 6.301 1.981e-06 6.1596 12.180
## 20 20 4.290 1.455 23 2.948 7.219e-03 1.2796 7.300
## 21 21 4.650 1.455 23 3.195 4.024e-03 1.6396 7.660
## 22 22 4.800 1.455 23 3.298 3.143e-03 1.7896 7.810
## 23 23 6.790 1.455 23 4.666 1.071e-04 3.7796 9.800
## 24 24 5.210 1.455 23 3.580 1.585e-03 2.1996 8.220
## Comparações múltiplas, contrastes de Tukey.
## Método de correção de p-valor: single-step.
g0 <- summary(glht(m2, linfct=mcp(faml="Tukey")),
test=adjusted(type="fdr"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = prdcm ~ rptc + faml, data = db, weights = n)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 2.1600 2.0581 1.05 0.91
## 3 - 1 == 0 1.7031 2.1147 0.81 0.96
## 4 - 1 == 0 -1.6300 2.0581 -0.79 0.96
## 5 - 1 == 0 -0.8200 2.0581 -0.40 0.98
## 6 - 1 == 0 -1.0600 2.0581 -0.52 0.98
## 7 - 1 == 0 1.2586 2.1147 0.60 0.98
## 8 - 1 == 0 1.4100 2.0581 0.69 0.98
## 9 - 1 == 0 -0.6700 2.0581 -0.33 0.98
## 10 - 1 == 0 -0.4500 2.0581 -0.22 0.98
## 11 - 1 == 0 -1.9253 2.1147 -0.91 0.93
## 12 - 1 == 0 -0.5800 2.0581 -0.28 0.98
## 13 - 1 == 0 4.8200 2.0581 2.34 0.43
## 14 - 1 == 0 0.2200 2.0581 0.11 0.98
## 15 - 1 == 0 -1.5700 2.0581 -0.76 0.97
## 16 - 1 == 0 -0.6100 2.0581 -0.30 0.98
## 17 - 1 == 0 0.6700 2.0581 0.33 0.98
## 18 - 1 == 0 0.6800 2.0581 0.33 0.98
## 19 - 1 == 0 3.5700 2.0581 1.73 0.76
## 20 - 1 == 0 -1.3100 2.0581 -0.64 0.98
## 21 - 1 == 0 -0.9500 2.0581 -0.46 0.98
## 22 - 1 == 0 -0.8000 2.0581 -0.39 0.98
## 23 - 1 == 0 1.1900 2.0581 0.58 0.98
## 24 - 1 == 0 -0.3900 2.0581 -0.19 0.98
## 3 - 2 == 0 -0.4569 2.1147 -0.22 0.98
## 4 - 2 == 0 -3.7900 2.0581 -1.84 0.68
## 5 - 2 == 0 -2.9800 2.0581 -1.45 0.87
## 6 - 2 == 0 -3.2200 2.0581 -1.56 0.84
## 7 - 2 == 0 -0.9014 2.1147 -0.43 0.98
## 8 - 2 == 0 -0.7500 2.0581 -0.36 0.98
## 9 - 2 == 0 -2.8300 2.0581 -1.38 0.87
## 10 - 2 == 0 -2.6100 2.0581 -1.27 0.91
## 11 - 2 == 0 -4.0853 2.1147 -1.93 0.59
## 12 - 2 == 0 -2.7400 2.0581 -1.33 0.87
## 13 - 2 == 0 2.6600 2.0581 1.29 0.89
## 14 - 2 == 0 -1.9400 2.0581 -0.94 0.93
## 15 - 2 == 0 -3.7300 2.0581 -1.81 0.69
## 16 - 2 == 0 -2.7700 2.0581 -1.35 0.87
## 17 - 2 == 0 -1.4900 2.0581 -0.72 0.97
## 18 - 2 == 0 -1.4800 2.0581 -0.72 0.97
## 19 - 2 == 0 1.4100 2.0581 0.69 0.98
## 20 - 2 == 0 -3.4700 2.0581 -1.69 0.79
## 21 - 2 == 0 -3.1100 2.0581 -1.51 0.87
## 22 - 2 == 0 -2.9600 2.0581 -1.44 0.87
## 23 - 2 == 0 -0.9700 2.0581 -0.47 0.98
## 24 - 2 == 0 -2.5500 2.0581 -1.24 0.91
## 4 - 3 == 0 -3.3331 2.1147 -1.58 0.84
## 5 - 3 == 0 -2.5231 2.1147 -1.19 0.91
## 6 - 3 == 0 -2.7631 2.1147 -1.31 0.88
## 7 - 3 == 0 -0.4444 2.1694 -0.20 0.98
## 8 - 3 == 0 -0.2931 2.1147 -0.14 0.98
## 9 - 3 == 0 -2.3731 2.1147 -1.12 0.91
## 10 - 3 == 0 -2.1531 2.1147 -1.02 0.91
## 11 - 3 == 0 -3.6284 2.1704 -1.67 0.79
## 12 - 3 == 0 -2.2831 2.1147 -1.08 0.91
## 13 - 3 == 0 3.1169 2.1147 1.47 0.87
## 14 - 3 == 0 -1.4831 2.1147 -0.70 0.98
## 15 - 3 == 0 -3.2731 2.1147 -1.55 0.85
## 16 - 3 == 0 -2.3131 2.1147 -1.09 0.91
## 17 - 3 == 0 -1.0331 2.1147 -0.49 0.98
## 18 - 3 == 0 -1.0231 2.1147 -0.48 0.98
## 19 - 3 == 0 1.8669 2.1147 0.88 0.93
## 20 - 3 == 0 -3.0131 2.1147 -1.42 0.87
## 21 - 3 == 0 -2.6531 2.1147 -1.25 0.91
## 22 - 3 == 0 -2.5031 2.1147 -1.18 0.91
## 23 - 3 == 0 -0.5131 2.1147 -0.24 0.98
## 24 - 3 == 0 -2.0931 2.1147 -0.99 0.91
## 5 - 4 == 0 0.8100 2.0581 0.39 0.98
## 6 - 4 == 0 0.5700 2.0581 0.28 0.98
## 7 - 4 == 0 2.8886 2.1147 1.37 0.87
## 8 - 4 == 0 3.0400 2.0581 1.48 0.87
## 9 - 4 == 0 0.9600 2.0581 0.47 0.98
## 10 - 4 == 0 1.1800 2.0581 0.57 0.98
## 11 - 4 == 0 -0.2953 2.1147 -0.14 0.98
## 12 - 4 == 0 1.0500 2.0581 0.51 0.98
## 13 - 4 == 0 6.4500 2.0581 3.13 0.35
## 14 - 4 == 0 1.8500 2.0581 0.90 0.93
## 15 - 4 == 0 0.0600 2.0581 0.03 0.99
## 16 - 4 == 0 1.0200 2.0581 0.50 0.98
## 17 - 4 == 0 2.3000 2.0581 1.12 0.91
## 18 - 4 == 0 2.3100 2.0581 1.12 0.91
## 19 - 4 == 0 5.2000 2.0581 2.53 0.35
## 20 - 4 == 0 0.3200 2.0581 0.16 0.98
## 21 - 4 == 0 0.6800 2.0581 0.33 0.98
## 22 - 4 == 0 0.8300 2.0581 0.40 0.98
## 23 - 4 == 0 2.8200 2.0581 1.37 0.87
## 24 - 4 == 0 1.2400 2.0581 0.60 0.98
## 6 - 5 == 0 -0.2400 2.0581 -0.12 0.98
## 7 - 5 == 0 2.0786 2.1147 0.98 0.91
## 8 - 5 == 0 2.2300 2.0581 1.08 0.91
## 9 - 5 == 0 0.1500 2.0581 0.07 0.98
## 10 - 5 == 0 0.3700 2.0581 0.18 0.98
## 11 - 5 == 0 -1.1053 2.1147 -0.52 0.98
## 12 - 5 == 0 0.2400 2.0581 0.12 0.98
## 13 - 5 == 0 5.6400 2.0581 2.74 0.35
## 14 - 5 == 0 1.0400 2.0581 0.51 0.98
## 15 - 5 == 0 -0.7500 2.0581 -0.36 0.98
## 16 - 5 == 0 0.2100 2.0581 0.10 0.98
## 17 - 5 == 0 1.4900 2.0581 0.72 0.97
## 18 - 5 == 0 1.5000 2.0581 0.73 0.97
## 19 - 5 == 0 4.3900 2.0581 2.13 0.54
## 20 - 5 == 0 -0.4900 2.0581 -0.24 0.98
## 21 - 5 == 0 -0.1300 2.0581 -0.06 0.98
## 22 - 5 == 0 0.0200 2.0581 0.01 1.00
## 23 - 5 == 0 2.0100 2.0581 0.98 0.91
## 24 - 5 == 0 0.4300 2.0581 0.21 0.98
## 7 - 6 == 0 2.3186 2.1147 1.10 0.91
## 8 - 6 == 0 2.4700 2.0581 1.20 0.91
## 9 - 6 == 0 0.3900 2.0581 0.19 0.98
## 10 - 6 == 0 0.6100 2.0581 0.30 0.98
## 11 - 6 == 0 -0.8653 2.1147 -0.41 0.98
## 12 - 6 == 0 0.4800 2.0581 0.23 0.98
## 13 - 6 == 0 5.8800 2.0581 2.86 0.35
## 14 - 6 == 0 1.2800 2.0581 0.62 0.98
## 15 - 6 == 0 -0.5100 2.0581 -0.25 0.98
## 16 - 6 == 0 0.4500 2.0581 0.22 0.98
## 17 - 6 == 0 1.7300 2.0581 0.84 0.94
## 18 - 6 == 0 1.7400 2.0581 0.85 0.94
## 19 - 6 == 0 4.6300 2.0581 2.25 0.49
## 20 - 6 == 0 -0.2500 2.0581 -0.12 0.98
## 21 - 6 == 0 0.1100 2.0581 0.05 0.99
## 22 - 6 == 0 0.2600 2.0581 0.13 0.98
## 23 - 6 == 0 2.2500 2.0581 1.09 0.91
## 24 - 6 == 0 0.6700 2.0581 0.33 0.98
## 8 - 7 == 0 0.1514 2.1147 0.07 0.98
## 9 - 7 == 0 -1.9286 2.1147 -0.91 0.93
## 10 - 7 == 0 -1.7086 2.1147 -0.81 0.96
## 11 - 7 == 0 -3.1839 2.1704 -1.47 0.87
## 12 - 7 == 0 -1.8386 2.1147 -0.87 0.93
## 13 - 7 == 0 3.5614 2.1147 1.68 0.79
## 14 - 7 == 0 -1.0386 2.1147 -0.49 0.98
## 15 - 7 == 0 -2.8286 2.1147 -1.34 0.87
## 16 - 7 == 0 -1.8686 2.1147 -0.88 0.93
## 17 - 7 == 0 -0.5886 2.1147 -0.28 0.98
## 18 - 7 == 0 -0.5786 2.1147 -0.27 0.98
## 19 - 7 == 0 2.3114 2.1147 1.09 0.91
## 20 - 7 == 0 -2.5686 2.1147 -1.21 0.91
## 21 - 7 == 0 -2.2086 2.1147 -1.04 0.91
## 22 - 7 == 0 -2.0586 2.1147 -0.97 0.91
## 23 - 7 == 0 -0.0686 2.1147 -0.03 0.99
## 24 - 7 == 0 -1.6486 2.1147 -0.78 0.96
## 9 - 8 == 0 -2.0800 2.0581 -1.01 0.91
## 10 - 8 == 0 -1.8600 2.0581 -0.90 0.93
## 11 - 8 == 0 -3.3353 2.1147 -1.58 0.84
## 12 - 8 == 0 -1.9900 2.0581 -0.97 0.91
## 13 - 8 == 0 3.4100 2.0581 1.66 0.79
## 14 - 8 == 0 -1.1900 2.0581 -0.58 0.98
## 15 - 8 == 0 -2.9800 2.0581 -1.45 0.87
## 16 - 8 == 0 -2.0200 2.0581 -0.98 0.91
## 17 - 8 == 0 -0.7400 2.0581 -0.36 0.98
## 18 - 8 == 0 -0.7300 2.0581 -0.35 0.98
## 19 - 8 == 0 2.1600 2.0581 1.05 0.91
## 20 - 8 == 0 -2.7200 2.0581 -1.32 0.87
## 21 - 8 == 0 -2.3600 2.0581 -1.15 0.91
## 22 - 8 == 0 -2.2100 2.0581 -1.07 0.91
## 23 - 8 == 0 -0.2200 2.0581 -0.11 0.98
## 24 - 8 == 0 -1.8000 2.0581 -0.87 0.93
## 10 - 9 == 0 0.2200 2.0581 0.11 0.98
## 11 - 9 == 0 -1.2553 2.1147 -0.59 0.98
## 12 - 9 == 0 0.0900 2.0581 0.04 0.99
## 13 - 9 == 0 5.4900 2.0581 2.67 0.35
## 14 - 9 == 0 0.8900 2.0581 0.43 0.98
## 15 - 9 == 0 -0.9000 2.0581 -0.44 0.98
## 16 - 9 == 0 0.0600 2.0581 0.03 0.99
## 17 - 9 == 0 1.3400 2.0581 0.65 0.98
## 18 - 9 == 0 1.3500 2.0581 0.66 0.98
## 19 - 9 == 0 4.2400 2.0581 2.06 0.55
## 20 - 9 == 0 -0.6400 2.0581 -0.31 0.98
## 21 - 9 == 0 -0.2800 2.0581 -0.14 0.98
## 22 - 9 == 0 -0.1300 2.0581 -0.06 0.98
## 23 - 9 == 0 1.8600 2.0581 0.90 0.93
## 24 - 9 == 0 0.2800 2.0581 0.14 0.98
## 11 - 10 == 0 -1.4753 2.1147 -0.70 0.98
## 12 - 10 == 0 -0.1300 2.0581 -0.06 0.98
## 13 - 10 == 0 5.2700 2.0581 2.56 0.35
## 14 - 10 == 0 0.6700 2.0581 0.33 0.98
## 15 - 10 == 0 -1.1200 2.0581 -0.54 0.98
## 16 - 10 == 0 -0.1600 2.0581 -0.08 0.98
## 17 - 10 == 0 1.1200 2.0581 0.54 0.98
## 18 - 10 == 0 1.1300 2.0581 0.55 0.98
## 19 - 10 == 0 4.0200 2.0581 1.95 0.59
## 20 - 10 == 0 -0.8600 2.0581 -0.42 0.98
## 21 - 10 == 0 -0.5000 2.0581 -0.24 0.98
## 22 - 10 == 0 -0.3500 2.0581 -0.17 0.98
## 23 - 10 == 0 1.6400 2.0581 0.80 0.96
## 24 - 10 == 0 0.0600 2.0581 0.03 0.99
## 12 - 11 == 0 1.3453 2.1147 0.64 0.98
## 13 - 11 == 0 6.7453 2.1147 3.19 0.35
## 14 - 11 == 0 2.1453 2.1147 1.01 0.91
## 15 - 11 == 0 0.3553 2.1147 0.17 0.98
## 16 - 11 == 0 1.3153 2.1147 0.62 0.98
## 17 - 11 == 0 2.5953 2.1147 1.23 0.91
## 18 - 11 == 0 2.6053 2.1147 1.23 0.91
## 19 - 11 == 0 5.4953 2.1147 2.60 0.35
## 20 - 11 == 0 0.6153 2.1147 0.29 0.98
## 21 - 11 == 0 0.9753 2.1147 0.46 0.98
## 22 - 11 == 0 1.1253 2.1147 0.53 0.98
## 23 - 11 == 0 3.1153 2.1147 1.47 0.87
## 24 - 11 == 0 1.5353 2.1147 0.73 0.97
## 13 - 12 == 0 5.4000 2.0581 2.62 0.35
## 14 - 12 == 0 0.8000 2.0581 0.39 0.98
## 15 - 12 == 0 -0.9900 2.0581 -0.48 0.98
## 16 - 12 == 0 -0.0300 2.0581 -0.01 1.00
## 17 - 12 == 0 1.2500 2.0581 0.61 0.98
## 18 - 12 == 0 1.2600 2.0581 0.61 0.98
## 19 - 12 == 0 4.1500 2.0581 2.02 0.55
## 20 - 12 == 0 -0.7300 2.0581 -0.35 0.98
## 21 - 12 == 0 -0.3700 2.0581 -0.18 0.98
## 22 - 12 == 0 -0.2200 2.0581 -0.11 0.98
## 23 - 12 == 0 1.7700 2.0581 0.86 0.93
## 24 - 12 == 0 0.1900 2.0581 0.09 0.98
## 14 - 13 == 0 -4.6000 2.0581 -2.24 0.49
## 15 - 13 == 0 -6.3900 2.0581 -3.10 0.35
## 16 - 13 == 0 -5.4300 2.0581 -2.64 0.35
## 17 - 13 == 0 -4.1500 2.0581 -2.02 0.55
## 18 - 13 == 0 -4.1400 2.0581 -2.01 0.55
## 19 - 13 == 0 -1.2500 2.0581 -0.61 0.98
## 20 - 13 == 0 -6.1300 2.0581 -2.98 0.35
## 21 - 13 == 0 -5.7700 2.0581 -2.80 0.35
## 22 - 13 == 0 -5.6200 2.0581 -2.73 0.35
## 23 - 13 == 0 -3.6300 2.0581 -1.76 0.74
## 24 - 13 == 0 -5.2100 2.0581 -2.53 0.35
## 15 - 14 == 0 -1.7900 2.0581 -0.87 0.93
## 16 - 14 == 0 -0.8300 2.0581 -0.40 0.98
## 17 - 14 == 0 0.4500 2.0581 0.22 0.98
## 18 - 14 == 0 0.4600 2.0581 0.22 0.98
## 19 - 14 == 0 3.3500 2.0581 1.63 0.81
## 20 - 14 == 0 -1.5300 2.0581 -0.74 0.97
## 21 - 14 == 0 -1.1700 2.0581 -0.57 0.98
## 22 - 14 == 0 -1.0200 2.0581 -0.50 0.98
## 23 - 14 == 0 0.9700 2.0581 0.47 0.98
## 24 - 14 == 0 -0.6100 2.0581 -0.30 0.98
## 16 - 15 == 0 0.9600 2.0581 0.47 0.98
## 17 - 15 == 0 2.2400 2.0581 1.09 0.91
## 18 - 15 == 0 2.2500 2.0581 1.09 0.91
## 19 - 15 == 0 5.1400 2.0581 2.50 0.35
## 20 - 15 == 0 0.2600 2.0581 0.13 0.98
## 21 - 15 == 0 0.6200 2.0581 0.30 0.98
## 22 - 15 == 0 0.7700 2.0581 0.37 0.98
## 23 - 15 == 0 2.7600 2.0581 1.34 0.87
## 24 - 15 == 0 1.1800 2.0581 0.57 0.98
## 17 - 16 == 0 1.2800 2.0581 0.62 0.98
## 18 - 16 == 0 1.2900 2.0581 0.63 0.98
## 19 - 16 == 0 4.1800 2.0581 2.03 0.55
## 20 - 16 == 0 -0.7000 2.0581 -0.34 0.98
## 21 - 16 == 0 -0.3400 2.0581 -0.17 0.98
## 22 - 16 == 0 -0.1900 2.0581 -0.09 0.98
## 23 - 16 == 0 1.8000 2.0581 0.87 0.93
## 24 - 16 == 0 0.2200 2.0581 0.11 0.98
## 18 - 17 == 0 0.0100 2.0581 0.00 1.00
## 19 - 17 == 0 2.9000 2.0581 1.41 0.87
## 20 - 17 == 0 -1.9800 2.0581 -0.96 0.91
## 21 - 17 == 0 -1.6200 2.0581 -0.79 0.96
## 22 - 17 == 0 -1.4700 2.0581 -0.71 0.97
## 23 - 17 == 0 0.5200 2.0581 0.25 0.98
## 24 - 17 == 0 -1.0600 2.0581 -0.52 0.98
## 19 - 18 == 0 2.8900 2.0581 1.40 0.87
## 20 - 18 == 0 -1.9900 2.0581 -0.97 0.91
## 21 - 18 == 0 -1.6300 2.0581 -0.79 0.96
## 22 - 18 == 0 -1.4800 2.0581 -0.72 0.97
## 23 - 18 == 0 0.5100 2.0581 0.25 0.98
## 24 - 18 == 0 -1.0700 2.0581 -0.52 0.98
## 20 - 19 == 0 -4.8800 2.0581 -2.37 0.43
## 21 - 19 == 0 -4.5200 2.0581 -2.20 0.50
## 22 - 19 == 0 -4.3700 2.0581 -2.12 0.54
## 23 - 19 == 0 -2.3800 2.0581 -1.16 0.91
## 24 - 19 == 0 -3.9600 2.0581 -1.92 0.59
## 21 - 20 == 0 0.3600 2.0581 0.17 0.98
## 22 - 20 == 0 0.5100 2.0581 0.25 0.98
## 23 - 20 == 0 2.5000 2.0581 1.21 0.91
## 24 - 20 == 0 0.9200 2.0581 0.45 0.98
## 22 - 21 == 0 0.1500 2.0581 0.07 0.98
## 23 - 21 == 0 2.1400 2.0581 1.04 0.91
## 24 - 21 == 0 0.5600 2.0581 0.27 0.98
## 23 - 22 == 0 1.9900 2.0581 0.97 0.91
## 24 - 22 == 0 0.4100 2.0581 0.20 0.98
## 24 - 23 == 0 -1.5800 2.0581 -0.77 0.97
## (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 12 13 14 15 16 17 18 19 20 21 22
## "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a" "a"
## 23 24
## "a" "a"
## Formatando média seguida de letras.
p0$let <- sprintf("%.2f %s", p0$estimate, l0$mcletters$Letters)
## p0$let <- sprintf("%5.2f %s", p0$estimate, l0$mcletters$Letters)
p0$let
## [1] "5.60 a" "7.76 a" "7.30 a" "3.97 a" "4.78 a" "4.54 a" "6.86 a" "7.01 a"
## [9] "4.93 a" "5.15 a" "3.67 a" "5.02 a" "10.42 a" "5.82 a" "4.03 a" "4.99 a"
## [17] "6.27 a" "6.28 a" "9.17 a" "4.29 a" "4.65 a" "4.80 a" "6.79 a" "5.21 a"
##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- "Produção"
ylab <- "Família de feijão"
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
p0 <- transform(p0, faml=reorder(faml, estimate))
levels(p0$faml)
## [1] "11" "4" "15" "20" "6" "21" "5" "22" "9" "16" "12" "10" "24" "1" "14" "17" "18"
## [18] "23" "7" "8" "3" "2" "19" "13"
p0 <- p0[order(p0$faml),]
xlim <- c(-1, 15)
segplot(faml~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)
})
## A amplitude dos intervalos reflete o número diferente de observações
## para cada nível de família.
##-----------------------------------------------------------------------------
## Experimento com mudas de ?? onde está se avaliando os substratos.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/mudas.txt",
header=TRUE, sep=";")
da$bloc <- factor(da$bloc)
str(da)
## 'data.frame': 1600 obs. of 14 variables:
## $ trat: Factor w/ 16 levels "BcanaTFiltro 3:2",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bloc: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ rep : int 1 2 3 4 5 6 7 8 9 10 ...
## $ viab: int 0 1 1 1 1 1 0 1 1 0 ...
## $ alt : num NA 9.5 17.5 19 21 21 NA 8 15 NA ...
## $ diam: num NA 2.42 4.1 6.09 8.48 5.72 NA 1.68 4.26 NA ...
## $ nfol: int NA 9 8 9 10 13 NA 7 9 NA ...
## $ afol: num NA NA 168 NA NA ...
## $ pff : num NA NA 2.58 NA NA NA NA 0.8 3.22 NA ...
## $ pfc : num NA NA 2.58 NA NA NA NA 0.46 2.66 NA ...
## $ pfr : num NA NA 4.44 NA NA NA NA 1.08 4 NA ...
## $ psf : num NA NA 0.84 NA NA NA NA 0.26 1.13 NA ...
## $ psc : num NA NA 0.82 NA NA NA NA 0.1 0.75 NA ...
## $ psr : num NA NA 1.13 NA NA NA NA 0.2 1.02 NA ...
##-----------------------------------------------------------------------------
## Seleciona só uma variável para não tomar espaço com o print dos
## resultados. Em situação real, não precisa fazer isso. Omitir linhas
## com NA. Cuidado ao usar na.omit() com muitas colunas de resposta,
## pois no caso as linhas serão eliminadas pela união de linhas com NA,
## ou seja, se houver NA para apenas uma variável na linha, a linha toda
## é removida, o que descarta valores úteis das demais.
da <- subset(da, select=c("trat","bloc","rep","psr"))
da <- droplevels(na.omit(da))
str(da)
## 'data.frame': 588 obs. of 4 variables:
## $ trat: Factor w/ 15 levels "BcanaTFiltro 3:2",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ bloc: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 2 2 2 ...
## $ rep : int 3 8 9 13 16 18 19 1 3 4 ...
## $ psr : num 1.13 0.2 1.02 0.46 0.52 0.44 0.27 0.4 0.81 1.02 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:1012] 1 2 4 5 6 7 10 11 12 14 ...
## .. ..- attr(*, "names")= chr [1:1012] "1" "2" "4" "5" ...
## Alto nível de desbalanceamento.
xtabs(~trat+bloc, da)
## bloc
## trat 1 2 3 4 5
## BcanaTFiltro 3:2 6 11 11 12 7
## CpinusCAC 7:3 4 13 10 10 9
## Plantmax 7 14 11 8 11
## TBCEstBov 1:1 0 6 8 8 11
## TBVermCAC 1:1:1 0 2 2 4 0
## TerraBCAC 1:1 2 7 7 7 9
## Terra Plantas 6 10 9 10 9
## Trimix 9 6 9 13 8
## TSAreiaEsterco 2:2:1 5 10 4 10 6
## VermCAC 1:1 3 12 10 12 7
## VermCAC 1:3 2 9 12 13 9
## VermCAC 3:1 3 9 8 9 9
## Vermiculita 5 14 8 10 12
## VermTSCAC 2:1:2 5 9 9 8 6
## VermTSMO 3:1:6 4 8 8 10 4
## Este é um caso em que começou-se com 20 plantas por parcela e algumas
## morreram. Para-se fazer a análise deve-se assumir que as perdas
## foram ao acaso (MAR: missing at random) e não podem de forma alguma
## estarem associadas aos níveis dos fatores pois enviesaria conclusões.
## psr: peso seco de raizes
xyplot(psr~trat|bloc, groups=rep, data=da, layout=c(NA,1),
scales=list(x=list(rot=90, alternating=1)))
##-----------------------------------------------------------------------------
## Fazendo a análise com as médias e os pesos.
db <- ddply(da, .(bloc, trat), summarise,
psrm=mean(psr), n=length(psr))
str(db)
## 'data.frame': 72 obs. of 4 variables:
## $ bloc: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trat: Factor w/ 15 levels "BcanaTFiltro 3:2",..: 1 2 3 6 7 8 9 10 11 12 ...
## $ psrm: num 0.8 0.383 0.577 1.13 0.637 ...
## $ n : int 6 4 7 2 6 9 5 3 2 3 ...
m0 <- lm(psrm~bloc+trat, db, weight=n)
par(mfrow=c(2,2)); plot(m0); layout(1)
anova(m0)
## Analysis of Variance Table
##
## Response: psrm
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 7.70 1.924 9.65 6.1e-06 ***
## trat 14 5.06 0.361 1.81 0.061 .
## Residuals 53 10.57 0.199
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Comparações múltiplas (embora não tenha necessidade, mas é usual
## reportar o valor das médias ajustadas).
X <- LSmatrix(m0, effect="trat")
rownames(X) <- levels(da$trat)
## A função apmc() engloba todos os comandos antes usados para chegar a
## tabela de comparações múltiplas.
## grid <- apmc(X=X, model=m0, focus="trat", test="bonferroni")
grid <- apmc(X=X, model=m0, focus="trat", test="fdr")
grid <- transform(grid, trat=reorder(trat, estimate))
grid <- grid[order(grid$trat),]
##-----------------------------------------------------------------------------
## Gráfico.
xlab <- "Peso seco de raízes"
ylab <- "Substrato"
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
xlim <- c(0.1, 1)
segplot(trat~lwr+upr, centers=estimate,
data=grid, draw=FALSE, xlim=xlim,
xlab=xlab, ylab=ylab, sub=sub,
letras=sprintf("%.2f %s", grid$estimate, grid$cld),
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)
})
## Porque o F e as comparações múltiplas às vezes discordam?