Análise de VCU de soja
##-----------------------------------------------------------------------------
## 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
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/data/grupoexperimentos.txt",
header=TRUE, sep="\t", na.string=".")
str(da)
## 'data.frame': 300 obs. of 11 variables:
## $ local: Factor w/ 4 levels "CPAO-DDOS","MARAC",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gen : Factor w/ 25 levels "BMX MAGNA RR",..: 15 15 15 16 16 16 25 25 25 9 ...
## $ trat : int 1 1 1 2 2 2 3 3 3 4 ...
## $ bl : int 1 2 3 1 2 3 1 2 3 1 ...
## $ df : int 34 35 35 44 44 44 36 36 36 37 ...
## $ dm : int 101 98 99 102 104 105 97 98 99 98 ...
## $ aca : num 1 2 1 1 1 1 1 1 1 1 ...
## $ ap : int 91 98 96 105 95 91 81 87 84 75 ...
## $ avag : int 13 16 15 15 12 13 14 16 14 14 ...
## $ rend : num 2444 2737 2948 2398 2579 ...
## $ p100 : num 11.2 11.3 11.2 10.2 11.2 ...
##-----------------------------------------------------------------------------
## Editar.
da <- transform(da,
bl=factor(bl),
blin=interaction(bl, local, drop=TRUE))
str(da)
## 'data.frame': 300 obs. of 12 variables:
## $ local: Factor w/ 4 levels "CPAO-DDOS","MARAC",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gen : Factor w/ 25 levels "BMX MAGNA RR",..: 15 15 15 16 16 16 25 25 25 9 ...
## $ trat : int 1 1 1 2 2 2 3 3 3 4 ...
## $ bl : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
## $ df : int 34 35 35 44 44 44 36 36 36 37 ...
## $ dm : int 101 98 99 102 104 105 97 98 99 98 ...
## $ aca : num 1 2 1 1 1 1 1 1 1 1 ...
## $ ap : int 91 98 96 105 95 91 81 87 84 75 ...
## $ avag : int 13 16 15 15 12 13 14 16 14 14 ...
## $ rend : num 2444 2737 2948 2398 2579 ...
## $ p100 : num 11.2 11.3 11.2 10.2 11.2 ...
## $ blin : Factor w/ 12 levels "1.CPAO-DDOS",..: 1 2 3 1 2 3 1 2 3 1 ...
sum(is.na(da$rend))
## [1] 11
da <- subset(da, !is.na(rend), select=c("local","gen","blin","rend"))
str(da)
## 'data.frame': 289 obs. of 4 variables:
## $ local: Factor w/ 4 levels "CPAO-DDOS","MARAC",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gen : Factor w/ 25 levels "BMX MAGNA RR",..: 15 15 15 16 16 25 25 25 9 9 ...
## $ blin : Factor w/ 12 levels "1.CPAO-DDOS",..: 1 2 3 1 2 1 2 3 1 2 ...
## $ rend : num 2444 2737 2948 2398 2579 ...
## Nos níveis dos fatores não pode ter nome com traço, dá problema na
## cld(). Trocar por espaço.
levels(da$gen) <- gsub("-", " ", levels(da$gen))
levels(da$local) <- gsub("-", " ", levels(da$local))
levels(da$gen)
## [1] "BMX MAGNA RR" "BMX POTÊNCIA RR" "BR 01 25656" "BR 02 04844"
## [5] "BR 02 22425" "BR 02 68661" "BR 02 72914" "BR 02 78838"
## [9] "BRS 239" "BRS 243 RR" "BRS 246 RR" "BRS 255 RR"
## [13] "BRS 268" "BRS 282" "BRS 284" "BRS 285"
## [17] "BRS 291 RR" "BRS 292 RR" "BRS 294 RR" "BRS Favorita RR"
## [21] "CD 202" "Embrapa 48" "M SOY 7908 RR" "NK 7059 RR"
## [25] "Vmax"
levels(da$local)
## [1] "CPAO DDOS" "MARAC" "NAV" "SIDROL"
##-----------------------------------------------------------------------------
## Layout.
x <- xtabs(~local+gen, data=da)
## dimnames(x) <- NULL
mosaicplot(x, off=10, las=4)
t(x)
## local
## gen CPAO DDOS MARAC NAV SIDROL
## BMX MAGNA RR 3 3 3 3
## BMX POTÊNCIA RR 3 3 3 3
## BR 01 25656 3 3 3 3
## BR 02 04844 3 2 3 3
## BR 02 22425 3 3 3 3
## BR 02 68661 3 3 3 3
## BR 02 72914 2 3 3 3
## BR 02 78838 3 3 3 3
## BRS 239 3 2 3 3
## BRS 243 RR 3 3 3 3
## BRS 246 RR 2 3 3 3
## BRS 255 RR 3 3 3 3
## BRS 268 3 3 2 3
## BRS 282 3 3 3 3
## BRS 284 3 3 2 3
## BRS 285 2 3 3 3
## BRS 291 RR 3 2 3 3
## BRS 292 RR 3 3 3 3
## BRS 294 RR 3 3 3 3
## BRS Favorita RR 3 3 3 3
## CD 202 3 3 2 3
## Embrapa 48 3 3 3 3
## M SOY 7908 RR 2 3 3 3
## NK 7059 RR 3 3 3 3
## Vmax 3 3 2 3
##-----------------------------------------------------------------------------
## Ver.
xyplot(rend~gen|local, data=da, as.table=TRUE)
##-----------------------------------------------------------------------------
## Ajuste do modelo que considera apenas o efeito de genótipo como
## fixo, os demais e interações são aleatórios.
## * fixo: gen;
## * aleatório: local, bl dentro de local e local interação com gen.
## da$y <- sqrt(da$rend)
## da$y <- log(da$rend)
da$y <- da$rend
m0 <- lmer(y~gen+(1|local)+(1|blin)+(1|local:gen),
data=da, REML=FALSE)
summary(m0)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: y ~ gen + (1 | local) + (1 | blin) + (1 | local:gen)
## Data: da
##
## AIC BIC logLik deviance df.resid
## 4442 4548 -2192 4384 260
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6264 -0.5557 -0.0478 0.5610 2.7147
##
## Random effects:
## Groups Name Variance Std.Dev.
## local:gen (Intercept) 130327 361.0
## blin (Intercept) 2330 48.3
## local (Intercept) 155226 394.0
## Residual 136231 369.1
## Number of obs: 289, groups: local:gen, 100; blin, 12; local, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3647.9 288.0 12.67
## genBMX POTÊNCIA RR 188.9 296.4 0.64
## genBR 01 25656 -25.6 296.4 -0.09
## genBR 02 04844 300.9 298.6 1.01
## genBR 02 22425 -405.3 296.4 -1.37
## genBR 02 68661 69.3 296.4 0.23
## genBR 02 72914 182.9 298.6 0.61
## genBR 02 78838 25.4 296.4 0.09
## genBRS 239 -224.9 298.6 -0.75
## genBRS 243 RR -278.7 296.4 -0.94
## genBRS 246 RR -279.8 298.6 -0.94
## genBRS 255 RR -117.8 296.4 -0.40
## genBRS 268 -334.5 298.6 -1.12
## genBRS 282 -567.9 296.4 -1.92
## genBRS 284 -124.1 298.6 -0.42
## genBRS 285 -231.2 298.6 -0.77
## genBRS 291 RR -552.3 298.6 -1.85
## genBRS 292 RR -382.4 296.4 -1.29
## genBRS 294 RR -237.6 296.4 -0.80
## genBRS Favorita RR 20.8 296.4 0.07
## genCD 202 -619.8 298.6 -2.08
## genEmbrapa 48 -539.6 296.4 -1.82
## genM SOY 7908 RR -125.8 298.6 -0.42
## genNK 7059 RR -181.8 296.4 -0.61
## genVmax -321.0 298.6 -1.07
##
## Correlation matrix not shown by default, as p = 25 > 20.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
## Abandona a interação.
m1 <- update(m0, .~.-(1|local:gen), data=da)
anova(m0, m1)
## Data: da
## Models:
## m1: y ~ gen + (1 | local) + (1 | blin)
## m0: y ~ gen + (1 | local) + (1 | blin) + (1 | local:gen)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 28 4498 4600 -2221 4442
## m0 29 4442 4549 -2192 4384 57.5 1 3.3e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Existe interação local:gen.
## Predição dos efeitos aleatórios.
lapply(ranef(m0), c)
## $`local:gen`
## $`local:gen`$`(Intercept)`
## [1] -12.112 -312.990 -100.398 -313.796 216.049 103.478 7.843 -151.404 106.546
## [10] -56.593 -119.331 -224.648 104.150 466.176 -188.231 -246.795 7.630 -211.662
## [19] -61.261 -366.083 282.578 206.647 -256.914 326.162 326.883 -137.756 -352.995
## [28] 204.368 604.739 621.789 247.846 742.952 318.692 35.043 -166.014 -99.038
## [37] -267.696 -345.011 -250.213 -30.507 -313.334 -236.468 264.362 -9.198 -147.150
## [46] -131.382 -400.463 307.601 -256.039 -325.879 165.536 483.546 353.434 -348.749
## [55] -597.197 -104.512 -656.107 -141.670 -132.629 254.426 373.047 711.019 -1.065
## [64] -503.473 -34.381 230.395 531.091 158.146 50.800 319.030 -769.016 43.068
## [73] -427.198 81.722 321.334 -15.668 182.440 -457.404 57.806 -240.641 -246.811
## [82] -94.688 -25.618 -8.960 -31.818 -154.677 -218.675 241.926 287.510 253.119
## [91] 329.734 -302.254 -210.846 19.659 194.203 617.820 150.748 376.511 -151.846
## [100] -322.338
##
##
## $blin
## $blin$`(Intercept)`
## [1] 0.5278 9.5558 -18.4518 -3.2469 10.9290 -9.8587 -46.1611 6.6619 45.9459
## [10] 0.1164 -13.4596 17.4413
##
##
## $local
## $local$`(Intercept)`
## [1] -557.5 -145.0 429.5 273.0
##-----------------------------------------------------------------------------
## Teste para os termos de efeito fixo.
anova(m0)
## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## gen 24 4527371 188640 1.38
## 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 ~ gen
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 = 33.2, df = 24, P(> X2) = 0.099
## Chi-quadrado da razão de verossimilhanças (não baseado em aproximação
## de função, melhor opção).
m1 <- update(m0, .~.-gen)
anova(m0, m1)
## Data: da
## Models:
## m1: y ~ (1 | local) + (1 | blin) + (1 | local:gen)
## m0: y ~ gen + (1 | local) + (1 | blin) + (1 | local:gen)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 5 4423 4441 -2206 4413
## m0 29 4442 4549 -2192 4384 28.8 24 0.23
##-----------------------------------------------------------------------------
## 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$local)
##-----------------------------------------------------------------------------
## Médias ajustadas.
## Formula só da parte de efeito fixo.
f <- nobars(formula(m0))
X <- LSmatrix(lm(f, da), effect=c("gen"))
rownames(X) <- levels(da$gen)
## Estimativas das médias.
summary(glht(m0, linfct=X), test=adjusted(type="none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = y ~ gen + (1 | local) + (1 | blin) + (1 | local:gen),
## data = da, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## BMX MAGNA RR == 0 3648 288 12.7 <2e-16 ***
## BMX POTÊNCIA RR == 0 3837 288 13.3 <2e-16 ***
## BR 01 25656 == 0 3622 288 12.6 <2e-16 ***
## BR 02 04844 == 0 3949 290 13.6 <2e-16 ***
## BR 02 22425 == 0 3243 288 11.3 <2e-16 ***
## BR 02 68661 == 0 3717 288 12.9 <2e-16 ***
## BR 02 72914 == 0 3831 290 13.2 <2e-16 ***
## BR 02 78838 == 0 3673 288 12.8 <2e-16 ***
## BRS 239 == 0 3423 290 11.8 <2e-16 ***
## BRS 243 RR == 0 3369 288 11.7 <2e-16 ***
## BRS 246 RR == 0 3368 290 11.6 <2e-16 ***
## BRS 255 RR == 0 3530 288 12.3 <2e-16 ***
## BRS 268 == 0 3313 290 11.4 <2e-16 ***
## BRS 282 == 0 3080 288 10.7 <2e-16 ***
## BRS 284 == 0 3524 290 12.1 <2e-16 ***
## BRS 285 == 0 3417 290 11.8 <2e-16 ***
## BRS 291 RR == 0 3096 290 10.7 <2e-16 ***
## BRS 292 RR == 0 3266 288 11.3 <2e-16 ***
## BRS 294 RR == 0 3410 288 11.8 <2e-16 ***
## BRS Favorita RR == 0 3669 288 12.7 <2e-16 ***
## CD 202 == 0 3028 290 10.4 <2e-16 ***
## Embrapa 48 == 0 3108 288 10.8 <2e-16 ***
## M SOY 7908 RR == 0 3522 290 12.1 <2e-16 ***
## NK 7059 RR == 0 3466 288 12.0 <2e-16 ***
## Vmax == 0 3327 290 11.5 <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] 300 25
g <- summary(glht(m0, linfct=Xc), test=adjusted(type="fdr")); g
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = y ~ gen + (1 | local) + (1 | blin) + (1 | local:gen),
## data = da, REML = FALSE)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## BMX MAGNA RR-BMX POTÊNCIA RR == 0 -188.90 296.43 -0.64 0.83
## BMX MAGNA RR-BR 01 25656 == 0 25.56 296.43 0.09 0.98
## BMX MAGNA RR-BR 02 04844 == 0 -300.87 298.63 -1.01 0.73
## BMX MAGNA RR-BR 02 22425 == 0 405.27 296.43 1.37 0.63
## BMX MAGNA RR-BR 02 68661 == 0 -69.29 296.43 -0.23 0.94
## BMX MAGNA RR-BR 02 72914 == 0 -182.88 298.63 -0.61 0.84
## BMX MAGNA RR-BR 02 78838 == 0 -25.39 296.43 -0.09 0.98
## BMX MAGNA RR-BRS 239 == 0 224.91 298.63 0.75 0.81
## BMX MAGNA RR-BRS 243 RR == 0 278.74 296.43 0.94 0.74
## BMX MAGNA RR-BRS 246 RR == 0 279.82 298.63 0.94 0.74
## BMX MAGNA RR-BRS 255 RR == 0 117.84 296.43 0.40 0.89
## BMX MAGNA RR-BRS 268 == 0 334.47 298.63 1.12 0.73
## BMX MAGNA RR-BRS 282 == 0 567.89 296.43 1.92 0.48
## BMX MAGNA RR-BRS 284 == 0 124.14 298.63 0.42 0.89
## BMX MAGNA RR-BRS 285 == 0 231.22 298.63 0.77 0.80
## BMX MAGNA RR-BRS 291 RR == 0 552.31 298.63 1.85 0.51
## BMX MAGNA RR-BRS 292 RR == 0 382.38 296.43 1.29 0.66
## BMX MAGNA RR-BRS 294 RR == 0 237.63 296.43 0.80 0.79
## BMX MAGNA RR-BRS Favorita RR == 0 -20.83 296.43 -0.07 0.98
## BMX MAGNA RR-CD 202 == 0 619.81 298.63 2.08 0.48
## BMX MAGNA RR-Embrapa 48 == 0 539.55 296.43 1.82 0.52
## BMX MAGNA RR-M SOY 7908 RR == 0 125.79 298.63 0.42 0.89
## BMX MAGNA RR-NK 7059 RR == 0 181.84 296.43 0.61 0.84
## BMX MAGNA RR-Vmax == 0 321.02 298.63 1.07 0.73
## BMX POTÊNCIA RR-BR 01 25656 == 0 214.46 296.43 0.72 0.82
## BMX POTÊNCIA RR-BR 02 04844 == 0 -111.97 298.63 -0.37 0.89
## BMX POTÊNCIA RR-BR 02 22425 == 0 594.17 296.43 2.00 0.48
## BMX POTÊNCIA RR-BR 02 68661 == 0 119.61 296.43 0.40 0.89
## BMX POTÊNCIA RR-BR 02 72914 == 0 6.02 298.63 0.02 0.99
## BMX POTÊNCIA RR-BR 02 78838 == 0 163.51 296.43 0.55 0.86
## BMX POTÊNCIA RR-BRS 239 == 0 413.81 298.63 1.39 0.63
## BMX POTÊNCIA RR-BRS 243 RR == 0 467.63 296.43 1.58 0.62
## BMX POTÊNCIA RR-BRS 246 RR == 0 468.71 298.63 1.57 0.62
## BMX POTÊNCIA RR-BRS 255 RR == 0 306.73 296.43 1.03 0.73
## BMX POTÊNCIA RR-BRS 268 == 0 523.37 298.63 1.75 0.54
## BMX POTÊNCIA RR-BRS 282 == 0 756.78 296.43 2.55 0.39
## BMX POTÊNCIA RR-BRS 284 == 0 313.04 298.63 1.05 0.73
## BMX POTÊNCIA RR-BRS 285 == 0 420.12 298.63 1.41 0.63
## BMX POTÊNCIA RR-BRS 291 RR == 0 741.21 298.63 2.48 0.39
## BMX POTÊNCIA RR-BRS 292 RR == 0 571.27 296.43 1.93 0.48
## BMX POTÊNCIA RR-BRS 294 RR == 0 426.52 296.43 1.44 0.63
## BMX POTÊNCIA RR-BRS Favorita RR == 0 168.07 296.43 0.57 0.86
## BMX POTÊNCIA RR-CD 202 == 0 808.70 298.63 2.71 0.38
## BMX POTÊNCIA RR-Embrapa 48 == 0 728.45 296.43 2.46 0.39
## BMX POTÊNCIA RR-M SOY 7908 RR == 0 314.69 298.63 1.05 0.73
## BMX POTÊNCIA RR-NK 7059 RR == 0 370.74 296.43 1.25 0.69
## BMX POTÊNCIA RR-Vmax == 0 509.92 298.63 1.71 0.55
## BR 01 25656-BR 02 04844 == 0 -326.43 298.63 -1.09 0.73
## BR 01 25656-BR 02 22425 == 0 379.71 296.43 1.28 0.66
## BR 01 25656-BR 02 68661 == 0 -94.85 296.43 -0.32 0.89
## BR 01 25656-BR 02 72914 == 0 -208.44 298.63 -0.70 0.82
## BR 01 25656-BR 02 78838 == 0 -50.95 296.43 -0.17 0.94
## BR 01 25656-BRS 239 == 0 199.35 298.63 0.67 0.82
## BR 01 25656-BRS 243 RR == 0 253.17 296.43 0.85 0.76
## BR 01 25656-BRS 246 RR == 0 254.25 298.63 0.85 0.76
## BR 01 25656-BRS 255 RR == 0 92.27 296.43 0.31 0.90
## BR 01 25656-BRS 268 == 0 308.91 298.63 1.03 0.73
## BR 01 25656-BRS 282 == 0 542.32 296.43 1.83 0.52
## BR 01 25656-BRS 284 == 0 98.58 298.63 0.33 0.89
## BR 01 25656-BRS 285 == 0 205.66 298.63 0.69 0.82
## BR 01 25656-BRS 291 RR == 0 526.75 298.63 1.76 0.54
## BR 01 25656-BRS 292 RR == 0 356.81 296.43 1.20 0.71
## BR 01 25656-BRS 294 RR == 0 212.06 296.43 0.72 0.82
## BR 01 25656-BRS Favorita RR == 0 -46.39 296.43 -0.16 0.94
## BR 01 25656-CD 202 == 0 594.24 298.63 1.99 0.48
## BR 01 25656-Embrapa 48 == 0 513.99 296.43 1.73 0.54
## BR 01 25656-M SOY 7908 RR == 0 100.23 298.63 0.34 0.89
## BR 01 25656-NK 7059 RR == 0 156.28 296.43 0.53 0.86
## BR 01 25656-Vmax == 0 295.46 298.63 0.99 0.73
## BR 02 04844-BR 02 22425 == 0 706.14 298.63 2.36 0.42
## BR 02 04844-BR 02 68661 == 0 231.58 298.63 0.78 0.80
## BR 02 04844-BR 02 72914 == 0 117.99 300.83 0.39 0.89
## BR 02 04844-BR 02 78838 == 0 275.49 298.63 0.92 0.75
## BR 02 04844-BRS 239 == 0 525.78 300.83 1.75 0.54
## BR 02 04844-BRS 243 RR == 0 579.61 298.63 1.94 0.48
## BR 02 04844-BRS 246 RR == 0 580.69 300.83 1.93 0.48
## BR 02 04844-BRS 255 RR == 0 418.71 298.63 1.40 0.63
## BR 02 04844-BRS 268 == 0 635.35 300.83 2.11 0.48
## BR 02 04844-BRS 282 == 0 868.76 298.63 2.91 0.37
## BR 02 04844-BRS 284 == 0 425.01 300.83 1.41 0.63
## BR 02 04844-BRS 285 == 0 532.09 300.83 1.77 0.54
## BR 02 04844-BRS 291 RR == 0 853.18 300.83 2.84 0.37
## BR 02 04844-BRS 292 RR == 0 683.25 298.63 2.29 0.44
## BR 02 04844-BRS 294 RR == 0 538.50 298.63 1.80 0.52
## BR 02 04844-BRS Favorita RR == 0 280.04 298.63 0.94 0.74
## BR 02 04844-CD 202 == 0 920.68 300.83 3.06 0.37
## BR 02 04844-Embrapa 48 == 0 840.43 298.63 2.81 0.37
## BR 02 04844-M SOY 7908 RR == 0 426.66 300.83 1.42 0.63
## BR 02 04844-NK 7059 RR == 0 482.71 298.63 1.62 0.60
## BR 02 04844-Vmax == 0 621.89 300.83 2.07 0.48
## BR 02 22425-BR 02 68661 == 0 -474.56 296.43 -1.60 0.61
## BR 02 22425-BR 02 72914 == 0 -588.15 298.63 -1.97 0.48
## BR 02 22425-BR 02 78838 == 0 -430.66 296.43 -1.45 0.63
## BR 02 22425-BRS 239 == 0 -180.36 298.63 -0.60 0.84
## BR 02 22425-BRS 243 RR == 0 -126.54 296.43 -0.43 0.89
## BR 02 22425-BRS 246 RR == 0 -125.46 298.63 -0.42 0.89
## BR 02 22425-BRS 255 RR == 0 -287.44 296.43 -0.97 0.74
## BR 02 22425-BRS 268 == 0 -70.80 298.63 -0.24 0.94
## BR 02 22425-BRS 282 == 0 162.61 296.43 0.55 0.86
## BR 02 22425-BRS 284 == 0 -281.13 298.63 -0.94 0.74
## BR 02 22425-BRS 285 == 0 -174.05 298.63 -0.58 0.86
## BR 02 22425-BRS 291 RR == 0 147.04 298.63 0.49 0.86
## BR 02 22425-BRS 292 RR == 0 -22.90 296.43 -0.08 0.98
## BR 02 22425-BRS 294 RR == 0 -167.65 296.43 -0.57 0.86
## BR 02 22425-BRS Favorita RR == 0 -426.10 296.43 -1.44 0.63
## BR 02 22425-CD 202 == 0 214.53 298.63 0.72 0.82
## BR 02 22425-Embrapa 48 == 0 134.28 296.43 0.45 0.88
## BR 02 22425-M SOY 7908 RR == 0 -279.48 298.63 -0.94 0.74
## BR 02 22425-NK 7059 RR == 0 -223.44 296.43 -0.75 0.81
## BR 02 22425-Vmax == 0 -84.25 298.63 -0.28 0.91
## BR 02 68661-BR 02 72914 == 0 -113.59 298.63 -0.38 0.89
## BR 02 68661-BR 02 78838 == 0 43.90 296.43 0.15 0.94
## BR 02 68661-BRS 239 == 0 294.20 298.63 0.99 0.73
## BR 02 68661-BRS 243 RR == 0 348.02 296.43 1.17 0.73
## BR 02 68661-BRS 246 RR == 0 349.10 298.63 1.17 0.73
## BR 02 68661-BRS 255 RR == 0 187.12 296.43 0.63 0.83
## BR 02 68661-BRS 268 == 0 403.76 298.63 1.35 0.63
## BR 02 68661-BRS 282 == 0 637.17 296.43 2.15 0.48
## BR 02 68661-BRS 284 == 0 193.43 298.63 0.65 0.83
## BR 02 68661-BRS 285 == 0 300.51 298.63 1.01 0.73
## BR 02 68661-BRS 291 RR == 0 621.60 298.63 2.08 0.48
## BR 02 68661-BRS 292 RR == 0 451.66 296.43 1.52 0.63
## BR 02 68661-BRS 294 RR == 0 306.91 296.43 1.04 0.73
## BR 02 68661-BRS Favorita RR == 0 48.46 296.43 0.16 0.94
## BR 02 68661-CD 202 == 0 689.09 298.63 2.31 0.44
## BR 02 68661-Embrapa 48 == 0 608.84 296.43 2.05 0.48
## BR 02 68661-M SOY 7908 RR == 0 195.08 298.63 0.65 0.83
## BR 02 68661-NK 7059 RR == 0 251.13 296.43 0.85 0.76
## BR 02 68661-Vmax == 0 390.31 298.63 1.31 0.66
## BR 02 72914-BR 02 78838 == 0 157.49 298.63 0.53 0.86
## BR 02 72914-BRS 239 == 0 407.79 300.83 1.36 0.63
## BR 02 72914-BRS 243 RR == 0 461.61 298.63 1.55 0.63
## BR 02 72914-BRS 246 RR == 0 462.69 300.83 1.54 0.63
## BR 02 72914-BRS 255 RR == 0 300.71 298.63 1.01 0.73
## BR 02 72914-BRS 268 == 0 517.35 300.83 1.72 0.55
## BR 02 72914-BRS 282 == 0 750.76 298.63 2.51 0.39
## BR 02 72914-BRS 284 == 0 307.02 300.83 1.02 0.73
## BR 02 72914-BRS 285 == 0 414.10 300.83 1.38 0.63
## BR 02 72914-BRS 291 RR == 0 735.19 300.83 2.44 0.39
## BR 02 72914-BRS 292 RR == 0 565.25 298.63 1.89 0.48
## BR 02 72914-BRS 294 RR == 0 420.50 298.63 1.41 0.63
## BR 02 72914-BRS Favorita RR == 0 162.05 298.63 0.54 0.86
## BR 02 72914-CD 202 == 0 802.68 300.83 2.67 0.38
## BR 02 72914-Embrapa 48 == 0 722.43 298.63 2.42 0.39
## BR 02 72914-M SOY 7908 RR == 0 308.67 300.76 1.03 0.73
## BR 02 72914-NK 7059 RR == 0 364.71 298.63 1.22 0.71
## BR 02 72914-Vmax == 0 503.90 300.83 1.68 0.56
## BR 02 78838-BRS 239 == 0 250.29 298.63 0.84 0.76
## BR 02 78838-BRS 243 RR == 0 304.12 296.43 1.03 0.73
## BR 02 78838-BRS 246 RR == 0 305.20 298.63 1.02 0.73
## BR 02 78838-BRS 255 RR == 0 143.22 296.43 0.48 0.86
## BR 02 78838-BRS 268 == 0 359.86 298.63 1.21 0.71
## BR 02 78838-BRS 282 == 0 593.27 296.43 2.00 0.48
## BR 02 78838-BRS 284 == 0 149.53 298.63 0.50 0.86
## BR 02 78838-BRS 285 == 0 256.61 298.63 0.86 0.76
## BR 02 78838-BRS 291 RR == 0 577.70 298.63 1.93 0.48
## BR 02 78838-BRS 292 RR == 0 407.76 296.43 1.38 0.63
## BR 02 78838-BRS 294 RR == 0 263.01 296.43 0.89 0.76
## BR 02 78838-BRS Favorita RR == 0 4.56 296.43 0.02 0.99
## BR 02 78838-CD 202 == 0 645.19 298.63 2.16 0.48
## BR 02 78838-Embrapa 48 == 0 564.94 296.43 1.91 0.48
## BR 02 78838-M SOY 7908 RR == 0 151.18 298.63 0.51 0.86
## BR 02 78838-NK 7059 RR == 0 207.22 296.43 0.70 0.82
## BR 02 78838-Vmax == 0 346.41 298.63 1.16 0.73
## BRS 239-BRS 243 RR == 0 53.83 298.63 0.18 0.94
## BRS 239-BRS 246 RR == 0 54.91 300.83 0.18 0.94
## BRS 239-BRS 255 RR == 0 -107.07 298.63 -0.36 0.89
## BRS 239-BRS 268 == 0 109.57 300.83 0.36 0.89
## BRS 239-BRS 282 == 0 342.98 298.63 1.15 0.73
## BRS 239-BRS 284 == 0 -100.77 300.83 -0.33 0.89
## BRS 239-BRS 285 == 0 6.31 300.83 0.02 0.99
## BRS 239-BRS 291 RR == 0 327.40 300.83 1.09 0.73
## BRS 239-BRS 292 RR == 0 157.47 298.63 0.53 0.86
## BRS 239-BRS 294 RR == 0 12.72 298.63 0.04 0.99
## BRS 239-BRS Favorita RR == 0 -245.74 298.63 -0.82 0.77
## BRS 239-CD 202 == 0 394.90 300.83 1.31 0.66
## BRS 239-Embrapa 48 == 0 314.65 298.63 1.05 0.73
## BRS 239-M SOY 7908 RR == 0 -99.12 300.83 -0.33 0.89
## BRS 239-NK 7059 RR == 0 -43.07 298.63 -0.14 0.94
## BRS 239-Vmax == 0 96.11 300.83 0.32 0.89
## BRS 243 RR-BRS 246 RR == 0 1.08 298.63 0.00 1.00
## BRS 243 RR-BRS 255 RR == 0 -160.90 296.43 -0.54 0.86
## BRS 243 RR-BRS 268 == 0 55.74 298.63 0.19 0.94
## BRS 243 RR-BRS 282 == 0 289.15 296.43 0.98 0.74
## BRS 243 RR-BRS 284 == 0 -154.59 298.63 -0.52 0.86
## BRS 243 RR-BRS 285 == 0 -47.51 298.63 -0.16 0.94
## BRS 243 RR-BRS 291 RR == 0 273.57 298.63 0.92 0.75
## BRS 243 RR-BRS 292 RR == 0 103.64 296.43 0.35 0.89
## BRS 243 RR-BRS 294 RR == 0 -41.11 296.43 -0.14 0.94
## BRS 243 RR-BRS Favorita RR == 0 -299.57 296.43 -1.01 0.73
## BRS 243 RR-CD 202 == 0 341.07 298.63 1.14 0.73
## BRS 243 RR-Embrapa 48 == 0 260.82 296.43 0.88 0.76
## BRS 243 RR-M SOY 7908 RR == 0 -152.94 298.63 -0.51 0.86
## BRS 243 RR-NK 7059 RR == 0 -96.90 296.43 -0.33 0.89
## BRS 243 RR-Vmax == 0 42.29 298.63 0.14 0.94
## BRS 246 RR-BRS 255 RR == 0 -161.98 298.63 -0.54 0.86
## BRS 246 RR-BRS 268 == 0 54.66 300.83 0.18 0.94
## BRS 246 RR-BRS 282 == 0 288.07 298.63 0.96 0.74
## BRS 246 RR-BRS 284 == 0 -155.67 300.83 -0.52 0.86
## BRS 246 RR-BRS 285 == 0 -48.59 300.76 -0.16 0.94
## BRS 246 RR-BRS 291 RR == 0 272.50 300.83 0.91 0.76
## BRS 246 RR-BRS 292 RR == 0 102.56 298.63 0.34 0.89
## BRS 246 RR-BRS 294 RR == 0 -42.19 298.63 -0.14 0.94
## BRS 246 RR-BRS Favorita RR == 0 -300.65 298.63 -1.01 0.73
## BRS 246 RR-CD 202 == 0 339.99 300.83 1.13 0.73
## BRS 246 RR-Embrapa 48 == 0 259.74 298.63 0.87 0.76
## BRS 246 RR-M SOY 7908 RR == 0 -154.02 300.83 -0.51 0.86
## BRS 246 RR-NK 7059 RR == 0 -97.98 298.63 -0.33 0.89
## BRS 246 RR-Vmax == 0 41.21 300.83 0.14 0.94
## BRS 255 RR-BRS 268 == 0 216.64 298.63 0.73 0.82
## BRS 255 RR-BRS 282 == 0 450.05 296.43 1.52 0.63
## BRS 255 RR-BRS 284 == 0 6.31 298.63 0.02 0.99
## BRS 255 RR-BRS 285 == 0 113.39 298.63 0.38 0.89
## BRS 255 RR-BRS 291 RR == 0 434.48 298.63 1.45 0.63
## BRS 255 RR-BRS 292 RR == 0 264.54 296.43 0.89 0.76
## BRS 255 RR-BRS 294 RR == 0 119.79 296.43 0.40 0.89
## BRS 255 RR-BRS Favorita RR == 0 -138.66 296.43 -0.47 0.87
## BRS 255 RR-CD 202 == 0 501.97 298.63 1.68 0.56
## BRS 255 RR-Embrapa 48 == 0 421.72 296.43 1.42 0.63
## BRS 255 RR-M SOY 7908 RR == 0 7.96 298.63 0.03 0.99
## BRS 255 RR-NK 7059 RR == 0 64.00 296.43 0.22 0.94
## BRS 255 RR-Vmax == 0 203.19 298.63 0.68 0.82
## BRS 268-BRS 282 == 0 233.41 298.63 0.78 0.80
## BRS 268-BRS 284 == 0 -210.33 300.83 -0.70 0.82
## BRS 268-BRS 285 == 0 -103.25 300.83 -0.34 0.89
## BRS 268-BRS 291 RR == 0 217.84 300.83 0.72 0.82
## BRS 268-BRS 292 RR == 0 47.90 298.63 0.16 0.94
## BRS 268-BRS 294 RR == 0 -96.85 298.63 -0.32 0.89
## BRS 268-BRS Favorita RR == 0 -355.30 298.63 -1.19 0.72
## BRS 268-CD 202 == 0 285.33 300.83 0.95 0.74
## BRS 268-Embrapa 48 == 0 205.08 298.63 0.69 0.82
## BRS 268-M SOY 7908 RR == 0 -208.68 300.83 -0.69 0.82
## BRS 268-NK 7059 RR == 0 -152.64 298.63 -0.51 0.86
## BRS 268-Vmax == 0 -13.45 300.83 -0.04 0.99
## BRS 282-BRS 284 == 0 -443.74 298.63 -1.49 0.63
## BRS 282-BRS 285 == 0 -336.66 298.63 -1.13 0.73
## BRS 282-BRS 291 RR == 0 -15.58 298.63 -0.05 0.99
## BRS 282-BRS 292 RR == 0 -185.51 296.43 -0.63 0.83
## BRS 282-BRS 294 RR == 0 -330.26 296.43 -1.11 0.73
## BRS 282-BRS Favorita RR == 0 -588.72 296.43 -1.99 0.48
## BRS 282-CD 202 == 0 51.92 298.63 0.17 0.94
## BRS 282-Embrapa 48 == 0 -28.33 296.43 -0.10 0.98
## BRS 282-M SOY 7908 RR == 0 -442.09 298.63 -1.48 0.63
## BRS 282-NK 7059 RR == 0 -386.05 296.43 -1.30 0.66
## BRS 282-Vmax == 0 -246.86 298.63 -0.83 0.77
## BRS 284-BRS 285 == 0 107.08 300.83 0.36 0.89
## BRS 284-BRS 291 RR == 0 428.17 300.83 1.42 0.63
## BRS 284-BRS 292 RR == 0 258.23 298.63 0.86 0.76
## BRS 284-BRS 294 RR == 0 113.48 298.63 0.38 0.89
## BRS 284-BRS Favorita RR == 0 -144.97 298.63 -0.49 0.86
## BRS 284-CD 202 == 0 495.66 300.76 1.65 0.58
## BRS 284-Embrapa 48 == 0 415.41 298.63 1.39 0.63
## BRS 284-M SOY 7908 RR == 0 1.65 300.83 0.01 1.00
## BRS 284-NK 7059 RR == 0 57.69 298.63 0.19 0.94
## BRS 284-Vmax == 0 196.88 300.76 0.65 0.83
## BRS 285-BRS 291 RR == 0 321.09 300.83 1.07 0.73
## BRS 285-BRS 292 RR == 0 151.15 298.63 0.51 0.86
## BRS 285-BRS 294 RR == 0 6.40 298.63 0.02 0.99
## BRS 285-BRS Favorita RR == 0 -252.05 298.63 -0.84 0.76
## BRS 285-CD 202 == 0 388.58 300.83 1.29 0.66
## BRS 285-Embrapa 48 == 0 308.33 298.63 1.03 0.73
## BRS 285-M SOY 7908 RR == 0 -105.43 300.83 -0.35 0.89
## BRS 285-NK 7059 RR == 0 -49.39 298.63 -0.17 0.94
## BRS 285-Vmax == 0 89.80 300.83 0.30 0.90
## BRS 291 RR-BRS 292 RR == 0 -169.93 298.63 -0.57 0.86
## BRS 291 RR-BRS 294 RR == 0 -314.68 298.63 -1.05 0.73
## BRS 291 RR-BRS Favorita RR == 0 -573.14 298.63 -1.92 0.48
## BRS 291 RR-CD 202 == 0 67.49 300.83 0.22 0.94
## BRS 291 RR-Embrapa 48 == 0 -12.76 298.63 -0.04 0.99
## BRS 291 RR-M SOY 7908 RR == 0 -426.52 300.83 -1.42 0.63
## BRS 291 RR-NK 7059 RR == 0 -370.47 298.63 -1.24 0.69
## BRS 291 RR-Vmax == 0 -231.29 300.83 -0.77 0.80
## BRS 292 RR-BRS 294 RR == 0 -144.75 296.43 -0.49 0.86
## BRS 292 RR-BRS Favorita RR == 0 -403.21 296.43 -1.36 0.63
## BRS 292 RR-CD 202 == 0 237.43 298.63 0.80 0.79
## BRS 292 RR-Embrapa 48 == 0 157.18 296.43 0.53 0.86
## BRS 292 RR-M SOY 7908 RR == 0 -256.58 298.63 -0.86 0.76
## BRS 292 RR-NK 7059 RR == 0 -200.54 296.43 -0.68 0.82
## BRS 292 RR-Vmax == 0 -61.35 298.63 -0.21 0.94
## BRS 294 RR-BRS Favorita RR == 0 -258.46 296.43 -0.87 0.76
## BRS 294 RR-CD 202 == 0 382.18 298.63 1.28 0.66
## BRS 294 RR-Embrapa 48 == 0 301.93 296.43 1.02 0.73
## BRS 294 RR-M SOY 7908 RR == 0 -111.83 298.63 -0.37 0.89
## BRS 294 RR-NK 7059 RR == 0 -55.79 296.43 -0.19 0.94
## BRS 294 RR-Vmax == 0 83.40 298.63 0.28 0.91
## BRS Favorita RR-CD 202 == 0 640.64 298.63 2.15 0.48
## BRS Favorita RR-Embrapa 48 == 0 560.38 296.43 1.89 0.48
## BRS Favorita RR-M SOY 7908 RR == 0 146.62 298.63 0.49 0.86
## BRS Favorita RR-NK 7059 RR == 0 202.67 296.43 0.68 0.82
## BRS Favorita RR-Vmax == 0 341.85 298.63 1.14 0.73
## CD 202-Embrapa 48 == 0 -80.25 298.63 -0.27 0.92
## CD 202-M SOY 7908 RR == 0 -494.01 300.83 -1.64 0.58
## CD 202-NK 7059 RR == 0 -437.97 298.63 -1.47 0.63
## CD 202-Vmax == 0 -298.78 300.76 -0.99 0.73
## Embrapa 48-M SOY 7908 RR == 0 -413.76 298.63 -1.39 0.63
## Embrapa 48-NK 7059 RR == 0 -357.72 296.43 -1.21 0.71
## Embrapa 48-Vmax == 0 -218.53 298.63 -0.73 0.82
## M SOY 7908 RR-NK 7059 RR == 0 56.05 298.63 0.19 0.94
## M SOY 7908 RR-Vmax == 0 195.23 300.83 0.65 0.83
## NK 7059 RR-Vmax == 0 139.19 298.63 0.47 0.87
## (Adjusted p values reported -- fdr method)
## Estimativas das médias com comparações.
## grid <- apmc(X=X, model=m0, focus="gen", test="bonferroni")
grid <- apmc(X=X, model=m0, focus="gen", test="fdr")
grid
## gen estimate lwr upr cld
## 1 BMX MAGNA RR 3648 3083 4212 a
## 2 BMX POTÊNCIA RR 3837 3272 4401 a
## 3 BR 01 25656 3622 3058 4187 a
## 4 BR 02 04844 3949 3380 4518 a
## 5 BR 02 22425 3243 2678 3807 a
## 6 BR 02 68661 3717 3153 4282 a
## 7 BR 02 72914 3831 3262 4400 a
## 8 BR 02 78838 3673 3109 4238 a
## 9 BRS 239 3423 2854 3992 a
## 10 BRS 243 RR 3369 2805 3934 a
## 11 BRS 246 RR 3368 2799 3937 a
## 12 BRS 255 RR 3530 2966 4095 a
## 13 BRS 268 3313 2745 3882 a
## 14 BRS 282 3080 2516 3644 a
## 15 BRS 284 3524 2955 4093 a
## 16 BRS 285 3417 2848 3986 a
## 17 BRS 291 RR 3096 2527 3664 a
## 18 BRS 292 RR 3266 2701 3830 a
## 19 BRS 294 RR 3410 2846 3975 a
## 20 BRS Favorita RR 3669 3104 4233 a
## 21 CD 202 3028 2459 3597 a
## 22 Embrapa 48 3108 2544 3673 a
## 23 M SOY 7908 RR 3522 2953 4091 a
## 24 NK 7059 RR 3466 2902 4031 a
## 25 Vmax 3327 2758 3896 a
##-----------------------------------------------------------------------------
## Gráfico.
segplot(reorder(gen, estimate)~lwr+upr, centers=estimate,
ylab="Cultivar de soja", xlab="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 local:gen, tratar como efeito fixo.
formula(m0)
## y ~ gen + (1 | local) + (1 | blin) + (1 | local:gen)
m0 <- lmer(y~local*gen+(1|blin),
data=da, REML=FALSE)
##-----------------------------------------------------------------------------
## Fazer o desdobramento dentro do primeiro local. Os demais seguem a
## mesma regra.
X <- LSmatrix(lm(nobars(formula(m0)), data=da), effect=c("local","gen"))
dim(X)
## [1] 100 100
L <- by(X, attr(X, "grid")$local, as.matrix)
L <- lapply(L, "rownames<-", levels(da$gen))
str(L)
## List of 4
## $ CPAO DDOS: num [1:25, 1:100] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:25] "BMX MAGNA RR" "BMX POTÊNCIA RR" "BR 01 25656" "BR 02 04844" ...
## .. ..$ : chr [1:100] "(Intercept)" "localMARAC" "localNAV" "localSIDROL" ...
## $ MARAC : num [1:25, 1:100] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:25] "BMX MAGNA RR" "BMX POTÊNCIA RR" "BR 01 25656" "BR 02 04844" ...
## .. ..$ : chr [1:100] "(Intercept)" "localMARAC" "localNAV" "localSIDROL" ...
## $ NAV : num [1:25, 1:100] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:25] "BMX MAGNA RR" "BMX POTÊNCIA RR" "BR 01 25656" "BR 02 04844" ...
## .. ..$ : chr [1:100] "(Intercept)" "localMARAC" "localNAV" "localSIDROL" ...
## $ SIDROL : num [1:25, 1:100] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:25] "BMX MAGNA RR" "BMX POTÊNCIA RR" "BR 01 25656" "BR 02 04844" ...
## .. ..$ : chr [1:100] "(Intercept)" "localMARAC" "localNAV" "localSIDROL" ...
grid <- lapply(L, apmc, model=m0, focus="gen", test="fdr")
str(grid)
## List of 4
## $ CPAO DDOS:'data.frame': 25 obs. of 5 variables:
## ..$ gen : Factor w/ 25 levels "BMX MAGNA RR",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..$ estimate: num [1:25] 3071 2854 2927 2965 2974 ...
## ..$ lwr : num [1:25] 2731 2514 2586 2625 2633 ...
## ..$ upr : num [1:25] 3412 3195 3267 3306 3314 ...
## ..$ cld : chr [1:25] "ab" "ab" "ab" "ab" ...
## $ MARAC :'data.frame': 25 obs. of 5 variables:
## ..$ gen : Factor w/ 25 levels "BMX MAGNA RR",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..$ estimate: num [1:25] 3316 3215 3752 4725 3935 ...
## ..$ lwr : num [1:25] 2976 2875 3412 4309 3595 ...
## ..$ upr : num [1:25] 3657 3556 4093 5142 4276 ...
## ..$ cld : chr [1:25] "cef" "efgi" "cdh" "a" ...
## $ NAV :'data.frame': 25 obs. of 5 variables:
## ..$ gen : Factor w/ 25 levels "BMX MAGNA RR",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..$ estimate: num [1:25] 4303 4920 4531 3910 2869 ...
## ..$ lwr : num [1:25] 3962 4580 4190 3570 2529 ...
## ..$ upr : num [1:25] 4643 5261 4871 4251 3209 ...
## ..$ cld : chr [1:25] "cd" "a" "ac" "bdf" ...
## $ SIDROL :'data.frame': 25 obs. of 5 variables:
## ..$ gen : Factor w/ 25 levels "BMX MAGNA RR",..: 1 2 3 4 5 6 7 8 9 10 ...
## ..$ estimate: num [1:25] 3901 4357 3280 4301 3193 ...
## ..$ lwr : num [1:25] 3561 4017 2939 3961 2852 ...
## ..$ upr : num [1:25] 4242 4698 3620 4642 3533 ...
## ..$ cld : chr [1:25] "dfg" "f" "ac" "bf" ...
grid <- ldply(grid); names(grid)[1] <- "local"
grid <- equallevels(grid, da)
str(grid)
## 'data.frame': 100 obs. of 6 variables:
## $ local : Factor w/ 4 levels "CPAO DDOS","MARAC",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gen : Factor w/ 25 levels "BMX MAGNA RR",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ estimate: num 3071 2854 2927 2965 2974 ...
## $ lwr : num 2731 2514 2586 2625 2633 ...
## $ upr : num 3412 3195 3267 3306 3314 ...
## $ cld : chr "ab" "ab" "ab" "ab" ...
##-----------------------------------------------------------------------------
## Gráfico.
segplot(gen~lwr+upr|local, centers=estimate,
ylab="Genivar 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 gen dentro dos locais 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 gen em cada local.
ngen <- apply(xtabs(~local+gen, da), 1, function(i) sum(i>0))
ngen
## CPAO DDOS MARAC NAV SIDROL
## 25 25 25 25
grid$locgen <- with(grid, interaction(local, gen))
## Ordem original dos níveis do fator locgen.
on <- levels(grid$locgen)
neworder <- with(grid, order(local, estimate))
orderin <- by(grid, INDICE=grid$local,
function(i){
as.character(i$locgen[order(i$estimate)])
})
## Uma cópia do fator com ondem diferente dos níveis.
grid$locgen2 <- factor(grid$locgen, levels=unlist(orderin))
segplot(locgen2~lwr+upr|local, centers=estimate,
ylab="Cultivar de soja", xlab="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)
