Walmes Zeviani
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "lme4", "plyr",
"aod")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra doBy multcomp lme4 plyr aod
## TRUE TRUE TRUE TRUE 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.0 (2014-04-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=pt_BR.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=pt_BR.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] aod_1.3 plyr_1.8 lme4_1.0-5 Matrix_1.1-3
## [5] multcomp_1.2-21 mvtnorm_0.9-9994 doBy_4.5-10 MASS_7.3-33
## [9] survival_2.37-7 latticeExtra_0.6-24 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.1 formatR_0.9 grid_3.1.0 minqa_1.2.2 nlme_3.1-110 stringr_0.6.2
## [7] tools_3.1.0
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" "AB062045"
## [5] "AB062048" "AB062104" "AB062138" "BRA052023"
## [9] "BRA052033" "BRA052034" "BRA052045" "BRSGO_Serra Dourada"
## [13] "BRSMG_Curinga" "BRS_Primavera" "BRS_Sertaneja" "CMG_1152"
## [17] "H1"
levels(da$ensaio)
## [1] "EA09VCU006" "EA09VCU008" "EA09VCU009" "EA09VCU017" "EA09VCU018" "EA09VCU022" "EA09VCU023"
## [8] "EA09VCU025" "EA09VCU026" "EA09VCU029" "EA09VCU030" "EA09VCU031" "EA09VCU036" "EA09VCU037"
## [15] "EA09VCU038" "EA09VCU039" "EA09VCU040" "EA09VCU050" "EA09VCU056"
##-----------------------------------------------------------------------------
## Layout.
x <- xtabs(~ensaio+cult, data=da)
## dimnames(x) <- NULL
x
## cult
## ensaio AB062008 AB062037 AB062041 AB062045 AB062048 AB062104 AB062138 BRA052023 BRA052033
## EA09VCU006 3 4 4 4 3 2 4 4 4
## EA09VCU008 4 4 4 4 4 4 4 4 4
## EA09VCU009 4 4 4 4 4 4 4 4 4
## EA09VCU017 4 4 4 4 4 4 4 4 4
## EA09VCU018 3 4 4 3 3 3 4 4 4
## EA09VCU022 4 4 4 4 4 4 4 4 3
## EA09VCU023 4 4 4 4 4 4 4 4 4
## EA09VCU025 4 4 4 4 4 4 4 4 4
## EA09VCU026 4 4 4 4 4 4 4 4 4
## EA09VCU029 4 4 4 4 4 4 4 4 4
## EA09VCU030 4 4 4 4 4 4 4 4 4
## EA09VCU031 4 4 4 4 4 4 4 4 4
## EA09VCU036 3 4 4 3 3 4 4 4 4
## EA09VCU037 4 4 4 4 4 4 4 4 4
## EA09VCU038 4 4 4 4 4 4 4 4 4
## EA09VCU039 4 4 4 4 4 4 4 4 4
## EA09VCU040 4 4 4 4 4 4 4 4 4
## EA09VCU050 4 4 4 4 4 4 4 4 4
## EA09VCU056 4 4 4 4 4 4 4 4 3
## cult
## ensaio BRA052034 BRA052045 BRSGO_Serra Dourada BRSMG_Curinga BRS_Primavera BRS_Sertaneja
## EA09VCU006 4 4 4 3 4 4
## EA09VCU008 4 4 4 4 4 4
## EA09VCU009 4 4 4 4 4 4
## EA09VCU017 4 4 4 4 4 4
## EA09VCU018 4 4 3 3 4 4
## EA09VCU022 4 4 4 4 4 4
## EA09VCU023 4 4 4 4 4 4
## EA09VCU025 4 4 4 4 4 4
## EA09VCU026 4 4 4 4 4 4
## EA09VCU029 4 4 4 4 4 4
## EA09VCU030 4 4 4 4 4 4
## EA09VCU031 4 4 4 4 4 4
## EA09VCU036 4 4 4 4 0 4
## EA09VCU037 4 4 4 4 4 4
## EA09VCU038 4 4 4 4 4 4
## EA09VCU039 4 4 4 4 4 4
## EA09VCU040 4 4 4 4 4 4
## EA09VCU050 4 4 4 4 4 4
## EA09VCU056 4 4 4 4 4 4
## cult
## ensaio CMG_1152 H1
## EA09VCU006 4 0
## EA09VCU008 4 0
## EA09VCU009 4 0
## EA09VCU017 3 0
## EA09VCU018 4 0
## EA09VCU022 4 3
## EA09VCU023 4 4
## EA09VCU025 4 4
## EA09VCU026 4 4
## EA09VCU029 4 4
## EA09VCU030 4 4
## EA09VCU031 4 4
## EA09VCU036 3 0
## EA09VCU037 4 0
## EA09VCU038 4 0
## EA09VCU039 4 0
## EA09VCU040 3 0
## EA09VCU050 4 0
## EA09VCU056 4 0
##-----------------------------------------------------------------------------
## Ver.
xyplot(prod~cult|ensaio, data=da)
##-----------------------------------------------------------------------------
## 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
## 8132 8239 -4045 8090
##
## 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 cBRA052033
## 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 0.502
## clBRA052045 -0.354 0.503 0.503 0.502 0.502 0.501 0.503 0.503 0.502
## cltBRSGO_SD -0.353 0.503 0.503 0.502 0.501 0.501 0.503 0.503 0.502
## cltBRSMG_Cr -0.353 0.502 0.502 0.501 0.500 0.500 0.502 0.502 0.501
## cltBRS_Prmv -0.349 0.496 0.496 0.495 0.494 0.494 0.496 0.496 0.495
## cltBRS_Srtn -0.354 0.503 0.503 0.502 0.502 0.501 0.503 0.503 0.502
## cltCMG_1152 -0.353 0.502 0.502 0.500 0.500 0.500 0.502 0.502 0.501
## cultH1 -0.255 0.363 0.363 0.362 0.362 0.362 0.363 0.363 0.362
## cBRA052034 cBRA05204 cBRSGD cBRSMG cBRS_P cBRS_S cCMG_1
## cltAB062037
## cltAB062041
## cltAB062045
## cltAB062048
## cltAB062104
## cltAB062138
## clBRA052023
## clBRA052033
## clBRA052034
## clBRA052045 0.503
## cltBRSGO_SD 0.503 0.503
## cltBRSMG_Cr 0.502 0.502 0.501
## cltBRS_Prmv 0.496 0.496 0.495 0.495
## cltBRS_Srtn 0.503 0.503 0.503 0.502 0.496
## cltCMG_1152 0.502 0.502 0.501 0.500 0.494 0.502
## cultH1 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 4.18480
## [10] -2.29390 2.88689 -0.63400 -1.93824 -2.47383 6.75154 -4.68854 0.86237 0.97966
## [19] -0.11547 0.30360 0.77902 0.13170 2.45682 -6.63295 4.00586 0.61846 -4.42697
## [28] -2.87491 -1.72847 2.31666 3.47303 2.73059 -4.81966 2.78984 0.95819 -0.09177
## [37] 1.29147 -2.89770 0.35472 5.24165 1.67010 -0.59983 0.66762 -4.43800 2.06264
## [46] 0.19727 -1.35551 1.17418 -1.43011 1.04700 -3.20342 -0.30287 0.56930 2.25830
## [55] -5.59013 -2.64498 1.08031 2.64248 2.66101 -3.64234 3.21152 5.11934 -0.70051
## [64] -0.70744 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 3.34359
## [82] 1.88810 1.74742 2.67323 3.00407 1.92374 -8.07050 -6.32087 -1.26288 2.91123
## [91] -0.61189 3.49278 -8.48065 6.20034 0.26076 9.39534 -10.75283 0.58836 2.74771
## [100] -3.20979 2.27157 2.30492 -0.25785 0.07463 6.07015 -7.90748 2.23157 2.59991
## [109] 0.48453 -3.84321 -6.05739 7.00991 -3.23794 0.02838 -1.72012 -1.43281 -4.40067
## [118] -1.89248 -1.19365 -2.25991 -1.17625 -2.30334 -2.48334 1.45476 0.54478 -0.99995
## [127] 4.11722 0.15665 1.92024 4.09915 2.52709 0.85527 -3.16430 1.90197 2.64858
## [136] -4.61689 -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 2.29243
## [154] -3.15407 2.56207 -1.80946 -3.50562 -4.19856 0.74123 -0.33622 -0.57205 6.36900
## [163] -2.24373 -0.78727 3.68269 0.19775 -0.31370 2.47426 2.07426 3.18364 1.28146
## [172] -0.69831 0.78709 -2.38608 -5.30012 -0.11018 -1.71046 -2.25318 -0.88865 2.12257
## [181] -1.48256 2.10977 -6.37378 4.03589 1.92111 3.26072 0.68570 5.35984 -5.21568
## [190] -3.37169 3.21174 -0.30206 2.48500 0.96640 5.13148 -0.68116 -3.55814 0.29236
## [199] -6.83386 3.96870 -0.88793 1.50170 -7.10415 -5.83591 -1.08633 1.04237 8.94752
## [208] 1.92426 1.48197 -3.69372 3.07148 2.17886 -7.33850 0.31260 5.50026 4.24638
## [217] 4.80210 -0.12514 1.98326 -1.32606 2.91053 -0.41463 4.30112 -0.34461 -1.20697
## [226] 4.65820 1.57812 -15.35474 1.20405 -10.40959 8.68346 -1.58484 4.21149 -9.03063
## [235] 3.47425 0.32169 -2.14663 0.25097 -4.05587 2.95149 -11.13702 7.85321 -0.63913
## [244] -3.74364 -2.61566 5.42754 2.11950 -6.95640 -6.97926 -1.53647 -6.78594 -4.08226
## [253] 7.55832 2.97906 5.36991 2.00565 10.53526 -2.71909 6.42736 -4.96990 -2.82673
## [262] 1.69795 -1.25100 -1.00577 -2.68915 -2.10799 -0.79120 -1.18733 -3.07100 -0.70909
## [271] -0.52350 -2.67091 -0.07579 -1.97587 -0.15410 7.28005 6.95356 2.32278 -4.74683
## [280] -3.47991 -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 1.77127
## [298] 6.00044 2.53874 8.10975 0.83639 1.98285 -1.95206 -0.71545 -2.27759 -4.24629
## [307] 1.31626 -2.88910 -6.42625 -0.22703
##
##
## $`ensaio:bl`
## $`ensaio:bl`$`(Intercept)`
## [1] -1.149052 -0.805353 0.896687 1.169605 -0.003725 -0.895021 -0.031755 1.231627 1.011008
## [10] 0.675104 -0.259515 -1.195947 -1.111428 0.519820 1.063905 -0.433864 1.006309 0.576912
## [19] -1.799872 -0.050145 -1.697355 2.087295 0.588831 -0.838513 -0.170985 -0.694252 0.114137
## [28] 0.949617 -0.308502 1.407372 -0.100469 -1.525828 0.415087 -0.131149 -0.587646 0.540098
## [37] -0.107367 0.377821 -0.165000 0.132901 -2.249656 1.110511 1.230161 -0.186453 -2.184552
## [46] -1.110153 1.048896 2.351851 -0.479350 -0.919155 -0.379656 1.619483 -0.544100 0.607737
## [55] 0.822489 -0.676701 -1.006698 0.122182 0.348434 0.349976 -0.708539 0.019243 0.064723
## [64] 0.816707 0.440223 -0.151110 -0.964315 0.501962 1.896316 -0.281803 0.205618 -1.992245
## [73] 0.050965 -1.284280 -0.033626 0.843520
##
##
## $ensaio
## $ensaio$`(Intercept)`
## [1] 4.319 11.623 8.903 1.484 -10.298 5.414 7.663 -20.359 9.125 9.201 -3.684 4.093
## [13] -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. Leia:
## https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
## Alternativas (em order de recomendação):
## * Teste de razão de verossimilhanças entre modelos aninhados (um com
## outro sem o termo fixo de interesse, usar REML=FALSE).
## * 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
wt$result
## $chi2
## chi2 df P
## 8.129e+01 1.600e+01 9.727e-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"))
qqmath(r)
qqmath(~r|da$ensaio)
qqmath(bec)
qqmath(beb)
qqmath(be)
##-----------------------------------------------------------------------------
## 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 <- desdobr.glht(X=X, m0=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)
## Error: rank of X = 310 < ncol(X) = 323
## Erro devido cultivares que não ocorrem em alguns ensaios.
with(da, nlevels(interaction(ensaio, cult, drop=TRUE)))
## [1] 310
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)
## 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
glht(m0, linfct=Xs)
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## EA09VCU006.AB062008 == 0 57.8
## EA09VCU006.AB062037 == 0 61.6
## EA09VCU006.AB062041 == 0 68.3
## EA09VCU006.AB062045 == 0 66.9
## EA09VCU006.AB062048 == 0 51.4
## EA09VCU006.AB062104 == 0 64.8
## EA09VCU006.AB062138 == 0 56.0
## EA09VCU006.BRA052023 == 0 65.6
## EA09VCU006.BRA052033 == 0 69.4
## EA09VCU006.BRA052034 == 0 61.7
## EA09VCU006.BRA052045 == 0 64.1
## EA09VCU006.BRSGO_Serra Dourada == 0 60.4
## EA09VCU006.BRSMG_Curinga == 0 53.1
## EA09VCU006.BRS_Primavera == 0 55.3
## EA09VCU006.BRS_Sertaneja == 0 69.3
## EA09VCU006.CMG_1152 == 0 49.9
##-----------------------------------------------------------------------------
## 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 EA09VCU026
## 16 16 16 16 16 17 17 17 17
## EA09VCU029 EA09VCU030 EA09VCU031 EA09VCU036 EA09VCU037 EA09VCU038 EA09VCU039 EA09VCU040 EA09VCU050
## 17 17 17 15 16 16 16 16 16
## EA09VCU056
## 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)
