Análise de VCU de arroz
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "lme4", "plyr",
"gridExtra", "aod", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra doBy multcomp lme4 plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## gridExtra aod wzRfun
## TRUE TRUE TRUE
#source("http://dl.dropboxusercontent.com/u/48140237/apc.R")
#source("http://dl.dropboxusercontent.com/u/48140237/equallevels.R")
#source("http://dl.dropboxusercontent.com/u/48140237/desdobrglht.R")
sessionInfo()
## R version 3.1.1 (2014-07-10)
## Platform: x86_64-pc-linux-gnu (64-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] grid methods splines stats graphics grDevices utils datasets
## [9] base
##
## other attached packages:
## [1] wzRfun_0.2 aod_1.3 gridExtra_0.9.1 plyr_1.8.1
## [5] lme4_1.1-6 Rcpp_0.11.2 Matrix_1.1-4 multcomp_1.3-3
## [9] TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-10 MASS_7.3-34
## [13] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [17] knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 minqa_1.2.3 nlme_3.1-117
## [5] RcppEigen_0.3.2.0.2 sandwich_2.3-0 stringr_0.6.2 tools_3.1.1
## [9] zoo_1.7-10
trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Leitura dos dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf2/VCU0910.csv",
header=TRUE, sep=";")
str(da)
## 'data.frame': 1244 obs. of 9 variables:
## $ par : int 101 102 103 104 105 106 107 108 109 110 ...
## $ rep : int 1 1 1 1 1 1 1 1 1 1 ...
## $ bl : int 1 1 1 1 1 1 1 1 1 1 ...
## $ trat : int 16 9 5 13 15 17 14 10 3 12 ...
## $ cult : Factor w/ 17 levels "AB062008","AB062037",..: 16 2 9 6 12 17 7 3 13 5 ...
## $ prod : int 1193 1287 1007 1320 1327 1320 1287 1260 1213 1167 ...
## $ ensaio: Factor w/ 19 levels "EA09VCU006","EA09VCU008",..: 8 8 8 8 8 8 8 8 8 8 ...
## $ estado: Factor w/ 7 levels "GO","MA","MT",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ local : Factor w/ 18 levels "ALTAMIRA","ANAPOLIS",..: 15 15 15 15 15 15 15 15 15 15 ...
##-----------------------------------------------------------------------------
## Editar.
da <- transform(da,
bl=factor(bl),
blin=interaction(bl, ensaio, drop=TRUE))
str(da)
## 'data.frame': 1244 obs. of 10 variables:
## $ par : int 101 102 103 104 105 106 107 108 109 110 ...
## $ rep : int 1 1 1 1 1 1 1 1 1 1 ...
## $ bl : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ trat : int 16 9 5 13 15 17 14 10 3 12 ...
## $ cult : Factor w/ 17 levels "AB062008","AB062037",..: 16 2 9 6 12 17 7 3 13 5 ...
## $ prod : int 1193 1287 1007 1320 1327 1320 1287 1260 1213 1167 ...
## $ ensaio: Factor w/ 19 levels "EA09VCU006","EA09VCU008",..: 8 8 8 8 8 8 8 8 8 8 ...
## $ estado: Factor w/ 7 levels "GO","MA","MT",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ local : Factor w/ 18 levels "ALTAMIRA","ANAPOLIS",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ blin : Factor w/ 76 levels "1.EA09VCU006",..: 29 29 29 29 29 29 29 29 29 29 ...
sum(is.na(da$prod))
## [1] 24
da <- na.omit(da)
str(da)
## 'data.frame': 1220 obs. of 10 variables:
## $ par : int 101 102 103 104 105 106 107 108 109 110 ...
## $ rep : int 1 1 1 1 1 1 1 1 1 1 ...
## $ bl : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ trat : int 16 9 5 13 15 17 14 10 3 12 ...
## $ cult : Factor w/ 17 levels "AB062008","AB062037",..: 16 2 9 6 12 17 7 3 13 5 ...
## $ prod : int 1193 1287 1007 1320 1327 1320 1287 1260 1213 1167 ...
## $ ensaio: Factor w/ 19 levels "EA09VCU006","EA09VCU008",..: 8 8 8 8 8 8 8 8 8 8 ...
## $ estado: Factor w/ 7 levels "GO","MA","MT",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ local : Factor w/ 18 levels "ALTAMIRA","ANAPOLIS",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ blin : Factor w/ 76 levels "1.EA09VCU006",..: 29 29 29 29 29 29 29 29 29 29 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:24] 143 484 490 493 502 517 522 523 528 752 ...
## .. ..- attr(*, "names")= chr [1:24] "143" "484" "490" "493" ...
## Nos níveis dos fatores não pode ter nome com traço, dá problema na
## cld(). Trocar por espaço.
levels(da$cult) <- gsub("-", " ", levels(da$cult))
levels(da$cult)
## [1] "AB062008" "AB062037" "AB062041"
## [4] "AB062045" "AB062048" "AB062104"
## [7] "AB062138" "BRA052023" "BRA052033"
## [10] "BRA052034" "BRA052045" "BRSGO_Serra Dourada"
## [13] "BRSMG_Curinga" "BRS_Primavera" "BRS_Sertaneja"
## [16] "CMG_1152" "H1"
levels(da$ensaio)
## [1] "EA09VCU006" "EA09VCU008" "EA09VCU009" "EA09VCU017" "EA09VCU018" "EA09VCU022"
## [7] "EA09VCU023" "EA09VCU025" "EA09VCU026" "EA09VCU029" "EA09VCU030" "EA09VCU031"
## [13] "EA09VCU036" "EA09VCU037" "EA09VCU038" "EA09VCU039" "EA09VCU040" "EA09VCU050"
## [19] "EA09VCU056"
##-----------------------------------------------------------------------------
## Layout.
x <- xtabs(~ensaio+cult, data=da)
## dimnames(x) <- NULL
mosaicplot(x, off=10, las=4)
x
## cult
## ensaio AB062008 AB062037 AB062041 AB062045 AB062048 AB062104 AB062138 BRA052023
## EA09VCU006 3 4 4 4 3 2 4 4
## EA09VCU008 4 4 4 4 4 4 4 4
## EA09VCU009 4 4 4 4 4 4 4 4
## EA09VCU017 4 4 4 4 4 4 4 4
## EA09VCU018 3 4 4 3 3 3 4 4
## EA09VCU022 4 4 4 4 4 4 4 4
## EA09VCU023 4 4 4 4 4 4 4 4
## EA09VCU025 4 4 4 4 4 4 4 4
## EA09VCU026 4 4 4 4 4 4 4 4
## EA09VCU029 4 4 4 4 4 4 4 4
## EA09VCU030 4 4 4 4 4 4 4 4
## EA09VCU031 4 4 4 4 4 4 4 4
## EA09VCU036 3 4 4 3 3 4 4 4
## EA09VCU037 4 4 4 4 4 4 4 4
## EA09VCU038 4 4 4 4 4 4 4 4
## EA09VCU039 4 4 4 4 4 4 4 4
## EA09VCU040 4 4 4 4 4 4 4 4
## EA09VCU050 4 4 4 4 4 4 4 4
## EA09VCU056 4 4 4 4 4 4 4 4
## cult
## ensaio BRA052033 BRA052034 BRA052045 BRSGO_Serra Dourada BRSMG_Curinga
## EA09VCU006 4 4 4 4 3
## EA09VCU008 4 4 4 4 4
## EA09VCU009 4 4 4 4 4
## EA09VCU017 4 4 4 4 4
## EA09VCU018 4 4 4 3 3
## EA09VCU022 3 4 4 4 4
## EA09VCU023 4 4 4 4 4
## EA09VCU025 4 4 4 4 4
## EA09VCU026 4 4 4 4 4
## EA09VCU029 4 4 4 4 4
## EA09VCU030 4 4 4 4 4
## EA09VCU031 4 4 4 4 4
## EA09VCU036 4 4 4 4 4
## EA09VCU037 4 4 4 4 4
## EA09VCU038 4 4 4 4 4
## EA09VCU039 4 4 4 4 4
## EA09VCU040 4 4 4 4 4
## EA09VCU050 4 4 4 4 4
## EA09VCU056 3 4 4 4 4
## cult
## ensaio BRS_Primavera BRS_Sertaneja CMG_1152 H1
## EA09VCU006 4 4 4 0
## EA09VCU008 4 4 4 0
## EA09VCU009 4 4 4 0
## EA09VCU017 4 4 3 0
## EA09VCU018 4 4 4 0
## EA09VCU022 4 4 4 3
## EA09VCU023 4 4 4 4
## EA09VCU025 4 4 4 4
## EA09VCU026 4 4 4 4
## EA09VCU029 4 4 4 4
## EA09VCU030 4 4 4 4
## EA09VCU031 4 4 4 4
## EA09VCU036 0 4 3 0
## EA09VCU037 4 4 4 0
## EA09VCU038 4 4 4 0
## EA09VCU039 4 4 4 0
## EA09VCU040 4 4 3 0
## EA09VCU050 4 4 4 0
## EA09VCU056 4 4 4 0
##-----------------------------------------------------------------------------
## Ver.
xyplot(prod~cult|ensaio, data=da, as.table=TRUE)
## xyplot(sqrt(prod)~cult|ensaio, data=da, as.table=TRUE)
##-----------------------------------------------------------------------------
## Ajuste do modelo que considera apenas o efeito de cultivar como
## fixo, os demais e interações são aleatórios.
## * fixo: cult;
## * aleatório: ensaio, bl dentro de ensaio e ensaio interação com cult.
## * Será análisado a raíz da produtividade porque nessa escala os
## pressupostos são melhor atendidos.
da$y <- sqrt(da$prod)
m0 <- lmer(y~cult+(1|ensaio)+(1|ensaio:bl)+(1|ensaio:cult),
data=da, REML=FALSE)
summary(m0)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: y ~ cult + (1 | ensaio) + (1 | ensaio:bl) + (1 | ensaio:cult)
## Data: da
##
## AIC BIC logLik deviance df.resid
## 8132 8239 -4045 8090 1199
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.204 -0.542 0.000 0.555 3.569
##
## Random effects:
## Groups Name Variance Std.Dev.
## ensaio:cult (Intercept) 21.52 4.64
## ensaio:bl (Intercept) 2.25 1.50
## ensaio (Intercept) 86.87 9.32
## Residual 28.28 5.32
## Number of obs: 1220, groups: ensaio:cult, 310; ensaio:bl, 76; ensaio, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 60.3910 2.4749 24.40
## cultAB062037 -2.7927 1.7402 -1.60
## cultAB062041 0.9710 1.7402 0.56
## cultAB062045 -4.9676 1.7437 -2.85
## cultAB062048 -5.8911 1.7456 -3.37
## cultAB062104 -2.7027 1.7469 -1.55
## cultAB062138 -4.3011 1.7402 -2.47
## cultBRA052023 -4.4031 1.7402 -2.53
## cultBRA052033 -0.9462 1.7438 -0.54
## cultBRA052034 0.0224 1.7402 0.01
## cultBRA052045 -4.4376 1.7402 -2.55
## cultBRSGO_Serra Dourada -3.4740 1.7420 -1.99
## cultBRSMG_Curinga -8.7869 1.7435 -5.04
## cultBRS_Primavera -6.1942 1.7654 -3.51
## cultBRS_Sertaneja -4.3709 1.7402 -2.51
## cultCMG_1152 -8.6274 1.7457 -4.94
## cultH1 -4.5380 2.4105 -1.88
##
## Correlation of Fixed Effects:
## (Intr) cAB06203 cAB062041 cAB062045 cAB062048 cAB06210 cAB06213 cBRA05202
## cltAB062037 -0.354
## cltAB062041 -0.354 0.503
## cltAB062045 -0.353 0.502 0.502
## cltAB062048 -0.353 0.502 0.502 0.501
## cltAB062104 -0.352 0.501 0.501 0.500 0.500
## cltAB062138 -0.354 0.503 0.503 0.502 0.502 0.501
## clBRA052023 -0.354 0.503 0.503 0.502 0.502 0.501 0.503
## clBRA052033 -0.353 0.502 0.502 0.501 0.501 0.500 0.502 0.502
## clBRA052034 -0.354 0.503 0.503 0.502 0.502 0.501 0.503 0.503
## clBRA052045 -0.354 0.503 0.503 0.502 0.502 0.501 0.503 0.503
## cltBRSGO_SD -0.353 0.503 0.503 0.502 0.501 0.501 0.503 0.503
## cltBRSMG_Cr -0.353 0.502 0.502 0.501 0.500 0.500 0.502 0.502
## cltBRS_Prmv -0.349 0.496 0.496 0.495 0.494 0.494 0.496 0.496
## cltBRS_Srtn -0.354 0.503 0.503 0.502 0.502 0.501 0.503 0.503
## cltCMG_1152 -0.353 0.502 0.502 0.500 0.500 0.500 0.502 0.502
## cultH1 -0.255 0.363 0.363 0.362 0.362 0.362 0.363 0.363
## cBRA052033 cBRA052034 cBRA05204 cBRSGD cBRSMG cBRS_P cBRS_S cCMG_1
## cltAB062037
## cltAB062041
## cltAB062045
## cltAB062048
## cltAB062104
## cltAB062138
## clBRA052023
## clBRA052033
## clBRA052034 0.502
## clBRA052045 0.502 0.503
## cltBRSGO_SD 0.502 0.503 0.503
## cltBRSMG_Cr 0.501 0.502 0.502 0.501
## cltBRS_Prmv 0.495 0.496 0.496 0.495 0.495
## cltBRS_Srtn 0.502 0.503 0.503 0.503 0.502 0.496
## cltCMG_1152 0.501 0.502 0.502 0.501 0.500 0.494 0.502
## cultH1 0.362 0.363 0.363 0.363 0.362 0.359 0.363 0.362
## Abandona a interação.
m1 <- update(m0, .~.-(1|ensaio:cult), data=da)
anova(m0, m1)
## Data: da
## Models:
## m1: y ~ cult + (1 | ensaio) + (1 | ensaio:bl)
## m0: y ~ cult + (1 | ensaio) + (1 | ensaio:bl) + (1 | ensaio:cult)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 20 8375 8477 -4167 8335
## m0 21 8132 8239 -4045 8090 245 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Existe interação ensaio:cult.
## Predição dos efeitos aleatórios.
lapply(ranef(m0), c)
## $`ensaio:cult`
## $`ensaio:cult`$`(Intercept)`
## [1] -4.81859 -0.27096 1.98192 5.38821 -5.19403 1.61502 -3.37540 3.94883
## [9] 4.18480 -2.29390 2.88689 -0.63400 -1.93823 -2.47383 6.75154 -4.68854
## [17] 0.86237 0.97966 -0.11547 0.30360 0.77902 0.13170 2.45682 -6.63295
## [25] 4.00586 0.61846 -4.42697 -2.87491 -1.72847 2.31666 3.47303 2.73059
## [33] -4.81966 2.78984 0.95819 -0.09177 1.29147 -2.89770 0.35472 5.24165
## [41] 1.67010 -0.59983 0.66762 -4.43800 2.06264 0.19727 -1.35551 1.17418
## [49] -1.43011 1.04700 -3.20342 -0.30287 0.56930 2.25830 -5.59013 -2.64498
## [57] 1.08031 2.64248 2.66101 -3.64234 3.21152 5.11934 -0.70051 -0.70744
## [65] 1.66575 3.44527 -1.86397 -5.82604 -2.22231 -0.83852 -1.83695 -1.21456
## [73] 2.34173 -2.93457 4.30108 -0.85744 3.30099 6.29461 -1.24005 -5.06580
## [81] 3.34359 1.88810 1.74742 2.67323 3.00407 1.92374 -8.07050 -6.32087
## [89] -1.26288 2.91123 -0.61189 3.49278 -8.48065 6.20034 0.26076 9.39534
## [97] -10.75283 0.58836 2.74771 -3.20979 2.27157 2.30492 -0.25785 0.07463
## [105] 6.07015 -7.90748 2.23157 2.59991 0.48453 -3.84321 -6.05739 7.00991
## [113] -3.23794 0.02838 -1.72012 -1.43281 -4.40067 -1.89248 -1.19365 -2.25991
## [121] -1.17625 -2.30334 -2.48334 1.45476 0.54478 -0.99995 4.11722 0.15665
## [129] 1.92024 4.09915 2.52709 0.85527 -3.16430 1.90197 2.64858 -4.61689
## [137] -6.88558 13.29573 0.33409 0.33396 5.03602 -6.14835 0.25196 -5.17811
## [145] -0.04193 -3.64725 -1.95385 9.23877 0.06139 3.10066 0.13661 -0.06022
## [153] 2.29243 -3.15407 2.56207 -1.80946 -3.50562 -4.19856 0.74123 -0.33622
## [161] -0.57205 6.36900 -2.24373 -0.78727 3.68269 0.19775 -0.31370 2.47426
## [169] 2.07426 3.18364 1.28146 -0.69831 0.78709 -2.38608 -5.30012 -0.11018
## [177] -1.71046 -2.25318 -0.88865 2.12257 -1.48256 2.10977 -6.37378 4.03589
## [185] 1.92111 3.26072 0.68570 5.35984 -5.21568 -3.37169 3.21174 -0.30206
## [193] 2.48500 0.96640 5.13148 -0.68116 -3.55814 0.29236 -6.83386 3.96870
## [201] -0.88793 1.50169 -7.10415 -5.83591 -1.08633 1.04237 8.94752 1.92426
## [209] 1.48197 -3.69372 3.07148 2.17886 -7.33850 0.31260 5.50026 4.24638
## [217] 4.80209 -0.12515 1.98326 -1.32606 2.91053 -0.41463 4.30112 -0.34461
## [225] -1.20697 4.65820 1.57812 -15.35474 1.20405 -10.40959 8.68346 -1.58484
## [233] 4.21149 -9.03063 3.47425 0.32169 -2.14663 0.25097 -4.05587 2.95149
## [241] -11.13702 7.85321 -0.63913 -3.74364 -2.61566 5.42754 2.11950 -6.95640
## [249] -6.97926 -1.53647 -6.78594 -4.08226 7.55832 2.97906 5.36991 2.00565
## [257] 10.53526 -2.71909 6.42736 -4.96990 -2.82673 1.69795 -1.25100 -1.00577
## [265] -2.68915 -2.10799 -0.79120 -1.18733 -3.07100 -0.70909 -0.52350 -2.67091
## [273] -0.07579 -1.97587 -0.15410 7.28004 6.95356 2.32278 -4.74683 -3.47991
## [281] -0.94631 3.45714 4.53313 2.97413 0.08927 -5.12064 -4.34695 -1.97361
## [289] 2.26570 3.65602 -4.53731 3.16644 2.25666 1.10752 -2.68629 -5.18388
## [297] 1.77127 6.00044 2.53874 8.10975 0.83639 1.98285 -1.95206 -0.71545
## [305] -2.27759 -4.24629 1.31626 -2.88910 -6.42625 -0.22703
##
##
## $`ensaio:bl`
## $`ensaio:bl`$`(Intercept)`
## [1] -1.149056 -0.805356 0.896690 1.169609 -0.003725 -0.895024 -0.031754 1.231631
## [9] 1.011012 0.675106 -0.259516 -1.195950 -1.111431 0.519822 1.063908 -0.433865
## [17] 1.006312 0.576914 -1.799879 -0.050145 -1.697360 2.087302 0.588833 -0.838515
## [25] -0.170986 -0.694254 0.114138 0.949620 -0.308503 1.407375 -0.100470 -1.525833
## [33] 0.415088 -0.131149 -0.587647 0.540100 -0.107367 0.377823 -0.165001 0.132902
## [41] -2.249663 1.110515 1.230165 -0.186453 -2.184558 -1.110156 1.048899 2.351858
## [49] -0.479352 -0.919159 -0.379658 1.619489 -0.544102 0.607739 0.822492 -0.676703
## [57] -1.006701 0.122182 0.348435 0.349976 -0.708541 0.019244 0.064723 0.816709
## [65] 0.440224 -0.151110 -0.964318 0.501963 1.896321 -0.281805 0.205619 -1.992251
## [73] 0.050965 -1.284284 -0.033626 0.843522
##
##
## $ensaio
## $ensaio$`(Intercept)`
## [1] 4.319 11.623 8.903 1.484 -10.298 5.414 7.663 -20.359 9.125 9.201
## [11] -3.684 4.093 -6.125 8.084 -7.184 7.416 -6.687 -6.644 -16.344
##-----------------------------------------------------------------------------
## Teste para os termos de efeito fixo.
anova(m0)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## cult 16 2299 144 5.08
## Não tem p-valor! É de propósito. Douglas Bates não concorda que o
## procedimento adequado para ser avaliar a estatística F seja por meio
## de ajustes no número de graus de liberdade do denominador. Para mais
## detalhes leia:
## https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
## Alternativas (em order de recomendação):
## 1. Teste de razão de verossimilhanças entre modelos aninhados (um com
## outro sem o termo fixo de interesse, usar REML=FALSE).
## 2. Teste de Wald (inferência aproximada, pacote aod).
V <- vcov(m0)
b <- fixef(m0)
nobars(formula(m0))
## y ~ cult
Terms <- which(attr(m0@pp$X, "assign")==1)
## Chi-quadrado de Wald (baseado na aproximação quadrática da
## verossimilhança).
wt <- wald.test(Sigma=V, b=b, Terms=Terms)
wt
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 81.3, df = 16, P(> X2) = 9.7e-11
## Chi-quadrado da razão de verossimilhanças (não baseado em aproximação
## de função, melhor opção).
m1 <- update(m0, .~.-cult)
anova(m0, m1)
## Data: da
## Models:
## m1: y ~ (1 | ensaio) + (1 | ensaio:bl) + (1 | ensaio:cult)
## m0: y ~ cult + (1 | ensaio) + (1 | ensaio:bl) + (1 | ensaio:cult)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 5 8172 8197 -4081 8162
## m0 21 8132 8239 -4045 8090 71.7 16 5e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
bec <- unlist(ranef(m0)[1])
beb <- unlist(ranef(m0)[2])
be <- unlist(ranef(m0)[3])
xyplot(r~f, type=c("p", "smooth"))
xyplot(sqrt(abs(r))~f, type=c("p", "smooth"))
grid.arrange(qqmath(r),
qqmath(bec),
qqmath(beb),
qqmath(be),
nrow=2)
qqmath(~r|da$ensaio)
##-----------------------------------------------------------------------------
## Médias ajustadas.
## Formula só da parte de efeito fixo.
f <- nobars(formula(m0))
X <- LSmatrix(lm(f, da), effect=c("cult"))
rownames(X) <- levels(da$cult)
## Estimativas das médias.
summary(glht(m0, linfct=X), test=adjusted(type="none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = y ~ cult + (1 | ensaio) + (1 | ensaio:bl) + (1 |
## ensaio:cult), data = da, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## AB062008 == 0 60.39 2.47 24.4 <2e-16 ***
## AB062037 == 0 57.60 2.47 23.3 <2e-16 ***
## AB062041 == 0 61.36 2.47 24.8 <2e-16 ***
## AB062045 == 0 55.42 2.47 22.4 <2e-16 ***
## AB062048 == 0 54.50 2.47 22.0 <2e-16 ***
## AB062104 == 0 57.69 2.48 23.3 <2e-16 ***
## AB062138 == 0 56.09 2.47 22.7 <2e-16 ***
## BRA052023 == 0 55.99 2.47 22.7 <2e-16 ***
## BRA052033 == 0 59.44 2.47 24.0 <2e-16 ***
## BRA052034 == 0 60.41 2.47 24.4 <2e-16 ***
## BRA052045 == 0 55.95 2.47 22.6 <2e-16 ***
## BRSGO_Serra Dourada == 0 56.92 2.47 23.0 <2e-16 ***
## BRSMG_Curinga == 0 51.60 2.47 20.9 <2e-16 ***
## BRS_Primavera == 0 54.20 2.49 21.8 <2e-16 ***
## BRS_Sertaneja == 0 56.02 2.47 22.7 <2e-16 ***
## CMG_1152 == 0 51.76 2.47 20.9 <2e-16 ***
## H1 == 0 55.85 2.98 18.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
## Contraste.
Xc <- apc(X)
dim(Xc)
## [1] 136 17
g <- summary(glht(m0, linfct=Xc), test=adjusted(type="fdr")); g
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = y ~ cult + (1 | ensaio) + (1 | ensaio:bl) + (1 |
## ensaio:cult), data = da, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## AB062008-AB062037 == 0 2.7927 1.7402 1.60 0.21391
## AB062008-AB062041 == 0 -0.9710 1.7402 -0.56 0.72413
## AB062008-AB062045 == 0 4.9676 1.7437 2.85 0.02131 *
## AB062008-AB062048 == 0 5.8911 1.7456 3.37 0.00558 **
## AB062008-AB062104 == 0 2.7027 1.7469 1.55 0.22848
## AB062008-AB062138 == 0 4.3011 1.7402 2.47 0.04254 *
## AB062008-BRA052023 == 0 4.4031 1.7402 2.53 0.04179 *
## AB062008-BRA052033 == 0 0.9462 1.7438 0.54 0.72413
## AB062008-BRA052034 == 0 -0.0224 1.7402 -0.01 0.98974
## AB062008-BRA052045 == 0 4.4376 1.7402 2.55 0.04179 *
## AB062008-BRSGO_Serra Dourada == 0 3.4740 1.7420 1.99 0.11351
## AB062008-BRSMG_Curinga == 0 8.7869 1.7435 5.04 1.6e-05 ***
## AB062008-BRS_Primavera == 0 6.1942 1.7654 3.51 0.00510 **
## AB062008-BRS_Sertaneja == 0 4.3709 1.7402 2.51 0.04181 *
## AB062008-CMG_1152 == 0 8.6274 1.7457 4.94 1.8e-05 ***
## AB062008-H1 == 0 4.5380 2.4105 1.88 0.13107
## AB062037-AB062041 == 0 -3.7637 1.7347 -2.17 0.08168 .
## AB062037-AB062045 == 0 2.1749 1.7384 1.25 0.35851
## AB062037-AB062048 == 0 3.0984 1.7402 1.78 0.15937
## AB062037-AB062104 == 0 -0.0901 1.7413 -0.05 0.98974
## AB062037-AB062138 == 0 1.5084 1.7347 0.87 0.53918
## AB062037-BRA052023 == 0 1.6104 1.7347 0.93 0.51959
## AB062037-BRA052033 == 0 -1.8465 1.7383 -1.06 0.47212
## AB062037-BRA052034 == 0 -2.8151 1.7347 -1.62 0.20926
## AB062037-BRA052045 == 0 1.6449 1.7347 0.95 0.51262
## AB062037-BRSGO_Serra Dourada == 0 0.6813 1.7365 0.39 0.80764
## AB062037-BRSMG_Curinga == 0 5.9941 1.7383 3.45 0.00548 **
## AB062037-BRS_Primavera == 0 3.4015 1.7602 1.93 0.12153
## AB062037-BRS_Sertaneja == 0 1.5782 1.7347 0.91 0.51959
## AB062037-CMG_1152 == 0 5.8347 1.7402 3.35 0.00572 **
## AB062037-H1 == 0 1.7453 2.4071 0.73 0.62454
## AB062041-AB062045 == 0 5.9386 1.7384 3.42 0.00555 **
## AB062041-AB062048 == 0 6.8621 1.7402 3.94 0.00109 **
## AB062041-AB062104 == 0 3.6737 1.7413 2.11 0.09302 .
## AB062041-AB062138 == 0 5.2721 1.7347 3.04 0.01344 *
## AB062041-BRA052023 == 0 5.3742 1.7347 3.10 0.01262 *
## AB062041-BRA052033 == 0 1.9172 1.7383 1.10 0.45344
## AB062041-BRA052034 == 0 0.9486 1.7347 0.55 0.72413
## AB062041-BRA052045 == 0 5.4087 1.7347 3.12 0.01238 *
## AB062041-BRSGO_Serra Dourada == 0 4.4450 1.7365 2.56 0.04179 *
## AB062041-BRSMG_Curinga == 0 9.7579 1.7383 5.61 2.4e-06 ***
## AB062041-BRS_Primavera == 0 7.1652 1.7602 4.07 0.00071 ***
## AB062041-BRS_Sertaneja == 0 5.3419 1.7347 3.08 0.01282 *
## AB062041-CMG_1152 == 0 9.5985 1.7402 5.52 2.4e-06 ***
## AB062041-H1 == 0 5.5091 2.4071 2.29 0.06261 .
## AB062045-AB062048 == 0 0.9235 1.7435 0.53 0.72413
## AB062045-AB062104 == 0 -2.2649 1.7450 -1.30 0.33448
## AB062045-AB062138 == 0 -0.6665 1.7384 -0.38 0.80840
## AB062045-BRA052023 == 0 -0.5645 1.7384 -0.32 0.84478
## AB062045-BRA052033 == 0 -4.0214 1.7420 -2.31 0.06068 .
## AB062045-BRA052034 == 0 -4.9900 1.7384 -2.87 0.02064 *
## AB062045-BRA052045 == 0 -0.5300 1.7384 -0.30 0.85473
## AB062045-BRSGO_Serra Dourada == 0 -1.4936 1.7400 -0.86 0.53922
## AB062045-BRSMG_Curinga == 0 3.8192 1.7420 2.19 0.07868 .
## AB062045-BRS_Primavera == 0 1.2266 1.7636 0.70 0.64268
## AB062045-BRS_Sertaneja == 0 -0.5967 1.7384 -0.34 0.83588
## AB062045-CMG_1152 == 0 3.6598 1.7438 2.10 0.09374 .
## AB062045-H1 == 0 -0.4296 2.4093 -0.18 0.95497
## AB062048-AB062104 == 0 -3.1884 1.7469 -1.83 0.14672
## AB062048-AB062138 == 0 -1.5900 1.7402 -0.91 0.51959
## AB062048-BRA052023 == 0 -1.4880 1.7402 -0.86 0.53922
## AB062048-BRA052033 == 0 -4.9449 1.7438 -2.84 0.02145 *
## AB062048-BRA052034 == 0 -5.9135 1.7402 -3.40 0.00555 **
## AB062048-BRA052045 == 0 -1.4535 1.7402 -0.84 0.54888
## AB062048-BRSGO_Serra Dourada == 0 -2.4171 1.7419 -1.39 0.29186
## AB062048-BRSMG_Curinga == 0 2.8958 1.7439 1.66 0.19650
## AB062048-BRS_Primavera == 0 0.3031 1.7654 0.17 0.95497
## AB062048-BRS_Sertaneja == 0 -1.5202 1.7402 -0.87 0.53918
## AB062048-CMG_1152 == 0 2.7363 1.7457 1.57 0.22525
## AB062048-H1 == 0 -1.3531 2.4105 -0.56 0.72413
## AB062104-AB062138 == 0 1.5984 1.7413 0.92 0.51959
## AB062104-BRA052023 == 0 1.7005 1.7413 0.98 0.50242
## AB062104-BRA052033 == 0 -1.7565 1.7449 -1.01 0.49312
## AB062104-BRA052034 == 0 -2.7251 1.7413 -1.56 0.22525
## AB062104-BRA052045 == 0 1.7350 1.7413 1.00 0.49312
## AB062104-BRSGO_Serra Dourada == 0 0.7714 1.7431 0.44 0.77217
## AB062104-BRSMG_Curinga == 0 6.0842 1.7450 3.49 0.00512 **
## AB062104-BRS_Primavera == 0 3.4915 1.7668 1.98 0.11449
## AB062104-BRS_Sertaneja == 0 1.6682 1.7413 0.96 0.51083
## AB062104-CMG_1152 == 0 5.9248 1.7468 3.39 0.00555 **
## AB062104-H1 == 0 1.8354 2.4113 0.76 0.60131
## AB062138-BRA052023 == 0 0.1021 1.7347 0.06 0.98974
## AB062138-BRA052033 == 0 -3.3549 1.7383 -1.93 0.12153
## AB062138-BRA052034 == 0 -4.3235 1.7347 -2.49 0.04181 *
## AB062138-BRA052045 == 0 0.1366 1.7347 0.08 0.98974
## AB062138-BRSGO_Serra Dourada == 0 -0.8271 1.7365 -0.48 0.75621
## AB062138-BRSMG_Curinga == 0 4.4858 1.7383 2.58 0.04179 *
## AB062138-BRS_Primavera == 0 1.8931 1.7602 1.08 0.46795
## AB062138-BRS_Sertaneja == 0 0.0698 1.7347 0.04 0.98974
## AB062138-CMG_1152 == 0 4.3264 1.7402 2.49 0.04181 *
## AB062138-H1 == 0 0.2369 2.4071 0.10 0.98974
## BRA052023-BRA052033 == 0 -3.4569 1.7383 -1.99 0.11351
## BRA052023-BRA052034 == 0 -4.4255 1.7347 -2.55 0.04179 *
## BRA052023-BRA052045 == 0 0.0345 1.7347 0.02 0.98974
## BRA052023-BRSGO_Serra Dourada == 0 -0.9291 1.7365 -0.54 0.72413
## BRA052023-BRSMG_Curinga == 0 4.3837 1.7383 2.52 0.04179 *
## BRA052023-BRS_Primavera == 0 1.7911 1.7602 1.02 0.49312
## BRA052023-BRS_Sertaneja == 0 -0.0323 1.7347 -0.02 0.98974
## BRA052023-CMG_1152 == 0 4.2243 1.7402 2.43 0.04595 *
## BRA052023-H1 == 0 0.1349 2.4071 0.06 0.98974
## BRA052033-BRA052034 == 0 -0.9686 1.7383 -0.56 0.72413
## BRA052033-BRA052045 == 0 3.4914 1.7383 2.01 0.11230
## BRA052033-BRSGO_Serra Dourada == 0 2.5278 1.7401 1.45 0.26184
## BRA052033-BRSMG_Curinga == 0 7.8407 1.7420 4.50 0.00013 ***
## BRA052033-BRS_Primavera == 0 5.2480 1.7638 2.98 0.01592 *
## BRA052033-BRS_Sertaneja == 0 3.4247 1.7383 1.97 0.11449
## BRA052033-CMG_1152 == 0 7.6812 1.7438 4.40 0.00018 ***
## BRA052033-H1 == 0 3.5918 2.4098 1.49 0.25013
## BRA052034-BRA052045 == 0 4.4600 1.7347 2.57 0.04179 *
## BRA052034-BRSGO_Serra Dourada == 0 3.4964 1.7365 2.01 0.11230
## BRA052034-BRSMG_Curinga == 0 8.8092 1.7383 5.07 1.6e-05 ***
## BRA052034-BRS_Primavera == 0 6.2166 1.7602 3.53 0.00510 **
## BRA052034-BRS_Sertaneja == 0 4.3933 1.7347 2.53 0.04179 *
## BRA052034-CMG_1152 == 0 8.6498 1.7402 4.97 1.8e-05 ***
## BRA052034-H1 == 0 4.5604 2.4071 1.89 0.12964
## BRA052045-BRSGO_Serra Dourada == 0 -0.9636 1.7365 -0.55 0.72413
## BRA052045-BRSMG_Curinga == 0 4.3492 1.7383 2.50 0.04181 *
## BRA052045-BRS_Primavera == 0 1.7566 1.7602 1.00 0.49312
## BRA052045-BRS_Sertaneja == 0 -0.0668 1.7347 -0.04 0.98974
## BRA052045-CMG_1152 == 0 4.1898 1.7402 2.41 0.04746 *
## BRA052045-H1 == 0 0.1004 2.4071 0.04 0.98974
## BRSGO_Serra Dourada-BRSMG_Curinga == 0 5.3128 1.7402 3.05 0.01340 *
## BRSGO_Serra Dourada-BRS_Primavera == 0 2.7202 1.7620 1.54 0.22848
## BRSGO_Serra Dourada-BRS_Sertaneja == 0 0.8969 1.7365 0.52 0.72877
## BRSGO_Serra Dourada-CMG_1152 == 0 5.1534 1.7420 2.96 0.01618 *
## BRSGO_Serra Dourada-H1 == 0 1.0640 2.4082 0.44 0.77217
## BRSMG_Curinga-BRS_Primavera == 0 -2.5927 1.7638 -1.47 0.25674
## BRSMG_Curinga-BRS_Sertaneja == 0 -4.4160 1.7383 -2.54 0.04179 *
## BRSMG_Curinga-CMG_1152 == 0 -0.1594 1.7438 -0.09 0.98974
## BRSMG_Curinga-H1 == 0 -4.2488 2.4094 -1.76 0.16283
## BRS_Primavera-BRS_Sertaneja == 0 -1.8233 1.7602 -1.04 0.48616
## BRS_Primavera-CMG_1152 == 0 2.4332 1.7654 1.38 0.29312
## BRS_Primavera-H1 == 0 -1.6562 2.4232 -0.68 0.64641
## BRS_Sertaneja-CMG_1152 == 0 4.2566 1.7402 2.45 0.04464 *
## BRS_Sertaneja-H1 == 0 0.1671 2.4071 0.07 0.98974
## CMG_1152-H1 == 0 -4.0894 2.4105 -1.70 0.18503
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
## Estimativas das médias com comparações.
## grid <- apmc(X=X, model=m0, focus="cult", test="bonferroni")
grid <- apmc(X=X, model=m0, focus="cult", test="fdr")
grid
## cult estimate lwr upr cld
## 1 AB062008 60.39 55.54 65.24 ab
## 2 AB062037 57.60 52.76 62.44 acd
## 3 AB062041 61.36 56.52 66.21 a
## 4 AB062045 55.42 50.58 60.27 ce
## 5 AB062048 54.50 49.65 59.35 de
## 6 AB062104 57.69 52.84 62.54 acd
## 7 AB062138 56.09 51.25 60.93 cd
## 8 BRA052023 55.99 51.14 60.83 cd
## 9 BRA052033 59.44 54.60 64.29 ac
## 10 BRA052034 60.41 55.57 65.26 ab
## 11 BRA052045 55.95 51.11 60.80 cd
## 12 BRSGO_Serra Dourada 56.92 52.07 61.76 bcd
## 13 BRSMG_Curinga 51.60 46.76 56.45 e
## 14 BRS_Primavera 54.20 49.32 59.08 de
## 15 BRS_Sertaneja 56.02 51.18 60.86 cd
## 16 CMG_1152 51.76 46.91 56.61 e
## 17 H1 55.85 50.01 61.70 ace
##-----------------------------------------------------------------------------
## Gráfico.
segplot(reorder(cult, estimate)~lwr+upr, centers=estimate,
ylab="Cultivar de arroz", xlab="raíz da Produtividade",
data=grid, draw=FALSE, cld=grid$cld,
panel=function(x, y, z, centers, subscripts, cld, ...){
panel.segplot(x, y, z, centers=centers,
subscripts=subscripts, ...)
panel.text(x=centers[subscripts], y=as.numeric(z)[subscripts],
labels=cld[subscripts], pos=4)
})
##-----------------------------------------------------------------------------
## Para desdobrar a interação ensaio:cult, tratar como efeito fixo.
formula(m0)
## y ~ cult + (1 | ensaio) + (1 | ensaio:bl) + (1 | ensaio:cult)
m0 <- lmer(y~ensaio*cult+(1|ensaio:bl),
data=da, REML=FALSE)
## fixed-effect model matrix is rank deficient so dropping 13 columns / coefficients
## Aviso devido cultivares que não ocorrem em alguns ensaios. Com isso a
## estrutura não é completamente cruzada.
length(fixef(m0))
## [1] 310
## Combinações previstas se fosse completamente cruzado.
with(da, nlevels(interaction(ensaio, cult, drop=TRUE)))
## [1] 310
## Combinações que de fato existem.
with(da, nlevels(interaction(ensaio, cult, drop=FALSE)))
## [1] 323
## Cria variável combinando níveis de ensaio com cult.
da$enscult <- with(da, interaction(ensaio, cult, drop=TRUE))
## O mesmo modelo escrito com uma variável combinada.
m0 <- lmer(y~enscult+(1|ensaio:bl),
data=da, REML=FALSE)
head(fixef(m0))
## (Intercept) enscultEA09VCU008.AB062008 enscultEA09VCU009.AB062008
## 57.773 15.462 5.175
## enscultEA09VCU017.AB062008 enscultEA09VCU018.AB062008 enscultEA09VCU022.AB062008
## 2.211 -5.306 12.509
##-----------------------------------------------------------------------------
## Fazer o desdobramento dentro do primeiro ensaio. Os demais seguem a
## mesma regra.
X <- LSmatrix(lm(nobars(formula(m0)), data=da), effect="enscult")
rownames(X) <- levels(da$enscult)
dim(X)
## [1] 310 310
## O nome do primeiro ensaio.
split <- levels(da$ensaio)[1]
split
## [1] "EA09VCU006"
## Seleciona as linhas correspondentes à esse ensaio.
Xs <- X[grep(split, rownames(X)),]
dim(Xs)
## [1] 16 310
summary(glht(m0, linfct=Xs), test=adjusted(type="fdr"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = y ~ enscult + (1 | ensaio:bl), data = da, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## EA09VCU006.AB062008 == 0 57.77 2.74 21.1 <2e-16 ***
## EA09VCU006.AB062037 == 0 61.59 2.38 25.8 <2e-16 ***
## EA09VCU006.AB062041 == 0 68.34 2.38 28.7 <2e-16 ***
## EA09VCU006.AB062045 == 0 66.93 2.38 28.1 <2e-16 ***
## EA09VCU006.AB062048 == 0 51.36 2.74 18.8 <2e-16 ***
## EA09VCU006.AB062104 == 0 64.78 3.34 19.4 <2e-16 ***
## EA09VCU006.AB062138 == 0 55.95 2.38 23.5 <2e-16 ***
## EA09VCU006.BRA052023 == 0 65.58 2.38 27.5 <2e-16 ***
## EA09VCU006.BRA052033 == 0 69.35 2.38 29.1 <2e-16 ***
## EA09VCU006.BRA052034 == 0 61.71 2.38 25.9 <2e-16 ***
## EA09VCU006.BRA052045 == 0 64.14 2.38 26.9 <2e-16 ***
## EA09VCU006.BRSGO_Serra Dourada == 0 60.42 2.38 25.4 <2e-16 ***
## EA09VCU006.BRSMG_Curinga == 0 53.13 2.74 19.4 <2e-16 ***
## EA09VCU006.BRS_Primavera == 0 55.26 2.38 23.2 <2e-16 ***
## EA09VCU006.BRS_Sertaneja == 0 69.34 2.38 29.1 <2e-16 ***
## EA09VCU006.CMG_1152 == 0 49.88 2.38 20.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
confint(glht(m0, linfct=Xs), calpha=univariate_calpha())
##
## Simultaneous Confidence Intervals
##
## Fit: lmer(formula = y ~ enscult + (1 | ensaio:bl), data = da, REML = FALSE)
##
## Quantile = 1.96
## 95% confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## EA09VCU006.AB062008 == 0 57.773 52.402 63.144
## EA09VCU006.AB062037 == 0 61.585 56.913 66.258
## EA09VCU006.AB062041 == 0 68.342 63.670 73.014
## EA09VCU006.AB062045 == 0 66.929 62.257 71.601
## EA09VCU006.AB062048 == 0 51.364 45.994 56.734
## EA09VCU006.AB062104 == 0 64.783 58.236 71.330
## EA09VCU006.AB062138 == 0 55.952 51.280 60.625
## EA09VCU006.BRA052023 == 0 65.581 60.909 70.254
## EA09VCU006.BRA052033 == 0 69.352 64.679 74.024
## EA09VCU006.BRA052034 == 0 61.713 57.040 66.385
## EA09VCU006.BRA052045 == 0 64.136 59.463 68.808
## EA09VCU006.BRSGO_Serra Dourada == 0 60.421 55.749 65.094
## EA09VCU006.BRSMG_Curinga == 0 53.129 47.758 58.500
## EA09VCU006.BRS_Primavera == 0 55.257 50.585 59.929
## EA09VCU006.BRS_Sertaneja == 0 69.337 64.665 74.009
## EA09VCU006.CMG_1152 == 0 49.881 45.209 54.554
##-----------------------------------------------------------------------------
## Comparações multiplas dentro do primeiro ensaio.
## Ordenar as linhas da matriz Xs pela estimativa das médias.
a <- confint(glht(m0, linfct=Xs), calpha=univariate_calpha())
Xs <- Xs[order(a$confint[,"Estimate"], decreasing=FALSE),]
Xc <- apc(Xs)
a <- confint(glht(m0, linfct=Xs), calpha=univariate_calpha()); a
##
## Simultaneous Confidence Intervals
##
## Fit: lmer(formula = y ~ enscult + (1 | ensaio:bl), data = da, REML = FALSE)
##
## Quantile = 1.96
## 95% confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## EA09VCU006.CMG_1152 == 0 49.881 45.209 54.554
## EA09VCU006.AB062048 == 0 51.364 45.994 56.734
## EA09VCU006.BRSMG_Curinga == 0 53.129 47.758 58.500
## EA09VCU006.BRS_Primavera == 0 55.257 50.585 59.929
## EA09VCU006.AB062138 == 0 55.952 51.280 60.625
## EA09VCU006.AB062008 == 0 57.773 52.402 63.144
## EA09VCU006.BRSGO_Serra Dourada == 0 60.421 55.749 65.094
## EA09VCU006.AB062037 == 0 61.585 56.913 66.258
## EA09VCU006.BRA052034 == 0 61.713 57.040 66.385
## EA09VCU006.BRA052045 == 0 64.136 59.463 68.808
## EA09VCU006.AB062104 == 0 64.783 58.236 71.330
## EA09VCU006.BRA052023 == 0 65.581 60.909 70.254
## EA09VCU006.AB062045 == 0 66.929 62.257 71.601
## EA09VCU006.AB062041 == 0 68.342 63.670 73.014
## EA09VCU006.BRS_Sertaneja == 0 69.337 64.665 74.009
## EA09VCU006.BRA052033 == 0 69.352 64.679 74.024
Xc <- apc(Xs)
s <- summary(glht(m0, linfct=Xc), test=adjusted(type="fdr"))
ret <- cld2(s)
ci <- cbind(data.frame(a$confint),
cld=ret$mcletters$Letters)
arrange(ci, -Estimate)
## Estimate lwr upr cld
## 1 69.35 64.68 74.02 a
## 2 69.34 64.66 74.01 a
## 3 68.34 63.67 73.01 ab
## 4 66.93 62.26 71.60 abc
## 5 65.58 60.91 70.25 ad
## 6 64.78 58.24 71.33 abce
## 7 64.14 59.46 68.81 ad
## 8 61.71 57.04 66.38 bdf
## 9 61.59 56.91 66.26 bdf
## 10 60.42 55.75 65.09 cdfg
## 11 57.77 52.40 63.14 deh
## 12 55.95 51.28 60.62 efh
## 13 55.26 50.58 59.93 fh
## 14 53.13 47.76 58.50 gh
## 15 51.36 45.99 56.73 h
## 16 49.88 45.21 54.55 h
##-----------------------------------------------------------------------------
## Obter todas as médias das cult em cada ensaio.
g <- confint(glht(m0, linfct=X), calpha=univariate_calpha())
grid <- attr(X, "grid")
## Quebrando o nome no ponto.
quebra <- do.call(rbind, strsplit(grid$enscult, split="\\."))
colnames(quebra) <- c("ensaio","cult")
grid <- cbind(quebra, grid, g$confint)
grid <- equallevels(grid, da)
str(grid)
## 'data.frame': 310 obs. of 6 variables:
## $ ensaio : Factor w/ 19 levels "EA09VCU006","EA09VCU008",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ cult : Factor w/ 17 levels "AB062008","AB062037",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ enscult : Factor w/ 310 levels "EA09VCU006.AB062008",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Estimate: num 57.8 73.2 62.9 60 52.5 ...
## $ lwr : num 52.4 68.6 58.3 55.3 47.1 ...
## $ upr : num 63.1 77.9 67.6 64.7 57.8 ...
##-----------------------------------------------------------------------------
## Gráfico.
segplot(cult~lwr+upr|ensaio, centers=Estimate,
ylab="Cultivar de arroz", xlab="Produtividade",
data=grid, draw=FALSE, strip=FALSE, strip.left=TRUE,
scales=list(
y=list(relation="free", rot=0)),
layout=c(2,NA))
##-----------------------------------------------------------------------------
## Truque para ordenar cultivares dentro dos ensaios no gráfico.
## Função que remove o texto antes do ponto separador.
yscale.components.dropend <- function(...){
ans <- yscale.components.default(...)
lab <- ans$left$labels$labels
ans$left$labels$labels <- gsub("^.*\\.", "", lab)
ans
}
## Número de cult em cada ensaio.
ncult <- apply(xtabs(~ensaio+cult, da), 1, function(i) sum(i>0))
ncult
## EA09VCU006 EA09VCU008 EA09VCU009 EA09VCU017 EA09VCU018 EA09VCU022 EA09VCU023 EA09VCU025
## 16 16 16 16 16 17 17 17
## EA09VCU026 EA09VCU029 EA09VCU030 EA09VCU031 EA09VCU036 EA09VCU037 EA09VCU038 EA09VCU039
## 17 17 17 17 15 16 16 16
## EA09VCU040 EA09VCU050 EA09VCU056
## 16 16 16
## Ordem original dos níveis do fator enscult.
on <- levels(grid$enscult)
neworder <- with(grid, order(ensaio, Estimate))
orderin <- by(grid, INDICE=grid$ensaio,
function(i){
as.character(i$enscult[order(i$Estimate)])
})
## Uma cópia do fator com ondem diferente dos níveis.
grid$enscult2 <- factor(grid$enscult, levels=unlist(orderin))
segplot(enscult2~lwr+upr|ensaio, centers=Estimate,
ylab="Cultivar de arroz", xlab="raíz da Produtividade",
data=grid, draw=FALSE,
scales=list(y=list(relation="free", rot=0, tck=0.5, cex=0.7)),
yscale.components=yscale.components.dropend,
between=list(y=0.2),
layout=c(2,NA), as.table=TRUE)
