Desenvolvimento ao longo do ciclo

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     = {Desenvolvimento ao longo do ciclo},
    doi       = {},
    journal   = {},
    month     = {},
    pages     = {},
    publisher = {},
    url       = {},
    year      = {2014}
}

Donwnload:

Definições da sessão

##-----------------------------------------------------------------------------
## 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", "nlme", "reshape", "plyr", "doBy",
         "multcomp", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=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); dev.off()
rm(list=c("colp", "mycol", "ps"))

##-----------------------------------------------------------------------------
## Para colocar as médias sob o centro dos segmentos nos gráficos.

panel.segplot.text <- function(z, centers, txt, ...){
    panel.segplot(z=z, centers=centers, ...)
    panel.text(y=as.integer(z), x=centers,
               label=txt, pos=3)
}

Importação dos dados

##-----------------------------------------------------------------------------
## Lendo arquivos de dados.

## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/abacaxiciclo.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", dec=",")
str(da)
## 'data.frame':    1152 obs. of  9 variables:
##  $ trat : int  11 12 13 14 21 22 23 24 31 32 ...
##  $ cober: Factor w/ 2 levels "com cobertura",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ gesso: Factor w/ 2 levels "com gesso","sem gesso": 2 2 2 2 2 2 2 2 2 2 ...
##  $ vari : Factor w/ 4 levels "HAVAI","IAC - FANTÁSTICO",..: 4 4 4 4 1 1 1 1 2 2 ...
##  $ bloc : int  1 2 3 4 1 2 3 4 1 2 ...
##  $ dia  : Factor w/ 18 levels "2012/02/28","2012/03/28",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ caule: num  41.2 34.8 34.9 35 33.6 ...
##  $ copa : num  24.6 18.8 21.2 24.2 18.4 13.2 10 12.8 17.2 15 ...
##  $ alt  : num  44.4 38.6 41.6 40.8 38.4 32.8 36 31.8 30.4 33 ...
## Converte variáveis para fator.
da <- transform(da,
                bloc=factor(bloc),
                gesso=factor(gesso, labels=c("4", "0")),
                cober=factor(cober, labels=c("Com cobertura",
                                        "Sem cobertura")),
                ue=interaction(cober, gesso, vari, bloc,
                    drop=TRUE, sep=":"))

## Passa dia para formato de data.
da$dia <- as.POSIXct(da$dia, format="%Y/%m/%d")

## Mostra a estrutura do objeto.
str(da)
## 'data.frame':    1152 obs. of  10 variables:
##  $ trat : int  11 12 13 14 21 22 23 24 31 32 ...
##  $ cober: Factor w/ 2 levels "Com cobertura",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ gesso: Factor w/ 2 levels "4","0": 2 2 2 2 2 2 2 2 2 2 ...
##  $ vari : Factor w/ 4 levels "HAVAI","IAC - FANTÁSTICO",..: 4 4 4 4 1 1 1 1 2 2 ...
##  $ bloc : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ dia  : POSIXct, format: "2012-02-28" "2012-02-28" ...
##  $ caule: num  41.2 34.8 34.9 35 33.6 ...
##  $ copa : num  24.6 18.8 21.2 24.2 18.4 13.2 10 12.8 17.2 15 ...
##  $ alt  : num  44.4 38.6 41.6 40.8 38.4 32.8 36 31.8 30.4 33 ...
##  $ ue   : Factor w/ 64 levels "com cobertura:com gesso:HAVAI:1",..: 16 32 48 64 4 20 36 52 8 24 ...

Análise exploratória

##-----------------------------------------------------------------------------
## Análise exploratória.

da$gescob <- with(da, interaction(gesso, cober))

p1 <- xyplot(caule~dia|gescob+vari,
             groups=bloc, data=da, type="b",
             ylab="Diâmetro do caule",
             xlab="Data no ciclo")
useOuterStrips(p1)

p1 <- xyplot(copa~dia|gescob+vari,
             groups=bloc, data=da, type="b",
             ylab="Diâmetro da copa",
             xlab="Data no ciclo")
useOuterStrips(p1)

p1 <- xyplot(alt~dia|gescob+vari,
             groups=bloc, data=da, type="b",
             ylab="Altura da planta",
             xlab="Data no ciclo")
useOuterStrips(p1)

## Número de curvas.
nlevels(da$ue)
## [1] 64

Análise exploratória sobre a última medida (fim do vegetativo)

##-----------------------------------------------------------------------------

db <- ddply(na.omit(da), .(vari), mutate, diaa=max(dia)-dia)
units(db$diaa) <- "days"
db$diaa <- as.numeric(db$diaa)
db <- subset(db, diaa==0)

p1 <- xyplot(caule~vari|gesso+cober,
             groups=bloc, data=db, type="b",
             ylab="Diametro do caule",
             xlab="Variedade")
useOuterStrips(p1)

p1 <- xyplot(copa~vari|gesso+cober,
             groups=bloc, data=db, type="b",
             ylab="Diâmetro da copa",
             xlab="Variedade")
useOuterStrips(p1)

p1 <- xyplot(alt~vari|gesso+cober,
             groups=bloc, data=db, type="b",
             ylab="Altura da planta",
             xlab="Variedade")
useOuterStrips(p1)

Diâmetro do colmo

##-----------------------------------------------------------------------------

m0 <- lm(caule~bloc+vari*cober*gesso, data=db)
par(mfrow=c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: caule
##                  Df  Sum Sq Mean Sq F value    Pr(>F)    
## bloc              3  226.84   75.61  2.9296   0.04369 *  
## vari              3 1053.35  351.12 13.6037 1.891e-06 ***
## cober             1   70.06   70.06  2.7143   0.10642    
## gesso             1   12.67   12.67  0.4910   0.48708    
## vari:cober        3   42.72   14.24  0.5518   0.64958    
## vari:gesso        3  105.13   35.04  1.3577   0.26781    
## cober:gesso       1    2.13    2.13  0.0826   0.77514    
## vari:cober:gesso  3   33.70   11.23  0.4353   0.72882    
## Residuals        45 1161.47   25.81                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
LSmeans(m0, effect="vari")
##   estimate       se df   t.stat      p.value             vari
## 1 64.65625 1.270102 45 50.90634 1.999032e-41            HAVAI
## 2 54.63125 1.270102 45 43.01328 3.374871e-38 IAC - FANTÁSTICO
## 3 55.54313 1.270102 45 43.73123 1.629831e-38         IMPERIAL
## 4 55.86188 1.270102 45 43.98220 1.267101e-38           PÉROLA
L <- LSmatrix(m0, effect="vari")
rownames(L) <- levels(da$vari)

grid <- apmc(model=m0, X=L, focus="vari", test="single-step")
grid$vari <- reorder(grid$vari, grid$estimate)
grid$txt <- with(grid, sprintf("%.2f %s", estimate, cld))
## dev.off()
## x11(w=7, h=3.5)
segplot(vari~lwr+upr, centers=estimate, data=grid, draw=FALSE,
        xlab=expression(Diâmetro~da~colmo~(mm)),
        ylab="Variedade",
        txt=grid$txt,
        panel=panel.segplot.text)

Diâmetro da copa

##-----------------------------------------------------------------------------

m0 <- lm(copa~bloc+vari*cober*gesso, data=db)
par(mfrow=c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: copa
##                  Df Sum Sq Mean Sq F value    Pr(>F)    
## bloc              3   77.6   25.85  0.3886 0.7617629    
## vari              3 6480.4 2160.15 32.4706 2.511e-11 ***
## cober             1  975.0  975.00 14.6559 0.0003962 ***
## gesso             1  546.4  546.39  8.2132 0.0063023 ** 
## vari:cober        3  172.7   57.58  0.8655 0.4659857    
## vari:gesso        3  611.3  203.78  3.0632 0.0375246 *  
## cober:gesso       1   54.4   54.39  0.8176 0.3707034    
## vari:cober:gesso  3   45.5   15.15  0.2278 0.8765824    
## Residuals        45 2993.7   66.53                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Efeito da cobertura.
LSmeans(m0, effect="cober")
##   estimate       se df   t.stat      p.value         cober
## 1 64.95000 1.441854 45 45.04617 4.424085e-39 Com cobertura
## 2 72.75625 1.441854 45 50.46020 2.950879e-41 Sem cobertura
L <- LSmatrix(m0, effect="cober")
rownames(L) <- levels(da$cober)

grid <- apmc(model=m0, X=L, focus="cober", test="single-step")
grid$txt <- with(grid, sprintf("%.2f %s", estimate, cld))
grid
##           cober estimate      lwr      upr cld     txt
## 1 Com cobertura 64.95000 62.04596 67.85404   b 64.95 b
## 2 Sem cobertura 72.75625 69.85221 75.66029   a 72.76 a
##-----------------------------------------------------------------------------
## Estudar o efeito de gesso dentro de cada variedade.

LSmeans(m0, effect=c("gesso", "vari"))
##   estimate       se df   t.stat      p.value gesso             vari
## 1   71.575 2.883708 45 24.82047 6.752746e-28     4            HAVAI
## 2   69.900 2.883708 45 24.23962 1.822820e-27     0            HAVAI
## 3   61.550 2.883708 45 21.34405 3.563067e-25     4 IAC - FANTÁSTICO
## 4   68.100 2.883708 45 23.61543 5.422779e-27     0 IAC - FANTÁSTICO
## 5   48.450 2.883708 45 16.80128 5.188941e-21     4         IMPERIAL
## 6   63.725 2.883708 45 22.09828 8.531203e-26     0         IMPERIAL
## 7   82.150 2.883708 45 28.48763 1.970303e-30     4           PÉROLA
## 8   85.375 2.883708 45 29.60598 3.795842e-31     0           PÉROLA
Xm <- LSmatrix(m0, effect=c("gesso", "vari"))
grid <- equallevels(attr(Xm, "grid"), db)

L <- by(Xm, INDICES=grid$vari, FUN=as.matrix)
L <- lapply(L, "rownames<-", levels(grid$gesso))
g2 <- lapply(L, apmc, model=m0, focus="gesso")
g2
## $HAVAI
##   gesso estimate      lwr      upr cld
## 1     4   71.575 65.76691 77.38309   a
## 2     0   69.900 64.09191 75.70809   a
## 
## $`IAC - FANTÁSTICO`
##   gesso estimate      lwr      upr cld
## 1     4    61.55 55.74191 67.35809   a
## 2     0    68.10 62.29191 73.90809   a
## 
## $IMPERIAL
##   gesso estimate      lwr      upr cld
## 1     4   48.450 42.64191 54.25809   b
## 2     0   63.725 57.91691 69.53309   a
## 
## $PÉROLA
##   gesso estimate      lwr      upr cld
## 1     4   82.150 76.34191 87.95809   a
## 2     0   85.375 79.56691 91.18309   a
g2 <- ldply(g2); names(g2)[1] <- "vari"
g2 <- equallevels(g2, db)
g2$txt <- with(g2, sprintf("%.2f %s", estimate, cld))
## segplot(vari~lwr+upr, centers=estimate, data=g2, draw=FALSE,
##         groups=g2$gesso, pch=g2$gesso, f=0.1, panel=panel.segplotBy,
##         xlab=expression(Diâmetro~da~copa~(cm)),
##         ylab="Variedade",
##         key=list(columns=2, type="o", divide=1,
##             title=expression(Gesso~(ton~ha^{-1})), cex.title=1.1,
##             lines=list(pch=1:2),
##             text=list(levels(g2$gesso))))

segplot(vari~lwr+upr, centers=estimate, data=g2, draw=FALSE,
        groups=g2$gesso, pch=g2$gesso, f=0.1, txt=g2$txt,
        panel=function (x, y, z, centers, subscripts, groups, f, txt, ...){
            d <- 2*((as.numeric(groups)-1)/(nlevels(groups)-1))-1
            z <- as.integer(z)+f*d
            panel.segplot(x, y, z, centers=centers,
                          subscripts=subscripts,
                          ...)
            panel.text(y=z, x=centers,
                       label=txt, pos=2+sign(d), cex=0.8)
        },
        xlab=expression(Diâmetro~da~copa~(cm)),
        ylab="Variedade",
        key=list(columns=2, type="o", divide=1,
            title=expression(Gesso~(ton~ha^{-1})), cex.title=1.1,
            lines=list(pch=1:2),
            text=list(levels(g2$gesso))))

Altura de plantas

##-----------------------------------------------------------------------------

m0 <- lm(alt~bloc+vari*cober*gesso, data=db)
par(mfrow=c(2,2)); plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: alt
##                  Df  Sum Sq Mean Sq F value    Pr(>F)    
## bloc              3  204.63   68.21  1.3448 0.2718197    
## vari              3 2659.06  886.35 17.4742 1.152e-07 ***
## cober             1  820.82  820.82 16.1822 0.0002173 ***
## gesso             1   79.21   79.21  1.5616 0.2178921    
## vari:cober        3  173.82   57.94  1.1423 0.3422829    
## vari:gesso        3  437.80  145.93  2.8770 0.0463948 *  
## cober:gesso       1    8.41    8.41  0.1658 0.6858021    
## vari:cober:gesso  3    7.85    2.62  0.0516 0.9843159    
## Residuals        45 2282.57   50.72                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Efeito da cobertura.
LSmeans(m0, effect="cober")
##   estimate       se df   t.stat      p.value         cober
## 1  54.0500 1.259014 45 42.93042 3.673354e-38 Com cobertura
## 2  61.2125 1.259014 45 48.61939 1.525529e-40 Sem cobertura
L <- LSmatrix(m0, effect="cober")
rownames(L) <- levels(da$cober)

grid <- apmc(model=m0, X=L, focus="cober", test="single-step")
grid$txt <- with(grid, sprintf("%.2f %s", estimate, cld))
grid
##           cober estimate      lwr      upr cld     txt
## 1 Com cobertura  54.0500 51.51422 56.58578   b 54.05 b
## 2 Sem cobertura  61.2125 58.67672 63.74828   a 61.21 a
##-----------------------------------------------------------------------------
## Estudar o efeito de gesso dentro de cada variedade.

LSmeans(m0, effect=c("gesso", "vari"))
##   estimate       se df   t.stat      p.value gesso             vari
## 1   67.050 2.518028 45 26.62798 3.479250e-29     4            HAVAI
## 2   64.225 2.518028 45 25.50607 2.145688e-28     0            HAVAI
## 3   56.475 2.518028 45 22.42826 4.623173e-26     4 IAC - FANTÁSTICO
## 4   56.900 2.518028 45 22.59705 3.389334e-26     0 IAC - FANTÁSTICO
## 5   42.425 2.518028 45 16.84850 4.651534e-21     4         IMPERIAL
## 6   53.425 2.518028 45 21.21700 4.551735e-25     0         IMPERIAL
## 7   60.125 2.518028 45 23.87781 3.419014e-27     4           PÉROLA
## 8   60.425 2.518028 45 23.99695 2.776917e-27     0           PÉROLA
Xm <- LSmatrix(m0, effect=c("gesso", "vari"))
grid <- equallevels(attr(Xm, "grid"), db)

L <- by(Xm, INDICES=grid$vari, FUN=as.matrix)
L <- lapply(L, "rownames<-", levels(grid$gesso))
g2 <- lapply(L, apmc, model=m0, focus="gesso")
g2
## $HAVAI
##   gesso estimate      lwr      upr cld
## 1     4   67.050 61.97843 72.12157   a
## 2     0   64.225 59.15343 69.29657   a
## 
## $`IAC - FANTÁSTICO`
##   gesso estimate      lwr      upr cld
## 1     4   56.475 51.40343 61.54657   a
## 2     0   56.900 51.82843 61.97157   a
## 
## $IMPERIAL
##   gesso estimate      lwr      upr cld
## 1     4   42.425 37.35343 47.49657   b
## 2     0   53.425 48.35343 58.49657   a
## 
## $PÉROLA
##   gesso estimate      lwr      upr cld
## 1     4   60.125 55.05343 65.19657   a
## 2     0   60.425 55.35343 65.49657   a
g2 <- ldply(g2); names(g2)[1] <- "vari"
g2 <- equallevels(g2, db)
g2$txt <- with(g2, sprintf("%.2f %s", estimate, cld))
## segplot(vari~lwr+upr, centers=estimate, data=g2, draw=FALSE,
##         groups=g2$gesso, pch=g2$gesso, f=0.1, panel=panel.segplotBy,
##         xlab=expression(Altura~de~planta~(cm)),
##         ylab="Variedade",
##         key=list(columns=2, type="o", divide=1,
##             title=expression(Gesso~(ton~ha^{-1})), cex.title=1.1,
##             lines=list(pch=1:2),
##             text=list(levels(g2$gesso))))

segplot(vari~lwr+upr, centers=estimate, data=g2, draw=FALSE,
        groups=g2$gesso, pch=g2$gesso, f=0.1, txt=g2$txt,
        panel=function (x, y, z, centers, subscripts, groups, f, txt, ...){
            d <- 2*((as.numeric(groups)-1)/(nlevels(groups)-1))-1
            z <- as.integer(z)+f*d
            panel.segplot(x, y, z, centers=centers,
                          subscripts=subscripts,
                          ...)
            panel.text(y=z, x=centers,
                       label=txt, pos=2+sign(d), cex=0.8)
        },
        xlab=expression(Altura~de~planta~(cm)),
        ylab="Variedade",
        key=list(columns=2, type="o", divide=1,
            title=expression(Gesso~(ton~ha^{-1})), cex.title=1.1,
            lines=list(pch=1:2),
            text=list(levels(g2$gesso))))

Versões

##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.

sessionInfo()
## R version 3.2.0 (2015-04-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
## 
## 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] methods   stats     graphics  grDevices utils     datasets  base     
## 
## other attached packages:
##  [1] wzRfun_0.5          multcomp_1.3-8      TH.data_1.0-4       mvtnorm_1.0-1      
##  [5] doBy_4.5-12         survival_2.38-1     plyr_1.8.1          reshape_0.8.5      
##  [9] nlme_3.1-120        latticeExtra_0.6-26 RColorBrewer_1.0-5  lattice_0.20-31    
## [13] rmarkdown_0.3.10    knitr_1.8          
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.11.3     splines_3.2.0   MASS_7.3-39     stringr_0.6.2   tools_3.2.0    
##  [6] grid_3.2.0      htmltools_0.2.6 yaml_2.1.13     digest_0.6.4    Matrix_1.2-0   
## [11] formatR_1.0     evaluate_0.5.5  sandwich_2.3-2  zoo_1.7-11
Sys.time()
## [1] "2015-05-24 17:02:13 BRT"