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á
# 2021-mai-18 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(lme4)
library(lmerTest)
library(emmeans)
library(tidyverse)
library(readxl)
# theme_set(theme_gray(base_family = "Times New Roman"))
Importação e preparo
#-----------------------------------------------------------------------
# Importação.
# Importação de tipagem das variáveis.
tb <- read_excel("braquiária raiz.xlsx",
skip = 1) %>%
mutate_at(c("gen", "blc", "disc", "cort"), "factor") %>%
mutate(enraz = factor(enraz, levels = c("I", "II", "III", "IV", "V", "VI")),
ue = interaction(blc, gen, disc))
str(tb)
## tibble [120 × 10] (S3: tbl_df/tbl/data.frame)
## $ gen : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ blc : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ disc : Factor w/ 2 levels "0","1": 1 1 2 2 1 1 2 2 1 1 ...
## $ cort : Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ dlm : num [1:120] 19.4 21.4 18.3 13.4 12.2 22.1 24.8 15.6 16.9 19.9 ...
## $ dms : num [1:120] 1.7 5.6 1.6 2.3 5.6 6.2 8.6 1.4 0.2 6.9 ...
## $ res : num [1:120] NA 26.2 NA 25 NA 45.8 NA 29.6 NA 26.8 ...
## $ drm : num [1:120] NA 53.4 NA 40.2 NA ...
## $ enraz: Factor w/ 6 levels "I","II","III",..: NA NA NA 2 NA NA NA 1 NA NA ...
## $ ue : Factor w/ 60 levels "1.1.0","2.1.0",..: 1 1 31 31 6 6 36 36 11 11 ...
# Lê a tabela com a legenda.
tb_leg <- read_excel("braquiária raiz.xlsx", sheet = 2) %>%
fill(variavel)
# Aplica os rótulos descritivos.
walk(
unique(tb_leg$variavel),
function(x) {
lv <- tb_leg[tb_leg$variavel == x, ]$codigo
lb <- tb_leg[tb_leg$variavel == x, ]$rotulo
tb[, x, drop = TRUE] <<-
factor(tb[, x, drop = TRUE],
levels = lv,
labels = lb)
}
)
#-----------------------------------------------------------------------
# Passa para inglês.
# str(tb)
# dput(levels(tb$cort))
# dput(levels(tb$disc))
tb <- tb %>%
mutate(
disc = fct_recode(
disc,
"No" = "Ausente",
"Yes" = "Presente"),
cort = fct_recode(
cort,
"First cut on July 24, 2020" = "Primeiro corte em 24 de julho/2020",
"Second cut on August 31, 2020" = "Segundo corte em 31 de agosto/2020"))
Dry leaf mass (g)
#-----------------------------------------------------------------------
# Análise exploratória.
ggplot(data = tb,
mapping = aes(x = gen, y = dlm, color = disc)) +
facet_wrap(facets = ~cort) +
geom_point() +
stat_summary(geom = "line", fun = "mean",
mapping = aes(group = disc)) +
# labs(x = "Genótipo",
# y = "Massa seca de folhas (g)",
# color = "Disco")
labs(x = "Genotype",
y = "Leaf dry mass (g)",
color = "Disc")
# Declaração do modelo apenas para avaliação dos pressupostos.
m0 <- lm(dlm ~ blc + gen * disc * cort, data = tb)
# Avaliação dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
#-----------------------------------------------------------------------
# Modelo de medidas repetidas.
mm0 <- lmer(dlm ~ blc + gen * disc * cort + (1 | ue), data = tb)
anova(mm0)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## blc 10.33 2.582 4 92 0.1026 0.981323
## gen 105.11 21.022 5 92 0.8353 0.527990
## disc 0.11 0.108 1 92 0.0043 0.947913
## cort 163.80 163.800 1 92 6.5082 0.012388 *
## gen:disc 91.26 18.253 5 92 0.7252 0.606222
## gen:cort 372.04 74.407 5 92 2.9564 0.016077 *
## disc:cort 176.18 176.176 1 92 7.0000 0.009587 **
## gen:disc:cort 83.16 16.631 5 92 0.6608 0.654065
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# 1º desdobramento.
# Comparações múltiplas.
emm <- emmeans(m0, specs = ~gen | cort) %>%
multcomp::cld(Letters = letters, adjust = "tukey") %>%
as.data.frame()
emm
## gen cort emmean SE df
## 1 Piatã First cut on July 24, 2020 16.68 1.586446 92
## 3 Paiguás First cut on July 24, 2020 18.10 1.586446 92
## 2 Decumbens First cut on July 24, 2020 19.08 1.586446 92
## 5 Xaraés First cut on July 24, 2020 19.74 1.586446 92
## 6 Marandu First cut on July 24, 2020 20.19 1.586446 92
## 4 Ruzizienses First cut on July 24, 2020 21.58 1.586446 92
## 12 Marandu Second cut on August 31, 2020 12.69 1.586446 92
## 10 Ruzizienses Second cut on August 31, 2020 14.90 1.586446 92
## 8 Decumbens Second cut on August 31, 2020 17.37 1.586446 92
## 7 Piatã Second cut on August 31, 2020 18.36 1.586446 92
## 9 Paiguás Second cut on August 31, 2020 18.71 1.586446 92
## 11 Xaraés Second cut on August 31, 2020 19.32 1.586446 92
## lower.CL upper.CL .group
## 1 12.414297 20.9457 a
## 3 13.834297 22.3657 a
## 2 14.814297 23.3457 a
## 5 15.474297 24.0057 a
## 6 15.924297 24.4557 a
## 4 17.314297 25.8457 a
## 12 8.424297 16.9557 a
## 10 10.634297 19.1657 ab
## 8 13.104297 21.6357 ab
## 7 14.094297 22.6257 ab
## 9 14.444297 22.9757 ab
## 11 15.054297 23.5857 b
# Ordena pelo 2º corte.
lv_ord <- emm %>%
filter(cort == levels(cort)[2]) %>%
arrange(emmean) %>%
pull(gen) %>%
as.character()
lv_ord
## [1] "Marandu" "Ruzizienses" "Decumbens" "Piatã"
## [5] "Paiguás" "Xaraés"
# Aplica a nova order dos nÃveis.
tb_emm <- emm %>%
mutate(gen = fct_relevel(gen, lv_ord))
# Gráfico com segmentos de erro.
ggplot(data = tb_emm,
mapping = aes(y = gen,
x = emmean,
xmin = lower.CL,
xmax = upper.CL,
label = sprintf("%0.1f%s", emmean, .group))) +
facet_wrap(facets = ~cort) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_text(nudge_y = 0.25, vjust = 0, size = 3) +
expand_limits(y = c(NA, nlevels(tb_emm$gen) + 1)) +
# labs(y = "Genótipos",
# x = "Massa seca de folhas (g)")
labs(y = "Genotype",
x = "Leaf dry mass (g)")
#-----------------------------------------------------------------------
# 2º desdobramento.
# Comparações múltiplas.
emm <- emmeans(m0, specs = ~disc | cort) %>%
multcomp::cld(Letters = letters, adjust = "tukey") %>%
as.data.frame()
emm
## disc cort emmean SE df lower.CL
## 1 No First cut on July 24, 2020 18.04667 0.9159352 92 15.96426
## 2 Yes First cut on July 24, 2020 20.41000 0.9159352 92 18.32759
## 4 Yes Second cut on August 31, 2020 15.65000 0.9159352 92 13.56759
## 3 No Second cut on August 31, 2020 18.13333 0.9159352 92 16.05093
## upper.CL .group
## 1 20.12907 a
## 2 22.49241 a
## 4 17.73241 a
## 3 20.21574 a
# Gráfico com segmentos de erro.
ggplot(data = emm,
mapping = aes(y = disc,
x = emmean,
xmin = lower.CL,
xmax = upper.CL,
label = sprintf("%0.1f%s", emmean, .group))) +
facet_wrap(facets = ~cort) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_text(nudge_y = 0.05, vjust = 0, size = 3) +
expand_limits(y = c(NA, nlevels(tb_emm$disc) + 1)) +
# labs(y = "Disco",
# x = "Massa seca de folhas (g)")
labs(x = "Leaf dry mass (g)",
y = "Disc")
#-----------------------------------------------------------------------
Dry mass of stems
#-----------------------------------------------------------------------
# Análise exploratória.
ggplot(data = tb,
mapping = aes(x = gen, y = dms, color = disc)) +
facet_wrap(facets = ~cort) +
geom_point() +
stat_summary(geom = "line", fun = "mean",
mapping = aes(group = disc)) +
# labs(x = "Genótipos",
# y = "Massa seca de ramos (g)",
# color = "Disco")
labs(x = "Genotypes",
y = "Stem dry mass (g)",
color = "Disc")
# NOTE: variância depende do corte.
# Declaração do modelo apenas para avaliação dos pressupostos.
m0 <- lm(dms ~ blc + gen * disc * cort, data = tb)
# Avaliação dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# MASS::boxcox(m0)
# abline(v = 1/3, col = "red")
# ATTENTION: Em função das variâncias serem distintas entre os cortes,
# há dois cursos de ação: 1) fazer a análise separado que é a abordagem
# mais simples; 2) considerar um modelo que estima variâncias diferentes
# por grupo.
#-----------------------------------------------------------------------
# Separado por corte.
fits <- tb %>%
group_split(cort) %>%
map(function(tbi) {
lm(dms ~ blc + gen * disc, data = tbi)
})
names(fits) <- levels(tb$cort)
# Quadros de anova.
map(fits, anova)
## $`First cut on July 24, 2020`
## Analysis of Variance Table
##
## Response: dms
## Df Sum Sq Mean Sq F value Pr(>F)
## blc 4 53.04 13.2606 0.7927 0.5363
## gen 5 17.99 3.5970 0.2150 0.9543
## disc 1 1.01 1.0140 0.0606 0.8067
## gen:disc 5 30.68 6.1356 0.3668 0.8686
## Residuals 44 736.01 16.7276
##
## $`Second cut on August 31, 2020`
## Analysis of Variance Table
##
## Response: dms
## Df Sum Sq Mean Sq F value Pr(>F)
## blc 4 8.512 2.128 0.6699 0.61635
## gen 5 262.673 52.535 16.5377 3.62e-09 ***
## disc 1 12.696 12.696 3.9967 0.05179 .
## gen:disc 5 3.826 0.765 0.2409 0.94215
## Residuals 44 139.773 3.177
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Efeito do genótipo.
# Comparações múltiplas.
emm <-
map(fits,
function(fit_i) {
emmeans(fit_i, specs = ~gen) %>%
multcomp::cld(Letters = letters, adjust = "tukey") %>%
as.data.frame()
}) %>%
bind_rows(.id = "cort") %>%
mutate(cort = factor(cort, levels(tb$cort)))
# Ordena pelo 2º corte.
lv_ord <- emm %>%
filter(cort == levels(cort)[2]) %>%
arrange(emmean) %>%
pull(gen) %>%
as.character()
lv_ord
## [1] "Marandu" "Xaraés" "Piatã" "Ruzizienses"
## [5] "Decumbens" "Paiguás"
# Aplica a nova order dos nÃveis.
tb_emm <- emm %>%
mutate(gen = fct_relevel(gen, lv_ord))
# Gráfico com segmentos de erro.
ggplot(data = tb_emm,
mapping = aes(y = gen,
x = emmean,
xmin = lower.CL,
xmax = upper.CL,
label = sprintf("%0.1f%s", emmean, .group))) +
facet_wrap(facets = ~cort) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_text(nudge_y = 0.25, vjust = 0, size = 3) +
expand_limits(y = c(NA, nlevels(tb_emm$gen) + 1)) +
# labs(y = "Genótipos",
# x = "Massa seca de ramos (g)")
labs(y = "Genotypes",
x = "Stem dry mass (g)")
#-----------------------------------------------------------------------
# Efeito do disco.
# Comparações múltiplas.
emm <-
map(fits,
function(fit_i) {
emmeans(fit_i, specs = ~disc) %>%
multcomp::cld(Letters = letters, adjust = "tukey") %>%
as.data.frame()
}) %>%
bind_rows(.id = "cort") %>%
mutate(cort = factor(cort, levels(tb$cort)))
# Gráfico com segmentos de erro.
ggplot(data = emm,
mapping = aes(y = disc,
x = emmean,
xmin = lower.CL,
xmax = upper.CL,
label = sprintf("%0.1f%s", emmean, .group))) +
facet_wrap(facets = ~cort) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_text(nudge_y = 0.25, vjust = 0, size = 3) +
expand_limits(y = c(NA, nlevels(tb_emm$disc) + 1)) +
# labs(y = "Disco",
# x = "Massa seca de ramos (g)")
labs(y = "Disc",
x = "Stem dry mass (g)")
#-----------------------------------------------------------------------
Dry root mass (g)
#-----------------------------------------------------------------------
# Análise exploratória.
tb_cort2 <- tb %>%
filter(cort == levels(cort)[2])
ggplot(data = tb_cort2,
mapping = aes(x = gen, y = drm, color = disc)) +
geom_point() +
stat_summary(geom = "line", fun = "mean",
mapping = aes(group = disc)) +
# labs(x = "Genótipo",
# y = "Massa seca de raÃzes (g)",
# color = "Disco")
labs(x = "Genotypes",
y = "Root dry mass (g)",
color = "Disc")
# Declaração do modelo.
m0 <- lm(drm ~ blc + gen * disc, data = tb)
# Avaliação dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Transformação Box-Cox.
MASS::boxcox(m0)
abline(v = 0, col = "red")
# Modelo com as variáveis transformadas.
m0 <- update(m0, formula = log10(.) ~ .)
# Avaliação dos pressupostos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: log10(drm)
## Df Sum Sq Mean Sq F value Pr(>F)
## blc 4 0.10348 0.025870 2.1882 0.085911 .
## gen 5 0.24623 0.049246 4.1654 0.003482 **
## disc 1 0.29976 0.299764 25.3554 8.577e-06 ***
## gen:disc 5 0.03954 0.007908 0.6689 0.649087
## Residuals 44 0.52019 0.011823
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Efeito de genótipo.
# Comparações múltiplas.
emm <- emmeans(m0, specs = ~gen) %>%
multcomp::cld(Letters = letters, adjust = "tukey") %>%
as.data.frame()
emm
## gen emmean SE df lower.CL upper.CL .group
## 3 Paiguás 1.763211 0.03438387 44 1.668495 1.857926 a
## 1 Piatã 1.793537 0.03438387 44 1.698821 1.888252 a
## 2 Decumbens 1.844737 0.03438387 44 1.750022 1.939452 ab
## 6 Marandu 1.846222 0.03438387 44 1.751507 1.940937 ab
## 5 Xaraés 1.857968 0.03438387 44 1.763253 1.952683 ab
## 4 Ruzizienses 1.967953 0.03438387 44 1.873238 2.062668 b
# Volta pra escala original.
emm <- emm %>%
mutate_at(c("emmean", "lower.CL", "upper.CL"), ~10^.)
emm
## gen emmean SE df lower.CL upper.CL .group
## 3 Paiguás 57.97097 0.03438387 44 46.61176 72.09841 a
## 1 Piatã 62.16366 0.03438387 44 49.98290 77.31284 a
## 2 Decumbens 69.94188 0.03438387 44 56.23700 86.98661 ab
## 6 Marandu 70.18140 0.03438387 44 56.42960 87.28451 ab
## 5 Xaraés 72.10539 0.03438387 44 57.97659 89.67737 ab
## 4 Ruzizienses 92.88655 0.03438387 44 74.68575 115.52286 b
# Ordena pelo 2º corte.
lv_ord <- emm %>%
arrange(emmean) %>%
pull(gen) %>%
as.character()
lv_ord
## [1] "Paiguás" "Piatã" "Decumbens" "Marandu"
## [5] "Xaraés" "Ruzizienses"
# Aplica a nova order dos nÃveis.
tb_emm <- emm %>%
mutate(gen = fct_relevel(gen, lv_ord))
# Gráfico com segmentos de erro.
ggplot(data = tb_emm,
mapping = aes(y = gen,
x = emmean,
xmin = lower.CL,
xmax = upper.CL,
label = sprintf("%0.2f%s", emmean, .group))) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_text(nudge_y = 0.25, vjust = 0, size = 3) +
expand_limits(y = c(NA, nlevels(tb_emm$gen) + 1)) +
# labs(y = "Genótipos",
# x = "Massa seca de raÃzes (g)")
labs(y = "Genotypes",
x = "Root dry mass (g)")
#-----------------------------------------------------------------------
# Efeito de disco.
# Comparações múltiplas.
emm <- emmeans(m0, specs = ~disc) %>%
multcomp::cld(Letters = letters, adjust = "tukey") %>%
as.data.frame()
emm
## disc emmean SE df lower.CL upper.CL .group
## 2 Yes 1.774922 0.01985154 44 1.728959 1.820884 a
## 1 No 1.916287 0.01985154 44 1.870325 1.962250 b
# Volta pra escala original.
emm <- emm %>%
mutate_at(c("emmean", "lower.CL", "upper.CL"), ~10^.)
emm
## disc emmean SE df lower.CL upper.CL .group
## 2 Yes 59.55546 0.01985154 44 53.57462 66.20398 a
## 1 No 82.46837 0.01985154 44 74.18650 91.67478 b
# Gráfico com segmentos de erro.
ggplot(data = emm,
mapping = aes(y = disc,
x = emmean,
xmin = lower.CL,
xmax = upper.CL,
label = sprintf("%0.1f%s", emmean, .group))) +
geom_errorbarh(height = 0.1) +
geom_point() +
geom_text(nudge_y = 0.05, vjust = 0, size = 3) +
expand_limits(y = c(NA, nlevels(tb_emm$disc) + 1)) +
# labs(y = "Disco",
# x = "Massa seca de raÃzes (g)")
labs(y = "Disc",
x = "Root dry mass (g)")
#-----------------------------------------------------------------------
Grau de enraizamento
#-----------------------------------------------------------------------
library(MASS)
# Contagem das ocorrências.
tb %>%
drop_na(enraz) %>%
xtabs(formula = ~gen + enraz) %>%
addmargins()
## enraz
## gen I II III IV V VI Sum
## Piatã 0 1 1 1 1 1 5
## Decumbens 1 3 0 0 0 1 5
## Paiguás 0 2 0 0 2 1 5
## Ruzizienses 0 1 1 0 0 3 5
## Xaraés 0 0 0 0 1 4 5
## Marandu 1 0 0 0 1 3 5
## Sum 2 7 2 1 5 13 30
# Regressão logÃstica ordinal.
m0 <- polr(enraz ~ blc + gen,
data = filter(tb, !is.na(enraz)),
Hess = TRUE)
m_null <- update(m0, formula = . ~ blc)
# class(m0)
# methods(class = "polr")
# Teste para o efeito dos genótipos.
anova(m0, m_null)
## Likelihood ratio tests of ordinal regression models
##
## Response: enraz
## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 blc 21 85.03433
## 2 blc + gen 16 74.59243 1 vs 2 5 10.4419 0.06363975
# Resumo do modelo.
summary(m0)
## Call:
## polr(formula = enraz ~ blc + gen, data = filter(tb, !is.na(enraz)),
## Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## blc2 -0.7086 1.217 -0.5822
## blc3 1.3993 1.190 1.1762
## blc4 0.8801 1.123 0.7835
## blc5 1.2265 1.152 1.0644
## genDecumbens -1.4766 1.243 -1.1880
## genPaiguás 0.1206 1.083 0.1114
## genRuzizienses 1.2288 1.221 1.0067
## genXaraés 2.8072 1.489 1.8847
## genMarandu 1.1159 1.250 0.8929
##
## Intercepts:
## Value Std. Error t value
## I|II -2.4031 1.2934 -1.8580
## II|III 0.1049 1.0956 0.0957
## III|IV 0.6110 1.1071 0.5519
## IV|V 0.8250 1.1162 0.7392
## V|VI 1.7473 1.1596 1.5067
##
## Residual Deviance: 74.59243
## AIC: 102.5924
# Comparações múltiplas.
emm <- emmeans(m0, specs = ~gen) %>%
multcomp::cld(Letters = letters, adjust = "tukey", alpha = 0.1)
emm
## gen emmean SE df asymp.LCL asymp.UCL .group
## Decumbens -1.094 0.930 Inf -3.541 1.35 a
## Piatã 0.382 0.789 Inf -1.694 2.46 ab
## Paiguás 0.503 0.751 Inf -1.472 2.48 ab
## Marandu 1.498 0.997 Inf -1.125 4.12 ab
## Ruzizienses 1.611 0.936 Inf -0.850 4.07 ab
## Xaraés 3.190 1.253 Inf -0.106 6.49 b
##
## Results are averaged over the levels of: blc
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 6 estimates
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.1
tb_emm <- emm %>%
as.data.frame() %>%
mutate(gen = fct_reorder(gen, emmean),
.group = str_trim(.group))
tb_emm
## gen emmean SE df asymp.LCL asymp.UCL .group
## 2 Decumbens -1.0941279 0.9299108 Inf -3.5407588 1.352503 a
## 1 Piatã 0.3824320 0.7892570 Inf -1.6941334 2.458997 ab
## 3 Paiguás 0.5030502 0.7507362 Inf -1.4721653 2.478266 ab
## 6 Marandu 1.4983096 0.9972526 Inf -1.1255002 4.122119 ab
## 4 Ruzizienses 1.6112543 0.9355168 Inf -0.8501263 4.072635 ab
## 5 Xaraés 3.1896723 1.2527798 Inf -0.1064395 6.485784 b
ggplot(data = tb_emm,
mapping = aes(y = gen,
x = emmean,
xmin = asymp.LCL,
xmax = asymp.UCL,
label = sprintf("%0.2f %s", emmean, .group))) +
geom_errorbarh(height = 0.2) +
geom_point() +
geom_text(nudge_y = 0.25, vjust = 0, size = 3) +
expand_limits(y = c(NA, nlevels(tb_emm$gen) + 1)) +
# labs(y = "Genótipos",
# x = "Grau de enraizamento")
labs(y = "Genotypes",
x = "Rooting score")
#-----------------------------------------------------------------------
# Probabilidade estimada (média nos blocos).
tb_pred <- tb %>%
distinct(gen, blc) %>%
complete()
tb_pred <-
predict(m0, newdata = tb_pred, type = "prob") %>%
as.data.frame() %>%
bind_cols(tb_pred, .) %>%
group_by(gen) %>%
summarise_at(vars(I:VI), mean)
# tb_pred
tb_pred <- tb_pred %>%
pivot_longer(cols = -gen,
names_to = "grau",
values_to = "prob") %>%
mutate(gen = factor(gen, levels = tb_emm$gen),
grau = factor(grau, levels = levels(tb$enraz)))
# tb_pred
ggplot(data = tb_pred,
mapping = aes(y = gen, x = grau, fill = prob)) +
geom_tile() +
scale_fill_distiller(direction = 1) +
geom_text(mapping = aes(label = sprintf("%0.3f", prob))) +
# labs(y = "Gentótipos",
# x = "Grau de enraizamento",
# fill = "Probabilidade\nestimada média")
labs(y = "Genotypes",
x = "Rooting score",
fill = "Estimated mean\nprobability")
#-----------------------------------------------------------------------