Descrição do experimento

Este experimento estudou o efeito do carvão enriquecido com fósforo no desenvolvimento de mudas de eucalipto. Foram estudados 3 fontes do biocarvão (biocar):

A quantidade de fósforo usada também foi um fator, com 10 níveis, praticamente igualmente espaçados na escala log base 2.

O experimento foi feito em delineamento de blocos casualizados completos em casa de vegetação. As unidades experimentais foram um vaso com uma muda de eucalipto em cada. As variáveis de planta medidas foram:

Para cada variável fez-se a análise de variância considerando o modelo estatístico para um experimento de delineamento de blocos com fatores em arranjo fatorial 3 x 10. Em caso de interação significativa (5%), o desdobramento estudou as diferencias entre as 3 fontes de biocarvão para uma dose fixa de fósforo por meio do teste de Tukey (5%). No caso de não haver interação, fez-se o estudo dos efeitos principais. Os pressupostos do modelo estatístico foram avaliados graficamente. Havendo fuga dos pressupostos, será empregada transformação aos dados para que haja atendimento.

#-----------------------------------------------------------------------
# Pacotes.

library(gdata)
library(lattice)
library(latticeExtra)
library(splines)
library(doBy)
library(multcomp)
library(agricolae)
library(plyr)

source("~/Dropbox/templates/_setup.R")

opts_chunk$set(fig.width = 7,
               fig.height = 5)

# url <-
#     c("https://raw.githubusercontent.com/walmes/wzRfun/master/R/panel.cbH.R",
#       "https://raw.githubusercontent.com/walmes/wzRfun/master/R/prepanel.cbH.R",
#       "https://raw.githubusercontent.com/walmes/wzRfun/master/R/panel.segplotBy.R")
# sapply(url, source)

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

tukeymeans <- function(m0) {
    L <- list()
    L[[1]] <- do.call(rbind,
                      lapply(levels(db$catep),
                             FUN = function(p) {
                                 df <- df.residual(m0)
                                 ms <- deviance(m0)/df
                                 da <- m0$model
                                 names(da)[1] <- "y"
                                 bp <- with(
                                     subset(da, catep == p),
                                     HSD.test(y = y,
                                              trt = biocar,
                                              MSerror = ms,
                                              DFerror = df))
                                 cbind(catep = p, bp$groups)
                             }))
    names(L[[1]])[2] <- "biocar"
    L[[2]] <- do.call(rbind,
                      lapply(levels(db$biocar),
                             FUN = function(p) {
                                 df <- df.residual(m0)
                                 ms <- deviance(m0)/df
                                 da <- m0$model
                                 names(da)[1] <- "y"
                                 bp <- with(
                                     subset(da, biocar == p),
                                     HSD.test(y = y,
                                              trt = catep,
                                              MSerror = ms,
                                              DFerror = df))
                                 bp <- cbind(biocar = p, bp$groups)
                                 return(bp)
                             }))
    names(L[[2]])[2] <- "catep"
    L <- lapply(L, FUN = function(d) {
        as.data.frame((lapply(d, function(x) {
            if (is.factor(x)){
                levels(x) <- trimws(levels(x))
            }
            return(x)
        })))
    })
    md <- merge(L[[1]], L[[2]],
                by = c("catep", "biocar", "means"),
                all = TRUE)
    md$M.x <- toupper(md$M.x)
    names(md)[4:5] <- c("Mbio", "Mcat")
    arrange(md, biocar, as.numeric(as.character(catep)))
}

panel.segplotBy <- function(x, y, z,
                            centers, subscripts,
                            groups, f, letters, digits, ...){
    d <- 2 * ((as.numeric(groups) - 1)/(nlevels(groups) - 1)) - 1
    z <- as.numeric(z) + f * d
    panel.segplot(x, y, z, centers = centers,
                  subscripts = subscripts, ...)
    panel.text(x = z, y = centers,
               labels = sprintf(paste0("%0.", digits, "f %s"),
                                centers, letters),
               cex = 0.8, srt = 90, adj = c(0.5, -0.5))
}

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

# system("file -b biocarvao.xlsx")
da <- read.xls("biocarvao.xlsx", sheet = 1,
               encoding = "latin1")

da <- da[, !sapply(da, FUN = is.logical)]

da <- subset(da, select = -c(6:7))
str(da)

# dput(gsub(x = iconv(tolower(names(da)), to = "ASCII//TRANSLIT"),
#           pattern = "\\..*$", replacement = ""))

names(da) <-
    c("biocar", "dosep", "tensao", "bloco",
      "umid", "altura", "diametro", "mfpa", "mfr", "mspa", "msr")

# help(within, h = "html")

# with(da, {
#     list(
#         unique(biocar),
#         unique(dosep),
#         unique(tensao),
#         unique(bloco))
# })

levels(da$biocar) <- c("S_BP", "SP_B", "SP")

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

db <- subset(da, tensao == 0,
             select = -c(umid, tensao))
str(db)

Altura de Mudas

db$catep <- factor(db$dosep)
db$bloco <- factor(db$bloco)
db$dosept <- log(db$dosep + 1, base = 2)

any(is.na(db$altura))
## [1] FALSE
xyplot(altura ~ dosept, groups = biocar, data = db,
       type = c("p", "a"))

m0 <- lm(altura ~ bloco + biocar * catep, data = db)

par(mfrow = c(2, 2))
plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: altura
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## bloco         2  162.9   81.43  0.9373 0.3975388    
## biocar        2 2654.9 1327.43 15.2786 4.676e-06 ***
## catep         9 1431.4  159.05  1.8306 0.0818713 .  
## biocar:catep 18 4847.8  269.32  3.0999 0.0005685 ***
## Residuals    58 5039.1   86.88                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = c("biocar", "catep"))
grid <- attr(L, "grid")
grid$biocar <- factor(grid$biocar, levels = levels(db$biocar))
grid$catep <- factor(grid$catep, levels = levels(db$catep))

grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = L),
                          calpha = univariate_calpha())$confint))
names(grid)[3] <- "fit"

xlab <- expression("Dose de P"~(mg~dm^{-3}~", "~log[2](P+1)))
ylab <- "Altura de mudas (cm)"
key <- list(type = "o", divide = 1, columns = 3,
            text = list(levels(grid$biocar)),
            lines = list(pch = c(1, 5, 8)))

grid <- merge(grid, tukeymeans(m0)[, -3], sort = FALSE)

segplot(catep ~ lwr + upr,
        centers = fit,
        data = grid,
        draw = FALSE, horizontal = FALSE,
        groups = biocar,
        pch = key$lines$pch,
        key = key, xlab = xlab, ylab = ylab,
        f = 0.15, panel = panel.segplotBy,
        letters = tolower(grid$Mbio), digits = 2)

grid
##    biocar catep      fit      lwr       upr Mbio Mcat
## 1    S_BP     0 66.00000 55.22776  76.77224    A  abc
## 2    SP_B     0 70.00000 59.22776  80.77224    A   ab
## 3      SP     0 80.00000 69.22776  90.77224    A    a
## 4    S_BP   1.5 70.33333 59.56109  81.10557    A  abc
## 5    SP_B   1.5 87.66667 76.89443  98.43891    A    a
## 6      SP   1.5 80.66667 69.89443  91.43891    A    a
## 7    S_BP     3 62.00000 51.22776  72.77224    B   bc
## 8    SP_B     3 78.66667 67.89443  89.43891   AB   ab
## 9      SP     3 88.33333 77.56109  99.10557    A    a
## 10   S_BP     6 66.00000 55.22776  76.77224    A  abc
## 11   SP_B     6 69.33333 58.56109  80.10557    A   ab
## 12     SP     6 79.33333 68.56109  90.10557    A    a
## 13   S_BP    12 81.33333 70.56109  92.10557    A  abc
## 14   SP_B    12 55.33333 44.56109  66.10557    B    b
## 15     SP    12 90.00000 79.22776 100.77224    A    a
## 16   S_BP    24 85.00000 74.22776  95.77224    A   ab
## 17   SP_B    24 73.66667 62.89443  84.43891    A   ab
## 18     SP    24 82.66667 71.89443  93.43891    A    a
## 19   S_BP    48 88.00000 77.22776  98.77224    A    a
## 20   SP_B    48 54.00000 43.22776  64.77224    B    b
## 21     SP    48 73.00000 62.22776  83.77224    A    a
## 22   S_BP    96 80.00000 69.22776  90.77224   AB  abc
## 23   SP_B    96 68.33333 57.56109  79.10557    B   ab
## 24     SP    96 88.00000 77.22776  98.77224    A    a
## 25   S_BP   192 77.33333 66.56109  88.10557    A  abc
## 26   SP_B   192 74.00000 63.22776  84.77224    A   ab
## 27     SP   192 84.00000 73.22776  94.77224    A    a
## 28   S_BP   348 58.33333 47.56109  69.10557    B    c
## 29   SP_B   348 65.66667 54.89443  76.43891   AB   ab
## 30     SP   348 80.00000 69.22776  90.77224    A    a

Diâmetro de Mudas

any(is.na(db$diametro))
## [1] FALSE
xyplot(diametro ~ dosept, groups = biocar, data = db,
       type = c("p", "a"))

m0 <- lm(diametro ~ bloco + biocar * catep, data = db)

par(mfrow = c(2, 2))
plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: diametro
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## bloco         2  1.0854  0.5427  3.0453 0.0552531 .  
## biocar        2 12.8382  6.4191 36.0193 6.781e-11 ***
## catep         9  8.9020  0.9891  5.5502 1.699e-05 ***
## biocar:catep 18  9.9412  0.5523  3.0991 0.0005701 ***
## Residuals    58 10.3364  0.1782                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = c("biocar", "catep"))
grid <- attr(L, "grid")
grid$biocar <- factor(grid$biocar, levels = levels(db$biocar))
grid$catep <- factor(grid$catep, levels = levels(db$catep))

grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = L),
                          calpha = univariate_calpha())$confint))
names(grid)[3] <- "fit"

xlab <- expression("Dose de P"~(mg~dm^{-3}~", "~log[2](P+1)))
ylab <- "Diâmetro de mudas (mm)"
key <- list(type = "o", divide = 1, columns = 3,
            text = list(levels(grid$biocar)),
            lines = list(pch = c(1, 5, 8)))

grid <- merge(grid, tukeymeans(m0)[, -3], sort = FALSE)

segplot(catep ~ lwr + upr,
        centers = fit,
        data = grid,
        draw = FALSE, horizontal = FALSE,
        groups = biocar,
        pch = key$lines$pch,
        key = key, xlab = xlab, ylab = ylab,
        f = 0.15, panel = panel.segplotBy,
        letters = tolower(grid$Mbio), digits = 2)

grid
##    biocar catep      fit      lwr      upr Mbio Mcat
## 1    S_BP     0 5.616667 5.128788 6.104545    B   bc
## 2    SP_B     0 6.203333 5.715455 6.691212   AB   ab
## 3      SP     0 6.820000 6.332121 7.307879    A    a
## 4    S_BP   1.5 5.920000 5.432121 6.407879    B   ab
## 5    SP_B   1.5 6.343333 5.855455 6.831212   AB    a
## 6      SP   1.5 6.920000 6.432121 7.407879    A    a
## 7    S_BP     3 6.426667 5.938788 6.914545    B   ab
## 8    SP_B     3 6.493333 6.005455 6.981212    B    a
## 9      SP     3 7.366667 6.878788 7.854545    A    a
## 10   S_BP     6 6.036667 5.548788 6.524545    A   ab
## 11   SP_B     6 6.290000 5.802121 6.777879    A    a
## 12     SP     6 6.543333 6.055455 7.031212    A    a
## 13   S_BP    12 6.683333 6.195455 7.171212    A   ab
## 14   SP_B    12 6.173333 5.685455 6.661212    A   ab
## 15     SP    12 6.966667 6.478788 7.454545    A    a
## 16   S_BP    24 6.963333 6.475455 7.451212    A    a
## 17   SP_B    24 6.163333 5.675455 6.651212    A   ab
## 18     SP    24 6.870000 6.382121 7.357879    A    a
## 19   S_BP    48 6.780000 6.292121 7.267879    A    a
## 20   SP_B    48 5.696667 5.208788 6.184545    B   ab
## 21     SP    48 6.873333 6.385455 7.361212    A    a
## 22   S_BP    96 6.976667 6.488788 7.464545    A    a
## 23   SP_B    96 5.650000 5.162121 6.137879    B   ab
## 24     SP    96 6.900000 6.412121 7.387879    A    a
## 25   S_BP   192 6.193333 5.705455 6.681212   AB   ab
## 26   SP_B   192 5.763333 5.275455 6.251212    B   ab
## 27     SP   192 6.723333 6.235455 7.211212    A    a
## 28   S_BP   348 4.756667 4.268788 5.244545    B    c
## 29   SP_B   348 5.130000 4.642121 5.617879    B    b
## 30     SP   348 6.873333 6.385455 7.361212    A    a

Massa Fresca de Parte Aérea

any(is.na(db$mfpa))
## [1] FALSE
xyplot(mfpa ~ dosept, groups = biocar, data = db,
       type = c("p", "a"))

m0 <- lm(mfpa ~ bloco + biocar * catep, data = db)

par(mfrow = c(2, 2))
plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: mfpa
##              Df Sum Sq Mean Sq  F value    Pr(>F)    
## bloco         2  465.3   232.7   5.1510  0.008727 ** 
## biocar        2 9397.6  4698.8 104.0247 < 2.2e-16 ***
## catep         9 1080.5   120.1   2.6579  0.011858 *  
## biocar:catep 18 3860.1   214.4   4.7476 2.919e-06 ***
## Residuals    58 2619.9    45.2                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = c("biocar", "catep"))
grid <- attr(L, "grid")
grid$biocar <- factor(grid$biocar, levels = levels(db$biocar))
grid$catep <- factor(grid$catep, levels = levels(db$catep))

grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = L),
                          calpha = univariate_calpha())$confint))
names(grid)[3] <- "fit"

xlab <- expression("Dose de P"~(mg~dm^{-3}~", "~log[2](P+1)))
ylab <- "Massa fresca de parte aérea (g)"
key <- list(type = "o", divide = 1, columns = 3,
            text = list(levels(grid$biocar)),
            lines = list(pch = c(1, 5, 8)))

grid <- merge(grid, tukeymeans(m0)[, -3], sort = FALSE)

segplot(catep ~ lwr + upr,
        centers = fit,
        data = grid,
        draw = FALSE, horizontal = FALSE,
        groups = biocar,
        pch = key$lines$pch,
        key = key, xlab = xlab, ylab = ylab,
        f = 0.15, panel = panel.segplotBy,
        letters = tolower(grid$Mbio), digits = 2)

grid
##    biocar catep      fit      lwr      upr Mbio Mcat
## 1    S_BP     0 45.46000 37.69275 53.22725    B   de
## 2    SP_B     0 54.07333 46.30608 61.84058    B    a
## 3      SP     0 79.44000 71.67275 87.20725    A    a
## 4    S_BP   1.5 51.50000 43.73275 59.26725    B  cde
## 5    SP_B   1.5 61.19667 53.42942 68.96392    B    a
## 6      SP   1.5 80.64333 72.87608 88.41058    A    a
## 7    S_BP     3 59.07000 51.30275 66.83725    B abcd
## 8    SP_B     3 60.85000 53.08275 68.61725    B    a
## 9      SP     3 82.23000 74.46275 89.99725    A    a
## 10   S_BP     6 55.71333 47.94608 63.48058    B bcde
## 11   SP_B     6 60.50333 52.73608 68.27058    B    a
## 12     SP     6 74.59667 66.82942 82.36392    A    a
## 13   S_BP    12 68.31333 60.54608 76.08058    A  abc
## 14   SP_B    12 50.33667 42.56942 58.10392    B    a
## 15     SP    12 80.97667 73.20942 88.74392    A    a
## 16   S_BP    24 69.02000 61.25275 76.78725   AB  abc
## 17   SP_B    24 57.87667 50.10942 65.64392    B    a
## 18     SP    24 78.52667 70.75942 86.29392    A    a
## 19   S_BP    48 73.71333 65.94608 81.48058    A   ab
## 20   SP_B    48 51.48667 43.71942 59.25392    B    a
## 21     SP    48 75.84000 68.07275 83.60725    A    a
## 22   S_BP    96 74.68667 66.91942 82.45392    A    a
## 23   SP_B    96 49.43667 41.66942 57.20392    B    a
## 24     SP    96 80.30667 72.53942 88.07392    A    a
## 25   S_BP   192 72.14000 64.37275 79.90725    A   ab
## 26   SP_B   192 57.20000 49.43275 64.96725    B    a
## 27     SP   192 78.08333 70.31608 85.85058    A    a
## 28   S_BP   348 40.04333 32.27608 47.81058    B    e
## 29   SP_B   348 52.60000 44.83275 60.36725    B    a
## 30     SP   348 83.61000 75.84275 91.37725    A    a

Massa Seca de Parte Aérea

any(is.na(db$mspa))
## [1] FALSE
xyplot(mspa ~ dosept, groups = biocar, data = db,
       type = c("p", "a"))

m0 <- lm(mspa ~ bloco + biocar * catep, data = db)

par(mfrow = c(2, 2))
plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: mspa
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## bloco         2   40.16   20.08  3.3547  0.041817 *  
## biocar        2 1024.98  512.49 85.6277 < 2.2e-16 ***
## catep         9  164.88   18.32  3.0609  0.004566 ** 
## biocar:catep 18  535.48   29.75  4.9705 1.504e-06 ***
## Residuals    58  347.13    5.99                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = c("biocar", "catep"))
grid <- attr(L, "grid")
grid$biocar <- factor(grid$biocar, levels = levels(db$biocar))
grid$catep <- factor(grid$catep, levels = levels(db$catep))

grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = L),
                          calpha = univariate_calpha())$confint))
names(grid)[3] <- "fit"

xlab <- expression("Dose de P"~(mg~dm^{-3}~", "~log[2](P+1)))
ylab <- "Massa seca de parte aérea (g)"
key <- list(type = "o", divide = 1, columns = 3,
            text = list(levels(grid$biocar)),
            lines = list(pch = c(1, 5, 8)))

grid <- merge(grid, tukeymeans(m0)[, -3], sort = FALSE)

segplot(catep ~ lwr + upr,
        centers = fit,
        data = grid,
        draw = FALSE, horizontal = FALSE,
        groups = biocar,
        pch = key$lines$pch,
        key = key, xlab = xlab, ylab = ylab,
        f = 0.15, panel = panel.segplotBy,
        letters = tolower(grid$Mbio), digits = 2)

grid
##    biocar catep      fit       lwr      upr Mbio Mcat
## 1    S_BP     0 10.39867  7.571333 13.22600    C   cd
## 2    SP_B     0 15.74867 12.921333 18.57600    B    a
## 3      SP     0 22.17867 19.351333 25.00600    A    a
## 4    S_BP   1.5 12.14867  9.321333 14.97600    B  bcd
## 5    SP_B   1.5 15.82533 12.998000 18.65267    B    a
## 6      SP   1.5 23.42867 20.601333 26.25600    A    a
## 7    S_BP     3 14.93867 12.111333 17.76600    B  abc
## 8    SP_B     3 17.09533 14.268000 19.92267    B    a
## 9      SP     3 22.43200 19.604667 25.25933    A    a
## 10   S_BP     6 14.83200 12.004667 17.65933    B  abc
## 11   SP_B     6 17.36533 14.538000 20.19267   AB    a
## 12     SP     6 21.65867 18.831333 24.48600    A    a
## 13   S_BP    12 18.88867 16.061333 21.71600    A    a
## 14   SP_B    12 13.51200 10.684667 16.33933    B    a
## 15     SP    12 23.68867 20.861333 26.51600    A    a
## 16   S_BP    24 20.07533 17.248000 22.90267   AB    a
## 17   SP_B    24 16.36200 13.534667 19.18933    B    a
## 18     SP    24 22.44533 19.618000 25.27267    A    a
## 19   S_BP    48 21.26867 18.441333 24.09600    A    a
## 20   SP_B    48 12.26533  9.438000 15.09267    B    a
## 21     SP    48 22.05533 19.228000 24.88267    A    a
## 22   S_BP    96 21.06533 18.238000 23.89267    A    a
## 23   SP_B    96 12.71533  9.888000 15.54267    B    a
## 24     SP    96 22.31533 19.488000 25.14267    A    a
## 25   S_BP   192 18.03867 15.211333 20.86600   AB   ab
## 26   SP_B   192 16.52533 13.698000 19.35267    B    a
## 27     SP   192 21.84867 19.021333 24.67600    A    a
## 28   S_BP   348  8.25200  5.424667 11.07933    B    d
## 29   SP_B   348 12.51533  9.688000 15.34267    B    a
## 30     SP   348 23.93200 21.104667 26.75933    A    a

Massa Fresca de Raízes

any(is.na(db$mfr))
## [1] FALSE
xyplot(mfr ~ dosept, groups = biocar, data = db,
       type = c("p", "a"))

m0 <- lm(mfr ~ bloco + biocar * catep, data = db)

par(mfrow = c(2, 2))
plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: mfr
##              Df  Sum Sq Mean Sq F value   Pr(>F)   
## bloco         2  304.11 152.055  3.0721 0.053931 . 
## biocar        2  552.00 276.000  5.5763 0.006095 **
## catep         9 1093.31 121.479  2.4543 0.019196 * 
## biocar:catep 18 1559.24  86.624  1.7502 0.055953 . 
## Residuals    58 2870.72  49.495                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = c("biocar", "catep"))
grid <- attr(L, "grid")
grid$biocar <- factor(grid$biocar, levels = levels(db$biocar))
grid$catep <- factor(grid$catep, levels = levels(db$catep))

grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = L),
                          calpha = univariate_calpha())$confint))
names(grid)[3] <- "fit"

xlab <- expression("Dose de P"~(mg~dm^{-3}~", "~log[2](P+1)))
ylab <- "Massa fresca de raízes (g)"
key <- list(type = "o", divide = 1, columns = 3,
            text = list(levels(grid$biocar)),
            lines = list(pch = c(1, 5, 8)))

grid <- merge(grid, tukeymeans(m0)[, -3], sort = FALSE)

segplot(catep ~ lwr + upr,
        centers = fit,
        data = grid,
        draw = FALSE, horizontal = FALSE,
        groups = biocar,
        pch = key$lines$pch,
        key = key, xlab = xlab, ylab = ylab,
        f = 0.15, panel = panel.segplotBy,
        letters = tolower(grid$Mbio), digits = 2)

grid
##    biocar catep      fit       lwr      upr Mbio Mcat
## 1    S_BP     0 13.61000  5.479376 21.74062    A    a
## 2    SP_B     0 18.56667 10.436042 26.69729    A    a
## 3      SP     0 24.81667 16.686042 32.94729    A    a
## 4    S_BP   1.5 18.01667  9.886042 26.14729    A    a
## 5    SP_B   1.5 22.39667 14.266042 30.52729    A    a
## 6      SP   1.5 29.63000 21.499376 37.76062    A    a
## 7    S_BP     3 14.21333  6.082709 22.34396    B    a
## 8    SP_B     3 14.18333  6.052709 22.31396    B    a
## 9      SP     3 29.18000 21.049376 37.31062    A    a
## 10   S_BP     6 15.03000  6.899376 23.16062    A    a
## 11   SP_B     6 18.47667 10.346042 26.60729    A    a
## 12     SP     6 17.27000  9.139376 25.40062    A    a
## 13   S_BP    12 25.19000 17.059376 33.32062    A    a
## 14   SP_B    12 22.45667 14.326042 30.58729    A    a
## 15     SP    12 27.92333 19.792709 36.05396    A    a
## 16   S_BP    24 18.98000 10.849376 27.11062   AB    a
## 17   SP_B    24 11.47000  3.339376 19.60062    B    a
## 18     SP    24 29.79333 21.662709 37.92396    A    a
## 19   S_BP    48 26.70000 18.569376 34.83062    A    a
## 20   SP_B    48 21.54000 13.409376 29.67062    A    a
## 21     SP    48 32.33000 24.199376 40.46062    A    a
## 22   S_BP    96 26.73000 18.599376 34.86062    A    a
## 23   SP_B    96 14.40333  6.272709 22.53396    A    a
## 24     SP    96 24.23667 16.106042 32.36729    A    a
## 25   S_BP   192 31.81333 23.682709 39.94396    A    a
## 26   SP_B   192 26.91000 18.779376 35.04062    A    a
## 27     SP   192 24.92000 16.789376 33.05062    A    a
## 28   S_BP   348 20.59333 12.462709 28.72396    A    a
## 29   SP_B   348 24.94333 16.812709 33.07396    A    a
## 30     SP   348 13.79667  5.666042 21.92729    A    a

Massa Seca de Raízes

any(is.na(db$msr))
## [1] FALSE
xyplot(msr ~ dosept, groups = biocar, data = db,
       type = c("p", "a"))

m0 <- lm(msr ~ bloco + biocar * catep, data = db)

par(mfrow = c(2, 2))
plot(m0); layout(1)

anova(m0)
## Analysis of Variance Table
## 
## Response: msr
##              Df Sum Sq Mean Sq F value   Pr(>F)   
## bloco         2  23.80 11.8982  1.4074 0.253010   
## biocar        2  62.21 31.1032  3.6791 0.031311 * 
## catep         9 248.22 27.5800  3.2624 0.002840 **
## biocar:catep 18 363.52 20.1954  2.3889 0.006473 **
## Residuals    58 490.33  8.4539                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
L <- LSmatrix(m0, effect = c("biocar", "catep"))
grid <- attr(L, "grid")
grid$biocar <- factor(grid$biocar, levels = levels(db$biocar))
grid$catep <- factor(grid$catep, levels = levels(db$catep))

grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = L),
                          calpha = univariate_calpha())$confint))
names(grid)[3] <- "fit"

xlab <- expression("Dose de P"~(mg~dm^{-3}~", "~log[2](P+1)))
ylab <- "Massa seca de raízes (g)"
key <- list(type = "o", divide = 1, columns = 3,
            text = list(levels(grid$biocar)),
            lines = list(pch = c(1, 5, 8)))

grid <- merge(grid, tukeymeans(m0)[, -3], sort = FALSE)

segplot(catep ~ lwr + upr,
        centers = fit,
        data = grid,
        draw = FALSE, horizontal = FALSE,
        groups = biocar,
        pch = key$lines$pch,
        key = key, xlab = xlab, ylab = ylab,
        f = 0.15, panel = panel.segplotBy,
        letters = tolower(grid$Mbio), digits = 2)

grid
##    biocar catep       fit       lwr      upr Mbio Mcat
## 1    S_BP     0  7.516667  4.156414 10.87692    A    c
## 2    SP_B     0 10.150000  6.789747 13.51025    A    a
## 3      SP     0 11.423333  8.063081 14.78359    A    a
## 4    S_BP   1.5  7.763333  4.403081 11.12359    B    c
## 5    SP_B   1.5  8.986667  5.626414 12.34692   AB    a
## 6      SP   1.5 14.466667 11.106414 17.82692    A    a
## 7    S_BP     3  7.410000  4.049747 10.77025    A    c
## 8    SP_B     3  7.503333  4.143081 10.86359    A    a
## 9      SP     3 11.516667  8.156414 14.87692    A    a
## 10   S_BP     6  7.850000  4.489747 11.21025    A    c
## 11   SP_B     6  8.616667  5.256414 11.97692    A    a
## 12     SP     6  9.366667  6.006414 12.72692    A    a
## 13   S_BP    12 11.633333  8.273081 14.99359    A  abc
## 14   SP_B    12  9.073333  5.713081 12.43359    A    a
## 15     SP    12 12.026667  8.666414 15.38692    A    a
## 16   S_BP    24  9.293333  5.933081 12.65359   AB   bc
## 17   SP_B    24  7.463333  4.103081 10.82359    B    a
## 18     SP    24 13.220000  9.859747 16.58025    A    a
## 19   S_BP    48 13.426667 10.066414 16.78692    A  abc
## 20   SP_B    48  9.460000  6.099747 12.82025    A    a
## 21     SP    48 11.053333  7.693081 14.41359    A    a
## 22   S_BP    96 16.480000 13.119747 19.84025    A   ab
## 23   SP_B    96  8.090000  4.729747 11.45025    B    a
## 24     SP    96 11.833333  8.473081 15.19359   AB    a
## 25   S_BP   192 17.836667 14.476414 21.19692    A    a
## 26   SP_B   192 13.640000 10.279747 17.00025   AB    a
## 27     SP   192 11.990000  8.629747 15.35025    B    a
## 28   S_BP   348  7.556667  4.196414 10.91692    A    c
## 29   SP_B   348 12.340000  8.979747 15.70025    A    a
## 30     SP   348  8.736667  5.376414 12.09692    A    a