Walmes Zeviani
Milson Evaldo Serafim
Material complementar ao artigo com título ??? submetido à revista ???, volume ???, número ???. A citação será incluída assim que for determinado o volume e paginação.
%% Citação LaTex.
@article{SerafimAbacaxi2014,
author = {},
title = {Crescimento vegetativo inicial de abacaxizeiro em
função da cultura de cobertura e aplicação de gesso},
doi = {},
journal = {},
month = {},
pages = {},
publisher = {},
url = {},
year = {2014}
}
Donwnload:
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
## Instruções para instalação do wzRfun em:
## browseURL("https://github.com/walmes/wzRfun")
pkg <- c("lattice", "latticeExtra",
"reshape", "plyr", "doBy", "multcomp",
"wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
## lattice latticeExtra reshape plyr doBy multcomp
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUE
##-----------------------------------------------------------------------------
## Configurações de cores da lattice.
colp <- brewer.pal(11, "Spectral")
colp <- colorRampPalette(colp, space="rgb")
mycol <- brewer.pal(8, "Dark2")
ps <- list(
box.dot=list(col="black", pch="|"),
box.rectangle=list(col="black", fill=c("gray70")),
box.umbrella=list(col="black", lty=1),
dot.symbol=list(col="black"),
dot.line=list(col="gray90", lty=1),
plot.symbol=list(col="black", cex=0.8),
plot.line=list(col="black", lty=1.2),
plot.polygon=list(col="gray80"),
superpose.line=list(col=mycol),
superpose.symbol=list(col=mycol),
superpose.polygon=list(col=mycol),
regions=list(col=colp(100)),
strip.background=list(col=c("gray90","gray50")),
strip.shingle=list(col=c("gray50","gray90"))
)
trellis.par.set(ps)
rm(list=c("colp", "mycol", "ps"))
##-----------------------------------------------------------------------------
## Funções pára essa sessão.
Tukey.contrasts <- function(pm, levels=NULL, complement=NULL){
comp <- combn(1:nrow(pm), 2)
M <- apply(comp, 2, function(i) pm[i[1],]-pm[i[2],])
if(!is.null(levels)){
lab <- t(outer(levels, levels, paste, sep="-"))
lab <- c(lab[lower.tri(lab)])
if(!is.null(complement)) lab <- paste(complement, lab)
colnames(M) <- lab
}
t(M)
}
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
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 stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.1 multcomp_1.3-3 TH.data_1.0-3 mvtnorm_0.9-9997
## [5] doBy_4.5-10 MASS_7.3-33 survival_2.37-7 plyr_1.8.1
## [9] reshape_0.8.5 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 lme4_1.1-6
## [5] Matrix_1.1-4 methods_3.1.1 minqa_1.2.3 nlme_3.1-117
## [9] Rcpp_0.11.2 RcppEigen_0.3.2.0.2 sandwich_2.3-0 stringr_0.6.2
## [13] tools_3.1.1 zoo_1.7-10
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/abacaxiplanta.txt"
path <- ifelse(Sys.info()["user"]=="walmes",
yes=basename(url), no=url)
## Importa os dados a partir do arquivo na web.
da <- read.table(file=path, header=TRUE, sep="\t")
## Trunca os nomes para 4 digitos.
names(da) <- substr(names(da), 1, 4)
## Converte variáveis para fator.
da <- transform(da, cobe=factor(cobe), bloc=factor(bloc), gess=factor(gess))
## Mostra a estrutura do objeto.
str(da)
## 'data.frame': 400 obs. of 9 variables:
## $ epoc: int 0 0 0 0 0 0 0 0 0 0 ...
## $ cult: Factor w/ 4 levels "IAC Fantástico",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ gess: Factor w/ 2 levels "0","4": NA NA NA NA NA NA NA NA NA NA ...
## $ cobe: Factor w/ 3 levels "","milheto","sem": 1 1 1 1 1 1 1 1 1 1 ...
## $ bloc: Factor w/ 4 levels "1","2","3","4": NA NA NA NA NA NA NA NA NA NA ...
## $ plan: int 1 2 3 4 5 1 2 3 4 5 ...
## $ altu: num 35 29 30.4 31 31 ...
## $ dcau: num 33.4 32.5 28.2 36.5 36.5 ...
## $ dros: int NA NA NA NA NA NA NA NA NA NA ...
## Ordenar as linhas.
da <- arrange(da, epoc, cult, gess, cobe, bloc, plan)
##-----------------------------------------------------------------------------
## Separar os valores conforme a época.
db <- dlply(da, .(epoc), .fun=droplevels)
str(db)
## List of 2
## $ 0 :'data.frame': 80 obs. of 9 variables:
## ..$ epoc: int [1:80] 0 0 0 0 0 0 0 0 0 0 ...
## ..$ cult: Factor w/ 4 levels "IAC Fantástico",..: 1 1 1 1 1 1 1 1 1 1 ...
## ..$ gess: Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
## ..$ cobe: Factor w/ 1 level "": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ bloc: Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
## ..$ plan: int [1:80] 1 1 1 1 2 2 2 2 3 3 ...
## ..$ altu: num [1:80] 14.5 14.5 15.8 15.2 15.8 ...
## ..$ dcau: num [1:80] 26.3 27.2 24 26.9 26.1 ...
## ..$ dros: int [1:80] NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "vars")= chr "epoc"
## $ 90:'data.frame': 320 obs. of 9 variables:
## ..$ epoc: int [1:320] 90 90 90 90 90 90 90 90 90 90 ...
## ..$ cult: Factor w/ 4 levels "IAC Fantástico",..: 1 1 1 1 1 1 1 1 1 1 ...
## ..$ gess: Factor w/ 2 levels "0","4": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ cobe: Factor w/ 2 levels "milheto","sem": 1 1 1 1 1 1 1 1 1 1 ...
## ..$ bloc: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 2 2 2 2 2 ...
## ..$ plan: int [1:320] 1 2 3 4 5 1 2 3 4 5 ...
## ..$ altu: num [1:320] 42 32 34 41 34 41 32 32 29 32 ...
## ..$ dcau: num [1:320] 30.8 25.5 31.4 32.4 27.6 ...
## ..$ dros: int [1:320] 39 41 40 50 28 47 48 32 34 38 ...
## ..- attr(*, "vars")= chr "epoc"
## - attr(*, "split_type")= chr "data.frame"
## - attr(*, "split_labels")='data.frame': 2 obs. of 1 variable:
## ..$ epoc: int [1:2] 0 90
##-----------------------------------------------------------------------------
## Obter o crescimento relativo em altura e diâmentro de cada cultivar
## no período de 90 dias. O crescimento relativo é dividir o valor das
## observações aos 90 dias pelo seu valor esperado no dia 0. O valor
## experado no dia zero, uma vez que os fatores experimentais ainda não
## foram aplicados, é função apenas da cultivar e do bloco. Como o
## experimento é balanceado, o valor médio não depende dos blocos, só
## das cultivares.
dc <- transform(
db[["90"]],
altur=(altu-predict(lm(altu~cult, db[["0"]]),
newdata=db[["90"]]))/90,
dcaur=(dcau-predict(lm(dcau~cult, db[["0"]]),
newdata=db[["90"]]))/90)
str(dc)
## 'data.frame': 320 obs. of 11 variables:
## $ epoc : int 90 90 90 90 90 90 90 90 90 90 ...
## $ cult : Factor w/ 4 levels "IAC Fantástico",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gess : Factor w/ 2 levels "0","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ cobe : Factor w/ 2 levels "milheto","sem": 1 1 1 1 1 1 1 1 1 1 ...
## $ bloc : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 2 2 2 2 2 ...
## $ plan : int 1 2 3 4 5 1 2 3 4 5 ...
## $ altu : num 42 32 34 41 34 41 32 32 29 32 ...
## $ dcau : num 30.8 25.5 31.4 32.4 27.6 ...
## $ dros : int 39 41 40 50 28 47 48 32 34 38 ...
## $ altur: num 0.298 0.187 0.209 0.287 0.209 ...
## $ dcaur: num 0.05023 -0.00799 0.05734 0.06845 0.01434 ...
rm("db")
##-----------------------------------------------------------------------------
## Análise exploratória.
p1 <-
xyplot(altu~gess|cult, groups=cobe, type=c("p","a"),
data=dc, jitter.x=TRUE, layout=c(NA,1))
## xyplot(altur~gess|cult, groups=cobe, type=c("p","a"),
## data=dc, jitter.x=TRUE, layout=c(NA,1))
p2 <-
xyplot(dcau~gess|cult, groups=cobe, type=c("p","a"),
data=dc, jitter.x=TRUE, layout=c(NA,1))
## xyplot(dcaur~gess|cult, groups=cobe, type=c("p","a"),
## data=dc, jitter.x=TRUE, layout=c(NA,1))
p3 <-
xyplot(dros~gess|cult, groups=cobe, type=c("p","a"),
data=dc, jitter.x=TRUE, layout=c(NA,1))
plot(p1, split=c(1,1,1,3), more=TRUE)
plot(p2, split=c(1,2,1,3), more=TRUE)
plot(p3, split=c(1,3,1,3), more=FALSE)
rm(list=paste0("p", 1:3))
##-----------------------------------------------------------------------------
## Estimação dos efeitos para o modelo correspondente ao fatorial tripo
## em blocos.
m0 <- lm(altur~bloc+cult*gess*cobe, data=dc)
##-----------------------------------------------------------------------------
## Diagnóstico para afastamento dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Sem qualquer afastamento dos pressupostos.
##-----------------------------------------------------------------------------
## Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: altur
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 3 0.031 0.0105 1.11 0.34459
## cult 3 0.189 0.0631 6.70 0.00022 ***
## gess 1 0.002 0.0023 0.24 0.62247
## cobe 1 0.051 0.0510 5.42 0.02061 *
## cult:gess 3 0.069 0.0228 2.42 0.06578 .
## cult:cobe 3 0.093 0.0311 3.30 0.02067 *
## gess:cobe 1 0.043 0.0429 4.56 0.03354 *
## cult:gess:cobe 3 0.022 0.0073 0.78 0.50842
## Residuals 301 2.834 0.0094
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Sem efeito de interação tripla. Duplas apenas.
##-----------------------------------------------------------------------------
## Valor do CV (que não tem utilidade alguma mas os colaboradores e
## revisores exigem). Neste o valor alto é devido a resposta ser uma
## fração, uma variável relativa. O valor alto de fato pouco importa e
## nada implica na análise haja visto que os pressupostos não apresentar
## nenhum afastamento.
100*mean(dc$altur)/sd(residuals(m0))
## [1] 183
##-----------------------------------------------------------------------------
## Abandono da interação tripla e a dupla não significativa.
m1 <- lm(altur~bloc+(cult+gess)*cobe, data=dc)
## Verifica se os termos abandoandos têm efeito nulo.
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: altur ~ bloc + (cult + gess) * cobe
## Model 2: altur ~ bloc + cult * gess * cobe
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 307 2.92
## 2 301 2.83 6 0.0904 1.6 0.15
## Quadro de anova para os termos restantes.
anova(m1)
## Analysis of Variance Table
##
## Response: altur
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 3 0.031 0.0105 1.10 0.34988
## cult 3 0.189 0.0631 6.62 0.00024 ***
## gess 1 0.002 0.0023 0.24 0.62449
## cobe 1 0.051 0.0510 5.35 0.02134 *
## cult:cobe 3 0.093 0.0311 3.26 0.02172 *
## gess:cobe 1 0.043 0.0429 4.51 0.03456 *
## Residuals 307 2.925 0.0095
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Predição dos valores médios para o modelo escolhido com intervalos de
## confiança.
pred <- with(dc,
expand.grid(
bloc="1",
cult=levels(cult),
gess=levels(gess),
cobe=levels(cobe)))
pred <- cbind(pred, predict(m1, newdata=pred, interval="confidence"))
str(pred)
## 'data.frame': 16 obs. of 7 variables:
## $ bloc: Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## $ cult: Factor w/ 4 levels "IAC Fantástico",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ gess: Factor w/ 2 levels "0","4": 1 1 1 1 2 2 2 2 1 1 ...
## $ cobe: Factor w/ 2 levels "milheto","sem": 1 1 1 1 1 1 1 1 2 2 ...
## $ fit : num 0.222 0.173 0.257 0.17 0.204 ...
## $ lwr : num 0.183 0.134 0.219 0.131 0.166 ...
## $ upr : num 0.261 0.212 0.296 0.209 0.243 ...
##-----------------------------------------------------------------------------
## Gráfico de intervalos.
segplot(cobe:gess~lwr+upr|cult, data=pred,
draw.bands=FALSE, centers=fit, layout=c(4,1),
xlab=expression(Cobertura:Gesso~(t~ha^{-1})),
ylab=expression(Taxa~de~crescimento~em~altura~(cm~dia^{-1})),
segments.fun=panel.arrows, ends="both",
angle=90, length=1, unit="mm", horizontal=FALSE,
scales=list(x=list(rot=90)),
panel=function(...){
panel.abline(v=1:8, lty=2, col="gray90", h=seq(25,55,5))
panel.segplot(...)
})
## Cultivar não interage com gesso e por isso, dentro de um mesmo nível
## de cobertura a mesma diferença em resposta se observa ao mudar o
## nível de gesso para as cultivares. Com cobertura de milheto, adição
## de gesso tende a diminuir enquanto que sem cobertura tende a aumentar
## (é significativo?). As cultivares respondem de forma diferenciada à
## cobertura onde para milheto tem-se Pérola é superior mas sem
## cobertura é igual a Smooth C.
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## Comparar níveis de cobertura para cada cultivar.
## Comparar níveis de gesso para cada cobertura.
## Comparar cultivares não faz sentido devido as diferenças morfológicas
## serem genéticas.
pm <- list()
for(i in levels(dc$cobe))
pm[[i]] <- LSmatrix(m1, effect=c("gess","cobe"), at=list(cobe=i))
for(i in levels(dc$cult))
pm[[i]] <- LSmatrix(m1, effect=c("cobe","cult"), at=list(cult=i))
K <- list()
K[[1]] <- Tukey.contrasts(pm[[1]], levels=levels(dc$gess),
complement="(conven)")
K[[2]] <- Tukey.contrasts(pm[[2]], levels=levels(dc$gess),
complement="(milhet)")
K[[3]] <- Tukey.contrasts(pm[[3]], levels=levels(dc$cobe),
complement="(Fantas)")
K[[4]] <- Tukey.contrasts(pm[[4]], levels=levels(dc$cobe),
complement="(Imperi)")
K[[5]] <- Tukey.contrasts(pm[[5]], levels=levels(dc$cobe),
complement="(Pérola)")
K[[6]] <- Tukey.contrasts(pm[[6]], levels=levels(dc$cobe),
complement="(Smooth)")
KK <- do.call(rbind, K)
##-----------------------------------------------------------------------------
## Teste para os constranstes construídos.
summary(glht(m1, linfct=KK))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = altur ~ bloc + (cult + gess) * cobe, data = dc)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## (conven) 0-4 == 0 0.0178 0.0154 1.15 0.8191
## (milhet) 0-4 == 0 -0.0285 0.0154 -1.85 0.3332
## (Fantas) milheto-sem == 0 0.0211 0.0218 0.97 0.9118
## (Imperi) milheto-sem == 0 -0.0200 0.0218 -0.92 0.9305
## (Pérola) milheto-sem == 0 0.0762 0.0218 3.49 0.0033 **
## (Smooth) milheto-sem == 0 0.0237 0.0218 1.08 0.8581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##-----------------------------------------------------------------------------
## Estimação dos efeitos para o modelo correspondente ao fatorial tripo
## em blocos.
m0 <- lm(dcaur~bloc+cult*gess*cobe, data=dc)
##-----------------------------------------------------------------------------
## Diagnóstico para afastamento dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
##-----------------------------------------------------------------------------
## Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: dcaur
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 3 0.008 0.0026 0.45 0.718
## cult 3 0.459 0.1529 26.06 5.2e-15 ***
## gess 1 0.001 0.0012 0.21 0.649
## cobe 1 0.009 0.0087 1.49 0.223
## cult:gess 3 0.063 0.0210 3.58 0.014 *
## cult:cobe 3 0.045 0.0149 2.54 0.057 .
## gess:cobe 1 0.013 0.0134 2.29 0.131
## cult:gess:cobe 3 0.006 0.0021 0.36 0.782
## Residuals 301 1.766 0.0059
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Valor do CV.
100*mean(dc$dcaur)/sd(residuals(m0))
## [1] 118.4
##-----------------------------------------------------------------------------
## Abandono da interação tripla e a dupla não significativa.
m1 <- lm(dcaur~bloc+cult*gess, data=dc)
## Verifica se os termos abandoandos têm efeito nulo.
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: dcaur ~ bloc + cult * gess
## Model 2: dcaur ~ bloc + cult * gess * cobe
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 309 1.84
## 2 301 1.77 8 0.0732 1.56 0.14
## Quadro de anova para os termos restantes.
anova(m1)
## Analysis of Variance Table
##
## Response: dcaur
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 3 0.008 0.0026 0.44 0.723
## cult 3 0.459 0.1529 25.68 7.3e-15 ***
## gess 1 0.001 0.0012 0.21 0.651
## cult:gess 3 0.063 0.0210 3.53 0.015 *
## Residuals 309 1.839 0.0060
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Predição dos valores médios para o modelo escolhido com intervalos de
## confiança.
pred <- with(dc,
expand.grid(
bloc="1",
cult=levels(cult),
gess=levels(gess)))
pred <- cbind(pred, predict(m1, newdata=pred, interval="confidence"))
str(pred)
## 'data.frame': 8 obs. of 6 variables:
## $ bloc: Factor w/ 1 level "1": 1 1 1 1 1 1 1 1
## $ cult: Factor w/ 4 levels "IAC Fantástico",..: 1 2 3 4 1 2 3 4
## $ gess: Factor w/ 2 levels "0","4": 1 1 1 1 2 2 2 2
## $ fit : num 0.0253 0.0984 0.1709 0.0886 0.0542 ...
## $ lwr : num -0.00284 0.07025 0.14273 0.06045 0.02609 ...
## $ upr : num 0.0535 0.1265 0.199 0.1167 0.0824 ...
##-----------------------------------------------------------------------------
## Gráfico de intervalos.
segplot(gess~lwr+upr|cult, data=pred,
draw.bands=FALSE, centers=fit, layout=c(4,1),
xlab=expression(Gesso~(t~ha^{-1})),
ylab=expression(Taxa~de~crescimento~em~diâmetro~(mm~dia^{-1})),
segments.fun=panel.arrows, ends="both",
angle=90, length=1, unit="mm", horizontal=FALSE,
scales=list(x=list(rot=0)),
panel=function(...){
panel.abline(v=1:2, lty=2, col="gray90", h=seq(3.3,3.9,0.1))
panel.segplot(...)
})
##-----------------------------------------------------------------------------
## Comparações múltiplas.
pm <- list()
for(i in levels(dc$cult))
pm[[i]] <- LSmatrix(m1, effect=c("gess","cult"), at=list(cult=i))
K <- list()
K[[1]] <- Tukey.contrasts(pm[[1]], levels=levels(dc$gess),
complement="(Fantas)")
K[[2]] <- Tukey.contrasts(pm[[2]], levels=levels(dc$gess),
complement="(Imperi)")
K[[3]] <- Tukey.contrasts(pm[[3]], levels=levels(dc$gess),
complement="(Pérola)")
K[[4]] <- Tukey.contrasts(pm[[4]], levels=levels(dc$gess),
complement="(Smooth)")
KK <- do.call(rbind, K)
##-----------------------------------------------------------------------------
## Teste para os constranstes construídos.
summary(glht(m1, linfct=KK))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = dcaur ~ bloc + cult * gess, data = dc)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## (Fantas) 0-4 == 0 -0.02893 0.01725 -1.68 0.327
## (Imperi) 0-4 == 0 0.00120 0.01725 0.07 1.000
## (Pérola) 0-4 == 0 0.04848 0.01725 2.81 0.021 *
## (Smooth) 0-4 == 0 -0.00512 0.01725 -0.30 0.997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##-----------------------------------------------------------------------------
## Estimação dos efeitos para o modelo correspondente ao fatorial tripo
## em blocos.
m0 <- lm(dros~bloc+cult*gess*cobe, data=dc)
##-----------------------------------------------------------------------------
## Diagnóstico para afastamento dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)
##-----------------------------------------------------------------------------
## Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: dros
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 3 507 169 1.09 0.352
## cult 3 4305 1435 9.29 6.9e-06 ***
## gess 1 170 170 1.10 0.296
## cobe 1 142 142 0.92 0.339
## cult:gess 3 1200 400 2.59 0.053 .
## cult:cobe 3 303 101 0.65 0.581
## gess:cobe 1 306 306 1.98 0.160
## cult:gess:cobe 3 282 94 0.61 0.611
## Residuals 301 46515 155
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Valor do CV.
100*mean(dc$dros)/sd(residuals(m0))
## [1] 329.7
##-----------------------------------------------------------------------------
## Abandono da interação tripla e a dupla não significativa.
m1 <- lm(dros~bloc+cult, data=dc)
## Verifica se os termos abandoandos têm efeito nulo.
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: dros ~ bloc + cult
## Model 2: dros ~ bloc + cult * gess * cobe
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 313 48917
## 2 301 46515 12 2402 1.3 0.22
## Quadro de anova para os termos restantes.
anova(m1)
## Analysis of Variance Table
##
## Response: dros
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 3 507 169 1.08 0.36
## cult 3 4305 1435 9.18 7.7e-06 ***
## Residuals 313 48917 156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Predição dos valores médios para o modelo escolhido com intervalos de
## confiança.
pred <- with(dc,
expand.grid(
bloc="1",
cult=levels(cult)))
pred <- cbind(pred, predict(m1, newdata=pred, interval="confidence"))
str(pred)
## 'data.frame': 4 obs. of 5 variables:
## $ bloc: Factor w/ 1 level "1": 1 1 1 1
## $ cult: Factor w/ 4 levels "IAC Fantástico",..: 1 2 3 4
## $ fit : num 41.4 34.8 39 44.9
## $ lwr : num 37.8 31.2 35.3 41.3
## $ upr : num 45 38.5 42.6 48.6
##-----------------------------------------------------------------------------
## Gráfico de intervalos.
segplot(cult~lwr+upr, data=pred,
draw.bands=FALSE, centers=fit,
ylab=expression(Cultivar),
xlab=expression(Diâmetro~da~roseta~(mm)),
segments.fun=panel.arrows, ends="both",
angle=90, length=1, unit="mm", horizontal=TRUE,
scales=list(x=list(rot=0)),
panel=function(...){
panel.abline(h=1:4, lty=2, col="gray90", v=seq(3.3,3.9,0.1))
panel.segplot(...)
})
##-----------------------------------------------------------------------------
## Comparações múltiplas.
pm <- LSmatrix(m1, effect=c("cult"))
K <- Tukey.contrasts(pm, levels=levels(dc$cult))
KK <- K
##-----------------------------------------------------------------------------
## Teste para os constranstes construídos.
summary(glht(m1, linfct=KK))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = dros ~ bloc + cult, data = dc)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## IAC Fantástico-Imperial == 0 6.55 1.98 3.31 0.0057 **
## IAC Fantástico-Pérola == 0 2.44 1.98 1.23 0.6062
## IAC Fantástico-Smooth Cayenne == 0 -3.52 1.98 -1.78 0.2833
## Imperial-Pérola == 0 -4.11 1.98 -2.08 0.1617
## Imperial-Smooth Cayenne == 0 -10.07 1.98 -5.10 <0.001 ***
## Pérola-Smooth Cayenne == 0 -5.96 1.98 -3.02 0.0147 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)