Análise de experimento

Componentes de produção de soja em consórcio com gramíneas

Luís Armando Zago Machado
Walmes Marques Zeviani

Definições da sessão

##-----------------------------------------------------------------------------
## Pacotes necessários.

pkg <- c("latticeExtra", "plyr", "doBy", "multcomp",
         "agricolae", "wzRfun")
sapply(pkg, require, character.only=TRUE)

## Remove objetos da mémoria para começar vazia.
rm(list=ls())

##-----------------------------------------------------------------------------
## Definições gráficas.

## Tábua de cores.
mycol <- brewer.pal(6, "Set1")

## Definições de preenchimentos, linhas, etc.
ps <- list(box.rectangle=list(col=1, fill=c("gray70")),
           box.umbrella=list(col=1, lty=1),
           dot.symbol=list(col=1, pch=19),
           dot.line=list(col="gray50", lty=3),
           plot.symbol=list(col=1, cex=1.2),
           plot.line=list(col=1),
           plot.polygon=list(col="gray95"),
           superpose.line=list(col=mycol),
           superpose.symbol=list(col=mycol, pch=c(1:5)),
           superpose.polygon=list(col=mycol),
           strip.background=list(col=c("gray80","gray50"))
           )
trellis.par.set(ps)
## show.settings()

Carregar os dados

Download dos dados: http://www.leg.ufpr.br/~walmes/data/
Download do script: script.R

##-----------------------------------------------------------------------------
## Ler.

dados <- read.table("1-consRendGraosSoja12-13.csv",
                    header=TRUE, sep=";", dec=",")

dados$Ano <- factor(dados$Ano, labels=c(2012,2013))
dados$Rep  <- factor(dados$Rep)
dados$Trat <- factor(dados$Trat,
                     levels=c("SS", "SA", "ST", "SPa",
                         "SPi", "SX", "SD", "SR"))
str(dados)
## 'data.frame':    112 obs. of  11 variables:
##  $ Ano    : Factor w/ 2 levels "2012","2013": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Trat   : Factor w/ 8 levels "SS","SA","ST",..: 1 1 1 1 1 1 1 2 2 2 ...
##  $ Rep    : Factor w/ 7 levels "1","2","3","4",..: 1 2 3 4 5 6 7 1 2 3 ...
##  $ RG     : num  1501 1541 1977 1702 2058 ...
##  $ Imp    : num  1.35 1.02 2.1 0.88 0.65 0.73 0.52 1.28 2.01 0.87 ...
##  $ TUmid  : num  9 9.4 9.6 9.4 9.7 9.8 9.5 9.6 9.7 10.1 ...
##  $ AltSoja: num  96.7 91.3 96.7 100 114.7 ...
##  $ NPD    : num  0 0.83 2.92 0.42 0 1.25 0 0 1.25 0.42 ...
##  $ NVagem : num  45.7 39.4 40.7 37.5 44.9 ...
##  $ NGVagem: num  2.13 1.92 1.67 1.67 1.53 1.71 1.79 2.25 2.01 2.02 ...
##  $ M100G  : num  10 10.6 9.8 10.2 10 ...
xtabs(~Ano+Trat, data=dados)
##       Trat
## Ano    SS SA ST SPa SPi SX SD SR
##   2012  7  7  7   7   7  7  7  7
##   2013  7  7  7   7   7  7  7  7

Análise exploratória preliminar

Legenda dos fatores e níveis:

Legenda das variáveis resposta:

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

xyplot(RG+Imp+TUmid+AltSoja+NPD+NVagem+NGVagem+M100G~Trat,
       groups=Ano, data=dados, outer=TRUE, as.table=TRUE,
       scales=list(y="free", x=list(rot=90)), jitter.x=TRUE,
       xlab="Consócios", ylab="Valores observados")

plot of chunk unnamed-chunk-4

Análise

Experimento instalado no delineamento de blocos casualizados, com 7 repetições, avaliado em dois anos agrícolas consecutivos.

Foi considerado o modelo linear gaussiano com efeito fixo de ano, efeito aleatório de bloco dentro de ano, efeito fixo de consórcio e de ano interação com consórcio.

\[ y_{ijk} = \mu+\alpha_i+b_{i(j)}+\tau_k+\gamma_{ik}+e_{i(jk)}, \]

em que \(y_{ijk}\) são os valores observados, \(\mu\) é uma constante desconhecida e arbitrária e correspondente à todas as observações, \(\alpha_i\) é o efeito do ano, \(b_{i(j)}\) é o efeito do bloco dentro do ano, \(\tau_k\) é o efeito do nível de consórcio, \(\gamma_{ik}\) é o efeito da interação entre ano e consórcio e \(e_{i(jk)}\) é o erro resídual. Letras gregas indicam os termos de efeito fixo e letras arábes os termos de efeito aleatório. Aos termos de efeito aleatório assume-se distribuição normal com média 0 e variância \(\sigma^2_b\) e \(\sigma^2_e\) respectivamente.

O modelo teve seus parâmetros estimados por mínimos quadrados. A avaliação dos pressupostos do modelo (normalidade dos erros e variância constante) se deu por inspeção gráfica baseada nos resíduos e valores ajustados. Ao se verificar afastamento dos pressupostos, buscou-se aplicar transformações da família Box-Cos aos dados sugeridas pelo valor de \(\lambda\) que maximizam a log-verossimilhança perfilhada. Para inferência sobre diferença entre médias, considerou-se o teste de Tukey para comparações múltiplas. O estudo foi feito em separado por quando verificado que a interação ano e consórcio foi significativa.

RG: rendimento de grãos (kg/ha);

##-----------------------------------------------------------------------------
## Ver.

lab <- expression("Rendimento de grãos"~(kg~ha^{-1}))

xyplot(RG~Trat, groups=Ano, data=dados,
       as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
       auto.key=list(columns=2),
       xlab="Consórcios", ylab=lab)

plot of chunk unnamed-chunk-5

##-----------------------------------------------------------------------------
## Ajuste.

m0 <- aov(RG~Ano+Error(Ano:Rep)+Ano*Trat, data=dados)
summary(m0)
## 
## Error: Ano:Rep
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Ano        1  405403  405403    1.17    0.3
## Residuals 12 4163961  346997               
## 
## Error: Within
##           Df  Sum Sq Mean Sq F value Pr(>F)  
## Trat       7 1306920  186703    1.97  0.069 .
## Ano:Trat   7 1029923  147132    1.55  0.161  
## Residuals 84 7964676   94818                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Para os termos Trat e Ano:Trat, são os mesmos F.
m0 <- lm(RG~Ano+Ano:Rep+Ano*Trat, data=dados)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

plot of chunk unnamed-chunk-6

## Sem afastamentos dos pressupostos que sejam sistemáticos
## (corrigíveis) ou prejudiciais à inferência.

##-----------------------------------------------------------------------------
## Médias ajustadas.

glr <- df.residual(m0)
qmr <- deviance(m0)/glr

## Teste de Tukey.
tk <- with(dados,
           HSD.test(y=RG, trt=Trat, DFerror=glr, MSerror=qmr))
tk$statistics["HSD"] ## DMS.
##     HSD
##   361.8
tk$groups
##   trt means  M
## 1 ST   2060  a
## 2 SS   1991 ab
## 3 SPa  1957 ab
## 4 SX   1955 ab
## 5 SA   1943 ab
## 6 SD   1837 ab
## 7 SPi  1823 ab
## 8 SR   1695  b
## Remove os espaços extras adicionados pela função.
tk$groups[,-2] <- sapply(tk$groups[,-2],
                         FUN=gsub, pattern="\\s*$", replacement="")

## Para obter os IC para a média (ajustada).
## LSmeans(m0, effect="Trat")
L <- LSmatrix(m0, effect="Trat")
rownames(L) <- levels(dados$Trat)

## grid <- apmc(X=L, model=m0, focus="Trat", test="fdr", level=0.05)
grid <- apmc(X=L, model=m0, focus="Trat",
             test="single-step", level=0.05)
## grid

grid <- merge(grid, tk$groups, by.x="Trat", by.y="trt")
grid$Trat <- reorder(grid$Trat, grid$estimate)
grid
##   Trat estimate  lwr  upr cld means  M
## 1   SA     1943 1780 2107  ab  1943 ab
## 2   SD     1837 1674 2001  ab  1837 ab
## 3  SPa     1957 1793 2120  ab  1957 ab
## 4  SPi     1823 1660 1987  ab  1823 ab
## 5   SR     1695 1531 1859   b  1695  b
## 6   SS     1991 1827 2154  ab  1991 ab
## 7   ST     2060 1897 2224   a  2060  a
## 8   SX     1955 1791 2119  ab  1955 ab
## cld veio da glht() e M veio da HSD.test().
## estimate é uma média ajustada de mínimos quadrados, means é média
## amostral.

Apesar do teste F ter sido não significativo, a aplicação do teste de Tukey indicou diferenças entre as médias extremas.

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

l <- levels(grid$Trat)
ylim <- c(0.5, length(l)+1)
sc <- list(y=list(at=seq_along(l), labels=l))

segplot(Trat~lwr+upr, data=grid, centers=estimate, draw=FALSE,
        ## ann=grid$cld,
        ann=grid$M, ## Teste de Tukey.
        ylim=ylim, scales=sc,
        xlab=lab, ylab="Consórcios",
        panel=function(z, centers, ann, subscripts, ...){
            panel.segplot(z=z, centers=centers,
                          subscripts=subscripts,
                          ...)
            ann <- sprintf("%0.1f %s", centers, ann)
            panel.text(x=centers, y=z, labels=ann, pos=3, cex=0.8)
        })

plot of chunk unnamed-chunk-7

Imp: teor de impurezas (%);

##-----------------------------------------------------------------------------
## Ver.

lab <- expression("Teor de impurezas"~("%"))

xyplot(Imp~Trat, groups=Ano, data=dados,
       as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
       auto.key=list(columns=2),
       xlab="Consórcios", ylab=lab)

plot of chunk unnamed-chunk-8

Efeito de ano bem expressivo, inclusive com diferença de variabilidade.

##-----------------------------------------------------------------------------
## Ajuste.

## Para os termos Trat e Ano:Trat, são os mesmos F.
m0 <- lm(Imp~Ano+Ano:Rep+Ano*Trat, data=dados)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0, which=1:3)
MASS::boxcox(m0)
abline(v=1/3, col=2)

plot of chunk unnamed-chunk-9

layout(1)

## Aplicando transformação.
m0 <- update(m0, formula=.^(1/3)~.)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

plot of chunk unnamed-chunk-9

## Sem afastamentos dos pressupostos que sejam sistemáticos
## (corrigíveis) ou prejudiciais à inferência.

## ANOVA com valores de F corretos.
summary(aov(Imp^(1/3)~Ano+Error(Ano:Rep)+Ano*Trat, data=dados))
## 
## Error: Ano:Rep
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## Ano        1  28.17   28.17     303 7.1e-10 ***
## Residuals 12   1.12    0.09                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Trat       7  0.396  0.0566    2.03   0.06 .
## Ano:Trat   7  0.183  0.0261    0.94   0.48  
## Residuals 84  2.340  0.0279                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Médias ajustadas.

glr <- df.residual(m0)
qmr <- deviance(m0)/glr

## Teste de Tukey.
tk <- HSD.test(y=m0$model[,1], trt=m0$model["Trat"],
               DFerror=glr, MSerror=qmr)
tk$statistics["HSD"] ## DMS.
##      HSD
##   0.1961
tk$groups
##   trt means  M
## 1 SR  1.597  a
## 2 SPi 1.511 ab
## 3 ST  1.467 ab
## 4 SD  1.466 ab
## 5 SA  1.456 ab
## 6 SPa 1.427 ab
## 7 SX  1.412 ab
## 8 SS  1.399  b
## Remove os espaços extras adicionados pela função.
tk$groups[,-2] <- sapply(tk$groups[,-2],
                         FUN=gsub, pattern="\\s*$", replacement="")

## Para obter os IC para a média (ajustada).
## LSmeans(m0, effect="Trat")
L <- LSmatrix(m0, effect="Trat")
rownames(L) <- levels(dados$Trat)

## grid <- apmc(X=L, model=m0, focus="Trat", test="fdr", level=0.05)
grid <- apmc(X=L, model=m0, focus="Trat",
             test="single-step", level=0.05)
## grid

grid <- merge(grid, tk$groups, by.x="Trat", by.y="trt")
grid$Trat <- reorder(grid$Trat, grid$estimate)
grid
##   Trat estimate   lwr   upr cld means  M
## 1   SA    1.456 1.368 1.545  ab 1.456 ab
## 2   SD    1.466 1.378 1.555  ab 1.466 ab
## 3  SPa    1.427 1.339 1.516  ab 1.427 ab
## 4  SPi    1.511 1.422 1.600  ab 1.511 ab
## 5   SR    1.597 1.509 1.686   a 1.597  a
## 6   SS    1.399 1.311 1.488   b 1.399  b
## 7   ST    1.467 1.378 1.556  ab 1.467 ab
## 8   SX    1.412 1.323 1.500  ab 1.412 ab
## Representar na escala dos dados.
grid[,c("estimate","lwr","upr")] <- grid[,c("estimate","lwr","upr")]^3

## cld veio da glht() e M veio da HSD.test().
## estimate é uma média ajustada de mínimos quadrados, means é média
## amostral.

Apesar do teste F ter sido não significativo, a aplicação do teste de Tukey indicou diferenças entre as médias extremas. No gráfico estão sendo representados os valores médios na escala original. Eles foram obtidos pela transformação das médias ajustadas considerando a função inversa usada na resposta dados para adequação dos pressupostos. Em função disso, os intervalos na escala transformada não apresentam a mesma amplitude e nem simetria.

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

l <- levels(grid$Trat)
ylim <- c(0.5, length(l)+1)
sc <- list(y=list(at=seq_along(l), labels=l))

segplot(Trat~lwr+upr, data=grid, centers=estimate, draw=FALSE,
        ## ann=grid$cld,
        ann=grid$M, ## Teste de Tukey.
        ylim=ylim, scales=sc,
        xlab=lab, ylab="Consórcios",
        panel=function(z, centers, ann, subscripts, ...){
            panel.segplot(z=z, centers=centers,
                          subscripts=subscripts,
                          ...)
            ann <- sprintf("%0.2f %s", centers, ann)
            panel.text(x=centers, y=z, labels=ann, pos=3, cex=0.8)
        })

plot of chunk unnamed-chunk-10

TUmid: teor de umidade (%);

##-----------------------------------------------------------------------------
## Ver.

lab <- expression("Teor de umidade"~("%"))

xyplot(TUmid~Trat, groups=Ano, data=dados,
       as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
       auto.key=list(columns=2),
       xlab="Consórcios", ylab=lab)

plot of chunk unnamed-chunk-11

Efeito de ano bem expressivo, inclusive com diferença de variabilidade.

##-----------------------------------------------------------------------------
## Ajuste.

## Para os termos Trat e Ano:Trat, são os mesmos F.
m0 <- lm(TUmid~Ano+Ano:Rep+Ano*Trat, data=dados)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0, which=1:3)
MASS::boxcox(m0, lambda=seq(-3, 1, by=0.1))
abline(v=-2, col=2)

plot of chunk unnamed-chunk-12

layout(1)

## Aplicando transformação.
m0 <- update(m0, formula=.^(-2)~.)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

plot of chunk unnamed-chunk-12

## Sem afastamentos dos pressupostos que sejam sistemáticos
## (corrigíveis) ou prejudiciais à inferência.

## ANOVA com valores de F corretos.
summary(aov(TUmid^(-2)~Ano+Error(Ano:Rep)+Ano*Trat, data=dados))
## 
## Error: Ano:Rep
##           Df   Sum Sq  Mean Sq F value  Pr(>F)    
## Ano        1 0.001076 0.001076     813 2.1e-12 ***
## Residuals 12 0.000016 0.000001                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df   Sum Sq  Mean Sq F value  Pr(>F)    
## Trat       7 9.16e-06 1.31e-06    4.84 0.00013 ***
## Ano:Trat   7 1.95e-06 2.79e-07    1.03 0.41486    
## Residuals 84 2.27e-05 2.71e-07                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Médias ajustadas.

glr <- df.residual(m0)
qmr <- deviance(m0)/glr

## Teste de Tukey.
tk <- HSD.test(y=m0$model[,1], trt=m0$model["Trat"],
               DFerror=glr, MSerror=qmr)
tk$statistics["HSD"] ## DMS.
##         HSD
##   0.0006111
tk$groups
##   trt    means  M
## 1 SS  0.008035  a
## 2 SX  0.007524 ab
## 3 SPa 0.007516 ab
## 4 SD  0.007459 ab
## 5 SA  0.007354  b
## 6 SPi 0.007281  b
## 7 ST  0.007090  b
## 8 SR  0.007071  b
## Remove os espaços extras adicionados pela função.
tk$groups[,-2] <- sapply(tk$groups[,-2],
                         FUN=gsub, pattern="\\s*$", replacement="")

## Para obter os IC para a média (ajustada).
## LSmeans(m0, effect="Trat")
L <- LSmatrix(m0, effect="Trat")
rownames(L) <- levels(dados$Trat)

## grid <- apmc(X=L, model=m0, focus="Trat", test="fdr", level=0.05)
grid <- apmc(X=L, model=m0, focus="Trat",
             test="single-step", level=0.05)
## grid

grid <- merge(grid, tk$groups, by.x="Trat", by.y="trt")
grid$Trat <- reorder(grid$Trat, -grid$estimate)
grid
##   Trat estimate      lwr      upr cld    means  M
## 1   SA 0.007354 0.007078 0.007630   b 0.007354  b
## 2   SD 0.007459 0.007183 0.007736  ab 0.007459 ab
## 3  SPa 0.007516 0.007239 0.007792  ab 0.007516 ab
## 4  SPi 0.007281 0.007004 0.007557   b 0.007281  b
## 5   SR 0.007071 0.006795 0.007348   b 0.007071  b
## 6   SS 0.008035 0.007759 0.008312   a 0.008035  a
## 7   ST 0.007090 0.006813 0.007366   b 0.007090  b
## 8   SX 0.007524 0.007248 0.007801  ab 0.007524 ab
## Representar na escala dos dados.
grid[,c("estimate","lwr","upr")] <-
    grid[,c("estimate","lwr","upr")]^(-1/2)

## cld veio da glht() e M veio da HSD.test().
## estimate é uma média ajustada de mínimos quadrados, means é média
## amostral.

No gráfico estão sendo representados os valores médios na escala original. Eles foram obtidos pela transformação das médias ajustadas considerando a função inversa usada na resposta dados para adequação dos pressupostos. Em função disso, os intervalos na escala transformada não apresentam a mesma amplitude e nem simetria.

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

l <- levels(grid$Trat)
ylim <- c(0.5, length(l)+1)
sc <- list(y=list(at=seq_along(l), labels=l))

segplot(Trat~lwr+upr, data=grid, centers=estimate, draw=FALSE,
        ## ann=grid$cld,
        ann=grid$M, ## Teste de Tukey.
        ylim=ylim, scales=sc,
        xlab=lab, ylab="Consórcios",
        panel=function(z, centers, ann, subscripts, ...){
            panel.segplot(z=z, centers=centers,
                          subscripts=subscripts,
                          ...)
            ann <- sprintf("%0.2f %s", centers, ann)
            panel.text(x=centers, y=z, labels=ann, pos=3, cex=0.8)
        })

plot of chunk unnamed-chunk-13

AltSoja: Altura das plantas de soja (cm);

##-----------------------------------------------------------------------------
## Ver.

lab <- expression("Altura das plantas de soja"~("cm"))

xyplot(AltSoja~Trat, groups=Ano, data=dados,
       as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
       auto.key=list(columns=2),
       xlab="Consórcios", ylab=lab)

plot of chunk unnamed-chunk-14

##-----------------------------------------------------------------------------
## Ajuste.

## Para os termos Trat e Ano:Trat, são os mesmos F.
m0 <- lm(AltSoja~Ano+Ano:Rep+Ano*Trat, data=dados)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

plot of chunk unnamed-chunk-15

## Sem afastamentos dos pressupostos que sejam sistemáticos
## (corrigíveis) ou prejudiciais à inferência.

## ANOVA com valores de F corretos.
summary(aov(AltSoja~Ano+Error(Ano:Rep)+Ano*Trat, data=dados))
## 
## Error: Ano:Rep
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## Ano        1   6107    6107    25.6 0.00028 ***
## Residuals 12   2868     239                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Trat       7     64    9.17    0.44   0.88
## Ano:Trat   7     77   11.04    0.53   0.81
## Residuals 84   1757   20.92
##-----------------------------------------------------------------------------
## Médias ajustadas.

glr <- df.residual(m0)
qmr <- deviance(m0)/glr

## Teste de Tukey.
tk <- HSD.test(y=m0$model[,1], trt=m0$model["Trat"],
               DFerror=glr, MSerror=qmr)
tk$statistics["HSD"] ## DMS.
##     HSD
##   5.373
tk$groups
##   trt means M
## 1 SPa 94.87 a
## 2 SPi 94.75 a
## 3 SS  94.69 a
## 4 SD  94.35 a
## 5 SA  93.99 a
## 6 SX  93.29 a
## 7 SR  93.26 a
## 8 ST  92.73 a
## Remove os espaços extras adicionados pela função.
tk$groups[,-2] <- sapply(tk$groups[,-2],
                         FUN=gsub, pattern="\\s*$", replacement="")

## Para obter os IC para a média (ajustada).
## LSmeans(m0, effect="Trat")
L <- LSmatrix(m0, effect="Trat")
rownames(L) <- levels(dados$Trat)

## grid <- apmc(X=L, model=m0, focus="Trat", test="fdr", level=0.05)
grid <- apmc(X=L, model=m0, focus="Trat",
             test="single-step", level=0.05)
## grid

grid <- merge(grid, tk$groups, by.x="Trat", by.y="trt")
grid$Trat <- reorder(grid$Trat, grid$estimate)
grid
##   Trat estimate   lwr   upr cld means M
## 1   SA    93.99 91.56 96.42   a 93.99 a
## 2   SD    94.35 91.92 96.78   a 94.35 a
## 3  SPa    94.87 92.44 97.30   a 94.87 a
## 4  SPi    94.75 92.32 97.18   a 94.75 a
## 5   SR    93.26 90.83 95.69   a 93.26 a
## 6   SS    94.69 92.26 97.12   a 94.69 a
## 7   ST    92.73 90.30 95.16   a 92.73 a
## 8   SX    93.29 90.86 95.72   a 93.29 a
## cld veio da glht() e M veio da HSD.test().
## estimate é uma média ajustada de mínimos quadrados, means é média
## amostral.
##-----------------------------------------------------------------------------

l <- levels(grid$Trat)
ylim <- c(0.5, length(l)+1)
sc <- list(y=list(at=seq_along(l), labels=l))

segplot(Trat~lwr+upr, data=grid, centers=estimate, draw=FALSE,
        ## ann=grid$cld,
        ann=grid$M, ## Teste de Tukey.
        ylim=ylim, scales=sc,
        xlab=lab, ylab="Consórcios",
        panel=function(z, centers, ann, subscripts, ...){
            panel.segplot(z=z, centers=centers,
                          subscripts=subscripts,
                          ...)
            ann <- sprintf("%0.2f %s", centers, ann)
            panel.text(x=centers, y=z, labels=ann, pos=3, cex=0.8)
        })

plot of chunk unnamed-chunk-16

NPD: número de plantas daninhas (plantas/m²);

##-----------------------------------------------------------------------------
## Ver.

lab <- expression("Número de plantas daninhas"~(plantas~m^{-2}))

xyplot(NPD~Trat, groups=Ano, data=dados,
       as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
       auto.key=list(columns=2),
       xlab="Consórcios", ylab=lab)

plot of chunk unnamed-chunk-17

## xyplot(log(NPD+0.5)~Trat, groups=Ano, data=dados,
##        as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
##        auto.key=list(columns=2),
##        xlab="Consórcios", ylab=lab)

Efeito de ano bem expressivo, inclusive com diferença de variabilidade. Ainda é importante mencionar que essa é uma variável proveniente de uma contagem, como média de contagens, portanto não se pode aplicar à media de contagens a distribuição Poisson, a menos que o número de unidades (plantas) avaliadas para obter tal média, seja informado. Apresenta valores baixos incluindo zeros, portanto não contínua como se supõe para uso da análise de variância. Por outro lado, certas transformações, se adequadamente aplicadas, podem permitir que os pressupostos sejam atendidos, não comprometendo a validade das inferências.

##-----------------------------------------------------------------------------
## Ajuste.

summary(dados$NPD)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.42    1.25    4.58    5.68   32.70
## Presença de zeros. Somar um para poder se fazer log, caso necessário.

## Para os termos Trat e Ano:Trat, são os mesmos F.
m0 <- lm((NPD+0.5)~Ano+Ano:Rep+Ano*Trat, data=dados)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0, which=1:3)
MASS::boxcox(m0, lambda=seq(-2, 2, by=0.1))
abline(v=0, col=2)

plot of chunk unnamed-chunk-18

layout(1)

## Aplicando transformação.
m0 <- update(m0, formula=log(.)~.)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

plot of chunk unnamed-chunk-18

## Sem afastamentos dos pressupostos que sejam sistemáticos
## (corrigíveis) ou prejudiciais à inferência.

## ANOVA com valores de F corretos.
summary(aov(log(NPD+0.5)~Ano+Error(Ano:Rep)+Ano*Trat, data=dados))
## 
## Error: Ano:Rep
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## Ano        1    107   107.2     117 1.5e-07 ***
## Residuals 12     11     0.9                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Trat       7    5.6   0.794    1.39   0.22
## Ano:Trat   7    2.7   0.385    0.67   0.69
## Residuals 84   48.0   0.572
##-----------------------------------------------------------------------------
## Médias ajustadas.

glr <- df.residual(m0)
qmr <- deviance(m0)/glr

## Teste de Tukey.
tk <- HSD.test(y=m0$model[,1], trt=m0$model["Trat"],
               DFerror=glr, MSerror=qmr)
tk$statistics["HSD"] ## DMS.
##      HSD
##   0.8886
tk$groups
##   trt  means M
## 1 SX  1.1369 a
## 2 SPi 1.1052 a
## 3 SS  1.0524 a
## 4 SPa 0.9343 a
## 5 ST  0.8153 a
## 6 SD  0.6713 a
## 7 SA  0.5818 a
## 8 SR  0.5394 a
## Remove os espaços extras adicionados pela função.
tk$groups[,-2] <- sapply(tk$groups[,-2],
                         FUN=gsub, pattern="\\s*$", replacement="")

## Para obter os IC para a média (ajustada).
## LSmeans(m0, effect="Trat")
L <- LSmatrix(m0, effect="Trat")
rownames(L) <- levels(dados$Trat)

## grid <- apmc(X=L, model=m0, focus="Trat", test="fdr", level=0.05)
grid <- apmc(X=L, model=m0, focus="Trat",
             test="single-step", level=0.05)
## grid

grid <- merge(grid, tk$groups, by.x="Trat", by.y="trt")
grid$Trat <- reorder(grid$Trat, grid$estimate)
grid
##   Trat estimate    lwr    upr cld  means M
## 1   SA   0.5818 0.1798 0.9837   a 0.5818 a
## 2   SD   0.6713 0.2693 1.0732   a 0.6713 a
## 3  SPa   0.9343 0.5323 1.3362   a 0.9343 a
## 4  SPi   1.1052 0.7033 1.5072   a 1.1052 a
## 5   SR   0.5394 0.1374 0.9413   a 0.5394 a
## 6   SS   1.0524 0.6505 1.4544   a 1.0524 a
## 7   ST   0.8153 0.4133 1.2173   a 0.8153 a
## 8   SX   1.1369 0.7350 1.5389   a 1.1369 a
## Representar na escala dos dados.
grid[,c("estimate","lwr","upr")] <-
    exp(grid[,c("estimate","lwr","upr")])-0.5

## cld veio da glht() e M veio da HSD.test().
## estimate é uma média ajustada de mínimos quadrados, means é média
## amostral.

No gráfico estão sendo representados os valores médios na escala original. Eles foram obtidos pela transformação das médias ajustadas considerando a função inversa usada na resposta dados para adequação dos pressupostos. Em função disso, os intervalos na escala transformada não apresentam a mesma amplitude e nem simetria.

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

l <- levels(grid$Trat)
ylim <- c(0.5, length(l)+1)
sc <- list(y=list(at=seq_along(l), labels=l))

segplot(Trat~lwr+upr, data=grid, centers=estimate, draw=FALSE,
        ## ann=grid$cld,
        ann=grid$M, ## Teste de Tukey.
        ylim=ylim, scales=sc,
        xlab=lab, ylab="Consórcios",
        panel=function(z, centers, ann, subscripts, ...){
            panel.segplot(z=z, centers=centers,
                          subscripts=subscripts,
                          ...)
            ann <- sprintf("%0.2f %s", centers, ann)
            panel.text(x=centers, y=z, labels=ann, pos=3, cex=0.8)
        })

plot of chunk unnamed-chunk-19

NVagem: número de vagens por planta;

##-----------------------------------------------------------------------------
## Ver.

lab <- expression("Número de vagens por planta")

xyplot(NVagem~Trat, groups=Ano, data=dados,
       as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
       auto.key=list(columns=2),
       xlab="Consórcios", ylab=lab)

plot of chunk unnamed-chunk-20

Efeito de ano bem expressivo, inclusive com diferença de variabilidade. Ainda é importante mencionar que essa é uma variável do tipo contagem então se fazem justos os mesmos comentários feitos para NPD. Os valores são altos (distantes do zero), portanto não contínua como se supõe para uso da análise de variância. Por outro lado, certas transformações, se adequadamente aplicadas, podem permitir que os pressupostos sejam atendidos, não comprometendo a validade das inferências. Em especial, pelo fato de ser uma contagem alta, espera-se menor afastamento dos pressupostos do que aquele para contagens baixas.

##-----------------------------------------------------------------------------
## Ajuste.

summary(dados$NVagem)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    25.8    37.3    42.3    44.6    50.5    79.7
## Sem presença de zeros.

## Para os termos Trat e Ano:Trat, são os mesmos F.
m0 <- lm(NVagem~Ano+Ano:Rep+Ano*Trat, data=dados)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0, which=1:3)
MASS::boxcox(m0, lambda=seq(-2, 2, by=0.1))
abline(v=0, col=2)

plot of chunk unnamed-chunk-21

layout(1)

## Aplicando transformação.
m0 <- update(m0, formula=log(.)~.)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

plot of chunk unnamed-chunk-21

## Sem afastamentos dos pressupostos que sejam sistemáticos
## (corrigíveis) ou prejudiciais à inferência.

## ANOVA com valores de F corretos.
summary(aov(log(NVagem)~Ano+Error(Ano:Rep)+Ano*Trat, data=dados))
## 
## Error: Ano:Rep
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## Ano        1   2.11   2.109    22.9 0.00044 ***
## Residuals 12   1.10   0.092                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Trat       7  0.338  0.0483    1.60   0.15  
## Ano:Trat   7  0.388  0.0555    1.84   0.09 .
## Residuals 84  2.533  0.0302                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Médias ajustadas.

glr <- df.residual(m0)
qmr <- deviance(m0)/glr

## Teste de Tukey.
tk <- HSD.test(y=m0$model[,1], trt=m0$model["Trat"],
               DFerror=glr, MSerror=qmr)
tk$statistics["HSD"] ## DMS.
##     HSD
##   0.204
tk$groups
##   trt means M
## 1 SA  3.835 a
## 2 ST  3.833 a
## 3 SS  3.828 a
## 4 SD  3.786 a
## 5 SR  3.747 a
## 6 SPa 3.720 a
## 7 SPi 3.712 a
## 8 SX  3.692 a
## Remove os espaços extras adicionados pela função.
tk$groups[,-2] <- sapply(tk$groups[,-2],
                         FUN=gsub, pattern="\\s*$", replacement="")

## Para obter os IC para a média (ajustada).
## LSmeans(m0, effect="Trat")
L <- LSmatrix(m0, effect="Trat")
rownames(L) <- levels(dados$Trat)

## grid <- apmc(X=L, model=m0, focus="Trat", test="fdr", level=0.05)
grid <- apmc(X=L, model=m0, focus="Trat",
             test="single-step", level=0.05)
## grid

grid <- merge(grid, tk$groups, by.x="Trat", by.y="trt")
grid$Trat <- reorder(grid$Trat, grid$estimate)
grid
##   Trat estimate   lwr   upr cld means M
## 1   SA    3.835 3.743 3.927   a 3.835 a
## 2   SD    3.786 3.694 3.878   a 3.786 a
## 3  SPa    3.720 3.628 3.813   a 3.720 a
## 4  SPi    3.712 3.620 3.804   a 3.712 a
## 5   SR    3.747 3.655 3.839   a 3.747 a
## 6   SS    3.828 3.735 3.920   a 3.828 a
## 7   ST    3.833 3.740 3.925   a 3.833 a
## 8   SX    3.692 3.600 3.785   a 3.692 a
## Representar na escala dos dados.
grid[,c("estimate","lwr","upr")] <-
    exp(grid[,c("estimate","lwr","upr")])

## cld veio da glht() e M veio da HSD.test().
## estimate é uma média ajustada de mínimos quadrados, means é média
## amostral.

No gráfico estão sendo representados os valores médios na escala original. Eles foram obtidos pela transformação das médias ajustadas considerando a função inversa usada na resposta dados para adequação dos pressupostos. Em função disso, os intervalos na escala transformada não apresentam a mesma amplitude e nem simetria.

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

l <- levels(grid$Trat)
ylim <- c(0.5, length(l)+1)
sc <- list(y=list(at=seq_along(l), labels=l))

segplot(Trat~lwr+upr, data=grid, centers=estimate, draw=FALSE,
        ## ann=grid$cld,
        ann=grid$M, ## Teste de Tukey.
        ylim=ylim, scales=sc,
        xlab=lab, ylab="Consórcios",
        panel=function(z, centers, ann, subscripts, ...){
            panel.segplot(z=z, centers=centers,
                          subscripts=subscripts,
                          ...)
            ann <- sprintf("%0.2f %s", centers, ann)
            panel.text(x=centers, y=z, labels=ann, pos=3, cex=0.8)
        })

plot of chunk unnamed-chunk-22

NGVagem: Número de grãos por vagem (Número de vagem/número de plantas de soja);

##-----------------------------------------------------------------------------
## Ver.

lab <- expression("Número de grãos por vagem")

xyplot(NGVagem~Trat, groups=Ano, data=dados,
       as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
       auto.key=list(columns=2),
       xlab="Consórcios", ylab=lab)

plot of chunk unnamed-chunk-23

Efeito de ano bem expressivo, inclusive com diferença de variabilidade. Ainda é importante mencionar que essa é uma variável do tipo contagem então se fazem justos os mesmos comentários feitos para NPD. Os valores são altos (distantes do zero), portanto não contínua como se supõe para uso da análise de variância. Por outro lado, certas transformações, se adequadamente aplicadas, podem permitir que os pressupostos sejam atendidos, não comprometendo a validade das inferências. Em especial, pelo fato de ser uma contagem alta, espera-se menor afastamento dos pressupostos do que aquele para contagens baixas.

##-----------------------------------------------------------------------------
## Ajuste.

summary(dados$NGVagem)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.57    1.27    1.73    1.90    2.12    5.73
## MASS::fractions(dados$NGVagem)

## Sem presença de zeros.

## A observação 88 foi removida por ser considerada outlier por medidas
## de influencia.

## Para os termos Trat e Ano:Trat, são os mesmos F.
m0 <- lm(NGVagem~Ano+Ano:Rep+Ano*Trat, data=dados[-88,])
## anova(m0)

## im <- influence.measures(m0)
## summary(im)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0, which=1:3)
## MASS::boxcox(m0, lambda=seq(-2, 2, by=0.1))
## abline(v=0.5, col=2)
layout(1)

plot of chunk unnamed-chunk-24

## Sem afastamentos dos pressupostos que sejam sistemáticos
## (corrigíveis) ou prejudiciais à inferência.

## ANOVA com valores de F corretos.
summary(aov(NGVagem~Ano+Error(Ano:Rep)+Ano*Trat, data=dados[-88,]))
## 
## Error: Ano:Rep
##           Df Sum Sq Mean Sq F value Pr(>F)
## Ano        1    0.2    0.16    0.03   0.87
## Trat       1    2.9    2.89    0.47   0.51
## Residuals 11   67.0    6.09               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## Trat       7   1.88   0.269    2.07   0.056 .  
## Ano:Trat   7   4.94   0.705    5.44 3.7e-05 ***
## Residuals 83  10.77   0.130                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Médias ajustadas.

## O desbalanceamento impede de usar Tukey.

## glr <- df.residual(m0)
## qmr <- deviance(m0)/glr
## 
## ## Teste de Tukey.
## tk <- HSD.test(y=m0$model[,1], trt=m0$model["Trat"],
##                DFerror=glr, MSerror=qmr)
## tk$statistics["HSD"] ## DMS.
## tk$groups
## 
## ## Remove os espaços extras adicionados pela função.
## tk$groups[,-2] <- sapply(tk$groups[,-2],
##                          FUN=gsub, pattern="\\s*$", replacement="")

## Para obter os IC para a média (ajustada).
## LSmeans(m0, effect="Trat")
L <- LSmatrix(m0, effect=c("Ano","Trat"))
rownames(L) <- apply(attr(L, "grid"), 1, paste, collapse="/")

## grid <- apmc(X=L, model=m0, focus="Trat", test="fdr", level=0.05)
## grid <- apmc(X=L, model=m0, focus="Trat",
##              test="single-step", level=0.05)
L <- by(data=L, INDICES=attr(L, "grid")$Ano, FUN=as.matrix)
L <- lapply(L, apmc, model=m0, focus="Trat",
            test="single-step", level=0.05)
grid <- ldply(L, .id="Ano")
grid$Trat <- factor(gsub(pattern="^.*/", x=grid$Trat, replacement=""),
                    levels(dados$Trat))
## grid

## grid <- merge(grid, tk$groups, by.x="Trat", by.y="trt")
grid$Trat <- reorder(grid$Trat, grid$estimate)
## grid

## cld veio da glht() e M veio da HSD.test().
## estimate é uma média ajustada de mínimos quadrados, means é média
## amostral.

No gráfico estão sendo representados os valores médios na escala original. Eles foram obtidos pela transformação das médias ajustadas considerando a função inversa usada na resposta dados para adequação dos pressupostos. Em função disso, os intervalos na escala transformada não apresentam a mesma amplitude e nem simetria.

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

l <- levels(grid$Trat)
ylim <- c(0, length(l)+1)
sc <- list(y=list(at=seq_along(l), labels=l))

l <- levels(grid$Ano)
pch <- c(1,19)
key <- list(columns=length(l), type="o", divide=1,
            text=list(l),
            lines=list(pch=pch, lty=1))

segplot(Trat~lwr+upr, groups=grid$Ano, data=grid,
        centers=estimate, draw=FALSE,
        pch=pch[as.integer(grid$Ano)],
        ann=grid$cld, f=0.1,
        ## ann=grid$M, ## Teste de Tukey.
        ylim=ylim, scales=sc, key=key,
        xlab=lab, ylab="Consórcios",
        panel=function(z, centers, ann, subscripts, groups, f, ...){
            d <- 2*((as.integer(groups)-1)/(nlevels(groups)-1))-1
            z <- as.numeric(z)+f*d
            panel.segplot(z, centers=centers,
                          subscripts=subscripts, ...)
            ann <- sprintf("%0.2f %s", centers, ann)
            panel.text(x=centers, y=z, labels=ann,
                       pos=sign(d)+2, cex=0.8)
        })

plot of chunk unnamed-chunk-25

No gráfico acima, as letras distinguem médias de consórcios para um mesmo nível de ano.

M100G: massa de 100 grãos (g);

##-----------------------------------------------------------------------------
## Ver.

lab <- expression("Massa de 100 graõs"~("g"))

xyplot(M100G~Trat, groups=Ano, data=dados,
       as.table=TRUE, type=c("p","a"), jitter.x=TRUE,
       auto.key=list(columns=2),
       xlab="Consórcios", ylab=lab)

plot of chunk unnamed-chunk-26

##-----------------------------------------------------------------------------
## Ajuste.

## Para os termos Trat e Ano:Trat, são os mesmos F.
m0 <- lm(M100G~Ano+Ano:Rep+Ano*Trat, data=dados)
## anova(m0)

## Diagnóstico dos pressupostos.
par(mfrow=c(2,2)); plot(m0); layout(1)

plot of chunk unnamed-chunk-27

## Sem afastamentos dos pressupostos que sejam sistemáticos
## (corrigíveis) ou prejudiciais à inferência.

## ANOVA com valores de F corretos.
summary(aov(M100G~Ano+Error(Ano:Rep)+Ano*Trat, data=dados))
## 
## Error: Ano:Rep
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## Ano        1  144.1   144.1     279 1.1e-09 ***
## Residuals 12    6.2     0.5                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Trat       7   2.34   0.334    1.58   0.15
## Ano:Trat   7   2.29   0.327    1.55   0.16
## Residuals 84  17.75   0.211
##-----------------------------------------------------------------------------
## Médias ajustadas.

glr <- df.residual(m0)
qmr <- deviance(m0)/glr

## Teste de Tukey.
tk <- HSD.test(y=m0$model[,1], trt=m0$model["Trat"],
               DFerror=glr, MSerror=qmr)
tk$statistics["HSD"] ## DMS.
##      HSD
##   0.5401
tk$groups
##   trt means M
## 1 SX  11.50 a
## 2 SPa 11.38 a
## 3 SA  11.28 a
## 4 SPi 11.28 a
## 5 SD  11.24 a
## 6 ST  11.23 a
## 7 SS  11.11 a
## 8 SR  10.99 a
## Remove os espaços extras adicionados pela função.
tk$groups[,-2] <- sapply(tk$groups[,-2],
                         FUN=gsub, pattern="\\s*$", replacement="")

## Para obter os IC para a média (ajustada).
## LSmeans(m0, effect="Trat")
L <- LSmatrix(m0, effect="Trat")
rownames(L) <- levels(dados$Trat)

## grid <- apmc(X=L, model=m0, focus="Trat", test="fdr", level=0.05)
grid <- apmc(X=L, model=m0, focus="Trat",
             test="single-step", level=0.05)
## grid

grid <- merge(grid, tk$groups, by.x="Trat", by.y="trt")
grid$Trat <- reorder(grid$Trat, grid$estimate)
grid
##   Trat estimate   lwr   upr cld means M
## 1   SA    11.28 11.03 11.52   a 11.28 a
## 2   SD    11.24 10.99 11.48   a 11.24 a
## 3  SPa    11.38 11.13 11.62   a 11.38 a
## 4  SPi    11.28 11.03 11.52   a 11.28 a
## 5   SR    10.99 10.74 11.23   a 10.99 a
## 6   SS    11.11 10.87 11.35   a 11.11 a
## 7   ST    11.23 10.99 11.48   a 11.23 a
## 8   SX    11.49 11.25 11.74   a 11.50 a
## cld veio da glht() e M veio da HSD.test().
## estimate é uma média ajustada de mínimos quadrados, means é média
## amostral.
##-----------------------------------------------------------------------------

l <- levels(grid$Trat)
ylim <- c(0.5, length(l)+1)
sc <- list(y=list(at=seq_along(l), labels=l))

segplot(Trat~lwr+upr, data=grid, centers=estimate, draw=FALSE,
        ## ann=grid$cld,
        ann=grid$M, ## Teste de Tukey.
        ylim=ylim, scales=sc,
        xlab=lab, ylab="Consórcios",
        panel=function(z, centers, ann, subscripts, ...){
            panel.segplot(z=z, centers=centers,
                          subscripts=subscripts,
                          ...)
            ann <- sprintf("%0.2f %s", centers, ann)
            panel.text(x=centers, y=z, labels=ann, pos=3, cex=0.8)
        })

plot of chunk unnamed-chunk-28

##-----------------------------------------------------------------------------
## Versões dos pacotes e data do documento.

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              
##  [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  base     
## 
## other attached packages:
##  [1] wzRfun_0.4          agricolae_1.2-1     multcomp_1.3-6     
##  [4] TH.data_1.0-3       mvtnorm_0.9-9997    doBy_4.5-10        
##  [7] MASS_7.3-39         survival_2.38-1     plyr_1.8.1         
## [10] latticeExtra_0.6-26 lattice_0.20-31     RColorBrewer_1.0-5 
## [13] rmarkdown_0.5.1     knitr_1.6          
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.11.0     spdep_0.5-88    formatR_1.0     nloptr_1.0.4   
##  [5] methods_3.2.0   LearnBayes_2.15 tools_3.2.0     boot_1.3-15    
##  [9] digest_0.6.4    lme4_1.1-7      evaluate_0.5.5  nlme_3.1-120   
## [13] Matrix_1.2-0    yaml_2.1.13     coda_0.16-1     stringr_0.6.2  
## [17] cluster_2.0.1   combinat_0.0-8  grid_3.2.0      sp_1.0-15      
## [21] minqa_1.2.3     klaR_0.6-12     deldir_0.0-22   htmltools_0.2.6
## [25] splines_3.2.0   sandwich_2.3-0  zoo_1.7-10
Sys.time()
## [1] "2015-05-29 19:08:35 BRT"