1 Definições da sessão

#-----------------------------------------------------------------------
#                                            Prof. Dr. Walmes M. Zeviani
#                                leg.ufpr.br/~walmes · github.com/walmes
#                                        walmes@ufpr.br · @walmeszeviani
#                      Laboratory of Statistics and Geoinformation (LEG)
#                Department of Statistics · Federal University of Paraná
#                                       2020-dez-10 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------

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

rm(list = objects())

# library(nlme)
library(multcomp)
library(emmeans)
library(tidyverse)
library(readxl)

# source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/pairwise.R")

2 Importação e preparo

#-----------------------------------------------------------------------
# Importação e preparo.

# Leitura da planilha.
tb <- read_xlsx("DADOS TCC- FINAL THIAGO.xlsx",
                sheet = 1,
                skip = 1,
                guess_max = 3000)
str(tb)
## tibble [1,944 × 15] (S3: tbl_df/tbl/data.frame)
##  $ supply  : num [1:1944] 1 1 1 1 1 1 1 1 1 1 ...
##  $ dose    : num [1:1944] 0 0 0 0 0.75 0.75 0.75 0.75 1.5 1.5 ...
##  $ rept    : num [1:1944] 1 2 3 4 1 2 3 4 1 2 ...
##  $ date_   : chr [1:1944] "43073" "43073" "43073" "43073" ...
##  $ date    : chr [1:1944] "2017-04-12" NA NA NA ...
##  $ height  : num [1:1944] NA 315 NA 306 295 340 280 330 NA NA ...
##  $ diameter: num [1:1944] NA 4.7 NA 3.9 3.9 4.6 5.1 4.7 NA NA ...
##  $ fbap    : num [1:1944] NA NA NA NA NA NA NA NA NA NA ...
##  $ dmap    : num [1:1944] NA NA NA NA NA NA NA NA NA NA ...
##  $ frm     : num [1:1944] NA NA NA NA NA NA NA NA NA NA ...
##  $ rdm     : num [1:1944] NA NA NA NA NA NA NA NA NA NA ...
##  $ p       : num [1:1944] NA NA NA NA NA NA NA NA NA NA ...
##  $ k       : num [1:1944] NA NA NA NA NA NA NA NA NA NA ...
##  $ ca      : num [1:1944] NA NA NA NA NA NA NA NA NA NA ...
##  $ mg      : num [1:1944] NA NA NA NA NA NA NA NA NA NA ...
# Elimina variável.
tb$date_ <- NULL

# Preenche a variável data e converte.
tb <- tb %>%
    fill(date, .direction = "down") %>%
    mutate(date = as.Date(date))

# Atribuí rótulos para o fator.
# supply_labels <- c("briquetes de biocarvão e fosfato",
#                    "mistura de biocarvão e fosfato",
#                    "fosfato sem biocarvão")
supply_labels <- c("Biquetes",
                   "Mistura",
                   "Fosfato")
tb <- tb %>%
    mutate(supply = factor(supply, labels = supply_labels))

#
#-----------------------------------------------------------------------
# Dados de crescimento.

# Tabela com dados ao longo do tempo.
tb_cres <- tb %>%
    select(1:diameter) %>%
    drop_na()

# Contagem de celas experimentais.
xtabs(~supply + dose, data = tb_cres)
##           dose
## supply      0 0.75 1.5  3  6 12 24 48 96
##   Biquetes 66   72  52 68 72 72 36 34 72
##   Mistura  72   72  54 72 36 72 54 54 52
##   Fosfato  72   55  72 72 54 54 52 72 54
ggplot(data = tb_cres,
       mapping = aes(x = date, y = height, group = rept)) +
    facet_grid(facets = supply ~ dose) +
    geom_point() +
    geom_line() +
    labs(x = "Date",
         y = "Height (mm)")

ggplot(data = filter(tb_cres, diameter < 20),
       mapping = aes(x = date, y = diameter, group = rept)) +
    facet_grid(facets = supply ~ dose) +
    geom_point() +
    geom_line() +
    labs(x = "Date",
         y = "Diameter (mm)")

#
#-----------------------------------------------------------------------
# Outras variáveis.

# Seleciona a última data onde as demais variáveis são medidas.
tb_vars <- tb %>%
    filter(date == max(date)) %>%
    select(-date)
tb_vars
## # A tibble: 108 x 13
##    supply  dose  rept height diameter  fbap  dmap   frm   rdm     p
##    <fct>  <dbl> <dbl>  <dbl>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 Bique…  0        1    600      8.2  57.8  51.7  5.95  5.25  0.8 
##  2 Bique…  0        2    503      9.4  54.0  28.0  8.2   3.26  0.6 
##  3 Bique…  0        3    468      6.7  45.8  21.8  8.25  3.71  0.96
##  4 Bique…  0        4    645      8    58.6  32.2  9.05  1.85  1.02
##  5 Bique…  0.75     1    545      7.3  58.6  21.6  6.3   1.59  0.72
##  6 Bique…  0.75     2    620      8.1  72.4  29.2  9.34  4.61  0.79
##  7 Bique…  0.75     3    546      9.5  69.2  36.7 12.0   5.19  0.85
##  8 Bique…  0.75     4    610     10.2  95.5  45.5 18.1   4.06  0.84
##  9 Bique…  1.5      1    670      9.2  82.4  37.9  6.48  2.44  0.76
## 10 Bique…  1.5      2     NA     NA    NA    NA   NA    NA    NA   
## # … with 98 more rows, and 3 more variables: k <dbl>, ca <dbl>,
## #   mg <dbl>
tb_vars %>%
    gather(key = "variable", value = "value", height:ca) %>%
    ggplot(mapping = aes(y = value,
                         x = log10(dose + 0.5),
                         group = 1)) +
    # facet_wrap(facets = ~variable + supply, scale = "free_y", ncol = 3) +
    facet_grid(facets = variable ~ supply, scales = "free_y") +
    scale_x_sqrt() +
    geom_point() +
    stat_summary(geom = "line", fun = "mean")

#

3 Análises de cada variável

#-----------------------------------------------------------------------
# Preparo para a análise estatística.

# Transforma variáveis.
tb_vars <- tb_vars %>%
    mutate(Dose = factor(dose),
           # Supply = factor(supply, labels = LETTERS[1:nlevels(supply)]),
           Supply = factor(supply),
           Block = factor(rept),
           doset = log10(dose + 0.5))

# Sequência de valores de dose.
doset_seq <- with(tb_vars, {
    er <- extendrange(doset)
    seq(er[1], er[2], length.out = 51)
})

# Função para termos polinomiais.
f <- function(x, d) {
    poly(x, degree = d, raw = FALSE)
}

# Para classificar o nível de significância.
pvalue_cut <- function(pvalue,
                       breaks = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                       labels = c("***", "**", "*", ".", "ns")) {
    as.character(cut(pvalue,
                     breaks = breaks,
                     labels = labels,
                     include.lowest = TRUE,
                     right = TRUE))
}

#

3.1 Fresh bulk of air part (g plant-1)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Fresh bulk of air part" ~ (g ~ plant^{-1}))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = fbap,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(fbap ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: fbap
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## Block        3  67009 22336.2 17.3719 2.319e-08 ***
## Supply       2    197    98.3  0.0764   0.92649    
## Dose         8  58692  7336.6  5.7060 1.771e-05 ***
## Supply:Dose 16  36057  2253.6  1.7527   0.05887 .  
## Residuals   64  82289  1285.8                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: fbap ~ Block + Supply + f(doset, 3) + Supply:f(doset, 3)
## Model 2: fbap ~ Block + Supply * Dose
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     79 99772                           
## 2     64 82289 15     17482 0.9065  0.561
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: fbap ~ Block + Supply + f(doset, 2) + Supply:f(doset, 2)
## Model 2: fbap ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     82 108736                           
## 2     64  82289 18     26447 1.1427  0.335
# Quadro de anova.
anova(m2)
## Analysis of Variance Table
## 
## Response: fbap
##                    Df Sum Sq Mean Sq F value    Pr(>F)    
## Block               3  67009 22336.2 16.8442 1.296e-08 ***
## Supply              2    197    98.3  0.0741   0.92862    
## f(doset, 2)         2  51221 25610.3 19.3133 1.340e-07 ***
## Supply:f(doset, 2)  4  17082  4270.5  3.2205   0.01656 *  
## Residuals          82 108736  1326.0                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = fbap ~ Block + Supply + f(doset, 2) + Supply:f(doset, 
##     2), data = tb_vars, contrasts = list(Block = contr.sum))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -82.403 -23.870  -2.816  20.288 113.243 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 103.122      6.471  15.937  < 2e-16 ***
## Block1                      -28.841      6.718  -4.293 4.81e-05 ***
## Block2                      -11.535      6.795  -1.697 0.093401 .  
## Block3                       -1.555      6.239  -0.249 0.803802    
## SupplyMistura                -1.520      9.335  -0.163 0.871057    
## SupplyFosfato                -1.134      9.153  -0.124 0.901720    
## f(doset, 2)1                116.154     66.226   1.754 0.083180 .  
## f(doset, 2)2               -247.083     63.836  -3.871 0.000217 ***
## SupplyMistura:f(doset, 2)1   48.636     95.618   0.509 0.612366    
## SupplyFosfato:f(doset, 2)1 -152.146     93.752  -1.623 0.108460    
## SupplyMistura:f(doset, 2)2  -87.636     93.792  -0.934 0.352861    
## SupplyFosfato:f(doset, 2)2  166.330     93.388   1.781 0.078605 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.41 on 82 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.5548, Adjusted R-squared:  0.4951 
## F-statistic:  9.29 on 11 and 82 DF,  p-value: 1.387e-10
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     Supply:poly(doset, degree = 2, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()
## NOTE: A nesting structure was detected in the fitted model:
##     doset %in% Supply
# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])))
tb_coef <- tb_coef %>%
    mutate(x_opt = -b1/(2 * b2),
           y_opt = b0 + b1 * x_opt + b2 * x_opt^2)
tb_coef
##     Supply        b0        b1 sig.b1        b2 sig.b2     x_opt
## 1 Biquetes  78.67522 101.66565    *** -50.05890    *** 1.0154603
## 2  Mistura  67.64016 138.72300    *** -67.81389    *** 1.0228214
## 3  Fosfato 102.39205  23.30043     ns -16.36047     ns 0.7120953
##      y_opt
## 1 130.2939
## 2 138.5846
## 3 110.6881
#

3.2 Dry mass of air part (g plant-1)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Dry mass of air part" ~ (g ~ plant^{-1}))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = dmap,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(dmap ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: dmap
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block        3 13954.0  4651.3 13.4647 6.518e-07 ***
## Supply       2    51.3    25.6  0.0742  0.928551    
## Dose         8 10028.7  1253.6  3.6289  0.001547 ** 
## Supply:Dose 16  8016.0   501.0  1.4503  0.147636    
## Residuals   64 22108.7   345.4                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: dmap ~ Block + Supply + f(doset, 3) + Supply:f(doset, 3)
## Model 2: dmap ~ Block + Supply * Dose
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     79 26535                           
## 2     64 22109 15    4426.2 0.8542 0.6159
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: dmap ~ Block + Supply + f(doset, 2) + Supply:f(doset, 2)
## Model 2: dmap ~ Block + Supply * Dose
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     82 28892                           
## 2     64 22109 18    6783.1 1.0909 0.3815
# Quadro de anova.
anova(m2)
## Analysis of Variance Table
## 
## Response: dmap
##                    Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block               3 13954.0  4651.3 13.2013 4.102e-07 ***
## Supply              2    51.3    25.6  0.0728   0.92988    
## f(doset, 2)         2  8348.5  4174.3 11.8472 3.021e-05 ***
## Supply:f(doset, 2)  4  2913.1   728.3  2.0670   0.09262 .  
## Residuals          82 28891.8   352.3                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = dmap ~ Block + Supply + f(doset, 2) + Supply:f(doset, 
##     2), data = tb_vars, contrasts = list(Block = contr.sum))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.257 -12.413  -0.716   7.469  59.088 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 46.7468     3.3354  14.015  < 2e-16 ***
## Block1                     -13.0510     3.4632  -3.769 0.000308 ***
## Block2                      -5.9818     3.5027  -1.708 0.091468 .  
## Block3                       0.2521     3.2158   0.078 0.937694    
## SupplyMistura                0.1057     4.8120   0.022 0.982528    
## SupplyFosfato                0.6474     4.7178   0.137 0.891192    
## f(doset, 2)1                43.6803    34.1371   1.280 0.204310    
## f(doset, 2)2               -77.2910    32.9055  -2.349 0.021237 *  
## SupplyMistura:f(doset, 2)1  30.2281    49.2881   0.613 0.541379    
## SupplyFosfato:f(doset, 2)1 -59.6739    48.3263  -1.235 0.220427    
## SupplyMistura:f(doset, 2)2 -68.2428    48.3468  -1.412 0.161873    
## SupplyFosfato:f(doset, 2)2  31.5009    48.1385   0.654 0.514697    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.77 on 82 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.4665, Adjusted R-squared:  0.395 
## F-statistic: 6.519 on 11 and 82 DF,  p-value: 1.099e-07
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     Supply:poly(doset, degree = 2, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()
## NOTE: A nesting structure was detected in the fitted model:
##     doset %in% Supply
# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])))
tb_coef <- tb_coef %>%
    mutate(x_opt = -b1/(2 * b2),
           y_opt = b0 + b1 * x_opt + b2 * x_opt^2)
tb_coef
##     Supply       b0       b1 sig.b1         b2 sig.b2     x_opt
## 1 Biquetes 38.26581 32.78839     ** -15.659120      * 1.0469421
## 2  Mistura 31.82977 60.61910    *** -29.485084    *** 1.0279621
## 3  Fosfato 47.12206 13.80489     ns  -9.277048     ns 0.7440346
##      y_opt
## 1 55.42958
## 2 62.98684
## 3 52.25772
#

3.3 Fresh root mass (g plant-1)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Fresh root mass" ~ (g ~ plant^{-1}))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = frm,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(frm ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: frm
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block        3  9607.6  3202.5 18.7917 7.469e-09 ***
## Supply       2   757.8   378.9  2.2232 0.1165601    
## Dose         8  5414.8   676.9  3.9716 0.0007207 ***
## Supply:Dose 16  4304.2   269.0  1.5785 0.1009492    
## Residuals   64 10907.1   170.4                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: frm ~ Block + Supply + f(doset, 3)
## Model 2: frm ~ Block + Supply * Dose
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     85 15408                           
## 2     64 10907 21    4501.4 1.2578  0.238
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: frm ~ Block + Supply + f(doset, 2)
## Model 2: frm ~ Block + Supply * Dose
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     86 15874                           
## 2     64 10907 22      4967 1.3248  0.191
# Quadro de anova.
anova(m2)
## Analysis of Variance Table
## 
## Response: frm
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block        3  9607.6  3202.5 17.3501 6.775e-09 ***
## Supply       2   757.8   378.9  2.0526    0.1346    
## f(doset, 2)  2  4751.9  2376.0 12.8721 1.288e-05 ***
## Residuals   86 15874.1   184.6                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = frm ~ Block + Supply + f(doset, 2), data = tb_vars, 
##     contrasts = list(Block = contr.sum))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.636  -9.359  -2.284   5.524  41.684 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     19.972      2.410   8.287 1.39e-12 ***
## Block1          -9.390      2.496  -3.762 0.000307 ***
## Block2          -5.454      2.535  -2.152 0.034216 *  
## Block3          -1.347      2.327  -0.579 0.564154    
## SupplyMistura    5.908      3.469   1.703 0.092180 .  
## SupplyFosfato    6.418      3.404   1.885 0.062741 .  
## f(doset, 2)1    52.552     14.461   3.634 0.000474 ***
## f(doset, 2)2   -48.345     14.402  -3.357 0.001176 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.59 on 86 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.4878, Adjusted R-squared:  0.4461 
## F-statistic:  11.7 on 7 and 86 DF,  p-value: 2.281e-10
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     poly(doset, degree = 2, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()

# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])))
tb_coef <- tb_coef %>%
    mutate(x_opt = -b1/(2 * b2),
           y_opt = b0 + b1 * x_opt + b2 * x_opt^2)
tb_coef
##     Supply       b0       b1 sig.b1        b2 sig.b2    x_opt    y_opt
## 1 Biquetes 11.80319 23.89532    *** -9.794581     ** 1.219823 26.37722
## 2  Mistura 17.71138 23.89532    *** -9.794581     ** 1.219823 32.28541
## 3  Fosfato 18.22133 23.89532    *** -9.794581     ** 1.219823 32.79536
#

3.4 Root dry mass (g plant-1)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Root dry mass" ~ (g ~ plant^{-1}))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = rdm,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(rdm ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: rdm
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## Block        3 650.93 216.977 18.6040 8.655e-09 ***
## Supply       2  75.30  37.652  3.2283   0.04616 *  
## Dose         8 491.98  61.497  5.2729 4.346e-05 ***
## Supply:Dose 16 396.66  24.792  2.1257   0.01754 *  
## Residuals   64 746.43  11.663                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: rdm ~ Block + Supply + f(doset, 3) + Supply:f(doset, 3)
## Model 2: rdm ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     79 967.68                           
## 2     64 746.43 15    221.25 1.2647 0.2507
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: rdm ~ Block + Supply + f(doset, 2) + Supply:f(doset, 2)
## Model 2: rdm ~ Block + Supply * Dose
##   Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
## 1     82 1107.28                              
## 2     64  746.43 18    360.86 1.7189 0.05882 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Quadro de anova.
anova(m2)
## Analysis of Variance Table
## 
## Response: rdm
##                    Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block               3  650.93 216.977 16.0683 2.644e-08 ***
## Supply              2   75.30  37.652  2.7883   0.06737 .  
## f(doset, 2)         2  403.59 201.796 14.9440 2.926e-06 ***
## Supply:f(doset, 2)  4  124.19  31.048  2.2993   0.06571 .  
## Residuals          82 1107.28  13.503                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = rdm ~ Block + Supply + f(doset, 2) + Supply:f(doset, 
##     2), data = tb_vars, contrasts = list(Block = contr.sum))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1237 -2.2291 -0.7072  1.5457 11.2116 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  6.45210    0.65297   9.881 1.28e-15 ***
## Block1                      -2.55274    0.67798  -3.765 0.000312 ***
## Block2                      -1.40991    0.68572  -2.056 0.042954 *  
## Block3                      -0.06108    0.62955  -0.097 0.922949    
## SupplyMistura                2.39804    0.94204   2.546 0.012782 *  
## SupplyFosfato                1.34668    0.92360   1.458 0.148642    
## f(doset, 2)1                11.89521    6.68295   1.780 0.078792 .  
## f(doset, 2)2               -13.89314    6.44185  -2.157 0.033957 *  
## SupplyMistura:f(doset, 2)1  12.58669    9.64904   1.304 0.195730    
## SupplyFosfato:f(doset, 2)1  -5.68326    9.46074  -0.601 0.549683    
## SupplyMistura:f(doset, 2)2 -12.55522    9.46477  -1.327 0.188350    
## SupplyFosfato:f(doset, 2)2   8.58115    9.42398   0.911 0.365193    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.675 on 82 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.5311, Adjusted R-squared:  0.4682 
## F-statistic: 8.442 on 11 and 82 DF,  p-value: 9.669e-10
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     Supply:poly(doset, degree = 2, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()
## NOTE: A nesting structure was detected in the fitted model:
##     doset %in% Supply
# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])))
tb_coef <- tb_coef %>%
    mutate(x_opt = -b1/(2 * b2),
           y_opt = b0 + b1 * x_opt + b2 * x_opt^2)
tb_coef
##     Supply       b0        b1 sig.b1        b2 sig.b2    x_opt
## 1 Biquetes 4.468681  6.436505     ** -2.814745      * 1.143355
## 2  Mistura 4.865819 12.499739    *** -5.358427    *** 1.166363
## 3  Fosfato 6.851576  2.684305     ns -1.076207     ns 1.247113
##       y_opt
## 1  8.148286
## 2 12.155434
## 3  8.525392
#

3.5 Final plant height (mm)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Final plant height" ~ (mm))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = height,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(height ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: height
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## Block        3 278982   92994 12.1529 2.256e-06 ***
## Supply       2  28988   14494  1.8942   0.15891    
## Dose         8 342136   42767  5.5890 2.352e-05 ***
## Supply:Dose 16 195001   12188  1.5927   0.09721 .  
## Residuals   63 482076    7652                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: height ~ Block + Supply + f(doset, 3)
## Model 2: height ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     84 720369                           
## 2     63 482076 21    238293 1.4829 0.1166
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: height ~ Block + Supply + f(doset, 2)
## Model 2: height ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq     F  Pr(>F)  
## 1     85 747720                             
## 2     63 482076 22    265644 1.578 0.08157 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Quadro de anova.
anova(m2)
## Analysis of Variance Table
## 
## Response: height
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## Block        3 278982   92994 10.5715 5.596e-06 ***
## Supply       2  28988   14494  1.6477    0.1986    
## f(doset, 2)  2 271493  135747 15.4315 1.917e-06 ***
## Residuals   85 747720    8797                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = height ~ Block + Supply + f(doset, 2), data = tb_vars, 
##     contrasts = list(Block = contr.sum))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -332.31  -37.50   13.25   48.44  279.04 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     646.33      16.96  38.113  < 2e-16 ***
## Block1          -68.69      17.28  -3.975 0.000147 ***
## Block2          -16.00      17.83  -0.897 0.372225    
## Block3           11.34      16.11   0.704 0.483593    
## SupplyMistura    31.25      24.13   1.295 0.198889    
## SupplyFosfato   -12.04      23.72  -0.508 0.613111    
## f(doset, 2)1    248.17     100.62   2.466 0.015661 *  
## f(doset, 2)2   -481.53      99.47  -4.841 5.72e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 93.79 on 85 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.4366, Adjusted R-squared:  0.3902 
## F-statistic:  9.41 on 7 and 85 DF,  p-value: 1.365e-08
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     poly(doset, degree = 2, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()

# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])))
tb_coef <- tb_coef %>%
    mutate(x_opt = -b1/(2 * b2),
           y_opt = b0 + b1 * x_opt + b2 * x_opt^2)
tb_coef
##     Supply       b0       b1 sig.b1       b2 sig.b2    x_opt    y_opt
## 1 Biquetes 596.2135 201.0597    *** -97.5586    *** 1.030456 699.8051
## 2  Mistura 627.4623 201.0597    *** -97.5586    *** 1.030456 731.0538
## 3  Fosfato 584.1779 201.0597    *** -97.5586    *** 1.030456 687.7694
#

3.6 Final plant diameter (mm)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Final plant diameter" ~ (mm))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = diameter,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(diameter ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: diameter
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block        3 250.979  83.660 20.7893 1.617e-09 ***
## Supply       2   8.898   4.449  1.1056    0.3373    
## Dose         8 202.204  25.276  6.2809 5.539e-06 ***
## Supply:Dose 16  87.574   5.473  1.3601    0.1908    
## Residuals   64 257.547   4.024                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: diameter ~ Block + Supply + f(doset, 3)
## Model 2: diameter ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     85 386.54                           
## 2     64 257.55 21    128.99 1.5263 0.1001
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: diameter ~ Block + Supply + f(doset, 2)
## Model 2: diameter ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     86 389.89                           
## 2     64 257.55 22    132.34 1.4949 0.1084
# Quadro de anova.
anova(m2)
## Analysis of Variance Table
## 
## Response: diameter
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## Block        3 250.98  83.660 18.4531 2.491e-09 ***
## Supply       2   8.90   4.449  0.9813     0.379    
## f(doset, 2)  2 157.43  78.717 17.3630 4.635e-07 ***
## Residuals   86 389.89   4.534                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = diameter ~ Block + Supply + f(doset, 2), data = tb_vars, 
##     contrasts = list(Block = contr.sum))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5712 -1.0720 -0.0642  1.3214  4.8848 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    10.27198    0.37770  27.196  < 2e-16 ***
## Block1         -1.69868    0.39121  -4.342 3.84e-05 ***
## Block2         -0.91286    0.39725  -2.298    0.024 *  
## Block3          0.11244    0.36476   0.308    0.759    
## SupplyMistura  -0.06053    0.54371  -0.111    0.912    
## SupplyFosfato   0.60728    0.53348   1.138    0.258    
## f(doset, 2)1    3.47431    2.26642   1.533    0.129    
## f(doset, 2)2  -12.65605    2.25714  -5.607 2.45e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.129 on 86 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared:  0.517,  Adjusted R-squared:  0.4777 
## F-statistic: 13.15 on 7 and 86 DF,  p-value: 2.101e-11
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     poly(doset, degree = 2, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()

# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])))
tb_coef <- tb_coef %>%
    mutate(x_opt = -b1/(2 * b2),
           y_opt = b0 + b1 * x_opt + b2 * x_opt^2)
tb_coef
##     Supply       b0       b1 sig.b1        b2 sig.b2     x_opt    y_opt
## 1 Biquetes 9.300711 4.875256    *** -2.564111    *** 0.9506717 11.61809
## 2  Mistura 9.240182 4.875256    *** -2.564111    *** 0.9506717 11.55757
## 3  Fosfato 9.907988 4.875256    *** -2.564111    *** 0.9506717 12.22537
#

3.7 Phosphor (g kg-1)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Phosphor" ~ (g ~ kg^{-1}))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = p,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(p ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: p
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block        3  0.5593 0.18643  5.8589  0.001355 ** 
## Supply       2  0.3964 0.19822  6.2293  0.003399 ** 
## Dose         8 11.2465 1.40582 44.1799 < 2.2e-16 ***
## Supply:Dose 16  0.4412 0.02757  0.8665  0.608389    
## Residuals   63  2.0047 0.03182                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: p ~ Block + Supply + f(doset, 3)
## Model 2: p ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     84 2.7647                           
## 2     63 2.0047 21   0.76003 1.1374 0.3364
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: p ~ Block + Supply + f(doset, 2)
## Model 2: p ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     85 2.8789                           
## 2     63 2.0047 22   0.87426 1.2489 0.2427
# Quadro de anova.
anova(m2)
## Analysis of Variance Table
## 
## Response: p
##             Df  Sum Sq Mean Sq  F value    Pr(>F)    
## Block        3  0.5593  0.1864   5.5043  0.001671 ** 
## Supply       2  0.3964  0.1982   5.8523  0.004157 ** 
## f(doset, 2)  2 10.8134  5.4067 159.6321 < 2.2e-16 ***
## Residuals   85  2.8789  0.0339                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = p ~ Block + Supply + f(doset, 2), data = tb_vars, 
##     contrasts = list(Block = contr.sum))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48074 -0.11388 -0.01588  0.09342  0.51220 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.20959    0.03265  37.046  < 2e-16 ***
## Block1         0.07253    0.03444   2.106  0.03818 *  
## Block2         0.03506    0.03444   1.018  0.31152    
## Block3        -0.02236    0.03158  -0.708  0.48096    
## SupplyMistura -0.11756    0.04700  -2.501  0.01429 *  
## SupplyFosfato -0.15192    0.04658  -3.261  0.00160 ** 
## f(doset, 2)1   3.47201    0.19800  17.535  < 2e-16 ***
## f(doset, 2)2   0.75977    0.19810   3.835  0.00024 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.184 on 85 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.8035, Adjusted R-squared:  0.7873 
## F-statistic: 49.64 on 7 and 85 DF,  p-value: < 2.2e-16
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     poly(doset, degree = 2, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()

# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])))
tb_coef <- tb_coef %>%
    mutate(x_opt = -b1/(2 * b2),
           y_opt = b0 + b1 * x_opt + b2 * x_opt^2)
tb_coef
##     Supply        b0        b1 sig.b1        b2 sig.b2      x_opt
## 1 Biquetes 0.8501586 0.2013533     ** 0.1539287    *** -0.6540475
## 2  Mistura 0.7326016 0.2013533     ** 0.1539287    *** -0.6540475
## 3  Fosfato 0.6982349 0.2013533     ** 0.1539287    *** -0.6540475
##       y_opt
## 1 0.7843113
## 2 0.6667542
## 3 0.6323875
#

3.8 Potassium (g kg-1)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Potassium" ~ (g ~ kg^{-1}))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = k,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(k ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: k
##             Df  Sum Sq Mean Sq F value   Pr(>F)   
## Block        3  21.469  7.1563  3.6908 0.016334 * 
## Supply       2  24.808 12.4042  6.3973 0.002955 **
## Dose         8  15.670  1.9587  1.0102 0.437641   
## Supply:Dose 16  33.422  2.0889  1.0773 0.394776   
## Residuals   63 122.155  1.9390                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 2))
m1 <- update(m0,
             formula = . ~ Block + Supply + f(doset, 1))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: k ~ Block + Supply + f(doset, 3)
## Model 2: k ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq     F Pr(>F)
## 1     84 164.62                          
## 2     63 122.16 21    42.468 1.043 0.4296
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: k ~ Block + Supply + f(doset, 2)
## Model 2: k ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     85 164.65                           
## 2     63 122.16 22    42.493 0.9961 0.4818
anova(m1, m0)
## Analysis of Variance Table
## 
## Model 1: k ~ Block + Supply + f(doset, 1)
## Model 2: k ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     86 169.41                           
## 2     63 122.16 23    47.257 1.0597 0.4124
# Quadro de anova.
anova(m1)
## Analysis of Variance Table
## 
## Response: k
##             Df  Sum Sq Mean Sq F value   Pr(>F)   
## Block        3  21.469  7.1563  3.6328 0.016034 * 
## Supply       2  24.808 12.4042  6.2968 0.002805 **
## f(doset, 1)  1   1.835  1.8348  0.9314 0.337201   
## Residuals   86 169.412  1.9699                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m1)
## 
## Call:
## lm(formula = k ~ Block + Supply + f(doset, 1), data = tb_vars, 
##     contrasts = list(Block = contr.sum))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7942 -0.8804 -0.2616  0.9671  3.9711 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     7.8888     0.2490  31.688  < 2e-16 ***
## Block1          0.8412     0.2626   3.203  0.00191 ** 
## Block2         -0.1034     0.2619  -0.395  0.69402    
## Block3         -0.1829     0.2408  -0.759  0.44970    
## SupplyMistura   0.9271     0.3584   2.587  0.01137 *  
## SupplyFosfato  -0.2928     0.3551  -0.824  0.41198    
## f(doset, 1)    -1.4570     1.5097  -0.965  0.33720    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.404 on 86 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.2212, Adjusted R-squared:  0.1668 
## F-statistic: 4.071 on 6 and 86 DF,  p-value: 0.001219
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     poly(doset, degree = 1, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()

# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)$"), 4])))
tb_coef
##                  Supply       b0         b1 sig.b1
## SupplyBiquetes Biquetes 8.054207 -0.1955607     ns
## SupplyMistura   Mistura 8.981338 -0.1955607     ns
## SupplyFosfato   Fosfato 7.761415 -0.1955607     ns
#

3.9 Calcium (g kg-1)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Calcium" ~ (g ~ kg^{-1}))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = ca,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(ca ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: ca
##             Df  Sum Sq  Mean Sq F value    Pr(>F)    
## Block        3 0.43338 0.144459  4.3332  0.007697 ** 
## Supply       2 0.12280 0.061400  1.8418  0.166971    
## Dose         8 1.31013 0.163767  4.9124 9.636e-05 ***
## Supply:Dose 16 1.16613 0.072883  2.1862  0.014549 *  
## Residuals   63 2.10026 0.033337                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: ca ~ Block + Supply + f(doset, 3) + Supply:f(doset, 3)
## Model 2: ca ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     78 2.9668                              
## 2     63 2.1003 15   0.86657 1.7329 0.06686 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: ca ~ Block + Supply + f(doset, 2) + Supply:f(doset, 2)
## Model 2: ca ~ Block + Supply * Dose
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     81 3.0771                              
## 2     63 2.1003 18   0.97685 1.6279 0.07998 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Quadro de anova.
anova(m2)
## Analysis of Variance Table
## 
## Response: ca
##                    Df  Sum Sq Mean Sq F value    Pr(>F)    
## Block               3 0.43338 0.14446  3.8027   0.01323 *  
## Supply              2 0.12280 0.06140  1.6163   0.20498    
## f(doset, 2)         2 0.98402 0.49201 12.9513 1.317e-05 ***
## Supply:f(doset, 2)  4 0.51540 0.12885  3.3918   0.01286 *  
## Residuals          81 3.07711 0.03799                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = ca ~ Block + Supply + f(doset, 2) + Supply:f(doset, 
##     2), data = tb_vars, contrasts = list(Block = contr.sum))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36120 -0.14695 -0.02167  0.11764  0.44035 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.76814    0.03464  22.178  < 2e-16 ***
## Block1                      0.07625    0.03662   2.082  0.04048 *  
## Block2                      0.04106    0.03649   1.125  0.26371    
## Block3                     -0.11297    0.03345  -3.377  0.00113 ** 
## SupplyMistura              -0.06854    0.04997  -1.372  0.17400    
## SupplyFosfato              -0.06251    0.04940  -1.265  0.20937    
## f(doset, 2)1                0.29366    0.35447   0.828  0.40984    
## f(doset, 2)2               -0.41733    0.34168  -1.221  0.22548    
## SupplyMistura:f(doset, 2)1  1.63878    0.51184   3.202  0.00195 ** 
## SupplyFosfato:f(doset, 2)1  0.20054    0.51112   0.392  0.69582    
## SupplyMistura:f(doset, 2)2 -0.40286    0.50201  -0.802  0.42462    
## SupplyFosfato:f(doset, 2)2  0.14394    0.51332   0.280  0.77988    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1949 on 81 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.4005, Adjusted R-squared:  0.3191 
## F-statistic: 4.919 on 11 and 81 DF,  p-value: 8.474e-06
# Modelo final.
mfinal <- update(m0,
                 formula = . ~ 0 + Supply + Block +
                     Supply:poly(doset, degree = 2, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()
## NOTE: A nesting structure was detected in the fitted model:
##     doset %in% Supply
# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])))
tb_coef <- tb_coef %>%
    mutate(x_opt = -b1/(2 * b2),
           y_opt = b0 + b1 * x_opt + b2 * x_opt^2)
tb_coef
##     Supply        b0        b1 sig.b1          b2 sig.b2    x_opt
## 1 Biquetes 0.7157813 0.1847989     ns -0.08455043     ns 1.092833
## 2  Mistura 0.4428806 0.5451040    *** -0.16616912      * 1.640209
## 3  Fosfato 0.6370733 0.1615739     ns -0.05538851     ns 1.458551
##       y_opt
## 1 0.8167584
## 2 0.8899228
## 3 0.7549051
#

3.10 Magnesium (g kg-1)

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

xlab <- expression(log[10](dose + 0.5))
ylab <- expression("Magnesium" ~ (g ~ kg^{-1}))
clab <- "Supply of P"

ggplot(data = tb_vars,
       mapping = aes(y = mg,
                     x = doset,
                     color = supply)) +
    geom_point() +
    # facet_wrap(facets = ~supply) +
    # stat_summary(fun = "mean", geom = "line") +
    geom_smooth(method = "lm",
                se = FALSE,
                formula = y ~ poly(x, degree = 3)) +
    labs(x = xlab, y = ylab, color = clab)

# Modelo saturado.
m0 <- lm(mg ~ Block + Supply * Dose,
         data = tb_vars,
         contrasts = list(Block = contr.sum))

# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)

layout(1)

# Quadro de anova.
anova(m0)
## Analysis of Variance Table
## 
## Response: mg
##             Df    Sum Sq    Mean Sq F value    Pr(>F)    
## Block        3 0.0008559 0.00028531  1.7550   0.16489    
## Supply       2 0.0003099 0.00015497  0.9533   0.39097    
## Dose         8 0.0147087 0.00183859 11.3096 9.018e-10 ***
## Supply:Dose 16 0.0047115 0.00029447  1.8113   0.04928 *  
## Residuals   63 0.0102419 0.00016257                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo com efeito de dose contínua.
m3 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 3))
m2 <- update(m0,
             formula = . ~ Block + Supply * f(doset, 2))

# Teste da falta de ajuste.
anova(m3, m0)
## Analysis of Variance Table
## 
## Model 1: mg ~ Block + Supply + f(doset, 3) + Supply:f(doset, 3)
## Model 2: mg ~ Block + Supply * Dose
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     78 0.013509                           
## 2     63 0.010242 15 0.0032672 1.3398 0.2064
anova(m2, m0)
## Analysis of Variance Table
## 
## Model 1: mg ~ Block + Supply + f(doset, 2) + Supply:f(doset, 2)
## Model 2: mg ~ Block + Supply * Dose
##   Res.Df      RSS Df Sum of Sq      F  Pr(>F)  
## 1     81 0.016580                              
## 2     63 0.010242 18  0.006338 2.1659 0.01275 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Quadro de anova.
anova(m3)
## Analysis of Variance Table
## 
## Response: mg
##                    Df    Sum Sq   Mean Sq F value    Pr(>F)    
## Block               3 0.0008559 0.0002853  1.6474   0.18529    
## Supply              2 0.0003099 0.0001550  0.8948   0.41285    
## f(doset, 3)         3 0.0137132 0.0045711 26.3929 6.931e-12 ***
## Supply:f(doset, 3)  6 0.0024399 0.0004066  2.3479   0.03887 *  
## Residuals          78 0.0135090 0.0001732                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimativas dos parâmetros.
summary(m2)
## 
## Call:
## lm(formula = mg ~ Block + Supply + f(doset, 2) + Supply:f(doset, 
##     2), data = tb_vars, contrasts = list(Block = contr.sum))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.046078 -0.008807 -0.001758  0.008084  0.037970 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.084687   0.002542  33.310  < 2e-16 ***
## Block1                      0.004229   0.002688   1.573  0.11955    
## Block2                      0.001482   0.002678   0.553  0.58150    
## Block3                     -0.001627   0.002456  -0.662  0.50958    
## SupplyMistura              -0.001306   0.003668  -0.356  0.72264    
## SupplyFosfato               0.002126   0.003626   0.586  0.55924    
## f(doset, 2)1                0.067543   0.026019   2.596  0.01120 *  
## f(doset, 2)2               -0.026038   0.025081  -1.038  0.30228    
## SupplyMistura:f(doset, 2)1  0.100905   0.037571   2.686  0.00878 ** 
## SupplyFosfato:f(doset, 2)1  0.041719   0.037518   1.112  0.26944    
## SupplyMistura:f(doset, 2)2  0.016663   0.036850   0.452  0.65234    
## SupplyFosfato:f(doset, 2)2  0.026744   0.037680   0.710  0.47989    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01431 on 81 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.4622, Adjusted R-squared:  0.3891 
## F-statistic: 6.328 on 11 and 81 DF,  p-value: 1.911e-07
# Modelo final.
mfinal <- update(m2,
                 formula = . ~ 0 + Supply + Block +
                     Supply:poly(doset, degree = 3, raw = TRUE))
# coef(mfinal)

# Determina as curvas médias com bandas de confiança.
grid <- emmeans(mfinal,
                specs = ~Supply + doset,
                at = list(doset = doset_seq)) %>%
    as.data.frame()
## NOTE: A nesting structure was detected in the fitted model:
##     doset %in% Supply
# Curvas ajustadas.
ggplot(data = grid,
       mapping = aes(y = emmean,
                     x = doset,
                     color = Supply)) +
    geom_ribbon(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.1, size = 0.25) +
    geom_line(size = 2) +
    labs(x = xlab, y = ylab, color = clab) +
    theme(legend.position = c(0.975, 0.025), ,
          legend.justification = c(1, 0))

# Termos da equação.
summ <- summary(mfinal)$coefficients
tb_coef <- data.frame(
    Supply = levels(tb_vars$Supply),
    b0 = c(summ[str_detect(names(coef(mfinal)), "^Supply.*\\D$"), 1]),
    b1 = c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 1]),
    "sig b1" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)1$"), 4])),
    b2 = c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 1]),
    "sig b2" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)2$"), 4])),
    b3 = c(summ[str_detect(names(coef(mfinal)), "\\)3$"), 1]),
    "sig b3" = pvalue_cut(c(summ[str_detect(names(coef(mfinal)), "\\)3$"), 4])))
tb_coef
##                  Supply         b0            b1 sig.b1          b2
## SupplyBiquetes Biquetes 0.07481793  0.0131487689     ns 0.006335583
## SupplyMistura   Mistura 0.05741045 -0.0059597393     ns 0.072545206
## SupplyFosfato   Fosfato 0.07262569  0.0001595822     ns 0.031463156
##                sig.b2           b3 sig.b3
## SupplyBiquetes     ns -0.004659424     ns
## SupplyMistura     *** -0.030277979    ***
## SupplyFosfato      ns -0.012656172     ns
#