library(doBy)
library(gridExtra)
library(tidyverse)
da <- read_csv2(file = "RCS_Dados.csv",
col_names = TRUE,
skip = 1)
da1 <- da %>% select(1:20)
da2 <- da %>% select(1:4, 21:32)
names(da2) <- names(da2) %>%
gsub(pattern = "_1", replacement = "")
da1$cam <- "C1"
da2$cam <- "C2"
da <- bind_rows(da1, da2)
# View(da)
# dput(names(da))
names(da) <- c("trat", "descr", "bloc", "prod", "pHH2O", "pHCaCl2",
"Pmehlich", "K", "Ca", "Mg", "Al", "MO", "CTC", "V",
"Zn", "Cu", "Fe", "Mn", "S", "B", "cam", "SB")
str(da)
## Classes 'tbl_df', 'tbl' and 'data.frame': 80 obs. of 22 variables:
## $ trat : chr "1" "1" "1" "1" ...
## $ descr : chr "Soja (PD)" "Soja (PD)" "Soja (PD)" "Soja (PD)" ...
## $ bloc : int 1 2 3 4 1 2 3 4 1 2 ...
## $ prod : num 61 59 61.4 62.3 77.8 82.2 77.7 72 82.6 79.6 ...
## $ pHH2O : num 6.5 6.5 6.1 6.2 6.5 6.5 6.5 6.4 6.5 6.5 ...
## $ pHCaCl2 : num 5.7 5.7 5.4 5.5 5.7 5.7 5.7 5.6 5.7 5.7 ...
## $ Pmehlich: num 23.9 18.5 18.5 19.6 15.3 10.7 11.6 11.6 20.8 11.1 ...
## $ K : num 104.3 115.2 80.9 88.8 73.9 ...
## $ Ca : num 4.8 4.2 3.6 3.4 4.1 4.4 3.9 4.1 4.8 4.9 ...
## $ Mg : num 1.7 1.5 1.3 1.3 1.5 1.6 1.5 1.5 1.8 1.8 ...
## $ Al : num 0 0 0 0 0 0 0 0 0 0 ...
## $ MO : num 32.2 31.7 31.8 30.4 38.4 36.7 37.4 36.4 45.2 42 ...
## $ CTC : num 9.7 8.8 9.5 8.4 8.5 9.1 8.2 9 9.8 9.8 ...
## $ V : num 70 68.6 53.9 58.6 67.7 67.6 66.1 63.4 68.5 69.5 ...
## $ Zn : num 10.2 9.8 8.8 7.6 6.5 7.5 7.3 8.7 9.2 6.5 ...
## $ Cu : num 1.4 1.6 2.3 2.1 1.1 1.3 1.3 1.6 1.3 1.2 ...
## $ Fe : num 63 60 69 69 81 79 74 85 71 82 ...
## $ Mn : num 49.6 35.1 29.3 21.6 31.1 32.5 28.9 29.1 41.3 29.5 ...
## $ S : num 20.2 17.7 16.1 19.1 9.8 10.2 9.1 10.3 12.3 10.9 ...
## $ B : num 0.5 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.6 0.4 ...
## $ cam : chr "C1" "C1" "C1" "C1" ...
## $ SB : num NA NA NA NA NA NA NA NA NA NA ...
db <- da %>%
gather(key = "var", value = "value", pHH2O:SB, -cam)
str(db)
## Classes 'tbl_df', 'tbl' and 'data.frame': 1360 obs. of 7 variables:
## $ trat : chr "1" "1" "1" "1" ...
## $ descr: chr "Soja (PD)" "Soja (PD)" "Soja (PD)" "Soja (PD)" ...
## $ bloc : int 1 2 3 4 1 2 3 4 1 2 ...
## $ prod : num 61 59 61.4 62.3 77.8 82.2 77.7 72 82.6 79.6 ...
## $ cam : chr "C1" "C1" "C1" "C1" ...
## $ var : chr "pHH2O" "pHH2O" "pHH2O" "pHH2O" ...
## $ value: num 6.5 6.5 6.1 6.2 6.5 6.5 6.5 6.4 6.5 6.5 ...
gg1 <- ggplot(data = db,
mapping = aes(x = value,
y = prod,
shape = cam,
color = trat)) +
geom_point() +
facet_wrap(facets = ~var, scales = "free_x")
gg1
gg2 <- gg1 +
geom_smooth(method = lm, se = FALSE)
gg2
ggplot(data = db,
mapping = aes(x = value,
y = prod,
color = cam)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
facet_grid(facets = trat ~ var, scales = "free_x") +
xlab("Valor da variável de solo") +
ylab("Produtividade") +
theme(legend.position = "top")
mp <- mean(da$prod)
# Função que retorna o modelo ajustado e os dados usados para o ajuste
# contendo a resposta descontada do efeito dos tratamentos para fazer os
# gráfico dos valores ajustados contra os observados (corrigidos para o
# efeito de tratamentos).
fit_fun <- function(variable = "K") {
# variable <- "K"
dc <- db %>%
filter(var == variable) %>%
na.omit()
cams <- dc %>%
pull(cam) %>%
unique() %>%
length()
if (cams > 1) {
# Quando o nutriente está em mais de uma camada.
dd <- dc %>% spread(key = "cam", value = "value")
m0 <- lm(prod ~ trat + C1 + C2, data = dd)
dd$resp <- mp + dd$prod -
fitted(update(m0, . ~ trat, data = m0$model))
dd <- dd %>% select(prod, trat, C1, C2, resp)
} else {
# Quando o nutriente está em apenas uma camada.
dd <- dc %>% rename(C = "value")
m0 <- lm(prod ~ trat + C, data = dd)
dd$resp <- mp + dd$prod -
fitted(update(m0, . ~ trat, data = m0$model))
dd <- dd %>% select(prod, trat, C, resp)
}
return(list(m0, dd))
}
# Para criar uma progressão aritmética de uma variável com amplitude
# extendida e retornar em uma matriz. A matriz é usada para obter os
# valores ajustados.
er <- function(data, variable = "C1", length.out = 30) {
mm <- extendrange(data[, variable])
matrix(data = seq(mm[1],
mm[2],
length.out = length.out),
ncol = 1,
dimnames = list(NULL,
variable))
}
# Cria a matriz para ser usada na predição. Essa matriz considera o
# efeito médio dos tratamentos (média ajustada para efeito de
# tratamentos) e os valores variando para o fator representado pela
# matriz `bind`.
pred_matrix <- function(model, bind, columns = 1:10) {
lsm <- LSmatrix(model)[, columns]
lsm <- outer(rep(1, nrow(bind)), lsm, FUN = "*")
X <- cbind(lsm, bind)
return(X)
}
# Retorna os valores ajustados com a banda de confiança.
prediction <- function(model, X) {
beta <- coef(model)
stopifnot(all(names(beta) == colnames(X)))
fit <- X %*% cbind(beta)
v <- vcov(model)
s <- sqrt(diag(X %*% v %*% t(X)))
qtl <- qt(p = 0.975, df = nrow(model$model) - length(beta))
mr <- qtl * s
lwr <- fit - mr
upr <- fit + mr
flu <- cbind(fit, lwr, upr)
colnames(flu) <- c("fit", "lwr", "upr")
return(flu)
}
plot_fit <- function(model, variable = "C1", xlab = "K") {
ds <- model$model
ds$resp <- mp + ds$prod -
fitted(update(model, . ~ trat, data = ds))
pred <- er(ds, variable)
X <- pred_matrix(model, pred)
ci <- prediction(model, X)
grid <- as.data.frame(cbind(pred, ci))
gg_pts <- ggplot(data = ds,
aes(y = resp,
x = get(variable))) +
geom_point() +
geom_hline(yintercept = mp, linetype = 3) +
ylab("Produtividade") +
xlab(xlab)
col <- ifelse(tail(summary(model)$coefficients[, 4], n = 1) < 0.05,
yes = "red",
no = "blue")
gg <- gg_pts +
geom_line(data = grid, mapping = aes(x = get(variable), y = fit), col = col) +
geom_line(data = grid, mapping = aes(x = get(variable), y = lwr), col = "gray") +
geom_line(data = grid, mapping = aes(x = get(variable), y = upr), col = "gray")
return(gg)
}
plot_list <- list()
pH H2O
#-----------------------------------------------------------------------
xlab <- expression(pH ~ H[2] * O)
obj <- fit_fun("pHH2O")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 15.9671 7.803e-09 ***
## C1 1 104.68 104.676 6.3986 0.01733 *
## C2 1 2.98 2.980 0.1822 0.67279
## Residuals 28 458.06 16.359
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.5944 -2.2389 0.0931 2.5653 5.4861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -29.031 35.114 -0.827 0.415124
## trat2 14.367 2.939 4.888 3.47e-05 ***
## trat3 16.678 2.873 5.804 2.73e-06 ***
## trat4A 14.003 2.873 4.873 3.61e-05 ***
## trat4B 13.022 2.903 4.485 0.000106 ***
## trat5A 17.042 2.939 5.798 2.78e-06 ***
## trat6A 15.928 2.873 5.543 5.61e-06 ***
## trat6B 18.667 2.939 6.350 6.12e-07 ***
## trat7 19.525 2.819 6.925 1.30e-07 ***
## trat8 -1.850 2.819 -0.656 0.516886
## C1 14.222 5.543 2.566 0.015720 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.987 on 29 degrees of freedom
## Multiple R-squared: 0.8419, Adjusted R-squared: 0.7874
## F-statistic: 15.45 on 10 and 29 DF, p-value: 4.402e-09
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(pHH2O = gg))
pH CaCl2
#-----------------------------------------------------------------------
xlab <- expression(pH ~ CaCl[2])
obj <- fit_fun("pHCaCl2")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 17.3645 3.026e-09 ***
## C1 1 131.78 131.779 8.7604 0.006203 **
## C2 1 12.74 12.739 0.8469 0.365297
## Residuals 28 421.19 15.043
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.3653 -2.1806 0.2145 2.6272 5.0347
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -31.690 31.268 -1.013 0.31921
## trat2 14.839 2.792 5.315 1.06e-05 ***
## trat3 17.269 2.750 6.281 7.39e-07 ***
## trat4A 13.764 2.792 4.930 3.08e-05 ***
## trat4B 13.139 2.792 4.706 5.74e-05 ***
## trat5A 17.514 2.792 6.273 7.55e-07 ***
## trat6A 15.689 2.792 5.619 4.55e-06 ***
## trat6B 18.308 2.861 6.399 5.37e-07 ***
## trat7 19.940 2.739 7.281 5.11e-08 ***
## trat8 -1.850 2.735 -0.676 0.50417
## C1 16.613 5.598 2.968 0.00596 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.868 on 29 degrees of freedom
## Multiple R-squared: 0.8512, Adjusted R-squared: 0.7999
## F-statistic: 16.59 on 10 and 29 DF, p-value: 1.904e-09
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(pHCaCl2 = gg))
P-mehlich
#-----------------------------------------------------------------------
xlab <- expression("P-Mehlich")
obj <- fit_fun("Pmehlich")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 13.6699 4.311e-08 ***
## C1 1 3.53 3.525 0.1845 0.6708
## C2 1 27.15 27.154 1.4211 0.2432
## Residuals 28 535.03 19.108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.4628 -1.6951 0.4033 1.8445 5.4021
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.3550 6.4162 9.095 5.43e-10 ***
## trat2 17.4993 3.8967 4.491 0.000104 ***
## trat3 18.8534 3.5798 5.267 1.21e-05 ***
## trat4A 16.5615 4.0983 4.041 0.000358 ***
## trat4B 15.8599 3.9838 3.981 0.000421 ***
## trat5A 20.2956 4.0741 4.982 2.67e-05 ***
## trat6A 18.5344 4.1722 4.442 0.000119 ***
## trat6B 21.5534 3.5798 6.021 1.50e-06 ***
## trat7 20.5753 3.9699 5.183 1.53e-05 ***
## trat8 -0.9465 3.7659 -0.251 0.803325
## C1 0.1277 0.2995 0.426 0.672941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.403 on 29 degrees of freedom
## Multiple R-squared: 0.8072, Adjusted R-squared: 0.7408
## F-statistic: 12.14 on 10 and 29 DF, p-value: 6.691e-08
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(Pmehlich = gg))
K
#-----------------------------------------------------------------------
xlab <- expression("Potássio")
obj <- fit_fun("K")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 15.1568 1.393e-08 ***
## C1 1 78.56 78.557 4.5583 0.04165 *
## C2 1 4.61 4.610 0.2675 0.60908
## Residuals 28 482.55 17.234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0509 -1.5108 0.6049 2.8665 5.7019
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.22192 7.99118 5.534 5.76e-06 ***
## trat2 22.35380 3.96571 5.637 4.33e-06 ***
## trat3 25.63613 4.53253 5.656 4.11e-06 ***
## trat4A 21.22301 3.94814 5.375 8.94e-06 ***
## trat4B 21.19026 4.13901 5.120 1.82e-05 ***
## trat5A 25.01164 3.96029 6.316 6.73e-07 ***
## trat6A 25.24663 4.66192 5.415 8.00e-06 ***
## trat6B 26.75251 3.99700 6.693 2.43e-07 ***
## trat7 27.51604 4.69620 5.859 2.35e-06 ***
## trat8 4.03814 3.97656 1.015 0.318
## C1 0.17167 0.07938 2.163 0.039 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.099 on 29 degrees of freedom
## Multiple R-squared: 0.833, Adjusted R-squared: 0.7754
## F-statistic: 14.46 on 10 and 29 DF, p-value: 9.409e-09
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(K = gg))
Ca
#-----------------------------------------------------------------------
xlab <- expression("Cálcio")
obj <- fit_fun("Ca")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 19.3748 8.585e-10 ***
## C1 1 147.28 147.282 10.9244 0.002605 **
## C2 1 40.94 40.938 3.0365 0.092386 .
## Residuals 28 377.49 13.482
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.248 -2.334 0.472 2.377 5.029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.443 7.288 5.275 1.18e-05 ***
## trat2 15.797 2.695 5.862 2.33e-06 ***
## trat3 15.992 2.766 5.782 2.90e-06 ***
## trat4A 15.284 2.686 5.690 3.74e-06 ***
## trat4B 14.097 2.695 5.231 1.33e-05 ***
## trat5A 17.910 2.715 6.597 3.14e-07 ***
## trat6A 14.259 2.855 4.995 2.58e-05 ***
## trat6B 16.866 2.955 5.708 3.56e-06 ***
## trat7 18.260 2.715 6.726 2.22e-07 ***
## trat8 -1.428 2.689 -0.531 0.59933
## C1 5.620 1.759 3.195 0.00336 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.799 on 29 degrees of freedom
## Multiple R-squared: 0.8565, Adjusted R-squared: 0.8071
## F-statistic: 17.31 on 10 and 29 DF, p-value: 1.15e-09
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(Ca = gg))
Mg
#-----------------------------------------------------------------------
xlab <- expression("Magnésio")
obj <- fit_fun("Mg")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 18.9761 1.093e-09 ***
## C1 1 134.79 134.792 9.7923 0.00407 **
## C2 1 45.50 45.498 3.3053 0.07977 .
## Residuals 28 385.42 13.765
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.950 -1.835 0.528 2.603 5.450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.177 8.117 4.580 8.14e-05 ***
## trat2 15.272 2.756 5.541 5.65e-06 ***
## trat3 15.234 2.887 5.277 1.18e-05 ***
## trat4A 14.606 2.739 5.332 1.01e-05 ***
## trat4B 13.981 2.739 5.104 1.90e-05 ***
## trat5A 17.947 2.756 6.512 3.95e-07 ***
## trat6A 14.074 2.935 4.796 4.47e-05 ***
## trat6B 16.705 3.046 5.484 6.61e-06 ***
## trat7 17.887 2.779 6.436 4.86e-07 ***
## trat8 -1.850 2.726 -0.679 0.50270
## C1 16.378 5.438 3.012 0.00534 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.855 on 29 degrees of freedom
## Multiple R-squared: 0.8523, Adjusted R-squared: 0.8013
## F-statistic: 16.73 on 10 and 29 DF, p-value: 1.729e-09
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(Mg = gg))
MO
#-----------------------------------------------------------------------
xlab <- expression("Matéria orgânica")
obj <- fit_fun("MO")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 15.5725 1.032e-08 ***
## C1 1 93.25 93.249 5.5592 0.0256 *
## C2 1 2.80 2.799 0.1669 0.6860
## Residuals 28 469.66 16.774
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0669 -2.0914 0.3871 2.3899 5.4351
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.2866 16.6908 1.275 0.2123
## trat2 9.3330 4.1376 2.256 0.0318 *
## trat3 4.3004 6.4355 0.668 0.5093
## trat4A 9.8926 3.6733 2.693 0.0116 *
## trat4B 8.1988 3.9698 2.065 0.0479 *
## trat5A 11.0650 4.4314 2.497 0.0185 *
## trat6A 7.1653 5.1253 1.398 0.1727
## trat6B 11.3069 4.8878 2.313 0.0280 *
## trat7 11.2578 4.4818 2.512 0.0178 *
## trat8 1.3877 3.1587 0.439 0.6637
## C1 1.2574 0.5256 2.392 0.0234 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.036 on 29 degrees of freedom
## Multiple R-squared: 0.838, Adjusted R-squared: 0.7821
## F-statistic: 15 on 10 and 29 DF, p-value: 6.171e-09
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(MO = gg))
CTC
#-----------------------------------------------------------------------
xlab <- expression("CTC")
obj <- fit_fun("CTC")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 15.9181 8.075e-09 ***
## C1 1 16.84 16.838 1.0261 0.31974
## C2 1 89.41 89.407 5.4485 0.02699 *
## Residuals 28 459.47 16.410
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C1, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C2, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.9728 -1.8226 0.4976 1.5459 6.3476
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.757 11.536 3.013 0.00532 **
## trat2 18.011 2.945 6.115 1.16e-06 ***
## trat3 17.623 2.879 6.122 1.14e-06 ***
## trat4A 16.300 2.896 5.628 4.44e-06 ***
## trat4B 16.470 2.961 5.562 5.33e-06 ***
## trat5A 18.618 2.881 6.462 4.53e-07 ***
## trat6A 15.759 2.953 5.337 9.96e-06 ***
## trat6B 19.766 2.906 6.802 1.81e-07 ***
## trat7 18.173 2.931 6.201 9.19e-07 ***
## trat8 0.218 3.008 0.072 0.94273
## C2 3.182 1.381 2.304 0.02856 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.061 on 29 degrees of freedom
## Multiple R-squared: 0.8361, Adjusted R-squared: 0.7795
## F-statistic: 14.79 on 10 and 29 DF, p-value: 7.279e-09
gg <- plot_fit(m1, "C2", xlab)
gg
plot_list <- append(plot_list, list(CTC = gg))
SB
#-----------------------------------------------------------------------
xlab <- expression("V%")
obj <- fit_fun("V")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 17.1363 3.517e-09 ***
## C1 1 136.66 136.656 8.9652 0.005699 **
## C2 1 2.25 2.253 0.1478 0.703541
## Residuals 28 426.80 15.243
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4719 -2.6466 0.2978 2.4802 5.3754
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.8263 8.8001 3.957 0.000449 ***
## trat2 15.0761 2.7599 5.463 7.02e-06 ***
## trat3 16.8424 2.7511 6.122 1.14e-06 ***
## trat4A 14.0218 2.7587 5.083 2.02e-05 ***
## trat4B 13.6463 2.7462 4.969 2.76e-05 ***
## trat5A 17.6887 2.7635 6.401 5.34e-07 ***
## trat6A 16.0508 2.7532 5.830 2.54e-06 ***
## trat6B 18.4406 2.8285 6.520 3.87e-07 ***
## trat7 19.7121 2.7205 7.246 5.60e-08 ***
## trat8 -1.1952 2.7284 -0.438 0.664586
## C1 0.4157 0.1368 3.039 0.004986 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.846 on 29 degrees of freedom
## Multiple R-squared: 0.8529, Adjusted R-squared: 0.8022
## F-statistic: 16.81 on 10 and 29 DF, p-value: 1.628e-09
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(V = gg))
S
#-----------------------------------------------------------------------
xlab <- expression("Enxofre")
obj <- fit_fun("S")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 14.5678 2.156e-08 ***
## C1 1 25.37 25.372 1.4150 0.2442
## C2 1 38.29 38.287 2.1353 0.1551
## Residuals 28 502.05 17.930
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . -C2, data = m0$model)
summary(m1)
##
## Call:
## lm(formula = prod ~ trat + C1, data = m0$model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0439 -1.4910 -0.0601 2.9415 6.7166
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.5511 6.8824 9.960 7.22e-11 ***
## trat2 12.9843 4.2888 3.028 0.005134 **
## trat3 15.6797 3.6903 4.249 0.000203 ***
## trat4A 13.1507 3.6214 3.631 0.001077 **
## trat4B 11.3156 4.2700 2.650 0.012896 *
## trat5A 17.3285 3.4381 5.040 2.27e-05 ***
## trat6A 14.8880 3.7105 4.012 0.000387 ***
## trat6B 17.2739 4.2950 4.022 0.000377 ***
## trat7 18.0019 3.3196 5.423 7.84e-06 ***
## trat8 -5.3970 4.3076 -1.253 0.220253
## C1 -0.4173 0.3576 -1.167 0.252751
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.317 on 29 degrees of freedom
## Multiple R-squared: 0.8147, Adjusted R-squared: 0.7509
## F-statistic: 12.75 on 10 and 29 DF, p-value: 3.896e-08
gg <- plot_fit(m1, "C1", xlab)
gg
plot_list <- append(plot_list, list(S = gg))
SB
#-----------------------------------------------------------------------
xlab <- expression("Soma de bases")
obj <- fit_fun("SB")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 19.175 6.119e-10 ***
## C 1 170.67 170.675 12.529 0.001373 **
## Residuals 29 395.04 13.622
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = prod ~ trat + C, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.6320 -2.1883 0.6485 2.0438 5.8317
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.487 4.997 8.902 8.61e-10 ***
## trat2 16.030 2.613 6.134 1.10e-06 ***
## trat3 18.804 2.617 7.185 6.57e-08 ***
## trat4A 15.660 2.611 5.998 1.60e-06 ***
## trat4B 16.561 2.657 6.234 8.41e-07 ***
## trat5A 18.588 2.615 7.108 8.04e-08 ***
## trat6A 18.642 2.635 7.074 8.79e-08 ***
## trat6B 20.683 2.610 7.924 9.71e-09 ***
## trat7 19.055 2.613 7.292 4.96e-08 ***
## trat8 -6.429 2.913 -2.207 0.03537 *
## C 4.697 1.327 3.540 0.00137 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.691 on 29 degrees of freedom
## Multiple R-squared: 0.8646, Adjusted R-squared: 0.8178
## F-statistic: 18.51 on 10 and 29 DF, p-value: 5.167e-10
gg <- plot_fit(m0, "C", xlab)
gg
# plot_list <- append(plot_list, list(SB = gg))
Zn
#-----------------------------------------------------------------------
xlab <- expression("Zinco")
obj <- fit_fun("Zn")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 13.9740 2.386e-08 ***
## C 1 23.63 23.631 1.2642 0.2701
## Residuals 29 542.08 18.692
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = prod ~ trat + C, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6762 -1.6756 0.2031 2.2952 5.5631
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 67.0721 5.8790 11.409 3.06e-12 ***
## trat2 15.4192 3.2047 4.811 4.28e-05 ***
## trat3 17.2725 3.1445 5.493 6.45e-06 ***
## trat4A 14.1922 3.2478 4.370 0.000145 ***
## trat4B 13.5165 3.2633 4.142 0.000272 ***
## trat5A 18.7359 3.0820 6.079 1.28e-06 ***
## trat6A 17.1305 3.0634 5.592 4.90e-06 ***
## trat6B 20.5973 3.0625 6.726 2.22e-07 ***
## trat7 19.9810 3.0839 6.479 4.32e-07 ***
## trat8 -4.2311 3.7190 -1.138 0.264559
## C -0.6755 0.6008 -1.124 0.270082
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.323 on 29 degrees of freedom
## Multiple R-squared: 0.8141, Adjusted R-squared: 0.75
## F-statistic: 12.7 on 10 and 29 DF, p-value: 4.071e-08
gg <- plot_fit(m0, "C", xlab)
gg
plot_list <- append(plot_list, list(Zn = gg))
Cu
#-----------------------------------------------------------------------
xlab <- expression("Cobre")
obj <- fit_fun("Cu")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 13.3920 3.83e-08 ***
## C 1 0.07 0.072 0.0037 0.952
## Residuals 29 565.64 19.505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = prod ~ trat + C, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.3861 -1.8093 0.4305 1.9055 5.4525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 61.01597 2.66846 22.866 < 2e-16 ***
## trat2 16.47418 3.15169 5.227 1.35e-05 ***
## trat3 18.07787 3.14408 5.750 3.17e-06 ***
## trat4A 15.39058 3.17392 4.849 3.86e-05 ***
## trat4B 14.83565 3.17760 4.669 6.36e-05 ***
## trat5A 19.13812 3.18140 6.016 1.53e-06 ***
## trat6A 17.34017 3.12708 5.545 5.58e-06 ***
## trat6B 20.77664 3.14649 6.603 3.09e-07 ***
## trat7 19.52746 3.12315 6.252 7.99e-07 ***
## trat8 -1.87950 3.16046 -0.595 0.557
## C -0.04917 0.80981 -0.061 0.952
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.416 on 29 degrees of freedom
## Multiple R-squared: 0.8061, Adjusted R-squared: 0.7392
## F-statistic: 12.05 on 10 and 29 DF, p-value: 7.273e-08
gg <- plot_fit(m0, "C", xlab)
gg
plot_list <- append(plot_list, list(Cu = gg))
Fe
#-----------------------------------------------------------------------
xlab <- expression("Ferro")
obj <- fit_fun("Fe")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 13.391 3.833e-08 ***
## C 1 0.04 0.039 0.002 0.9648
## Residuals 29 565.67 19.506
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = prod ~ trat + C, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.3988 -1.8224 0.4146 1.9551 5.4500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.652324 6.518998 9.304 3.30e-10 ***
## trat2 16.439405 3.407465 4.825 4.13e-05 ***
## trat3 18.070747 3.191547 5.662 4.04e-06 ***
## trat4A 15.413508 3.133660 4.919 3.18e-05 ***
## trat4B 14.746719 3.345062 4.409 0.000131 ***
## trat5A 19.159329 3.142811 6.096 1.22e-06 ***
## trat6A 17.336418 3.137887 5.525 5.91e-06 ***
## trat6B 20.772837 3.182190 6.528 3.78e-07 ***
## trat7 19.504105 3.158149 6.176 9.84e-07 ***
## trat8 -1.895968 3.289709 -0.576 0.568837
## C 0.004179 0.094001 0.044 0.964846
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.417 on 29 degrees of freedom
## Multiple R-squared: 0.806, Adjusted R-squared: 0.7392
## F-statistic: 12.05 on 10 and 29 DF, p-value: 7.279e-08
gg <- plot_fit(m0, "C", xlab)
gg
plot_list <- append(plot_list, list(Fe = gg))
Mn
#-----------------------------------------------------------------------
xlab <- expression("Manganês")
obj <- fit_fun("Mn")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 13.4068 3.784e-08 ***
## C 1 0.70 0.698 0.0358 0.8512
## Residuals 29 565.01 19.483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = prod ~ trat + C, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.6180 -1.7507 0.4018 2.0072 5.4970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.03519 5.19319 11.560 2.23e-12 ***
## trat2 16.59187 3.15867 5.253 1.26e-05 ***
## trat3 18.14397 3.12979 5.797 2.78e-06 ***
## trat4A 15.63302 3.30898 4.724 5.45e-05 ***
## trat4B 14.92337 3.18848 4.680 6.16e-05 ***
## trat5A 19.12382 3.13285 6.104 1.20e-06 ***
## trat6A 17.48124 3.19724 5.468 6.92e-06 ***
## trat6B 20.72913 3.14353 6.594 3.16e-07 ***
## trat7 19.56437 3.12808 6.254 7.94e-07 ***
## trat8 -1.55143 3.49710 -0.444 0.661
## C 0.02625 0.13867 0.189 0.851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.414 on 29 degrees of freedom
## Multiple R-squared: 0.8063, Adjusted R-squared: 0.7395
## F-statistic: 12.07 on 10 and 29 DF, p-value: 7.165e-08
gg <- plot_fit(m0, "C", xlab)
gg
plot_list <- append(plot_list, list(Mn = gg))
B
#-----------------------------------------------------------------------
xlab <- expression("Boro")
obj <- fit_fun("B")
m0 <- obj[[1]]
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 9 2350.88 261.208 14.2142 1.972e-08 ***
## C 1 32.79 32.791 1.7844 0.192
## Residuals 29 532.92 18.377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = prod ~ trat + C, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.140 -1.535 0.150 1.830 6.085
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 53.502 5.956 8.983 7.09e-10 ***
## trat2 16.063 3.049 5.269 1.20e-05 ***
## trat3 17.227 3.101 5.555 5.43e-06 ***
## trat4A 13.678 3.301 4.144 0.000271 ***
## trat4B 13.927 3.101 4.491 0.000104 ***
## trat5A 17.865 3.186 5.608 4.70e-06 ***
## trat6A 16.040 3.186 5.035 2.30e-05 ***
## trat6B 19.490 3.186 6.118 1.15e-06 ***
## trat7 18.652 3.101 6.015 1.53e-06 ***
## trat8 -1.413 3.049 -0.464 0.646407
## C 17.465 13.075 1.336 0.192002
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.287 on 29 degrees of freedom
## Multiple R-squared: 0.8173, Adjusted R-squared: 0.7543
## F-statistic: 12.97 on 10 and 29 DF, p-value: 3.225e-08
gg <- plot_fit(m0, "C", xlab)
gg
plot_list <- append(plot_list, list(B = gg))
Gráfico final
# plot_list <- lapply(plot_list,
# FUN = function(gg) {
# gg + theme(legend.position = "none")
# })
do.call(what = function(...) grid.arrange(..., ncol = 5),
args = plot_list)
