Delineamento Estrela em Experimento de Nutrição

Descrição do experimento

TODO
TODO
TODO
TODO
TODO
TODO
TODO

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

library(gdata)
library(lattice)
library(latticeExtra)
library(gridExtra)
library(plyr)
library(doBy)
library(multcomp)

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

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

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

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

# 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 Dados_Elvis.xlsx")
da <- read.xls("Dados_Elvis.xlsx", sheet = 1,
               encoding = "latin1")

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

da <- subset(da, select = -c(1:2))
str(da)

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

names(da) <-
    c("trt", "var", "N", "P", "K", "Ca", "Mg", "blc",
      "dc", "hmm", "hd")
da$blc <- factor(da$blc)
da$trt <- relevel(da$trt, ref = "Tmin")

knitr::kable(da)

trt	var	N	P	K	Ca	Mg	blc	dc	hmm	hd
N_50	N	50	1500	100	0.8	0.4	1	4.50	160	35.56
N_50	N	50	1500	100	0.8	0.4	2	5.60	215	38.39
N_50	N	50	1500	100	0.8	0.4	3	5.72	250	43.71
N_50	N	50	1500	100	0.8	0.4	4	5.48	190	34.67
N_150	N	150	1500	100	0.8	0.4	1	4.60	265	57.61
N_150	N	150	1500	100	0.8	0.4	2	8.20	300	36.59
N_150	N	150	1500	100	0.8	0.4	3	4.30	260	60.47
N_150	N	150	1500	100	0.8	0.4	4	2.60	200	76.92
N_200	N	200	1500	100	0.8	0.4	1	5.58	215	38.53
N_200	N	200	1500	100	0.8	0.4	2	5.35	230	42.99
N_200	N	200	1500	100	0.8	0.4	3	7.11	170	23.91
N_200	N	200	1500	100	0.8	0.4	4	3.65	270	73.97
P_150	P	100	150	100	0.8	0.4	1	5.03	230	45.73
P_150	P	100	150	100	0.8	0.4	2	5.91	250	42.30
P_150	P	100	150	100	0.8	0.4	3	3.80	160	42.11
P_150	P	100	150	100	0.8	0.4	4	4.84	240	49.59
P_450	P	100	450	100	0.8	0.4	1	5.31	285	53.67
P_450	P	100	450	100	0.8	0.4	2	7.23	280	38.73
P_450	P	100	450	100	0.8	0.4	3	6.00	285	47.50
P_450	P	100	450	100	0.8	0.4	4	5.14	190	36.96
P_600	P	100	600	100	0.8	0.4	1	4.59	280	61.00
P_600	P	100	600	100	0.8	0.4	2	5.17	255	49.32
P_600	P	100	600	100	0.8	0.4	3	5.22	220	42.15
P_600	P	100	600	100	0.8	0.4	4	6.18	280	45.31
K_50	K	100	1500	50	0.8	0.4	1	4.80	260	54.17
K_50	K	100	1500	50	0.8	0.4	2	5.75	280	48.70
K_50	K	100	1500	50	0.8	0.4	3	6.60	280	42.42
K_50	K	100	1500	50	0.8	0.4	4	6.80	220	32.35
K_150	K	100	1500	150	0.8	0.4	1	8.53	310	36.34
K_150	K	100	1500	150	0.8	0.4	2	5.80	235	40.52
K_150	K	100	1500	150	0.8	0.4	3	5.60	265	47.32
K_150	K	100	1500	150	0.8	0.4	4	4.70	200	42.55
K_200	K	100	1500	200	0.8	0.4	1	3.58	175	48.88
K_200	K	100	1500	200	0.8	0.4	2	5.48	255	46.53
K_200	K	100	1500	200	0.8	0.4	3	6.16	240	38.96
K_200	K	100	1500	200	0.8	0.4	4	6.83	285	41.73
Ca_0.4	Ca	100	1500	100	0.4	0.4	1	4.48	210	46.88
Ca_0.4	Ca	100	1500	100	0.4	0.4	2	5.02	225	44.82
Ca_0.4	Ca	100	1500	100	0.4	0.4	3	3.96	170	42.93
Ca_0.4	Ca	100	1500	100	0.4	0.4	4	5.60	260	46.43
Ca_1.6	Ca	100	1500	100	1.6	0.4	1	6.80	340	50.00
Ca_1.6	Ca	100	1500	100	1.6	0.4	2	4.40	225	51.14
Ca_1.6	Ca	100	1500	100	1.6	0.4	3	4.70	155	32.98
Ca_1.6	Ca	100	1500	100	1.6	0.4	4	7.73	290	37.52
Ca_3.2	Ca	100	1500	100	3.2	0.4	1	5.60	160	28.57
Ca_3.2	Ca	100	1500	100	3.2	0.4	2	5.50	190	34.55
Ca_3.2	Ca	100	1500	100	3.2	0.4	3	6.72	250	37.20
Ca_3.2	Ca	100	1500	100	3.2	0.4	4	5.50	235	42.73
Mg_0.2	Mg	100	1500	100	0.8	0.2	1	5.43	190	34.99
Mg_0.2	Mg	100	1500	100	0.8	0.2	2	4.60	225	48.91
Mg_0.2	Mg	100	1500	100	0.8	0.2	3	3.45	190	55.07
Mg_0.2	Mg	100	1500	100	0.8	0.2	4	4.77	300	62.89
Mg_0.8	Mg	100	1500	100	0.8	0.8	1	5.16	280	54.26
Mg_0.8	Mg	100	1500	100	0.8	0.8	2	3.50	200	57.14
Mg_0.8	Mg	100	1500	100	0.8	0.8	3	3.30	205	62.12
Mg_0.8	Mg	100	1500	100	0.8	0.8	4	4.34	215	49.54
Mg_1.6	Mg	100	1500	100	0.8	1.6	1	4.70	200	42.55
Mg_1.6	Mg	100	1500	100	0.8	1.6	2	6.87	300	43.67
Mg_1.6	Mg	100	1500	100	0.8	1.6	3	5.62	295	52.49
Mg_1.6	Mg	100	1500	100	0.8	1.6	4	4.82	270	56.02
Tref	Tref	100	1500	100	0.8	0.4	1	4.60	160	34.78
Tref	Tref	100	1500	100	0.8	0.4	2	5.02	230	45.82
Tref	Tref	100	1500	100	0.8	0.4	3	4.20	200	47.62
Tref	Tref	100	1500	100	0.8	0.4	4	5.00	225	45.00
Tmax	Tmax	200	3000	200	3.2	1.6	1	9.86	260	26.37
Tmax	Tmax	200	3000	200	3.2	1.6	2	4.73	270	57.08
Tmax	Tmax	200	3000	200	3.2	1.6	3	5.30	250	47.17
Tmax	Tmax	200	3000	200	3.2	1.6	4	5.15	240	46.60
Tmin	Tmin	0	0	0	0.0	0.0	1	5.29	220	41.59
Tmin	Tmin	0	0	0	0.0	0.0	2	4.00	240	60.00
Tmin	Tmin	0	0	0	0.0	0.0	3	4.30	200	46.51
Tmin	Tmin	0	0	0	0.0	0.0	4	4.90	230	46.94

# Delineamento estrela com testemunhas positiva, negativa e ponto
# central.
xtabs(~N + P, data = da)

##      P
## N      0 150 450 600 1500 3000
##   0    4   0   0   0    0    0
##   50   0   0   0   0    4    0
##   100  0   4   4   4   40    0
##   150  0   0   0   0    4    0
##   200  0   0   0   0    4    4

xtabs(~trt, data = da)

## trt
##   Tmin Ca_0.4 Ca_1.6 Ca_3.2  K_150  K_200   K_50 Mg_0.2 Mg_0.8 
##      4      4      4      4      4      4      4      4      4 
## Mg_1.6  N_150  N_200   N_50  P_150  P_450  P_600   Tmax   Tref 
##      4      4      4      4      4      4      4      4      4

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

# Indica quem são as testemunhas.
da$ttm <- as.character(da$var) %in% tail(levels(da$var), n = 3)

md <- c(subset(da, trt == "Tref",
               select = c(N, P, K, Ca, Mg))[1, ])

# Padroniza os níveis.
da <- transform(da,
                N = N/md$N,
                P = P/md$P,
                K = K/md$K,
                Ca = Ca/md$Ca,
                Mg = Mg/md$Mg)

# Representa a parta fatorial.
da$Tmin <- as.integer(da$trt == "Tmin")
da$Tmax <- as.integer(da$trt == "Tmax")
da$fat <- with(da, 1 - Tmin - Tmax)
str(da)

## 'data.frame':    72 obs. of  15 variables:
##  $ trt : Factor w/ 18 levels "Tmin","Ca_0.4",..: 13 13 13 13 11 11 11 11 12 12 ...
##  $ var : Factor w/ 8 levels "Ca","K","Mg",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ N   : num  0.5 0.5 0.5 0.5 1.5 1.5 1.5 1.5 2 2 ...
##  $ P   : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ K   : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ Ca  : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ Mg  : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ blc : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ dc  : num  4.5 5.6 5.72 5.48 4.6 8.2 4.3 2.6 5.58 5.35 ...
##  $ hmm : num  160 215 250 190 265 300 260 200 215 230 ...
##  $ hd  : num  35.6 38.4 43.7 34.7 57.6 ...
##  $ ttm : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ Tmin: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Tmax: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ fat : num  1 1 1 1 1 1 1 1 1 1 ...

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

# Modelo para efeitos de segundo grau.
form <- ~blc +
    (N + I(N^2)) +
    (P + I(P^2)) +
    (K + I(K^2)) +
    (Ca + I(Ca^2)) +
    (Mg + I(Mg^2))
mm <- model.matrix(form, data = da)

bl <- which(attr(mm, "assign") == 1)

Diâmetro do coleto (DC)

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

y <- "dc"
ylab <- "Diametro do coleto (mm)"
da$y <- da[, y]
grid.arrange(ncol = 2,
             xyplot(y ~ N, data = subset(da, var == "N" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ P, data = subset(da, var == "P" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ K, data = subset(da, var == "K" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ Ca, data = subset(da, var == "Ca" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ Mg, data = subset(da, var == "Mg" | ttm),
                    ylab = ylab, groups = ttm))

#-----------------------------------------------------------------------
# Ajuste do modelo com níveis categóricos.

m0 <- lm(y ~ blc + trt, data = da)
par(mfrow = c(2, 2))
plot(m0); layout(1)

# X <- LSmatrix(m0, effect = "trt")
X <- LSmatrix(m0, effect = "trt",
              at = list(trt = c("Tmin", "Tref", "Tmax")))
grid <- attr(X, "grid")
grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = X),
                          calpha = univariate_calpha())$confint))
grid$trt <- factor(grid$trt, levels = grid$trt)

# Estimativas das testemunhas e ponto central.
grid

##    trt Estimate      lwr      upr
## 1 Tmin   4.6225 3.354992 5.890008
## 2 Tref   4.7050 3.437492 5.972508
## 3 Tmax   6.2600 4.992492 7.527508

# Contraste entre as testemunhas.
summary(glht(m0, linfct = apc(X)))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ blc + trt, data = da)
## 
## Linear Hypotheses:
##                Estimate Std. Error t value Pr(>|t|)
## Tmin-Tref == 0  -0.0825     0.8929  -0.092    0.995
## Tmin-Tmax == 0  -1.6375     0.8929  -1.834    0.169
## Tref-Tmax == 0  -1.5550     0.8929  -1.742    0.200
## (Adjusted p values reported -- single-step method)

#-----------------------------------------------------------------------
# Ajuste do modelo de regressão quadrática em cada nutriente.

m1 <- lm(y ~ 0 + mm, data = da)

# Falta de ajuste.
anova(m1, m0)

## Analysis of Variance Table
## 
## Model 1: y ~ 0 + mm
## Model 2: y ~ blc + trt
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     58 95.884                           
## 2     51 81.318  7    14.567 1.3051  0.267

# Resíduos do modelo.
par(mfrow = c(2, 2))
plot(m1); layout(1)

# Quadro de estimativas dos parâmetros.
summary(m1)

## 
## Call:
## lm(formula = y ~ 0 + mm, data = da)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5135 -0.6983 -0.1754  0.5429  3.5066 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## mm(Intercept)  4.92484    0.68163   7.225 1.23e-09 ***
## mmblc2        -0.01722    0.42859  -0.040    0.968    
## mmblc3        -0.35444    0.42859  -0.827    0.412    
## mmblc4        -0.24500    0.42859  -0.572    0.570    
## mmN            0.05912    2.15123   0.027    0.978    
## mmI(N^2)       0.06111    0.82639   0.074    0.941    
## mmP            0.36725    1.26535   0.290    0.773    
## mmI(P^2)      -0.44673    0.88531  -0.505    0.616    
## mmK           -1.51634    2.15123  -0.705    0.484    
## mmI(K^2)       0.74566    0.82639   0.902    0.371    
## mmCa           1.70875    1.16318   1.469    0.147    
## mmI(Ca^2)     -0.28968    0.24289  -1.193    0.238    
## mmMg          -0.47928    1.16318  -0.412    0.682    
## mmI(Mg^2)      0.11783    0.24289   0.485    0.629    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.286 on 58 degrees of freedom
## Multiple R-squared:  0.9553, Adjusted R-squared:  0.9445 
## F-statistic: 88.47 on 14 and 58 DF,  p-value: < 2.2e-16

#-----------------------------------------------------------------------
# Curvas preditas.

pred0 <- data.frame(N = 1, P = 1, K = 1, Ca = 1, Mg = 1,
                    blc = factor("1", levels(da$blc)))
n <- 30
pred0 <- pred0[rep(1, 30), ]

L <- vector(mode = "list", length = 5)
names(L) <- names(pred0)[1:length(L)]
for (k in names(L)) {
    pred1 <- pred0
    s <- with(da, seq(min(da[, k]), max(da[, k]), length.out = 30))
    pred1[, k] <- s
    X <- model.matrix(form, data = pred1)
    X[, bl] <- X[, bl] * 0 + 1/(length(bl) + 1)
    pred1 <- cbind(pred1,
                   as.data.frame(
                       confint(glht(m1, linfct = unique(X)),
                               calpha = univariate_calpha())$confint))
    pred1 <- pred1[, c(k, "Estimate", "lwr", "upr")]
    names(pred1)[1] <- "x"
    pred1$nutr <- k
    L[[k]] <- pred1
}

L <- do.call(rbind, L)
L$nutr <- factor(L$nutr, levels = names(pred0)[1:5])
str(L)

## 'data.frame':    150 obs. of  5 variables:
##  $ x       : num  0 0.069 0.138 0.207 0.276 ...
##  $ Estimate: num  4.98 4.98 4.99 4.99 5 ...
##  $ lwr     : num  2.26 2.54 2.81 3.06 3.29 ...
##  $ upr     : num  7.7 7.43 7.17 6.93 6.71 ...
##  $ nutr    : Factor w/ 5 levels "N","P","K","Ca",..: 1 1 1 1 1 1 1 1 1 1 ...

xyplot(Estimate ~ x | nutr,
       data = L, type = "l", #ylim = ylim,
       ylab = ylab,
       xlab = "Teor dos nutrientes (codificado)",
       scales = list(x = list(relation = "free")),
       as.table = TRUE,
       uy = L$upr, ly = L$lwr, cty = "bands",
       prepanel = prepanel.cbH,
       panel = panel.cbH) +
    layer(panel.arrows(c(0, 1, max(x)), grid$lwr,
                       c(0, 1, max(x)), grid$upr,
                       length = 0.025, code = 3, angle = 90)) +
    layer(panel.points(c(0, 1, max(x)), grid$Estimate))

Altura das mudas (H)

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

y <- "hmm"
ylab <- "Altura das mudas (cm)"
da$y <- da[, y]
grid.arrange(ncol = 2,
             xyplot(y ~ N, data = subset(da, var == "N" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ P, data = subset(da, var == "P" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ K, data = subset(da, var == "K" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ Ca, data = subset(da, var == "Ca" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ Mg, data = subset(da, var == "Mg" | ttm),
                    ylab = ylab, groups = ttm))

#-----------------------------------------------------------------------
# Ajuste do modelo com níveis categóricos.

m0 <- lm(y ~ blc + trt, data = da)
par(mfrow = c(2, 2))
plot(m0); layout(1)

# X <- LSmatrix(m0, effect = "trt")
X <- LSmatrix(m0, effect = "trt",
              at = list(trt = c("Tmin", "Tref", "Tmax")))
grid <- attr(X, "grid")
grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = X),
                          calpha = univariate_calpha())$confint))
grid$trt <- factor(grid$trt, levels = grid$trt)

# Estimativas das testemunhas e ponto central.
grid

##    trt Estimate      lwr      upr
## 1 Tmin   222.50 180.0043 264.9957
## 2 Tref   203.75 161.2543 246.2457
## 3 Tmax   255.00 212.5043 297.4957

# Contraste entre as testemunhas.
summary(glht(m0, linfct = apc(X)))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ blc + trt, data = da)
## 
## Linear Hypotheses:
##                Estimate Std. Error t value Pr(>|t|)
## Tmin-Tref == 0    18.75      29.94   0.626    0.806
## Tmin-Tmax == 0   -32.50      29.94  -1.086    0.527
## Tref-Tmax == 0   -51.25      29.94  -1.712    0.210
## (Adjusted p values reported -- single-step method)

#-----------------------------------------------------------------------
# Ajuste do modelo de regressão quadrática em cada nutriente.

m1 <- lm(y ~ 0 + mm, data = da)

# Falta de ajuste.
anova(m1, m0)

## Analysis of Variance Table
## 
## Model 1: y ~ 0 + mm
## Model 2: y ~ blc + trt
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     58 104451                           
## 2     51  91406  7     13045 1.0398 0.4156

# Resíduos do modelo.
par(mfrow = c(2, 2))
plot(m1); layout(1)

# Quadro de estimativas dos parâmetros.
summary(m1)

## 
## Call:
## lm(formula = y ~ 0 + mm, data = da)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -81.027 -23.022  -1.911  30.052  95.362 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## mm(Intercept)  218.2234    22.4972   9.700 9.48e-14 ***
## mmblc2          11.3889    14.1456   0.805   0.4240    
## mmblc3          -8.6111    14.1456  -0.609   0.5451    
## mmblc4           7.7778    14.1456   0.550   0.5845    
## mmN            112.8227    71.0018   1.589   0.1175    
## mmI(N^2)       -37.7826    27.2750  -1.385   0.1713    
## mmP            -14.7439    41.7630  -0.353   0.7253    
## mmI(P^2)        -0.8653    29.2199  -0.030   0.9765    
## mmK           -122.1773    71.0018  -1.721   0.0906 .  
## mmI(K^2)        47.2174    27.2750   1.731   0.0887 .  
## mmCa            50.3986    38.3911   1.313   0.1944    
## mmI(Ca^2)      -11.2920     8.0165  -1.409   0.1643    
## mmMg           -20.2593    38.3911  -0.528   0.5997    
## mmI(Mg^2)        6.5735     8.0165   0.820   0.4156    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 42.44 on 58 degrees of freedom
## Multiple R-squared:  0.9748, Adjusted R-squared:  0.9687 
## F-statistic: 159.9 on 14 and 58 DF,  p-value: < 2.2e-16

#-----------------------------------------------------------------------
# Curvas preditas.

pred0 <- data.frame(N = 1, P = 1, K = 1, Ca = 1, Mg = 1,
                    blc = factor("1", levels(da$blc)))
n <- 30
pred0 <- pred0[rep(1, 30), ]

L <- vector(mode = "list", length = 5)
names(L) <- names(pred0)[1:length(L)]
for (k in names(L)) {
    pred1 <- pred0
    s <- with(da, seq(min(da[, k]), max(da[, k]), length.out = 30))
    pred1[, k] <- s
    X <- model.matrix(form, data = pred1)
    X[, bl] <- X[, bl] * 0 + 1/(length(bl) + 1)
    pred1 <- cbind(pred1,
                   as.data.frame(
                       confint(glht(m1, linfct = unique(X)),
                               calpha = univariate_calpha())$confint))
    pred1 <- pred1[, c(k, "Estimate", "lwr", "upr")]
    names(pred1)[1] <- "x"
    pred1$nutr <- k
    L[[k]] <- pred1
}

L <- do.call(rbind, L)
L$nutr <- factor(L$nutr, levels = names(pred0)[1:5])
str(L)

## 'data.frame':    150 obs. of  5 variables:
##  $ x       : num  0 0.069 0.138 0.207 0.276 ...
##  $ Estimate: num  156 163 171 177 184 ...
##  $ lwr     : num  65.9 82.7 98.6 113.5 127.6 ...
##  $ upr     : num  246 244 243 241 240 ...
##  $ nutr    : Factor w/ 5 levels "N","P","K","Ca",..: 1 1 1 1 1 1 1 1 1 1 ...

xyplot(Estimate ~ x | nutr,
       data = L, type = "l", #ylim = ylim,
       ylab = ylab,
       xlab = "Teor dos nutrientes (codificado)",
       scales = list(x = list(relation = "free")),
       as.table = TRUE,
       uy = L$upr, ly = L$lwr, cty = "bands",
       prepanel = prepanel.cbH,
       panel = panel.cbH) +
    layer(panel.arrows(c(0, 1, max(x)), grid$lwr,
                       c(0, 1, max(x)), grid$upr,
                       length = 0.025, code = 3, angle = 90)) +
    layer(panel.points(c(0, 1, max(x)), grid$Estimate))

Relação entre a altura da parte aérea e o diâmetro do coleto (H/D)

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

y <- "hd"
ylab <- expression("H/D"~(cm~cm^{-1}))
da$y <- da[, y]
grid.arrange(ncol = 2,
             xyplot(y ~ N, data = subset(da, var == "N" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ P, data = subset(da, var == "P" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ K, data = subset(da, var == "K" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ Ca, data = subset(da, var == "Ca" | ttm),
                    ylab = ylab, groups = ttm),
             xyplot(y ~ Mg, data = subset(da, var == "Mg" | ttm),
                    ylab = ylab, groups = ttm))

#-----------------------------------------------------------------------
# Ajuste do modelo com níveis categóricos.

m0 <- lm(y ~ blc + trt, data = da)
par(mfrow = c(2, 2))
plot(m0); layout(1)

# X <- LSmatrix(m0, effect = "trt")
X <- LSmatrix(m0, effect = "trt",
              at = list(trt = c("Tmin", "Tref", "Tmax")))
grid <- attr(X, "grid")
grid <- cbind(grid,
              as.data.frame(
                  confint(glht(m0, linfct = X),
                          calpha = univariate_calpha())$confint))
grid$trt <- factor(grid$trt, levels = grid$trt)

# Estimativas das testemunhas e ponto central.
grid

##    trt Estimate     lwr     upr
## 1 Tmin   48.760 39.2105 58.3095
## 2 Tref   43.305 33.7555 52.8545
## 3 Tmax   44.305 34.7555 53.8545

# Contraste entre as testemunhas.
summary(glht(m0, linfct = apc(X)))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ blc + trt, data = da)
## 
## Linear Hypotheses:
##                Estimate Std. Error t value Pr(>|t|)
## Tmin-Tref == 0    5.455      6.727   0.811    0.698
## Tmin-Tmax == 0    4.455      6.727   0.662    0.786
## Tref-Tmax == 0   -1.000      6.727  -0.149    0.988
## (Adjusted p values reported -- single-step method)

#-----------------------------------------------------------------------
# Ajuste do modelo de regressão quadrática em cada nutriente.

m1 <- lm(y ~ 0 + mm, data = da)

# Falta de ajuste.
anova(m1, m0)

## Analysis of Variance Table
## 
## Model 1: y ~ 0 + mm
## Model 2: y ~ blc + trt
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     58 5731.0                           
## 2     51 4615.8  7    1115.2 1.7603  0.116

# Resíduos do modelo.
par(mfrow = c(2, 2))
plot(m1); layout(1)

# Quadro de estimativas dos parâmetros.
summary(m1)

## 
## Call:
## lm(formula = y ~ 0 + mm, data = da)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.880  -5.355  -1.146   4.017  25.462 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## mm(Intercept)  44.90797    5.26970   8.522 8.22e-12 ***
## mmblc2          1.98444    3.31344   0.599    0.552    
## mmblc3          1.17556    3.31344   0.355    0.724    
## mmblc4          4.23556    3.31344   1.278    0.206    
## mmN            21.66052   16.63135   1.302    0.198    
## mmI(N^2)       -6.53575    6.38887  -1.023    0.311    
## mmP            -4.60948    9.78251  -0.471    0.639    
## mmI(P^2)        3.29669    6.84442   0.482    0.632    
## mmK           -14.14062   16.63135  -0.850    0.399    
## mmI(K^2)        3.86561    6.38887   0.605    0.548    
## mmCa           -5.61006    8.99268  -0.624    0.535    
## mmI(Ca^2)       0.42667    1.87778   0.227    0.821    
## mmMg            1.33086    8.99268   0.148    0.883    
## mmI(Mg^2)      -0.03079    1.87778  -0.016    0.987    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.94 on 58 degrees of freedom
## Multiple R-squared:  0.9637, Adjusted R-squared:  0.955 
## F-statistic:   110 on 14 and 58 DF,  p-value: < 2.2e-16

#-----------------------------------------------------------------------
# Curvas preditas.

pred0 <- data.frame(N = 1, P = 1, K = 1, Ca = 1, Mg = 1,
                    blc = factor("1", levels(da$blc)))
n <- 30
pred0 <- pred0[rep(1, 30), ]

L <- vector(mode = "list", length = 5)
names(L) <- names(pred0)[1:length(L)]
for (k in names(L)) {
    pred1 <- pred0
    s <- with(da, seq(min(da[, k]), max(da[, k]), length.out = 30))
    pred1[, k] <- s
    X <- model.matrix(form, data = pred1)
    X[, bl] <- X[, bl] * 0 + 1/(length(bl) + 1)
    pred1 <- cbind(pred1,
                   as.data.frame(
                       confint(glht(m1, linfct = unique(X)),
                               calpha = univariate_calpha())$confint))
    pred1 <- pred1[, c(k, "Estimate", "lwr", "upr")]
    names(pred1)[1] <- "x"
    pred1$nutr <- k
    L[[k]] <- pred1
}

L <- do.call(rbind, L)
L$nutr <- factor(L$nutr, levels = names(pred0)[1:5])
str(L)

## 'data.frame':    150 obs. of  5 variables:
##  $ x       : num  0 0.069 0.138 0.207 0.276 ...
##  $ Estimate: num  31.3 32.7 34.1 35.5 36.8 ...
##  $ lwr     : num  10.3 13.9 17.3 20.5 23.6 ...
##  $ upr     : num  52.3 51.6 51 50.5 50 ...
##  $ nutr    : Factor w/ 5 levels "N","P","K","Ca",..: 1 1 1 1 1 1 1 1 1 1 ...

xyplot(Estimate ~ x | nutr,
       data = L, type = "l", #ylim = ylim,
       ylab = ylab,
       xlab = "Teor dos nutrientes (codificado)",
       scales = list(x = list(relation = "free")),
       as.table = TRUE,
       uy = L$upr, ly = L$lwr, cty = "bands",
       prepanel = prepanel.cbH,
       panel = panel.cbH) +
    layer(panel.arrows(c(0, 1, max(x)), grid$lwr,
                       c(0, 1, max(x)), grid$upr,
                       length = 0.025, code = 3, angle = 90)) +
    layer(panel.points(c(0, 1, max(x)), grid$Estimate))

Informações da sessão

Sys.time()

## [1] "2016-05-21 09:10:24 BRT"

sessionInfo()

## R version 3.3.0 (2016-05-03)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.4 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] methods   stats     graphics  grDevices utils     datasets 
## [7] base     
## 
## other attached packages:
##  [1] multcomp_1.4-4      TH.data_1.0-7       MASS_7.3-45        
##  [4] survival_2.39-4     mvtnorm_1.0-5       doBy_4.5-15        
##  [7] plyr_1.8.3          gridExtra_2.2.1     latticeExtra_0.6-28
## [10] RColorBrewer_1.1-2  lattice_0.20-33     gdata_2.17.0       
## [13] rmarkdown_0.9.6     knitr_1.12.3       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.4      magrittr_1.5     splines_3.3.0   
##  [4] highr_0.5.1      stringr_1.0.0    tools_3.3.0     
##  [7] grid_3.3.0       gtable_0.2.0     htmltools_0.3.5 
## [10] gtools_3.5.0     yaml_2.1.13      digest_0.6.9    
## [13] Matrix_1.2-6     formatR_1.3      codetools_0.2-14
## [16] evaluate_0.9     sandwich_2.3-4   stringi_1.0-1   
## [19] zoo_1.7-13