Adubação nitrogenada em trigo

Código

# message: false
#-----------------------------------------------------------------------
# Carrega pacotes.

library(car)
library(multcomp)
library(emmeans)
library(tidyverse)

Importação e preparo dos dados

Código

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

tb <- readxl::read_excel("Planilha Doses N Trigo x Rendimento.xlsx",
                         sheet = 1)
str(tb)

tibble [25 × 9] (S3: tbl_df/tbl/data.frame)
 $ Trat.        kg N.ha-1
        : chr [1:25] "0" "0" "0" "0" ...
 $ Rep.                            : num [1:25] 1 2 3 4 NA 1 2 3 4 NA ...
 $ Perfilhos.planta-1              : num [1:25] 25 26 23 30 26 ...
 $ Comprimento Espiga (cm)         : num [1:25] 7.85 8.4 6.8 7.7 7.69 ...
 $ Espigueta.Espiga-1              : num [1:25] 24.3 25.2 25.7 23.7 24.7 ...
 $ Massa de mil grãos (g)          : num [1:25] 41.8 36.9 36.4 42.5 39.4 ...
 $ Massa do hectolitro (kg.hL-1)   : num [1:25] 8.15 7.37 7.28 8.2 7.75 7.77 7.81 7.88 7.74 7.8 ...
 $ Produtividade kg.ha-1           : num [1:25] 2.89 2 1.89 3.11 2.47 ...
 $ Produtividade média das Parcelas: num [1:25] 5.21 3.61 3.41 5.6 4.46 ...

Código

# Legenda.
leg <- c("nitr" = "Dose de nitrogênio (kg/ha)",
         "bloc" = "Bloco",
         "perf" = "Perfilhos por planta",
         "comp" = "Comprimento de espiga (cm)",
         "espi" = "Espigueta por espiga",
         "mmil" = "Massa de mil grãos (g)",
         "mhec" = "Massa do hectolitro (kg/hL)",
         "prod" = "Produtividade (kg/ha)",
         "prmp" = "Produtividade média da parcela")

names(tb) <- names(leg)
str(tb)

tibble [25 × 9] (S3: tbl_df/tbl/data.frame)
 $ nitr: chr [1:25] "0" "0" "0" "0" ...
 $ bloc: num [1:25] 1 2 3 4 NA 1 2 3 4 NA ...
 $ perf: num [1:25] 25 26 23 30 26 ...
 $ comp: num [1:25] 7.85 8.4 6.8 7.7 7.69 ...
 $ espi: num [1:25] 24.3 25.2 25.7 23.7 24.7 ...
 $ mmil: num [1:25] 41.8 36.9 36.4 42.5 39.4 ...
 $ mhec: num [1:25] 8.15 7.37 7.28 8.2 7.75 7.77 7.81 7.88 7.74 7.8 ...
 $ prod: num [1:25] 2.89 2 1.89 3.11 2.47 ...
 $ prmp: num [1:25] 5.21 3.61 3.41 5.6 4.46 ...

Código

tb <- tb |>
    filter(nitr != "Média") |>
    mutate(nitr = as.integer(nitr),
           bloc = factor(bloc),
           prod = prod * 1000)

# IMPORTANT: define o tipo de contraste para `bloc`.
contrasts(tb$bloc) <- contr.Sum(unique(tb$bloc))

Análise exploratória

Código

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

# Troca o nome e valor no vetor.
nm <- names(leg) |>
    setNames(unname(leg))

tb |>
    select(-prmp) |>
    pivot_longer(cols = -(1:2),
                 names_to = "var",
                 values_to = "val") |>
    mutate(var = fct_recode(var, !!!nm)) |>
    ggplot(aes(x = nitr, y = val)) +
    facet_wrap(~var, scales = "free_y") +
    geom_point() +
    # geom_point(mapping = aes(color = bloc)) +
    stat_summary(geom = "line", fun = "mean",
                 lty = 2,
                 linewidth = 0.5,
                 color = "gray20") +
    geom_smooth(method = "lm",
                formula = y ~ poly(x, degree = 2),
                se = FALSE) +
    labs(x = leg["nitr"],
         y = "Valores",
         color = "Bloco")

Comprimento e espiga

Código

#-----------------------------------------------------------------------
# Comprimento da espiga.

# Ajuste do modelo saturado aos dados.
# m0 <- lm(comp ~ bloc + poly(nitr, degree = 4), data = tb)
m0 <- lm(comp ~ bloc + factor(nitr), data = tb)

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

Código

layout(1)

# Quadro de anova.
anova(m0)

Analysis of Variance Table

Response: comp
             Df Sum Sq Mean Sq F value Pr(>F)
bloc          3 1.4919 0.49729  1.1781 0.3590
factor(nitr)  4 1.8549 0.46373  1.0986 0.4011
Residuals    12 5.0653 0.42211               

Código

tb_means <-
    emmeans(m0, spec = ~nitr) |>
    cld(Letters = letters, sort = FALSE)
tb_means

 nitr emmean    SE df lower.CL upper.CL .group
    0   7.69 0.325 12     6.98     8.40  a    
   45   7.53 0.325 12     6.82     8.23  a    
   90   8.31 0.325 12     7.60     9.02  a    
  135   8.13 0.325 12     7.42     8.84  a    
  180   8.18 0.325 12     7.47     8.89  a    

Results are averaged over the levels of: bloc 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Código

# Gráfico com as médias.
tb_means |>
    as.data.frame() |>
    mutate_at(".group", trimws) |>
    ggplot(data = _, aes(x = nitr, y = emmean)) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0.2) +
    geom_text(aes(label = sprintf("%0.2f", emmean)),
              hjust = 0,
              nudge_x = 3) +
    geom_point(data = tb,
               mapping = aes(x = nitr, y = comp),
               position = position_nudge(x = -2),
               pch = 1,
               inherit.aes = FALSE) +
    expand_limits(x = c(NA, 200)) +
    labs(x = leg["nitr"],
         y = leg["comp"])

Os resultados da análise de variância sugerem que não há diferenças significativas associadas as doses de nitrogênio. As diferentes doses de nitrogênio não apresentaram efeitos significativos (valor p > 0,05). Portanto, com base nesses dados, não podemos afirmar que a variação no comprimento da espiga seja atribuída as doses de nitrogênio aplicadas.

Mesmo não sendo recomendado pelo teste F ter sido não significativo, apresenta-se o quadro de médias com o teste de Tukey. Sugere-se que a apresentação das médias não contenha as letras do teste de Tukey visto que não há diferenças pelo teste F e que também não é apropriado aplicar comparações de médias para fatores quantiativos.

Perfilhos por planta

Código

#-----------------------------------------------------------------------
# Perfilhos por planta.

# Ajuste do modelo saturado aos dados.
# m0 <- lm(perf ~ bloc + poly(nitr, degree = 4), data = tb)
m0 <- lm(perf ~ bloc + factor(nitr), data = tb)

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

Código

layout(1)

# Quadro de anova.
anova(m0)

Analysis of Variance Table

Response: perf
             Df  Sum Sq Mean Sq F value Pr(>F)
bloc          3   8.163  2.7211  0.2310 0.8730
factor(nitr)  4  50.620 12.6549  1.0742 0.4115
Residuals    12 141.363 11.7803               

Código

tb_means <-
    emmeans(m0, spec = ~nitr) |>
    cld(Letters = letters, sort = FALSE)
tb_means

 nitr emmean   SE df lower.CL upper.CL .group
    0   26.0 1.72 12     22.3     29.7  a    
   45   25.8 1.72 12     22.1     29.5  a    
   90   23.2 1.72 12     19.5     27.0  a    
  135   28.0 1.72 12     24.3     31.7  a    
  180   24.5 1.72 12     20.8     28.2  a    

Results are averaged over the levels of: bloc 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Código

# Gráfico com as médias.
tb_means |>
    as.data.frame() |>
    mutate_at(".group", trimws) |>
    ggplot(data = _, aes(x = nitr, y = emmean)) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0.2) +
    geom_text(aes(label = sprintf("%0.1f", emmean)),
              hjust = 0,
              nudge_x = 3) +
    geom_point(data = tb,
               mapping = aes(x = nitr, y = perf),
               position = position_nudge(x = -2),
               pch = 1,
               inherit.aes = FALSE) +
    expand_limits(x = c(NA, 200)) +
    labs(x = leg["nitr"],
         y = leg["perf"])

Os resultados da análise de variância sugerem que não há diferenças significativas associadas as doses de nitrogênio. As diferentes doses de nitrogênio não apresentaram efeitos significativos (valor p > 0,05). Portanto, com base nesses dados, não podemos afirmar que a variação no número de perfilhos por planta seja atribuída as doses de nitrogênio aplicadas.

Mesmo não sendo recomendado pelo teste F ter sido não significativo, apresenta-se o quadro de médias com o teste de Tukey. Sugere-se que a apresentação das médias não contenha as letras do teste de Tukey visto que não há diferenças pelo teste F e que também não é apropriado aplicar comparações de médias para fatores quantiativos.

Massa de mil grãos

Código

#-----------------------------------------------------------------------
# Massa de mil grãos.

# Ajuste do modelo saturado aos dados.
# m0 <- lm(mmil ~ bloc + poly(nitr, degree = 4), data = tb)
m0 <- lm(mmil ~ bloc + factor(nitr), data = tb)

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

Código

layout(1)

# Quadro de anova.
anova(m0)

Analysis of Variance Table

Response: mmil
             Df Sum Sq Mean Sq F value Pr(>F)
bloc          3  3.878  1.2928  0.1988 0.8952
factor(nitr)  4 19.901  4.9753  0.7649 0.5680
Residuals    12 78.055  6.5046               

Código

tb_means <-
    emmeans(m0, spec = ~nitr) |>
    cld(Letters = letters, sort = FALSE)
tb_means

 nitr emmean   SE df lower.CL upper.CL .group
    0   39.4 1.28 12     36.6     42.2  a    
   45   39.0 1.28 12     36.2     41.8  a    
   90   36.5 1.28 12     33.7     39.3  a    
  135   38.1 1.28 12     35.4     40.9  a    
  180   38.4 1.28 12     35.7     41.2  a    

Results are averaged over the levels of: bloc 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Código

# Gráfico com as médias.
tb_means |>
    as.data.frame() |>
    mutate_at(".group", trimws) |>
    ggplot(data = _, aes(x = nitr, y = emmean)) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0.2) +
    geom_text(aes(label = sprintf("%0.1f %s", emmean, .group)),
              hjust = 0,
              nudge_x = 3) +
    geom_point(data = tb,
               mapping = aes(x = nitr, y = mmil),
               position = position_nudge(x = -2),
               pch = 1,
               inherit.aes = FALSE) +
    expand_limits(x = c(NA, 200)) +
    labs(x = leg["nitr"],
         y = leg["mmil"])

Os resultados da análise de variância sugerem que não há diferenças significativas associadas as doses de nitrogênio. As diferentes doses de nitrogênio não apresentaram efeitos significativos (valor p > 0,05). Portanto, com base nesses dados, não podemos afirmar que a variação na massa de mil grãos seja atribuída as doses de nitrogênio aplicadas.

Mesmo não sendo recomendado pelo teste F ter sido não significativo, apresenta-se o quadro de médias com o teste de Tukey. Sugere-se que a apresentação das médias não contenha as letras do teste de Tukey visto que não há diferenças pelo teste F e que também não é apropriado aplicar comparações de médias para fatores quantiativos.

Massa do hectolitro

Código

#-----------------------------------------------------------------------
# Massa do hectolitro.

# Ajuste do modelo saturado aos dados.
# m0 <- lm(mhec ~ bloc + poly(nitr, degree = 4), data = tb)
m0 <- lm(mhec ~ bloc + factor(nitr), data = tb)

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

Código

layout(1)

# Quadro de anova.
anova(m0)

Analysis of Variance Table

Response: mhec
             Df  Sum Sq  Mean Sq F value Pr(>F)
bloc          3 0.20513 0.068377  0.3930 0.7603
factor(nitr)  4 0.66790 0.166974  0.9597 0.4641
Residuals    12 2.08778 0.173981               

Código

tb_means <-
    emmeans(m0, spec = ~nitr) |>
    cld(Letters = letters, sort = FALSE)
tb_means

 nitr emmean    SE df lower.CL upper.CL .group
    0   7.75 0.209 12     7.30     8.20  a    
   45   7.80 0.209 12     7.35     8.25  a    
   90   7.30 0.209 12     6.84     7.75  a    
  135   7.63 0.209 12     7.17     8.08  a    
  180   7.75 0.209 12     7.29     8.20  a    

Results are averaged over the levels of: bloc 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Código

# Gráfico com as médias.
tb_means |>
    as.data.frame() |>
    mutate_at(".group", trimws) |>
    ggplot(data = _, aes(x = nitr, y = emmean)) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0.2) +
    geom_text(aes(label = sprintf("%0.2f %s", emmean, .group)),
              hjust = 0,
              nudge_x = 3) +
    geom_point(data = tb,
               mapping = aes(x = nitr, y = mhec),
               position = position_nudge(x = -2),
               pch = 1,
               inherit.aes = FALSE) +
    expand_limits(x = c(NA, 200)) +
    labs(x = leg["nitr"],
         y = leg["mhec"])

Os resultados da análise de variância sugerem que não há diferenças significativas associadas as doses de nitrogênio. As diferentes doses de nitrogênio não apresentaram efeitos significativos (valor p > 0,05). Portanto, com base nesses dados, não podemos afirmar que a variação na massa do hectolitro seja atribuída as doses de nitrogênio aplicadas.

Mesmo não sendo recomendado pelo teste F ter sido não significativo, apresenta-se o quadro de médias com o teste de Tukey. Sugere-se que a apresentação das médias não contenha as letras do teste de Tukey visto que não há diferenças pelo teste F e que também não é apropriado aplicar comparações de médias para fatores quantiativos.

Produtividade

Código

#-----------------------------------------------------------------------
# Produtividade.

# Ajuste do modelo saturado aos dados.
# m0 <- lm(prod ~ bloc + poly(nitr, degree = 4), data = tb)
m0 <- lm(prod ~ bloc + factor(nitr), data = tb)

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

Código

layout(1)

# Quadro de anova.
anova(m0)

Analysis of Variance Table

Response: prod
             Df  Sum Sq Mean Sq F value Pr(>F)
bloc          3  504499  168166  1.3741 0.2978
factor(nitr)  4  391348   97837  0.7994 0.5483
Residuals    12 1468571  122381               

Código

tb_means <-
    emmeans(m0, spec = ~nitr) |>
    cld(Letters = letters, sort = FALSE)
tb_means

 nitr emmean  SE df lower.CL upper.CL .group
    0   2474 175 12     2093     2855  a    
   45   2736 175 12     2355     3118  a    
   90   2506 175 12     2125     2887  a    
  135   2822 175 12     2440     3203  a    
  180   2521 175 12     2140     2902  a    

Results are averaged over the levels of: bloc 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Código

# Gráfico com as médias.
tb_means |>
    as.data.frame() |>
    mutate_at(".group", trimws) |>
    ggplot(data = _, aes(x = nitr, y = emmean)) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0.2) +
    geom_text(aes(label = sprintf("%0.0f %s", emmean, .group)),
              hjust = 0,
              nudge_x = 3) +
    geom_point(data = tb,
               mapping = aes(x = nitr, y = prod),
               position = position_nudge(x = -2),
               pch = 1,
               inherit.aes = FALSE) +
    expand_limits(x = c(NA, 200)) +
    labs(x = leg["nitr"],
         y = leg["prod"])

Os resultados da análise de variância sugerem que não há diferenças significativas associadas as doses de nitrogênio. As diferentes doses de nitrogênio não apresentaram efeitos significativos (valor p > 0,05). Portanto, com base nesses dados, não podemos afirmar que a variação na produtividade seja atribuída as doses de nitrogênio aplicadas.

Mesmo não sendo recomendado pelo teste F ter sido não significativo, apresenta-se o quadro de médias com o teste de Tukey. Sugere-se que a apresentação das médias não contenha as letras do teste de Tukey visto que não há diferenças pelo teste F e que também não é apropriado aplicar comparações de médias para fatores quantiativos.

Espiguetas por planta

Código

#-----------------------------------------------------------------------
# Espiguetas por planta.

# Ajuste do modelo saturado aos dados.
# m0 <- lm(espi ~ bloc + poly(nitr, degree = 4), data = tb)
m0 <- lm(espi ~ bloc + factor(nitr), data = tb)

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)

Código

layout(1)

# Quadro de anova.
anova(m0)

Analysis of Variance Table

Response: espi
             Df Sum Sq Mean Sq F value    Pr(>F)    
bloc          3  2.532  0.8439  0.6805    0.5807    
factor(nitr)  4 88.079 22.0198 17.7556 5.592e-05 ***
Residuals    12 14.882  1.2402                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Código

tb_means <-
    emmeans(m0, spec = ~nitr) |>
    cld(Letters = letters, sort = FALSE, reverse = TRUE)
tb_means

 nitr emmean    SE df lower.CL upper.CL .group
    0   24.7 0.557 12     23.5     25.9    c  
   45   28.6 0.557 12     27.4     29.8   b   
   90   29.3 0.557 12     28.0     30.5  ab   
  135   31.2 0.557 12     30.0     32.4  a    
  180   28.5 0.557 12     27.2     29.7   b   

Results are averaged over the levels of: bloc 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 5 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 

Código

# Gráfico com as médias.
tb_means |>
    as.data.frame() |>
    mutate_at(".group", trimws) |>
    ggplot(data = _, aes(x = nitr, y = emmean)) +
    geom_point() +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                  width = 0.2) +
    geom_text(aes(label = sprintf("%0.1f", emmean)),
              hjust = 0,
              nudge_x = 3) +
    geom_point(data = tb,
               mapping = aes(x = nitr, y = espi),
               position = position_nudge(x = -2),
               pch = 1,
               inherit.aes = FALSE) +
    expand_limits(x = c(NA, 200)) +
    labs(x = leg["nitr"],
         y = leg["espi"])

Os resultados da análise de variância sugerem que há diferenças significativas associadas as doses de nitrogênio. As diferentes doses de nitrogênio apresentaram efeitos significativos (valor p < 0,05). Portanto, com base nesses dados, podemos afirmar que a variação no número de espiguetas por por planta seja atribuída as doses de nitrogênio aplicadas.

Apenas para conhecimento das médias dos tratamentos, aplicou-se o teste de Tukey. Sugere-se que a apresentação das médias não contenha as letras do teste de Tukey visto que não é apropriado aplicar comparações de médias para fatores quantiativos. A seguir é feita a análise de regressão para o efeito das doses.

Código

# Ajuste com modelo de regressão.
m1 <- lm(espi ~ bloc + poly(nitr, degree = 2), data = tb)

# Teste da falta de ajuste.
anova(m1, m0)

Analysis of Variance Table

Model 1: espi ~ bloc + poly(nitr, degree = 2)
Model 2: espi ~ bloc + factor(nitr)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     14 21.763                           
2     12 14.882  2     6.881 2.7742 0.1022

Código

# Exame dos pressupostos.
par(mfrow = c(2, 2))
plot(m1)

Código

layout(1)

# Quadro de anova.
anova(m1)

Analysis of Variance Table

Response: espi
                       Df Sum Sq Mean Sq F value    Pr(>F)    
bloc                    3  2.532   0.844  0.5429    0.6608    
poly(nitr, degree = 2)  2 81.198  40.599 26.1173 1.885e-05 ***
Residuals              14 21.763   1.554                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Código

m1 <- lm(espi ~ bloc + nitr + I(nitr^2),
         data = tb,
         contrasts = list(bloc = "contr.sum"))
summary(m1)

Call:
lm(formula = espi ~ bloc + nitr + I(nitr^2), data = tb, contrasts = list(bloc = "contr.sum"))

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1399 -0.8061 -0.0205  0.8900  1.5051 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.472e+01  5.867e-01  42.128 3.78e-16 ***
bloc1       -4.680e-01  4.829e-01  -0.969 0.348901    
bloc2        5.340e-01  4.829e-01   1.106 0.287428    
bloc3       -3.200e-02  4.829e-01  -0.066 0.948101    
nitr         9.807e-02  1.544e-02   6.350 1.80e-05 ***
I(nitr^2)   -4.205e-04  8.228e-05  -5.110 0.000159 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.247 on 14 degrees of freedom
Multiple R-squared:  0.7937,    Adjusted R-squared:   0.72 
F-statistic: 10.77 on 5 and 14 DF,  p-value: 0.0002101

Código

# Coeficientes da esquação estimada.
coef(m1)[c("(Intercept)", "nitr", "I(nitr^2)")] |>
    setNames(c("b0", "b1", "b2"))

           b0            b1            b2 
24.7161428571  0.0980714286 -0.0004204586 

Código

tb_means <-
    emmeans(m1, spec = ~nitr,
            at = list(nitr = seq(0, 180, by = 5)))
tb_means |>
    as.data.frame() |>
    filter(nitr %in% unique(tb$nitr))

 nitr emmean    SE df lower.CL upper.CL
    0   24.7 0.587 14     23.5     26.0
   45   28.3 0.380 14     27.5     29.1
   90   30.1 0.434 14     29.2     31.1
  135   30.3 0.380 14     29.5     31.1
  180   28.7 0.587 14     27.5     30.0

Results are averaged over the levels of: bloc 
Confidence level used: 0.95 

Código

# Gráfico com as médias.
tb_means |>
    as.data.frame() |>
    ggplot(data = _, aes(x = nitr, y = emmean)) +
    geom_line() +
    geom_ribbon(mapping = aes(x = nitr,
                              ymin = lower.CL,
                              ymax = upper.CL),
                alpha = 0.2) +
    annotate("text",
             x = 100,
             y = 25,
             parse = TRUE,
             label = sprintf("y == %0.2f + %0.4f * x - %0.6f * x^2",
                             coef(m1)["(Intercept)"],
                             coef(m1)["nitr"],
                             abs(coef(m1)["I(nitr^2)"]))) +
    geom_point(data = tb,
               mapping = aes(x = nitr, y = espi),
               pch = 1,
               inherit.aes = FALSE) +
    # expand_limits(x = c(NA, 200)) +
    labs(x = leg["nitr"],
         y = leg["espi"])

Código

# Dose para o valor máximo de espigueta por planta.
nitr_max <- -coef(m1)["nitr"]/(2 * coef(m1)["I(nitr^2)"])
nitr_max

    nitr 
116.6244 

Código

# Número de espiguetas por planta.
emmeans(m1, spec = ~nitr,
        at = list(nitr = nitr_max)) |>
    as.data.frame()

 nitr emmean    SE df lower.CL upper.CL
  117   30.4 0.409 14     29.6     31.3

Results are averaged over the levels of: bloc 
Confidence level used: 0.95 

O modelo de segundo grau ajustado apresentou a seguinte equação estimada \[ \hat{y} = 24.72 + 0.0981 x - 0.0004205 x^2 \quad (R^2 = 79,4\%) \] em que \(\hat{y}\) é o número de espiguetas por planta estimado por uma função de segundo grau da dose de nitrogênio \(x\). O teste maior ordem, foi significativo à 5%.

O modelo ajustado aos dados com as bandas de confiança de 95% foi apresentado acima.

A dose de nitrogênio que maximiza o número de espiguetas por plantas é 116.6 kg/ha. O número de espigutas médio nessa condição é 30.4.

Resumo dos resultados

Código

#-----------------------------------------------------------------------
# Resultados resumidos.

v <- c("perf", "comp", "espi", "mmil", "mhec", "prod")

# model <- lm(espi ~ bloc + factor(nitr), data = tb)
# model <- lm(espi ~ bloc + nitr + I(nitr^2), data = tb)
# rm(model)

# Como arredondar.
mag_to_digits <- function(x) {
    mag <- ceiling(log10(x)) |>
        min()
    if (mag > 3) {
        digits <- 0
    } else if (mag >= 0) {
        digits <- 3 - mag
    } else {
        digits <- abs(mag) + 2
    }
    return(digits)
}

# Para resumir o quadro de anova.
mean_square <- function(model) {
    aov <- anova(model) |>
        as.data.frame()
    ms <- aov[, "Mean Sq"]
    pv <- aov[, "Pr(>F)"] |>
        cut(c(0, 0.001, 0.01, 0.05, 1),
            labels = c("***", "**", "*", "ns")) |>
        as.character()
    pv <- ifelse(is.na(pv), "", pv)
    digits <- mag_to_digits(ms)
    fmt <- paste("%0.", digits, "f %s", sep = "")
    tb_aov <-
        data.frame(sprintf(fmt, ms, pv)) |>
        setNames(formula(model)[[2]] |> as.character())
    tb_aov$gl <- aov[, "Df"]
    tb_aov$fonte <- rownames(aov)
    tb_aov <- tb_aov |>
        dplyr::relocate(fonte, gl)
    rownames(tb_aov) <- NULL
    return(tb_aov)
}

# Para retornar a tabela de médias com teste de Tukey.
tukey_means <- function(model, spec = ~nitr) {
    tbm <- emmeans(model, spec = spec) |>
        cld(Letters = letters, sort = FALSE, reverse = TRUE) |>
        as.data.frame() |>
        dplyr::select(1, emmean, .group)
    digits <- mag_to_digits(diff(range(tbm$emmean)))
    # digits <- mag_to_digits(tbm$emmean)
    fmt <- paste("%0.", digits, "f %s", sep = "")
    tbm$mean <- sprintf(fmt, tbm$emmean, trimws(tbm$.group))
    tbm <- tbm |>
        dplyr::select(1, "mean")
    names(tbm)[2] <- formula(model)[[2]] |> as.character()
    return(tbm)
}

# model <- lm(espi ~ bloc + nitr + I(nitr^2), data = tb)
# second_order <- function(model) {
#     coef(model)
#     summary(model)$coefficients[model$assign != 1, ]
# }

# Para retornar os coeficientes do modelo de segundo grau.

model_fits <-
    sapply(v,
           simplify = FALSE,
           function(y) {
               f <- as.formula(sprintf("%s ~ bloc + factor(nitr)", y))
               m0 <- lm(f, data = tb)
               return(m0)
           })

cap <- "Resumo do quadro de análise de variância contendo as fontes de variação, graus de liberdade e quadrados médios com indicação da significância para a estatística F. ns: não significativo; .: significativo a 10%; \\*: signifivativo a 5%; \\*\\*: significativo a 1%; \\*\\*\\*: significativo a 0,1%."
tb_aov <-
    sapply(model_fits,
           simplify = FALSE,
           function(model) {
               mean_square(model)
           }) |>
    Reduce(f = dplyr::left_join,
           x = _) |>
    suppressMessages() |>
    mutate(fonte = c("Bloco", "Dose", "Resíduo"))
knitr::kable(tb_aov, caption = cap)

Resumo do quadro de análise de variância contendo as fontes de variação, graus de liberdade e quadrados médios com indicação da significância para a estatística F. ns: não significativo; .: significativo a 10%; *: signifivativo a 5%; **: significativo a 1%; ***: significativo a 0,1%.

fonte	gl	perf	comp	espi	mmil	mhec	prod
Bloco	3	2.72 ns	0.497 ns	0.844 ns	1.29 ns	0.068 ns	168166 ns
Dose	4	12.65 ns	0.464 ns	22.020 ***	4.98 ns	0.167 ns	97837 ns
Resíduo	12	11.78	0.422	1.240	6.50	0.174	122381

Código

cap <- "Resumo das médias com teste de Tukey. Médias seguidas de mesma letra não diferem entre si ao nível de 5% de significância pelo teste de Tukey."
tb_means <-
    sapply(model_fits,
           simplify = FALSE,
           function(model) {
               tukey_means(model)
           }) |>
    Reduce(f = dplyr::left_join,
           x = _) |>
    suppressMessages()
knitr::kable(tb_means, caption = cap)

Resumo das médias com teste de Tukey. Médias seguidas de mesma letra não diferem entre si ao nível de 5% de significância pelo teste de Tukey.

nitr	perf	comp	espi	mmil	mhec	prod
0	26.00 a	7.687 a	24.71 c	39.38 a	7.750 a	2474 a
45	25.80 a	7.525 a	28.59 b	39.00 a	7.800 a	2736 a
90	23.25 a	8.312 a	29.25 ab	36.49 a	7.298 a	2506 a
135	28.00 a	8.130 a	31.17 a	38.14 a	7.628 a	2821 a
180	24.50 a	8.177 a	28.46 b	38.44 a	7.749 a	2521 a

Código

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

ATENÇÃO: Não é recomendado usar teste de Tukey para fatores qualitativos. Como a maioria das respostas não apresentou efeito do nitrogênio, a tabela com as médias serve apenas como descritivo dos resultados. Recomenda-se omitir as letras do teste de Tukey.