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

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)

plot of chunk unnamed-chunk-3

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)

plot of chunk unnamed-chunk-3

## 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"))

plot of chunk unnamed-chunk-4

xyplot(sqrt(abs(r))~f, type=c("p", "smooth"))

plot of chunk unnamed-chunk-4

grid.arrange(qqmath(r),
             qqmath(bec),
             qqmath(beb),
             qqmath(be),
             nrow=2)

plot of chunk unnamed-chunk-4

qqmath(~r|da$ensaio)

plot of chunk unnamed-chunk-4

##-----------------------------------------------------------------------------
## 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)
        })

plot of chunk unnamed-chunk-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))

plot of chunk unnamed-chunk-5

##-----------------------------------------------------------------------------
## 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)