Upland rice yield as affected by row spacing, cultivar and water stress
Alexandre Bryan Heinemann
Adriano Stephan Nascente
Walmes Marques Zeviani
Luís Fernando Stone
Paulo Cesar Sentelhas
Data and analyses description
Four field experiments were conducted during the growing seasons of 2009 and 2010 in the state of Goiás and Tocantins, Central Brazil. Three experiments were carried out in 2009 in Formoso (lat. 12.3ºS and long. 49.8ºW), Porangatu (lat. 13.8ºS and long. 50.0ºW) and Santo Antônio de Goiás (lat. 16.5ºS and long. 49.3ºW) and one experiment in 2010 in Santo Antônio de Goiás.
The experiments were arranged as a randomized block design in a split-plot scheme with rowing spacing as the main plots and cultivar treatments as the subplots, totaling 12 treatments. Each treatment was replicated four times for each sowing date, totaling 192 evaluated subplots. For each experiment rice crop was sown at 0.25, 0.35, 0.45 and 0.55 m between rows, keeping the same seed density per area (230 seeds m-2).
The sowing dates were 12/12/2009, 16/11/2009, 13/11/2009 and 16/12/2010 respectively for Formoso, Porangatu, Santo Antônio de Goiás – 2009 and Santo Antônio de Goiás – 2010. Irrigation was applied only in the 2010 experiment in Santo Antônio de Goiás with a sprinkler system.
Yield and yield components (panicle per m2, 100-grain weight, percentage of sterile spikelets and number of grains per panicle) were measured from an area of 6 m2 in the center of each subplot, with each sample consisted of four rows of 0.5 m. Grain yield was expressed at 13% of moisture content.
Yield and yield components data were analyzed using linear mixed effects models. Block and plot were considered as a random effect and row spacing and cultivar as fixed effect. The statistical model parameters were estimated by maximum-likelihood. The Wald’s test was applied to determine the significance of fixed terms effect. All inferences were applied for 5% of significance. Statistical analyses were done in the R software using the nlme R package.
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sess?o.
pkg <- c("lattice", "latticeExtra", "gridExtra", "nlme",
"doBy", "multcomp", "plyr", "wzRfun")
sapply(pkg, require, character.only=TRUE)## lattice latticeExtra gridExtra nlme doBy multcomp
## TRUE TRUE TRUE TRUE TRUE TRUE
## plyr wzRfun
## TRUE TRUE## Instruções para instalação do pacote wzRfun em:
## https://github.com/walmes/wzRfun
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] splines grid stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.4 plyr_1.8.1 multcomp_1.3-6 TH.data_1.0-3
## [5] mvtnorm_0.9-9997 doBy_4.5-10 MASS_7.3-34 survival_2.37-7
## [9] nlme_3.1-117 gridExtra_0.9.1 latticeExtra_0.6-26 RColorBrewer_1.0-5
## [13] lattice_0.20-29 rmarkdown_0.2.68 knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 htmltools_0.2.6 lme4_1.1-7
## [6] Matrix_1.1-4 methods_3.1.1 minqa_1.2.3 nloptr_1.0.4 Rcpp_0.11.0
## [11] sandwich_2.3-0 stringr_0.6.2 tools_3.1.1 zoo_1.7-10##-----------------------------------------------------------------------------
## Usados nos gráficos.
## col <- gray.colors(3)[c(1, 1.3, 1.6)]
## col <- gray.colors(3)
col <- gray.colors(3)[c(1,1,1)]
cult.lab <- c("BRSMGCuringa", "Douradão", "BRSPrimavera")
la.lab <- c("Porangatu", "Formoso", "SGA - 2010", "SGA - 2009")
ef.names <- c("constante", "ambiente", "espaçaamento linear",
"espaçamento quadrático", "cultivar", "amb:esp lin",
"amb:esp qua", "amb:cult", "esp lin:cult", "esp quad:cult",
"amb:esp lin:cult",
"amb:esp qua:cult")
##-----------------------------------------------------------------------------
## Ler.
## Se não tiver os dados, fazer o download do arquivo.
## url <- "http://www.leg.ufpr.br/~walmes/data/grupoexpcov.txt"
## download.file(url=url, destfile=basename(url))
gec <- read.table("grupoexpcov.txt", header=TRUE, sep="\t")
str(gec)## 'data.frame': 192 obs. of 21 variables:
## $ local : int 1 1 1 1 1 1 1 1 1 1 ...
## $ ano : int 2009 2009 2009 2009 2009 2009 2009 2009 2009 2009 ...
## $ parc : int 102 104 109 112 202 206 209 210 301 305 ...
## $ bl : int 1 1 1 1 2 2 2 2 3 3 ...
## $ cult : Factor w/ 3 levels "C","D","P": 1 1 1 1 1 1 1 1 1 1 ...
## $ esp : num 0.45 0.35 0.55 0.25 0.45 0.35 0.55 0.25 0.45 0.35 ...
## $ rend : int 935 928 1928 1795 861 1499 1653 1821 1784 2272 ...
## $ ester : num 11.2 26.7 33.5 34.7 11.8 ...
## $ npan : num 296 326 329 480 304 ...
## $ p100 : num NA 2.8 2.6 2.7 3.5 3.3 NA 2.9 3.5 2 ...
## $ ng_pan: num 124.5 87.7 118.9 75.4 154.8 ...
## $ prect : num 811 811 811 811 811 811 811 811 811 811 ...
## $ precv : num 569 569 569 569 569 ...
## $ precr : num 128 128 128 128 128 ...
## $ gdet : num 1897 1897 1897 1897 1897 ...
## $ gdev : num 968 968 968 968 968 968 968 968 968 968 ...
## $ gder : num 437 437 437 437 437 437 437 437 437 437 ...
## $ rad : num 2112 2112 2112 2112 2112 ...
## $ txdesv: num 0.000671 0.000671 0.000671 0.000671 0.000671 ...
## $ tdesr : num 8e-04 8e-04 8e-04 8e-04 8e-04 ...
## $ local2: Factor w/ 4 levels "CNPAF09","CNPAF10",..: 1 1 1 1 1 1 1 1 1 1 ...## summary(gec)
## head(gec)
## Cria variáveis fatores.
gec <- transform(gec,
local=factor(local),
ano=factor(ano),
bl=factor(bl),
la=interaction(local, ano, drop=TRUE),
espac=ordered(esp))
## Troca ordem dos níveis.
## levels(gec$la)
gec$la <- factor(gec$la, levels=c("1.2010","1.2009","2.2009","3.2009"))
gec$la.bl <- interaction(gec$la, gec$bl)
## Espaçamento centrado em 0.4.
gec$esp0 <- gec$esp
gec$esp <- gec$esp-0.4
##-----------------------------------------------------------------------------
## Estrutura.
names(gec)[c(1:6,20:21)] ## Fatores experimentais.## [1] "local" "ano" "parc" "bl" "cult" "esp" "tdesr" "local2"names(gec)[7:10] ## Respostas observadas.## [1] "rend" "ester" "npan" "p100"names(gec)[11:19] ## Covariáveis auxiliares observadas.## [1] "ng_pan" "prect" "precv" "precr" "gdet" "gdev" "gder" "rad" "txdesv"##-----------------------------------------------------------------------------
## Informações.
## Em cada local tem-se um experimento em parcela subdividida em que a
## parcela recebe os níveis de espaçamento e a subparcela recebe os
## níveis de cultivar, delineamento em blocos.
## Fazer análise via modelo de parcela subdividida com estimação REML ou
## ML bloco e parcela aleatórios, local, cultivar e espaçaamento fixo
## espaçaamento efeito contínuo usando polinômiosRendimento de grãos
##-----------------------------------------------------------------------------
## Ver o rendimento de grãos.
xyplot(rend~esp|la, groups=cult,
data=subset(gec, !is.na(rend)), as.table=TRUE,
auto.key=TRUE, jitter.x=TRUE,
## scales=list(y=list(log=10)),
type=c("p","a"),
par.settings=list(superpose.symbol=list(pch=19)),
xlab="Espaçamento-40 (m)",
ylab=expression("Rendimento de graõs"~(kg~ha^{-1})))
xyplot(rend^0.33~esp|la, groups=cult,
data=subset(gec, !is.na(rend)), as.table=TRUE,
auto.key=TRUE, jitter.x=TRUE,
type=c("p","a"),
par.settings=list(superpose.symbol=list(pch=19)),
xlab="Espaçamento-40 (m)",
ylab=expression("Rendimento de graõs"~(kg~ha^{-1})))
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m0 <- lme(rend~(la*(esp+I(esp^2))*cult),
random=~1|la.bl/espac,
data=gec, na.action=na.omit,
method="ML")
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r))
## Embora esteja considerando só a parte de efeito fixo serve de
## indicação.
MASS::boxcox(lm(formula(m0), data=gec))
abline(v=1/3, col=2)
##-----------------------------------------------------------------------------
## Com transformação.
## Ajuste do modelo com a resposta transformada.
m0 <- update(m0, rend^(1/3)~.)
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r))
## Quadro de testes F de Wald para os termos de efeito fixo.
anova(m0)## numDF denDF F-value p-value
## (Intercept) 1 79 6959 <.0001
## la 3 12 367 <.0001
## esp 1 40 0 0.5122
## I(esp^2) 1 40 2 0.1876
## cult 2 79 25 <.0001
## la:esp 3 40 0 0.9224
## la:I(esp^2) 3 40 4 0.0176
## la:cult 6 79 22 <.0001
## esp:cult 2 79 3 0.0732
## I(esp^2):cult 2 79 0 0.6352
## la:esp:cult 6 79 1 0.2620
## la:I(esp^2):cult 6 79 1 0.5749## Quadro de anova com nomes modificados.
## Teste de Wald para os fatores fixos do modelo.
aux <- anova(m0)
rownames(aux) <- ef.names
aux## numDF denDF F-value p-value
## constante 1 79 6959 <.0001
## ambiente 3 12 367 <.0001
## espaçaamento linear 1 40 0 0.5122
## espaçamento quadrático 1 40 2 0.1876
## cultivar 2 79 25 <.0001
## amb:esp lin 3 40 0 0.9224
## amb:esp qua 3 40 4 0.0176
## amb:cult 6 79 22 <.0001
## esp lin:cult 2 79 3 0.0732
## esp quad:cult 2 79 0 0.6352
## amb:esp lin:cult 6 79 1 0.2620
## amb:esp qua:cult 6 79 1 0.5749## Em modelos linares mistos não se tem um quadro de análise de
## variância como se tem para modelos de efeito fixo. O "dispositivo" de
## teste de hipóteses é o teste de Wald. Em modelos lineares de efeito
## fixo fazer o teste de Wald e fazer a análise de variância retornam o
## mesmo valor da estatística F e mesmo p-valor. Para modelos mistos,
## como não se tem anova, usa-se o teste de Wald. A função anova() do R
## é uma função método que só retorna quadro de anova para objetos de
## classe lm. Para objetos de classe glm, por exemplo, ela retorna um
## quadro de deviance e para objetos de classe lme, retorna o quadro de
## testes de Wald.
## Estimativas dos parâmetros e medidas de ajuste.
summary(m0)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## 576 697.6 -249
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 0.2427
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 0.4838 0.9608
##
## Fixed effects: rend^(1/3) ~ la + esp + I(esp^2) + cult + la:esp + la:I(esp^2) + la:cult + esp:cult + I(esp^2):cult + la:esp:cult + la:I(esp^2):cult
## Value Std.Error DF t-value p-value
## (Intercept) 18.16 0.58 79 31.515 0.0000
## la1.2009 -7.07 0.77 12 -9.231 0.0000
## la2.2009 -13.74 0.81 12 -16.892 0.0000
## la3.2009 -8.59 0.77 12 -11.213 0.0000
## esp -2.86 2.95 40 -0.969 0.3384
## I(esp^2) -79.85 33.95 40 -2.352 0.0237
## cultD -3.01 0.74 79 -4.065 0.0001
## cultP -1.94 0.76 79 -2.572 0.0120
## la1.2009:esp 1.70 4.01 40 0.425 0.6732
## la2.2009:esp 3.33 4.03 40 0.825 0.4145
## la3.2009:esp 5.78 4.01 40 1.441 0.1572
## la1.2009:I(esp^2) 111.32 45.55 40 2.444 0.0190
## la2.2009:I(esp^2) 61.91 47.14 40 1.313 0.1965
## la3.2009:I(esp^2) 23.69 45.55 40 0.520 0.6058
## la1.2009:cultD 5.20 0.96 79 5.411 0.0000
## la2.2009:cultD 2.67 1.00 79 2.668 0.0093
## la3.2009:cultD 7.32 0.96 79 7.619 0.0000
## la1.2009:cultP 1.94 0.97 79 1.992 0.0499
## la2.2009:cultP 0.81 1.07 79 0.752 0.4541
## la3.2009:cultP 4.71 1.01 79 4.665 0.0000
## esp:cultD 5.48 3.90 79 1.405 0.1638
## esp:cultP 4.17 5.08 79 0.820 0.4144
## I(esp^2):cultD 93.79 44.82 79 2.093 0.0396
## I(esp^2):cultP 40.80 53.52 79 0.762 0.4481
## la1.2009:esp:cultD -2.52 5.19 79 -0.485 0.6292
## la2.2009:esp:cultD -5.37 5.31 79 -1.010 0.3155
## la3.2009:esp:cultD -9.58 5.19 79 -1.844 0.0689
## la1.2009:esp:cultP -5.32 6.13 79 -0.868 0.3879
## la2.2009:esp:cultP -6.62 6.43 79 -1.030 0.3063
## la3.2009:esp:cultP -14.41 6.14 79 -2.347 0.0214
## la1.2009:I(esp^2):cultD -92.64 58.99 79 -1.570 0.1203
## la2.2009:I(esp^2):cultD -86.85 60.80 79 -1.429 0.1571
## la3.2009:I(esp^2):cultD -102.18 58.99 79 -1.732 0.0871
## la1.2009:I(esp^2):cultP -37.50 65.84 79 -0.570 0.5706
## la2.2009:I(esp^2):cultP 5.64 70.09 79 0.080 0.9361
## la3.2009:I(esp^2):cultP -57.54 66.94 79 -0.860 0.3926
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.45677 -0.49480 0.05747 0.53768 2.52253
##
## Number of Observations: 167
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64tTable <- as.data.frame(summary(m0)$tTable)
## tTable[tTable[,5]<0.05,]
## fixef(m0)
## tTable
##-----------------------------------------------------------------------------
## Modelo com termos não significativos abandonados.
m1 <- update(m0, .~cult*(la+esp)+la*I(esp^2))
anova(m0, m1)## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 39 576.0 697.6 -249.0
## m1 2 22 558.5 627.1 -257.3 1 vs 2 16.47 0.4909anova(m1)## numDF denDF F-value p-value
## (Intercept) 1 93 7535 <.0001
## cult 2 93 22 <.0001
## la 3 12 399 <.0001
## esp 1 43 1 0.4056
## I(esp^2) 1 43 1 0.2354
## cult:la 6 93 22 <.0001
## cult:esp 2 93 3 0.0729
## la:I(esp^2) 3 43 5 0.0076summary(m1)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## 558.5 627.1 -257.3
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 0.2556
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 0.4347 1.033
##
## Fixed effects: rend^(1/3) ~ cult + la + esp + I(esp^2) + cult:la + cult:esp + la:I(esp^2)
## Value Std.Error DF t-value p-value
## (Intercept) 17.60 0.466 93 37.76 0.0000
## cultD -1.74 0.447 93 -3.89 0.0002
## cultP -1.44 0.513 93 -2.81 0.0061
## la1.2009 -6.53 0.620 12 -10.53 0.0000
## la2.2009 -13.38 0.648 12 -20.66 0.0000
## la3.2009 -7.93 0.622 12 -12.76 0.0000
## esp -0.11 1.361 43 -0.08 0.9374
## I(esp^2) -41.38 22.796 43 -1.82 0.0765
## cultD:la1.2009 3.95 0.592 93 6.67 0.0000
## cultP:la1.2009 1.48 0.643 93 2.30 0.0237
## cultD:la2.2009 1.51 0.608 93 2.49 0.0146
## cultP:la2.2009 0.89 0.683 93 1.31 0.1936
## cultD:la3.2009 5.95 0.592 93 10.06 0.0000
## cultP:la3.2009 3.99 0.653 93 6.10 0.0000
## cultD:esp 0.86 1.805 93 0.48 0.6356
## cultP:esp -3.44 1.908 93 -1.80 0.0748
## la1.2009:I(esp^2) 74.33 30.060 43 2.47 0.0174
## la2.2009:I(esp^2) 37.76 31.062 43 1.22 0.2308
## la3.2009:I(esp^2) -22.59 30.294 43 -0.75 0.4599
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.48027 -0.53285 0.04548 0.57352 2.52598
##
## Number of Observations: 167
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64##-----------------------------------------------------------------------------
## Predição.
predict.lme <- function(lme.model, newdata){
## F <- model.matrix(terms(m0), newdata)
F <- model.matrix(formula(lme.model)[-2], newdata)
U <- chol(vcov(lme.model))
haty <- F%*%fixef(lme.model)
se <- sqrt(apply(F%*%t(U), 1, function(x) sum(x^2)))
me <- outer(se,
## qt(c(fit=.5, lwr=.025, upr=.975),
## df=length(residuals(m0))-length(fixef(m0))),
qnorm(c(fit=.5, lwr=.025, upr=.975)))
## aux <- apply(limt, 2, function(x) x+haty)
aux <- sweep(me, MARGIN=1, STATS=haty, FUN="+")
pred <- cbind(pred, aux)
return(pred)
}
pred <- expand.grid(la=levels(gec$la),
cult=levels(gec$cult),
esp=seq(0.25, 0.55, by=0.025)-0.4)
## pred <- predict.lme(m0, pred)
pred <- predict.lme(m1, pred)
pred$esp <- pred$esp+0.4
str(pred)## 'data.frame': 156 obs. of 6 variables:
## $ la : Factor w/ 4 levels "1.2010","1.2009",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ cult: Factor w/ 3 levels "C","D","P": 1 1 1 1 2 2 2 2 3 3 ...
## $ esp : num 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ fit : num 16.68 11.83 4.15 8.25 14.81 ...
## $ lwr : num 15.81 11.03 3.33 7.44 13.92 ...
## $ upr : num 17.55 12.63 4.97 9.05 15.71 ...##-----------------------------------------------------------------------------
## Gráfico.
ylab <- expression("Cubinc root of grain yield "*(kg~ha^{-1}))
xlab <- expression("Row spacing (m)")
## tiff(file="Figure.tiff", width=17.15, height=17.15,
## units="cm", res=1200, pointsize=10, compression="lzw")
p1 <- xyplot(rend^(1/3)~esp0|la, data=gec,
groups=cult, pch=1:3,
ylab=ylab, xlab=xlab, layout=c(4,1),
## ylim=extendrange(gec$rend, f=0.1),
col=col,
strip=strip.custom(bg="gray90", factor.levels=la.lab),
key=list(
## x=0.55, y=0.95, corner=c(0,1),
x=0.025, y=0.05, corner=c(0,0),
lines=list(pch=1:3, lty=1:3, type="b", col=col, size=6),
text=list(lab=cult.lab, cex=0.8)))
p2 <- xyplot(fit~esp|la,
data=pred, groups=cult,
type="l", col=col, lty=1:3,
ly=pred$lwr, uy=pred$upr, cty="bands",
alpha=0.4,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)## x11(width=8, height=4)
p1+as.layer(p2, under=TRUE)
Esterelidade
##-----------------------------------------------------------------------------
## Ver esterelidade.
xyplot(ester~esp|la, groups=cult,
data=subset(gec, !is.na(ester)), as.table=TRUE,
auto.key=TRUE, jitter.x=TRUE,
type=c("p","a"),
par.settings=list(superpose.symbol=list(pch=19)),
xlab="Espaçamento-40 (m)",
ylab="Esterelidade (%)")
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m0 <- lme(ester~(la*(esp+I(esp^2))*cult),
random=~1|la.bl/espac,
data=gec, na.action=na.omit,
method="ML")
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r))
## Embora esteja considerando só a parte de efeito fixo serve de
## indicação.
MASS::boxcox(lm(formula(m0), data=gec))
abline(v=1, col=2)
## Quadro de testes F de Wald para os termos de efeito fixo.
anova(m0)## numDF denDF F-value p-value
## (Intercept) 1 79 2970.5 <.0001
## la 3 12 354.1 <.0001
## esp 1 40 3.7 0.0632
## I(esp^2) 1 40 6.2 0.0172
## cult 2 79 46.7 <.0001
## la:esp 3 40 2.3 0.0938
## la:I(esp^2) 3 40 1.1 0.3436
## la:cult 6 79 45.5 <.0001
## esp:cult 2 79 1.3 0.2900
## I(esp^2):cult 2 79 1.7 0.1880
## la:esp:cult 6 79 0.5 0.7854
## la:I(esp^2):cult 6 79 1.1 0.3843## Quadro de anova com nomes modificados.
## Teste de Wald para os fatores fixos do modelo.
aux <- anova(m0)
rownames(aux) <- ef.names
aux## numDF denDF F-value p-value
## constante 1 79 2970.5 <.0001
## ambiente 3 12 354.1 <.0001
## espaçaamento linear 1 40 3.7 0.0632
## espaçamento quadrático 1 40 6.2 0.0172
## cultivar 2 79 46.7 <.0001
## amb:esp lin 3 40 2.3 0.0938
## amb:esp qua 3 40 1.1 0.3436
## amb:cult 6 79 45.5 <.0001
## esp lin:cult 2 79 1.3 0.2900
## esp quad:cult 2 79 1.7 0.1880
## amb:esp lin:cult 6 79 0.5 0.7854
## amb:esp qua:cult 6 79 1.1 0.3843## Estimativas dos parâmetros e medidas de ajuste.
summary(m0)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## 1167 1288 -544.4
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 1.742
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 0.9282 6.056
##
## Fixed effects: ester ~ (la * (esp + I(esp^2)) * cult)
## Value Std.Error DF t-value p-value
## (Intercept) 10.2 2.9 79 3.457 0.0009
## la1.2009 16.0 4.2 12 3.840 0.0024
## la2.2009 47.4 4.5 12 10.564 0.0000
## la3.2009 21.5 4.2 12 5.168 0.0002
## esp 2.2 15.5 40 0.145 0.8851
## I(esp^2) 162.5 173.0 40 0.940 0.3531
## cultD 6.2 3.9 79 1.602 0.1132
## cultP 27.2 4.0 79 6.772 0.0000
## la1.2009:esp -61.6 22.7 40 -2.713 0.0098
## la2.2009:esp -37.1 26.7 40 -1.389 0.1725
## la3.2009:esp 3.4 22.7 40 0.149 0.8823
## la1.2009:I(esp^2) 289.2 249.8 40 1.158 0.2537
## la2.2009:I(esp^2) -107.5 284.4 40 -0.378 0.7073
## la3.2009:I(esp^2) -162.6 249.8 40 -0.651 0.5186
## la1.2009:cultD 5.7 5.7 79 1.002 0.3195
## la2.2009:cultD 25.0 6.0 79 4.206 0.0001
## la3.2009:cultD -16.9 5.5 79 -3.095 0.0027
## la1.2009:cultP -34.1 6.0 79 -5.705 0.0000
## la2.2009:cultP 7.1 5.8 79 1.222 0.2252
## la3.2009:cultP -32.0 5.6 79 -5.739 0.0000
## esp:cultD 9.4 21.6 79 0.436 0.6641
## esp:cultP 3.9 22.6 79 0.174 0.8620
## I(esp^2):cultD 75.6 241.8 79 0.313 0.7553
## I(esp^2):cultP 259.4 252.1 79 1.029 0.3065
## la1.2009:esp:cultD 7.8 36.5 79 0.214 0.8312
## la2.2009:esp:cultD 1.8 34.3 79 0.052 0.9589
## la3.2009:esp:cultD -14.6 31.2 79 -0.467 0.6418
## la1.2009:esp:cultP 43.1 32.6 79 1.320 0.1906
## la2.2009:esp:cultP 31.1 35.2 79 0.882 0.3805
## la3.2009:esp:cultP -3.0 33.4 79 -0.090 0.9286
## la1.2009:I(esp^2):cultD -844.2 386.3 79 -2.185 0.0318
## la2.2009:I(esp^2):cultD -94.1 377.9 79 -0.249 0.8041
## la3.2009:I(esp^2):cultD -193.0 345.6 79 -0.558 0.5782
## la1.2009:I(esp^2):cultP -359.0 368.4 79 -0.974 0.3328
## la2.2009:I(esp^2):cultP -324.3 381.6 79 -0.850 0.3980
## la3.2009:I(esp^2):cultP -91.7 362.8 79 -0.253 0.8012
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.13850 -0.58463 -0.05225 0.53240 2.75533
##
## Number of Observations: 167
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64tTable <- as.data.frame(summary(m0)$tTable)
## tTable[tTable[,5]<0.05,]
## fixef(m0)
## tTable
##-----------------------------------------------------------------------------
## Modelo com termos não significativos abandonados.
m1 <- update(m0, .~la*(cult+esp)+I(esp^2))
anova(m0, m1)## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 39 1167 1288 -544.4
## m1 2 20 1152 1215 -556.2 1 vs 2 23.53 0.215anova(m1)## numDF denDF F-value p-value
## (Intercept) 1 95 3383 <.0001
## la 3 12 403 <.0001
## cult 2 95 46 <.0001
## esp 1 43 3 0.1022
## I(esp^2) 1 43 5 0.0287
## la:cult 6 95 44 <.0001
## la:esp 3 43 3 0.0276summary(m1)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## 1152 1215 -556.2
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 1.582
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 0.7067 6.577
##
## Fixed effects: ester ~ la + cult + esp + I(esp^2) + la:cult + la:esp
## Value Std.Error DF t-value p-value
## (Intercept) 10.60 2.05 95 5.163 0.0000
## la1.2009 19.07 2.78 12 6.872 0.0000
## la2.2009 46.38 3.06 12 15.141 0.0000
## la3.2009 19.50 2.78 12 7.025 0.0000
## cultD 7.14 2.45 95 2.911 0.0045
## cultP 30.45 2.54 95 11.972 0.0000
## esp 5.61 9.32 43 0.602 0.5500
## I(esp^2) 126.78 55.15 43 2.299 0.0265
## la1.2009:cultD -3.09 3.79 95 -0.816 0.4168
## la2.2009:cultD 23.45 3.80 95 6.176 0.0000
## la3.2009:cultD -19.40 3.50 95 -5.541 0.0000
## la1.2009:cultP -38.02 3.72 95 -10.212 0.0000
## la2.2009:cultP 3.02 3.82 95 0.791 0.4310
## la3.2009:cultP -33.39 3.63 95 -9.210 0.0000
## la1.2009:esp -38.75 14.17 43 -2.734 0.0090
## la2.2009:esp -23.56 13.87 43 -1.699 0.0966
## la3.2009:esp -1.81 13.50 43 -0.134 0.8940
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.51826 -0.58702 -0.07987 0.54066 2.80970
##
## Number of Observations: 167
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64##-----------------------------------------------------------------------------
## Predição.
pred <- expand.grid(la=levels(gec$la),
cult=levels(gec$cult),
esp=seq(0.25, 0.55, by=0.025)-0.4)
## pred <- predict.lme(m0, pred)
pred <- predict.lme(m1, pred)
pred$esp <- pred$esp+0.4
str(pred)## 'data.frame': 156 obs. of 6 variables:
## $ la : Factor w/ 4 levels "1.2010","1.2009",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ cult: Factor w/ 3 levels "C","D","P": 1 1 1 1 2 2 2 2 3 3 ...
## $ esp : num 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ fit : num 12.6 37.5 62.5 32.4 19.8 ...
## $ lwr : num 8.08 32.53 57.2 27.8 15.22 ...
## $ upr : num 17.2 42.5 67.9 37 24.3 ...##-----------------------------------------------------------------------------
## Gráfico.
ylab <- expression("Sterile spikelets (%)")
xlab <- expression("Row spacing (m)")
p1 <- xyplot(ester~esp0|la, data=gec,
groups=cult, pch=1:3,
ylab=ylab, xlab=xlab, layout=c(4,1),
## ylim=extendrange(gec$rend, f=0.1),
col=col,
strip=strip.custom(bg="gray90", factor.levels=la.lab),
key=list(
## x=0.55, y=0.95, corner=c(0,1),
x=0.025, y=0.9, corner=c(0,1),
lines=list(pch=1:3, lty=1:3, type="b", col=col, size=6),
text=list(lab=cult.lab, cex=0.8)))
p2 <- xyplot(fit~esp|la,
data=pred, groups=cult,
type="l", col=col, lty=1:3,
ly=pred$lwr, uy=pred$upr, cty="bands",
alpha=0.4,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)## x11(width=8, height=4)
p1+as.layer(p2, under=TRUE)
Número de panículas
##-----------------------------------------------------------------------------
## Ver número de panículas.
xyplot(npan~esp|la, groups=cult,
data=subset(gec, !is.na(npan)), as.table=TRUE,
auto.key=TRUE, jitter.x=TRUE,
## scales=list(y=list(log=10)),
type=c("p","a"),
par.settings=list(superpose.symbol=list(pch=19)),
xlab="Espaçamento-40 (m)",
ylab="Número de panículas")
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m0 <- lme(npan~(la*(esp+I(esp^2))*cult),
random=~1|la.bl/espac,
data=gec, na.action=na.omit,
method="ML")
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r))
## Embora esteja considerando só a parte de efeito fixo serve de
## indicação.
MASS::boxcox(lm(formula(m0), data=gec))
abline(v=0.5, col=2)
##-----------------------------------------------------------------------------
## Com transformação.
## Ajuste do modelo com a resposta transformada.
m0 <- update(m0, sqrt(npan)~.)
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r))
## Quadro de testes F de Wald para os termos de efeito fixo.
anova(m0)## numDF denDF F-value p-value
## (Intercept) 1 100 17931 <.0001
## la 3 12 226 <.0001
## esp 1 40 166 <.0001
## I(esp^2) 1 40 19 0.0001
## cult 2 100 19 <.0001
## la:esp 3 40 1 0.3291
## la:I(esp^2) 3 40 2 0.1334
## la:cult 6 100 10 <.0001
## esp:cult 2 100 1 0.4978
## I(esp^2):cult 2 100 0 0.8275
## la:esp:cult 6 100 1 0.7315
## la:I(esp^2):cult 6 100 2 0.1386## Quadro de anova com nomes modificados.
## Teste de Wald para os fatores fixos do modelo.
aux <- anova(m0)
rownames(aux) <- ef.names
aux## numDF denDF F-value p-value
## constante 1 100 17931 <.0001
## ambiente 3 12 226 <.0001
## espaçaamento linear 1 40 166 <.0001
## espaçamento quadrático 1 40 19 0.0001
## cultivar 2 100 19 <.0001
## amb:esp lin 3 40 1 0.3291
## amb:esp qua 3 40 2 0.1334
## amb:cult 6 100 10 <.0001
## esp lin:cult 2 100 1 0.4978
## esp quad:cult 2 100 0 0.8275
## amb:esp lin:cult 6 100 1 0.7315
## amb:esp qua:cult 6 100 2 0.1386## Estimativas dos parâmetros e medidas de ajuste.
summary(m0)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## 721.5 847.7 -321.7
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 4.694e-05
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 0.3351 1.3
##
## Fixed effects: sqrt(npan) ~ la + esp + I(esp^2) + cult + la:esp + la:I(esp^2) + la:cult + esp:cult + I(esp^2):cult + la:esp:cult + la:I(esp^2):cult
## Value Std.Error DF t-value p-value
## (Intercept) 21.29 0.60 100 35.63 0.0000
## la1.2009 -3.33 0.84 12 -3.94 0.0019
## la2.2009 -11.46 0.84 12 -13.56 0.0000
## la3.2009 -6.92 0.84 12 -8.19 0.0000
## esp -17.66 3.34 40 -5.29 0.0000
## I(esp^2) 9.58 37.32 40 0.26 0.7987
## cultD -4.59 0.82 100 -5.61 0.0000
## cultP -3.70 0.82 100 -4.53 0.0000
## la1.2009:esp 7.42 4.72 40 1.57 0.1236
## la2.2009:esp 1.49 4.72 40 0.32 0.7542
## la3.2009:esp 2.88 4.72 40 0.61 0.5451
## la1.2009:I(esp^2) -2.95 52.78 40 -0.06 0.9557
## la2.2009:I(esp^2) 78.49 52.78 40 1.49 0.1448
## la3.2009:I(esp^2) 44.03 52.78 40 0.83 0.4091
## la1.2009:cultD 0.62 1.18 100 0.52 0.6015
## la2.2009:cultD 4.70 1.16 100 4.06 0.0001
## la3.2009:cultD 6.30 1.16 100 5.45 0.0000
## la1.2009:cultP 1.28 1.16 100 1.11 0.2712
## la2.2009:cultP 3.12 1.16 100 2.69 0.0083
## la3.2009:cultP 4.53 1.16 100 3.92 0.0002
## esp:cultD 3.31 4.57 100 0.72 0.4713
## esp:cultP 5.62 4.57 100 1.23 0.2219
## I(esp^2):cultD 99.93 51.10 100 1.96 0.0533
## I(esp^2):cultP 10.03 51.10 100 0.20 0.8448
## la1.2009:esp:cultD -2.30 6.96 100 -0.33 0.7413
## la2.2009:esp:cultD 3.32 6.46 100 0.51 0.6083
## la3.2009:esp:cultD -2.74 6.46 100 -0.42 0.6722
## la1.2009:esp:cultP -8.06 6.46 100 -1.25 0.2156
## la2.2009:esp:cultP -0.59 6.46 100 -0.09 0.9271
## la3.2009:esp:cultP -8.74 6.46 100 -1.35 0.1794
## la1.2009:I(esp^2):cultD -29.76 75.89 100 -0.39 0.6958
## la2.2009:I(esp^2):cultD -127.19 72.27 100 -1.76 0.0815
## la3.2009:I(esp^2):cultD -172.70 72.27 100 -2.39 0.0187
## la1.2009:I(esp^2):cultP 15.51 72.27 100 0.21 0.8306
## la2.2009:I(esp^2):cultP 35.93 72.27 100 0.50 0.6202
## la3.2009:I(esp^2):cultP -42.93 72.27 100 -0.59 0.5538
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.65950 -0.51789 -0.02355 0.59178 2.12918
##
## Number of Observations: 188
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64tTable <- as.data.frame(summary(m0)$tTable)
## tTable[tTable[,5]<0.05,]
## fixef(m0)
## tTable
##-----------------------------------------------------------------------------
## Modelo com termos não significativos abandonados.
m1 <- update(m0, .~la*cult+esp+I(esp^2))
anova(m0, m1)## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 39 721.5 847.7 -321.7
## m1 2 17 705.6 760.7 -335.8 1 vs 2 28.15 0.1707anova(m1)## numDF denDF F-value p-value
## (Intercept) 1 116 17341 <.0001
## la 3 12 218 <.0001
## cult 2 116 19 <.0001
## esp 1 46 162 <.0001
## I(esp^2) 1 46 18 1e-04
## la:cult 6 116 10 <.0001summary(m1)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## 705.6 760.7 -335.8
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 5.635e-05
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 0.3899 1.394
##
## Fixed effects: sqrt(npan) ~ la + cult + esp + I(esp^2) + la:cult
## Value Std.Error DF t-value p-value
## (Intercept) 20.80 0.404 116 51.52 0.0000
## la1.2009 -3.37 0.532 12 -6.33 0.0000
## la2.2009 -10.47 0.532 12 -19.69 0.0000
## la3.2009 -6.37 0.532 12 -11.97 0.0000
## cultD -3.34 0.512 116 -6.52 0.0000
## cultP -3.58 0.512 116 -6.98 0.0000
## esp -13.47 1.053 46 -12.79 0.0000
## I(esp^2) 48.35 11.735 46 4.12 0.0002
## la1.2009:cultD 0.12 0.756 116 0.16 0.8702
## la2.2009:cultD 3.11 0.724 116 4.29 0.0000
## la3.2009:cultD 4.14 0.724 116 5.72 0.0000
## la1.2009:cultP 1.47 0.724 116 2.03 0.0442
## la2.2009:cultP 3.57 0.724 116 4.92 0.0000
## la3.2009:cultP 3.99 0.724 116 5.51 0.0000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.21492 -0.50518 0.08308 0.71095 2.05879
##
## Number of Observations: 188
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64##-----------------------------------------------------------------------------
## Predição.
pred <- expand.grid(la=levels(gec$la),
cult=levels(gec$cult),
esp=seq(0.25, 0.55, by=0.025)-0.4)
## pred <- predict.lme(m0, pred)
pred <- predict.lme(m1, pred)
pred$esp <- pred$esp+0.4
str(pred)## 'data.frame': 156 obs. of 6 variables:
## $ la : Factor w/ 4 levels "1.2010","1.2009",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ cult: Factor w/ 3 levels "C","D","P": 1 1 1 1 2 2 2 2 3 3 ...
## $ esp : num 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ fit : num 23.9 20.5 13.4 17.5 20.6 ...
## $ lwr : num 23.1 19.7 12.6 16.7 19.8 ...
## $ upr : num 24.7 21.3 14.2 18.3 21.4 ...##-----------------------------------------------------------------------------
## Gráfico.
## ylab <- expression("Panicle " *(number~m^{-2}))
ylab <- expression("Square root of panicle " *(number~m^{-2}))
xlab <- expression("Row spacing (m)")
p1 <- xyplot(sqrt(npan)~esp0|la, data=gec,
groups=cult, pch=1:3,
ylab=ylab, xlab=xlab, layout=c(4,1),
## ylim=extendrange(gec$rend, f=0.1),
col=col,
strip=strip.custom(bg="gray90", factor.levels=la.lab),
key=list(
## x=0.55, y=0.95, corner=c(0,1),
x=0.025, y=0.05, corner=c(0,0),
lines=list(pch=1:3, lty=1:3, type="b", col=col, size=6),
text=list(lab=cult.lab, cex=0.8)))
p2 <- xyplot(fit~esp|la,
data=pred, groups=cult,
type="l", col=col, lty=1:3,
ly=pred$lwr, uy=pred$upr, cty="bands",
alpha=0.4,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)## x11(width=8, height=4)
p1+as.layer(p2, under=TRUE)
Peso de 100 grãos
##-----------------------------------------------------------------------------
## Ver peso de 100 grãos.
xyplot(p100~esp|la, groups=cult,
data=subset(gec, !is.na(p100)), as.table=TRUE,
auto.key=TRUE, jitter.x=TRUE,
## scales=list(y=list(log=10)),
type=c("p","a"),
par.settings=list(superpose.symbol=list(pch=19)),
xlab="Espaçamento-40 (m)",
ylab="Número de panículas")
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m0 <- lme(p100~(la*(esp+I(esp^2))*cult),
random=~1|la.bl/espac,
data=gec, na.action=na.omit,
method="ML")
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r))
## Embora esteja considerando só a parte de efeito fixo serve de
## indicação.
## MASS::boxcox(lm(formula(m0), data=gec))
## abline(v=0.5, col=2)
## Quadro de testes F de Wald para os termos de efeito fixo.
anova(m0)## numDF denDF F-value p-value
## (Intercept) 1 76 5089 <.0001
## la 3 12 65 <.0001
## esp 1 40 0 0.8467
## I(esp^2) 1 40 0 0.7516
## cult 2 76 104 <.0001
## la:esp 3 40 0 0.9043
## la:I(esp^2) 3 40 1 0.2425
## la:cult 6 76 19 <.0001
## esp:cult 2 76 0 0.7518
## I(esp^2):cult 2 76 2 0.0969
## la:esp:cult 6 76 0 0.9379
## la:I(esp^2):cult 6 76 1 0.2096## Quadro de anova com nomes modificados.
## Teste de Wald para os fatores fixos do modelo.
aux <- anova(m0)
rownames(aux) <- ef.names
aux## numDF denDF F-value p-value
## constante 1 76 5089 <.0001
## ambiente 3 12 65 <.0001
## espaçaamento linear 1 40 0 0.8467
## espaçamento quadrático 1 40 0 0.7516
## cultivar 2 76 104 <.0001
## amb:esp lin 3 40 0 0.9043
## amb:esp qua 3 40 1 0.2425
## amb:cult 6 76 19 <.0001
## esp lin:cult 2 76 0 0.7518
## esp quad:cult 2 76 2 0.0969
## amb:esp lin:cult 6 76 0 0.9379
## amb:esp qua:cult 6 76 1 0.2096## Estimativas dos parâmetros e medidas de ajuste.
summary(m0)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## 116.9 237.8 -19.44
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 0.08839
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 0.08449 0.2509
##
## Fixed effects: p100 ~ (la * (esp + I(esp^2)) * cult)
## Value Std.Error DF t-value p-value
## (Intercept) 2.564 0.130 76 19.734 0.0000
## la1.2009 0.415 0.190 12 2.186 0.0494
## la2.2009 -0.506 0.207 12 -2.439 0.0312
## la3.2009 -0.296 0.189 12 -1.564 0.1438
## esp -0.125 0.670 40 -0.187 0.8529
## I(esp^2) -0.625 7.491 40 -0.083 0.9339
## cultD 0.495 0.161 76 3.082 0.0029
## cultP -0.459 0.161 76 -2.858 0.0055
## la1.2009:esp 0.875 1.087 40 0.805 0.4259
## la2.2009:esp 0.470 0.996 40 0.472 0.6394
## la3.2009:esp -0.001 0.951 40 -0.001 0.9988
## la1.2009:I(esp^2) -9.328 11.672 40 -0.799 0.4289
## la2.2009:I(esp^2) -19.992 11.611 40 -1.722 0.0928
## la3.2009:I(esp^2) -1.285 10.786 40 -0.119 0.9058
## la1.2009:cultD 0.031 0.243 76 0.128 0.8985
## la2.2009:cultD -0.848 0.252 76 -3.363 0.0012
## la3.2009:cultD 0.476 0.232 76 2.052 0.0436
## la1.2009:cultP 0.345 0.232 76 1.484 0.1420
## la2.2009:cultP 0.954 0.403 76 2.367 0.0205
## la3.2009:cultP 0.410 0.232 76 1.769 0.0809
## esp:cultD -0.425 0.898 76 -0.473 0.6374
## esp:cultP 0.300 0.898 76 0.334 0.7392
## I(esp^2):cultD 6.875 10.040 76 0.685 0.4956
## I(esp^2):cultP -1.250 10.040 76 -0.125 0.9012
## la1.2009:esp:cultD -1.170 1.428 76 -0.820 0.4149
## la2.2009:esp:cultD 0.280 1.335 76 0.210 0.8342
## la3.2009:esp:cultD 0.726 1.273 76 0.571 0.5698
## la1.2009:esp:cultP -1.393 1.438 76 -0.969 0.3356
## la2.2009:esp:cultP -0.034 1.516 76 -0.022 0.9822
## la3.2009:esp:cultP -0.124 1.273 76 -0.097 0.9229
## la1.2009:I(esp^2):cultD 30.506 15.609 76 1.954 0.0543
## la2.2009:I(esp^2):cultD 8.959 15.335 76 0.584 0.5608
## la3.2009:I(esp^2):cultD -10.590 14.343 76 -0.738 0.4626
## la1.2009:I(esp^2):cultP 20.339 15.356 76 1.324 0.1893
## la2.2009:I(esp^2):cultP -4.061 21.231 76 -0.191 0.8488
## la3.2009:I(esp^2):cultP 5.660 14.343 76 0.395 0.6942
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.29209 -0.41599 0.01943 0.42158 3.84639
##
## Number of Observations: 164
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64tTable <- as.data.frame(summary(m0)$tTable)
## tTable[tTable[,5]<0.05,]
## fixef(m0)
## tTable
##-----------------------------------------------------------------------------
## Modelo com termos não significativos abandonados.
m1 <- update(m0, .~la*cult)
anova(m0, m1)## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 39 116.88 237.8 -19.44
## m1 2 15 97.87 144.4 -33.94 1 vs 2 28.99 0.2204anova(m1)## numDF denDF F-value p-value
## (Intercept) 1 92 7163 <.0001
## la 3 12 89 <.0001
## cult 2 92 102 <.0001
## la:cult 6 92 17 <.0001summary(m1)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## 97.87 144.4 -33.94
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 0.05498
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 0.107 0.276
##
## Fixed effects: p100 ~ la + cult + la:cult
## Value Std.Error DF t-value p-value
## (Intercept) 2.5562 0.0820 92 31.174 0.0000
## la1.2009 0.3270 0.1240 12 2.638 0.0217
## la2.2009 -0.8031 0.1243 12 -6.459 0.0000
## la3.2009 -0.3160 0.1176 12 -2.687 0.0198
## cultD 0.5812 0.1013 92 5.735 0.0000
## cultP -0.4750 0.1013 92 -4.687 0.0000
## la1.2009:cultD 0.3910 0.1563 92 2.502 0.0141
## la2.2009:cultD -0.6896 0.1533 92 -4.499 0.0000
## la3.2009:cultD 0.3472 0.1447 92 2.400 0.0184
## la1.2009:cultP 0.5608 0.1528 92 3.670 0.0004
## la2.2009:cultP 0.7767 0.1880 92 4.132 0.0001
## la3.2009:cultP 0.4847 0.1447 92 3.351 0.0012
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.26909 -0.46353 0.04392 0.47047 3.21854
##
## Number of Observations: 164
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64##-----------------------------------------------------------------------------
## Predição.
pred <- expand.grid(la=levels(gec$la),
cult=levels(gec$cult),
esp=seq(0.25, 0.55, by=0.025)-0.4)
## pred <- predict.lme(m0, pred)
pred <- predict.lme(m1, pred)
pred$esp <- pred$esp+0.4
str(pred)## 'data.frame': 156 obs. of 6 variables:
## $ la : Factor w/ 4 levels "1.2010","1.2009",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ cult: Factor w/ 3 levels "C","D","P": 1 1 1 1 2 2 2 2 3 3 ...
## $ esp : num 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ fit : num 2.56 2.88 1.75 2.24 3.14 ...
## $ lwr : num 2.4 2.71 1.58 2.08 2.98 ...
## $ upr : num 2.71 3.06 1.93 2.4 3.29 ...##-----------------------------------------------------------------------------
## Gráfico.
## ylab <- expression("Panicle " *(number~m^{-2}))
ylab <- expression("100 grain weight (g)") # 100 unidades
xlab <- expression("Row spacing (m)")
p1 <- xyplot(p100~esp0|la, data=gec,
groups=cult, pch=1:3,
ylab=ylab, xlab=xlab, layout=c(4,1),
## ylim=extendrange(gec$rend, f=0.1),
col=col,
strip=strip.custom(bg="gray90", factor.levels=la.lab),
key=list(
## x=0.55, y=0.95, corner=c(0,1),
x=0.025, y=0.05, corner=c(0,0),
lines=list(pch=1:3, lty=1:3, type="b", col=col, size=6),
text=list(lab=cult.lab, cex=0.8)))
p2 <- xyplot(fit~esp|la,
data=pred, groups=cult,
type="l", col=col, lty=1:3,
ly=pred$lwr, uy=pred$upr, cty="bands",
alpha=0.4,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)## x11(width=8, height=4)
p1+as.layer(p2, under=TRUE)
Rendimento de número de graõs por panícula
##-----------------------------------------------------------------------------
## Ver redimento de número de grãos.
xyplot(ng_pan~esp|la, groups=cult,
data=subset(gec, !is.na(p100)), as.table=TRUE,
auto.key=TRUE, jitter.x=TRUE,
## scales=list(y=list(log=10)),
type=c("p","a"),
par.settings=list(superpose.symbol=list(pch=19)),
xlab="Espaçamento-40 (m)",
ylab="Rendimento de número de grãos")
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m0 <- lme(ng_pan~(la*(esp+I(esp^2))*cult),
random=~1|la.bl/espac,
data=gec, na.action=na.omit,
method="ML")
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r))
## Embora esteja considerando só a parte de efeito fixo serve de
## indicação.
MASS::boxcox(lm(formula(m0), data=gec))
abline(v=0.25, col=2)
##-----------------------------------------------------------------------------
## Com transformação.
## Ajuste do modelo com a resposta transformada.
m0 <- update(m0, ng_pan^(1/4)~.)
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r))
## Quadro de testes F de Wald para os termos de efeito fixo.
anova(m0)## numDF denDF F-value p-value
## (Intercept) 1 97 61449 <.0001
## la 3 12 1844 <.0001
## esp 1 40 18 0.0001
## I(esp^2) 1 40 6 0.0175
## cult 2 97 12 <.0001
## la:esp 3 40 8 0.0002
## la:I(esp^2) 3 40 0 0.9520
## la:cult 6 97 21 <.0001
## esp:cult 2 97 1 0.4402
## I(esp^2):cult 2 97 0 0.9920
## la:esp:cult 6 97 2 0.0991
## la:I(esp^2):cult 6 97 1 0.4137## Quadro de anova com nomes modificados.
## Teste de Wald para os fatores fixos do modelo.
aux <- anova(m0)
rownames(aux) <- ef.names
aux## numDF denDF F-value p-value
## constante 1 97 61449 <.0001
## ambiente 3 12 1844 <.0001
## espaçaamento linear 1 40 18 0.0001
## espaçamento quadrático 1 40 6 0.0175
## cultivar 2 97 12 <.0001
## amb:esp lin 3 40 8 0.0002
## amb:esp qua 3 40 0 0.9520
## amb:cult 6 97 21 <.0001
## esp lin:cult 2 97 1 0.4402
## esp quad:cult 2 97 0 0.9920
## amb:esp lin:cult 6 97 2 0.0991
## amb:esp qua:cult 6 97 1 0.4137## Estimativas dos parâmetros e medidas de ajuste.
summary(m0)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## -127.9 -2.321 103
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 0.006975
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 3.164e-06 0.1385
##
## Fixed effects: ng_pan^(1/4) ~ la + esp + I(esp^2) + cult + la:esp + la:I(esp^2) + la:cult + esp:cult + I(esp^2):cult + la:esp:cult + la:I(esp^2):cult
## Value Std.Error DF t-value p-value
## (Intercept) 3.466 0.062 97 55.99 0.0000
## la1.2009 -0.256 0.088 12 -2.92 0.0128
## la2.2009 -1.853 0.088 12 -21.17 0.0000
## la3.2009 -0.095 0.088 12 -1.09 0.2975
## esp 0.288 0.345 40 0.83 0.4089
## I(esp^2) -1.438 3.859 40 -0.37 0.7114
## cultD -0.083 0.087 97 -0.95 0.3447
## cultP 0.369 0.087 97 4.23 0.0001
## la1.2009:esp 0.666 0.488 40 1.36 0.1803
## la2.2009:esp -0.350 0.488 40 -0.72 0.4771
## la3.2009:esp -0.359 0.488 40 -0.74 0.4663
## la1.2009:I(esp^2) -2.297 5.457 40 -0.42 0.6761
## la2.2009:I(esp^2) -0.084 5.457 40 -0.02 0.9878
## la3.2009:I(esp^2) -3.603 5.457 40 -0.66 0.5129
## la1.2009:cultD 0.037 0.124 97 0.30 0.7671
## la2.2009:cultD -0.028 0.124 97 -0.23 0.8215
## la3.2009:cultD 0.016 0.124 97 0.13 0.8955
## la1.2009:cultP 0.096 0.133 97 0.73 0.4690
## la2.2009:cultP -0.737 0.124 97 -5.97 0.0000
## la3.2009:cultP -0.496 0.124 97 -4.01 0.0001
## esp:cultD 1.166 0.488 97 2.39 0.0188
## esp:cultP 0.191 0.488 97 0.39 0.6971
## I(esp^2):cultD 1.180 5.457 97 0.22 0.8292
## I(esp^2):cultP -3.784 5.457 97 -0.69 0.4898
## la1.2009:esp:cultD -1.465 0.726 97 -2.02 0.0464
## la2.2009:esp:cultD -0.334 0.690 97 -0.48 0.6298
## la3.2009:esp:cultD -1.679 0.690 97 -2.43 0.0168
## la1.2009:esp:cultP 0.283 0.721 97 0.39 0.6951
## la2.2009:esp:cultP 0.295 0.690 97 0.43 0.6696
## la3.2009:esp:cultP -0.507 0.690 97 -0.73 0.4647
## la1.2009:I(esp^2):cultD 0.186 7.945 97 0.02 0.9814
## la2.2009:I(esp^2):cultD -7.023 7.718 97 -0.91 0.3651
## la3.2009:I(esp^2):cultD 1.259 7.718 97 0.16 0.8708
## la1.2009:I(esp^2):cultP -3.252 8.159 97 -0.40 0.6911
## la2.2009:I(esp^2):cultP 6.010 7.718 97 0.78 0.4380
## la3.2009:I(esp^2):cultP 9.078 7.718 97 1.18 0.2424
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.44398 -0.71408 0.01505 0.65317 2.74192
##
## Number of Observations: 185
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64tTable <- as.data.frame(summary(m0)$tTable)
## tTable[tTable[,5]<0.05,]
## fixef(m0)
## tTable
##-----------------------------------------------------------------------------
## Modelo com termos não significativos abandonados.
m1 <- update(m0, .~la*(cult+esp+I(esp^2)))
anova(m0, m1)## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 39 -127.9 -2.32 102.96
## m1 2 23 -137.9 -63.86 91.96 1 vs 2 21.99 0.1436anova(m1)## numDF denDF F-value p-value
## (Intercept) 1 113 62030 <.0001
## la 3 12 1864 <.0001
## cult 2 113 11 <.0001
## esp 1 40 17 0.0002
## I(esp^2) 1 40 7 0.0147
## la:cult 6 113 20 <.0001
## la:esp 3 40 8 0.0003
## la:I(esp^2) 3 40 0 0.8686summary(m1)## Linear mixed-effects model fit by maximum likelihood
## Data: gec
## AIC BIC logLik
## -137.9 -63.86 91.96
##
## Random effects:
## Formula: ~1 | la.bl
## (Intercept)
## StdDev: 1.135e-06
##
## Formula: ~1 | espac %in% la.bl
## (Intercept) Residual
## StdDev: 2.991e-06 0.1472
##
## Fixed effects: ng_pan^(1/4) ~ la + cult + esp + I(esp^2) + la:cult + la:esp + la:I(esp^2)
## Value Std.Error DF t-value p-value
## (Intercept) 3.477 0.048 113 72.35 0.0000
## la1.2009 -0.253 0.069 12 -3.67 0.0032
## la2.2009 -1.849 0.068 12 -27.21 0.0000
## la3.2009 -0.138 0.068 12 -2.04 0.0644
## cultD -0.068 0.055 113 -1.24 0.2184
## cultP 0.322 0.055 113 5.85 0.0000
## esp 0.740 0.201 40 3.68 0.0007
## I(esp^2) -2.306 2.250 40 -1.03 0.3115
## la1.2009:cultD 0.043 0.080 113 0.54 0.5930
## la2.2009:cultD -0.116 0.078 113 -1.49 0.1403
## la3.2009:cultD 0.032 0.078 113 0.41 0.6820
## la1.2009:cultP 0.048 0.082 113 0.59 0.5585
## la2.2009:cultP -0.662 0.078 113 -8.50 0.0000
## la3.2009:cultP -0.383 0.078 113 -4.91 0.0000
## la1.2009:esp 0.241 0.299 40 0.80 0.4256
## la2.2009:esp -0.363 0.285 40 -1.28 0.2093
## la3.2009:esp -1.088 0.285 40 -3.82 0.0005
## la1.2009:I(esp^2) -2.512 3.326 40 -0.76 0.4546
## la2.2009:I(esp^2) -0.422 3.181 40 -0.13 0.8951
## la3.2009:I(esp^2) -0.158 3.181 40 -0.05 0.9607
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.0311736 -0.6831712 0.0006591 0.5795065 2.8652428
##
## Number of Observations: 185
## Number of Groups:
## la.bl espac %in% la.bl
## 16 64##-----------------------------------------------------------------------------
## Predição.
pred <- expand.grid(la=levels(gec$la),
cult=levels(gec$cult),
esp=seq(0.25, 0.55, by=0.025)-0.4)
## pred <- predict.lme(m0, pred)
pred <- predict.lme(m1, pred)
pred$esp <- pred$esp+0.4
str(pred)## 'data.frame': 156 obs. of 6 variables:
## $ la : Factor w/ 4 levels "1.2010","1.2009",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ cult: Factor w/ 3 levels "C","D","P": 1 1 1 1 2 2 2 2 3 3 ...
## $ esp : num 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
## $ fit : num 3.31 2.97 1.51 3.33 3.25 ...
## $ lwr : num 3.21 2.86 1.41 3.23 3.15 ...
## $ upr : num 3.41 3.07 1.61 3.44 3.35 ...##-----------------------------------------------------------------------------
## Gráfico.
ylab <- expression("Grain number per panicle")
xlab <- expression("Row spacing (m)")
p1 <- xyplot(ng_pan^0.25~esp0|la, data=gec,
groups=cult, pch=1:3,
ylab=ylab, xlab=xlab, layout=c(4,1),
## ylim=extendrange(gec$rend, f=0.1),
col=col,
strip=strip.custom(bg="gray90", factor.levels=la.lab),
key=list(
## x=0.55, y=0.95, corner=c(0,1),
x=0.025, y=0.05, corner=c(0,0),
lines=list(pch=1:3, lty=1:3, type="b", col=col, size=6),
text=list(lab=cult.lab, cex=0.8)))
p2 <- xyplot(fit~esp|la,
data=pred, groups=cult,
type="l", col=col, lty=1:3,
ly=pred$lwr, uy=pred$upr, cty="bands",
alpha=0.4,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)## x11(width=8, height=4)
p1+as.layer(p2, under=TRUE)