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-23 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(multcomp)
library(emmeans)
library(tidyverse)
# source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/pairwise.R")
Importação e preparo
#-----------------------------------------------------------------------
# Importação.
# Leitura.
tb <- read_csv2("dados_reproducao.csv", col_types = cols())
## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
attr(tb, "spec") <- NULL
str(tb)
## tibble [40 × 14] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ teste: chr [1:40] "A" "A" "A" "A" ...
## $ trat : chr [1:40] "Abamectina" "Abamectina" "Abamectina" "Abamectina" ...
## $ bloco: chr [1:40] "A" "B" "C" "D" ...
## $ alt : num [1:40] 139.4 139.3 57.3 136.3 53.1 ...
## $ dia : num [1:40] 8 6.9 7.2 11.3 6 6.6 7.3 7.6 7.3 8.4 ...
## $ msd : num [1:40] 36.1 21 10.7 46.8 10.3 20.2 28.9 27.2 28.8 27.1 ...
## $ msr : num [1:40] 2.66 1.23 1.1 2.56 1.37 ...
## $ mtt : num [1:40] 38.8 22.2 11.8 49.4 11.7 ...
## $ vol : num [1:40] 30 21 19 27 21 31 30 19 19 35 ...
## $ ig : num [1:40] 2 3 5 1 1 3 2 1 9 4 ...
## $ nes : num [1:40] 7392 19168 28500 14324 1680 ...
## $ ner : num [1:40] 14024 6395 12024 10023 842 ...
## $ ntt : num [1:40] 21416 25563 40524 24347 2522 ...
## $ fre : num [1:40] 4.28 5.11 8.1 4.87 0.5 1.7 1.93 0.65 7.21 5.65 ...
# Converte variáveis para fator.
tb <- tb %>%
mutate_at(c("trat", "bloco", "teste"), "factor") %>%
rename_at(c("trat", "bloco", "teste"), "str_to_title")
tb
## # A tibble: 40 x 14
## Teste Trat Bloco alt dia msd msr mtt vol ig nes
## <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A Abamecti… A 139. 8 36.1 2.66 38.8 30 2 7392
## 2 A Abamecti… B 139. 6.9 21 1.23 22.2 21 3 19168
## 3 A Abamecti… C 57.3 7.2 10.7 1.10 11.8 19 5 28500
## 4 A Abamecti… D 136. 11.3 46.8 2.56 49.4 27 1 14324
## 5 A Carbofur… A 53.1 6 10.3 1.37 11.7 21 1 1680
## 6 A Carbofur… B 126. 6.6 20.2 1.47 21.7 31 3 8240
## 7 A Carbofur… C 102. 7.3 28.9 2.66 31.6 30 2 8359
## 8 A Carbofur… D 141. 7.6 27.2 1.12 28.3 19 1 3024
## 9 A Flex A 115. 7.3 28.8 0.826 29.6 19 9 31342
## 10 A Flex B 96.8 8.4 27.1 2.56 29.7 35 4 18816
## # … with 30 more rows, and 3 more variables: ner <dbl>, ntt <dbl>,
## # fre <dbl>
# Inspeção das frequências.
tb %>%
xtabs(formula = ~Trat + Teste + Bloco) %>%
ftable()
## Bloco A B C D
## Trat Teste
## Abamectina A 1 1 1 1
## B 1 1 1 1
## Carbofurano A 1 1 1 1
## B 1 1 1 1
## Control A 1 1 1 1
## B 1 1 1 1
## Flex A 1 1 1 1
## B 1 1 1 1
## Plus A 1 1 1 1
## B 1 1 1 1
Análise exploratória
#-----------------------------------------------------------------------
# Gráficos.
tb %>%
gather(key = "variable", value = "value", 4:(last_col())) %>%
ggplot(data = .,
mapping = aes(x = Trat, y = value,
color = Teste, group = Teste)) +
facet_wrap(facets = ~variable, scale = "free_y") +
geom_point() +
stat_summary(fun = "mean", geom = "line") +
labs(x = "Treatment",
y = "Response",
color = "Test") +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
#
Fator de reprodução
#-----------------------------------------------------------------------
# Ajuste do modelo com fatores qualitativos (modelo maximal).
# Cria a variável resposta (nome padrão para fácil reuso).
tb$y <- tb$fre %>%
log10()
leg <- list(Trat = "Treatments",
Test = "Test",
resp = expression(log[10] ~ "of reproduction factor"))
# Gráfico.
gg0 <- ggplot(data = tb,
mapping = aes(x = Trat, y = y,
color = Teste, group = Teste)) +
geom_point() +
labs(x = leg$Trat, y = leg$resp, color = leg$Test) +
theme(legend.position = c(0.025, 0.025),
legend.justification = c(0, 0))
gg0 +
stat_summary(fun = "mean", geom = "line")
# Declara o modelo maximal.
m0 <- lm(y ~ Bloco + Trat * Teste, data = tb)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de ANOVA.
anova(m0)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## Bloco 3 0.0717 0.02390 1.2152 0.323190
## Trat 4 4.2513 1.06283 54.0410 1.689e-12 ***
## Teste 1 0.2515 0.25152 12.7889 0.001342 **
## Trat:Teste 4 0.1213 0.03032 1.5416 0.218418
## Residuals 27 0.5310 0.01967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Comparações múltiplas para Trat.
# Teste de Tukey.
agricolae::HSD.test(m0, trt = "Trat")[["groups"]] %>%
rownames_to_column(var = "Trat") %>%
rename("cld" = "groups", "mean" = "y")
## Trat mean cld
## 1 Control 0.989725562 a
## 2 Flex 0.724150022 b
## 3 Plus 0.597161117 b
## 4 Abamectina 0.589119630 b
## 5 Carbofurano -0.005140142 c
# Correção de p-valor conforme o teste de Tukey.
pred <- emmeans(m0, specs = ~Trat) %>%
multcomp::cld(Letters = letters) %>%
as.data.frame() %>%
rename("cld" = ".group") %>%
mutate(cld = str_trim(cld))
## NOTE: Results may be misleading due to involvement in interactions
pred
## Trat emmean SE df lower.CL upper.CL cld
## 1 Carbofurano -0.005140142 0.0495822 27 -0.1068744 0.09659413 a
## 2 Abamectina 0.589119630 0.0495822 27 0.4873854 0.69085390 b
## 3 Plus 0.597161117 0.0495822 27 0.4954268 0.69889539 b
## 4 Flex 0.724150022 0.0495822 27 0.6224157 0.82588429 b
## 5 Control 0.989725562 0.0495822 27 0.8879913 1.09145983 c
# Gráfico com os ICs.
ggplot(data = pred,
mapping = aes(y = reorder(Trat, emmean),
x = emmean,
xmin = lower.CL,
xmax = upper.CL)) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f %s", emmean, cld)),
vjust = 0, nudge_y = 0.05) +
labs(y = leg$Trat, x = leg$resp)
#-----------------------------------------------------------------------
# Comparações múltiplas para Test.
# Teste de Tukey.
agricolae::HSD.test(m0, trt = "Teste")[["groups"]] %>%
rownames_to_column(var = "Teste") %>%
rename("cld" = "groups", "mean" = "y")
## Teste mean cld
## 1 A 0.6583003 a
## 2 B 0.4997061 b
# Correção de p-valor conforme o teste de Tukey.
pred <- emmeans(m0, specs = ~Teste) %>%
multcomp::cld(Letters = letters) %>%
as.data.frame() %>%
rename("cld" = ".group") %>%
mutate(cld = str_trim(cld))
## NOTE: Results may be misleading due to involvement in interactions
pred
## Teste emmean SE df lower.CL upper.CL cld
## 1 B 0.4997061 0.03135854 27 0.4353637 0.5640485 a
## 2 A 0.6583003 0.03135854 27 0.5939579 0.7226427 b
# Gráfico com os ICs.
ggplot(data = pred,
mapping = aes(y = reorder(Teste, emmean),
x = emmean,
xmin = lower.CL,
xmax = upper.CL)) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f %s", emmean, cld)),
vjust = 0, nudge_y = 0.05) +
labs(y = leg$Test, x = leg$resp)
#
Matéria seca de raízes
#-----------------------------------------------------------------------
# Ajuste do modelo com fatores qualitativos (modelo maximal).
# Cria a variável resposta (nome padrão para fácil reuso).
tb$y <- tb$msr# %>%
# log10()
leg <- list(Trat = "Treatments",
Test = "Test",
resp = "Root dry matter")
# Gráfico.
gg0 <- ggplot(data = tb,
mapping = aes(x = Trat, y = y,
color = Teste, group = Teste)) +
geom_point() +
labs(x = leg$Trat, y = leg$resp, color = leg$Test) +
theme(legend.position = c(0.025, 0.025),
legend.justification = c(0, 0))
gg0 +
stat_summary(fun = "mean", geom = "line")
# Declara o modelo maximal.
m0 <- lm(y ~ Bloco + Trat * Teste, data = tb)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de ANOVA.
anova(m0)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## Bloco 3 6.0186 2.00621 2.3705 0.09265 .
## Trat 4 2.6752 0.66879 0.7902 0.54175
## Teste 1 0.6168 0.61678 0.7288 0.40079
## Trat:Teste 4 4.0974 1.02434 1.2103 0.32939
## Residuals 27 22.8510 0.84633
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Comparações múltiplas para Trat.
# Teste de Tukey.
agricolae::HSD.test(m0, trt = "Trat")[["groups"]] %>%
rownames_to_column(var = "Trat") %>%
rename("cld" = "groups", "mean" = "y")
## Trat mean cld
## 1 Carbofurano 2.273500 a
## 2 Plus 2.004625 a
## 3 Abamectina 1.913375 a
## 4 Flex 1.735875 a
## 5 Control 1.502750 a
# Correção de p-valor conforme o teste de Tukey.
pred <- emmeans(m0, specs = ~Trat) %>%
multcomp::cld(Letters = letters) %>%
as.data.frame() %>%
rename("cld" = ".group") %>%
mutate(cld = str_trim(cld))
## NOTE: Results may be misleading due to involvement in interactions
pred
## Trat emmean SE df lower.CL upper.CL cld
## 1 Control 1.502750 0.325256 27 0.8353799 2.170120 a
## 2 Flex 1.735875 0.325256 27 1.0685049 2.403245 a
## 3 Abamectina 1.913375 0.325256 27 1.2460049 2.580745 a
## 4 Plus 2.004625 0.325256 27 1.3372549 2.671995 a
## 5 Carbofurano 2.273500 0.325256 27 1.6061299 2.940870 a
# Gráfico com os ICs.
ggplot(data = pred,
mapping = aes(y = reorder(Trat, emmean),
x = emmean,
xmin = lower.CL,
xmax = upper.CL)) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f %s", emmean, cld)),
vjust = 0, nudge_y = 0.05) +
labs(y = leg$Trat, x = leg$resp)
#-----------------------------------------------------------------------
# Comparações múltiplas para Test.
# Teste de Tukey.
agricolae::HSD.test(m0, trt = "Teste")[["groups"]] %>%
rownames_to_column(var = "Teste") %>%
rename("cld" = "groups", "mean" = "y")
## Teste mean cld
## 1 B 2.01020 a
## 2 A 1.76185 a
# Correção de p-valor conforme o teste de Tukey.
pred <- emmeans(m0, specs = ~Teste) %>%
multcomp::cld(Letters = letters) %>%
as.data.frame() %>%
rename("cld" = ".group") %>%
mutate(cld = str_trim(cld))
## NOTE: Results may be misleading due to involvement in interactions
pred
## Teste emmean SE df lower.CL upper.CL cld
## 1 A 1.76185 0.2057099 27 1.339768 2.183932 a
## 2 B 2.01020 0.2057099 27 1.588118 2.432282 a
# Gráfico com os ICs.
ggplot(data = pred,
mapping = aes(y = reorder(Teste, emmean),
x = emmean,
xmin = lower.CL,
xmax = upper.CL)) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f %s", emmean, cld)),
vjust = 0, nudge_y = 0.05) +
labs(y = leg$Test, x = leg$resp)
#