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-21 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
# library(multcomp)
# library(nlme)
# library(drc)
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_penetracao.csv", col_types = cols())
## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
attr(tb, "spec") <- NULL
str(tb)
## tibble [48 × 14] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ prod : chr [1:48] "Flex" "Flex" "Flex" "Flex" ...
## $ conc : num [1:48] 0 0 0 0.25 0.25 0.25 0.75 0.75 0.75 1.5 ...
## $ bloco: chr [1:48] "A" "B" "C" "A" ...
## $ alt : num [1:48] 27.9 28.8 26.4 30.8 30.8 30.8 36.4 32.4 24.5 30.6 ...
## $ dia : num [1:48] 0.41 0.44 0.45 0.41 0.38 0.41 0.38 0.31 0.44 0.31 ...
## $ msd : num [1:48] 1.03 0.86 1.01 0.94 0.88 0.79 1.2 0.96 0.76 1.11 ...
## $ msr : num [1:48] 1.03 1.25 1.21 1.26 1.08 ...
## $ mst : num [1:48] 2.06 2.11 2.22 2.2 1.96 ...
## $ ngal : num [1:48] 258 171 196 178 139 193 99 110 138 52 ...
## $ gal : num [1:48] 5 5 5 5 5 5 4 5 5 3 ...
## $ j2s : num [1:48] 745 500 916 228 375 464 88 34 95 16 ...
## $ ovo : num [1:48] 624 827 571 608 523 284 92 228 253 104 ...
## $ ntt : num [1:48] 1369 1327 1487 836 898 ...
## $ fre : num [1:48] 13.69 13.27 14.87 8.36 8.98 ...
# Converte variáveis para fator.
tb <- tb %>%
mutate_at(c("prod", "bloco"), "factor") %>%
rename_at(c("prod", "bloco"), "str_to_title") %>%
mutate(Conc = factor(conc))
tb
## # A tibble: 48 x 15
## Prod conc Bloco alt dia msd msr mst ngal gal j2s
## <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Flex 0 A 27.9 0.41 1.03 1.03 2.06 258 5 745
## 2 Flex 0 B 28.8 0.44 0.86 1.25 2.11 171 5 500
## 3 Flex 0 C 26.4 0.45 1.01 1.21 2.22 196 5 916
## 4 Flex 0.25 A 30.8 0.41 0.94 1.26 2.2 178 5 228
## 5 Flex 0.25 B 30.8 0.38 0.88 1.08 1.96 139 5 375
## 6 Flex 0.25 C 30.8 0.41 0.79 1.54 2.33 193 5 464
## 7 Flex 0.75 A 36.4 0.38 1.2 1.09 2.29 99 4 88
## 8 Flex 0.75 B 32.4 0.31 0.96 1.70 2.66 110 5 34
## 9 Flex 0.75 C 24.5 0.44 0.76 0.884 1.64 138 5 95
## 10 Flex 1.5 A 30.6 0.31 1.11 1.06 2.17 52 3 16
## # … with 38 more rows, and 4 more variables: ovo <dbl>, ntt <dbl>,
## # fre <dbl>, Conc <fct>
# Inspeção das frequências.
tb %>%
xtabs(formula = ~Prod + Conc + Bloco) %>%
ftable()
## Bloco A B C
## Prod Conc
## Flex 0 2 2 2
## 0.25 1 1 1
## 0.5 1 1 1
## 0.75 1 1 1
## 1 1 1 1
## 1.5 1 1 1
## 2 1 1 1
## Plus 0 2 2 2
## 0.25 1 1 1
## 0.5 1 1 1
## 0.75 1 1 1
## 1 1 1 1
## 1.5 1 1 1
## 2 1 1 1
Análise exploratória
#-----------------------------------------------------------------------
# Gráficos.
tb %>%
gather(key = "variable", value = "value", 4:(last_col() - 1)) %>%
ggplot(data = .,
mapping = aes(x = conc, y = value, color = Prod)) +
facet_wrap(facets = ~variable, scale = "free_y") +
geom_point() +
stat_summary(fun = "mean", geom = "line") +
labs(x = "Concentration",
y = "Response")
# {fre, j2s, ntt, ovo}: tem função média com decaimento tipo exponencial
# e relação média variância. Talvez tenha que se analisar a respostra
# transformada para o atendimento dos pressupostos.
# As demais variáveis apresentam maior variação e nenhuma tendência
# específica que sugira a adoção de algum modelo não linear. Talvez
# empregar polinômios seja o suficiente para tirar conclusões.
tb %>%
select(1:3, fre, j2s, ntt, ovo) %>%
gather(key = "variable", value = "value", 4:last_col()) %>%
ggplot(data = .,
mapping = aes(x = conc, y = value, color = Prod)) +
facet_wrap(facets = ~variable, scale = "free_y") +
scale_y_log10() +
geom_point() +
stat_summary(fun = "mean", geom = "line") +
labs(x = "Concentration",
y = "Response (log)")
# Como esperado, na escala log(), o comportamente se torna linear e com
# variância homogênea.
# ngal: galhas.
# fre: fator de reprodução.
# Uma versão de nome curto para a função `poly()`.
f <- function(x, d = 1, ...) {
poly(x = x, degree = d, raw = TRUE, ...)
}
# Para sinalizar o nível de significância.
signif_cut <- function(pvalue,
breaks = c(0, 0.001, 0.01, 0.05, 0.1, 1),
labels = c("***", "**", "*", ".", "ns")) {
cut(x = pvalue,
breaks = breaks,
labels = labels,
include.lowest = TRUE,
right = TRUE)
}
# Sequência de valores de concentração para uso na predição.
conc_seq <- with(tb, seq(0, 1.1 * max(conc), length.out = 37))
Número de galhas
#-----------------------------------------------------------------------
# Ajuste do modelo com fatores qualitativos (modelo maximal).
# Cria a variável resposta (nome padrão para fácil reuso).
tb$y <- tb$ngal %>%
log10()
# Gráfico.
gg0 <- ggplot(data = tb,
mapping = aes(x = conc, y = y, color = Prod)) +
geom_point() +
labs(x = "Concentration",
y = "Response",
color = "Product") +
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 + Prod * Conc, 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 2 0.06035 0.030175 2.2479 0.122044
## Prod 1 0.14051 0.140515 10.4675 0.002822 **
## Conc 6 1.83770 0.306283 22.8163 2.848e-10 ***
## Prod:Conc 6 0.16476 0.027460 2.0456 0.088132 .
## Residuals 32 0.42956 0.013424
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Ajuste do modelo reduzido com `conc` quantitativo.
# Modelo com polinômio de 2 grau.
m2 <- update(m0, formula = . ~ Bloco + Prod * f(conc, 2))
# Modelo com polinômio de 1 grau.
m1 <- update(m0, formula = . ~ Bloco + Prod * f(conc, 1))
# Teste da falta de ajuste.
anova(m2, m0)
## Analysis of Variance Table
##
## Model 1: y ~ Bloco + Prod + f(conc, 2) + Prod:f(conc, 2)
## Model 2: y ~ Bloco + Prod * Conc
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 40 0.56702
## 2 32 0.42956 8 0.13745 1.2799 0.2884
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: y ~ Bloco + Prod + f(conc, 1) + Prod:f(conc, 1)
## Model 2: y ~ Bloco + Prod * Conc
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 42 0.62237
## 2 32 0.42956 10 0.19281 1.4363 0.2095
# Quadro de anova do modelo mais parcimonioso.
anova(m1)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## Bloco 2 0.06035 0.03018 2.0364 0.143190
## Prod 1 0.14051 0.14051 9.4825 0.003649 **
## f(conc, 1) 1 1.75464 1.75464 118.4107 8.483e-14 ***
## Prod:f(conc, 1) 1 0.05501 0.05501 3.7121 0.060804 .
## Residuals 42 0.62237 0.01482
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo final.
m_fit <- m1
#-----------------------------------------------------------------------
# Para extrair as equações ajustadas.
# Modelo para obter as equações.
m_expr <- update(m0,
formula = . ~ 0 + Prod/f(conc, 1) + Bloco,
contrasts = list(Bloco = "contr.sum"))
# Confere se o modelo é o mesmo.
all.equal(deviance(m_fit), deviance(m_expr))
## [1] TRUE
# Manipula as estimativas.
coef_stats <- cbind(summary(m_expr)[["coefficients"]][, -3],
rename(as.data.frame(confint(m_expr)),
Lower = 1,
Upper = 2))
coef_stats <- coef_stats %>%
rownames_to_column(var = "Param") %>%
filter(str_detect(Param, "Bloco", negate = TRUE)) %>%
mutate(signif = signif_cut(`Pr(>|t|)`),
`Pr(>|t|)` = NULL)
coef_stats
## Param Estimate Std. Error Lower Upper
## 1 ProdFlex 2.2773839 0.03720070 2.2023098 2.3524579
## 2 ProdPlus 2.2446075 0.03720070 2.1695335 2.3196816
## 3 ProdFlex:f(conc, 1) -0.2337409 0.03691343 -0.3082352 -0.1592465
## 4 ProdPlus:f(conc, 1) -0.3343201 0.03691343 -0.4088144 -0.2598258
## signif
## 1 ***
## 2 ***
## 3 ***
## 4 ***
# Separa os coeficientes para reportar a equação.
sapply(levels(tb$Prod),
simplify = FALSE,
function(p) {
coef_stats %>%
filter(str_detect(Param, p)) %>%
mutate(Param = str_remove(Param, "Prod"))
})
## $Flex
## Param Estimate Std. Error Lower Upper signif
## 1 Flex 2.2773839 0.03720070 2.2023098 2.3524579 ***
## 2 Flex:f(conc, 1) -0.2337409 0.03691343 -0.3082352 -0.1592465 ***
##
## $Plus
## Param Estimate Std. Error Lower Upper signif
## 1 Plus 2.2446075 0.03720070 2.1695335 2.3196816 ***
## 2 Plus:f(conc, 1) -0.3343201 0.03691343 -0.4088144 -0.2598258 ***
#-----------------------------------------------------------------------
# Curvas ajustadas com as bandas.
# Valores preditos com intervalo de confiança para a média.
pred <- emmeans(m_fit,
specs = ~Prod + conc,
at = list(conc = conc_seq)) %>%
as.data.frame()
# Gráfico.
gg0 +
geom_ribbon(data = pred,
mapping = aes(x = conc,
ymin = lower.CL,
ymax = upper.CL,
color = Prod),
size = 0.25, lty = 2,
fill = "gray30", alpha = 0.1, inherit.aes = FALSE) +
geom_line(data = pred,
mapping = aes(x = conc,
y = emmean,
color = Prod),
inherit.aes = FALSE)
#
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()
# Gráfico.
gg0 <- ggplot(data = tb,
mapping = aes(x = conc, y = y, color = Prod)) +
geom_point() +
labs(x = "Concentration",
y = "Response",
color = "Product") +
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 + Prod * Conc, 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 2 0.0005 0.00025 0.0103 0.9898
## Prod 1 0.0099 0.00988 0.4138 0.5246
## Conc 6 9.0770 1.51284 63.3534 2.295e-16 ***
## Prod:Conc 6 0.1109 0.01849 0.7742 0.5960
## Residuals 32 0.7641 0.02388
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Ajuste do modelo reduzido com `conc` quantitativo.
# Modelo com polinômio de 2 grau.
m2 <- update(m0, formula = . ~ Bloco + f(conc, 2))
# Modelo com polinômio de 1 grau.
m1 <- update(m0, formula = . ~ Bloco + f(conc, 1))
# Teste da falta de ajuste.
anova(m2, m0)
## Analysis of Variance Table
##
## Model 1: y ~ Bloco + f(conc, 2)
## Model 2: y ~ Bloco + Prod * Conc
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 43 0.94873
## 2 32 0.76414 11 0.18459 0.7028 0.7266
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: y ~ Bloco + f(conc, 1)
## Model 2: y ~ Bloco + Prod * Conc
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 44 1.23661
## 2 32 0.76414 12 0.47247 1.6488 0.1271
# Quadro de anova do modelo mais parcimonioso.
anova(m1)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## Bloco 2 0.0005 0.0002 0.0087 0.9913
## f(conc, 1) 1 8.7253 8.7253 310.4574 <2e-16 ***
## Residuals 44 1.2366 0.0281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo final.
m_fit <- m1
#-----------------------------------------------------------------------
# Para extrair as equações ajustadas.
# Modelo para obter as equações.
m_expr <- update(m0,
formula = . ~ f(conc, 1) + Bloco,
contrasts = list(Bloco = "contr.sum"))
# Confere se o modelo é o mesmo.
all.equal(deviance(m_fit), deviance(m_expr))
## [1] TRUE
# Manipula as estimativas.
coef_stats <- cbind(summary(m_expr)[["coefficients"]][, -3],
rename(as.data.frame(confint(m_expr)),
Lower = 1,
Upper = 2))
coef_stats <- coef_stats %>%
rownames_to_column(var = "Param") %>%
filter(str_detect(Param, "Bloco", negate = TRUE)) %>%
mutate(signif = signif_cut(`Pr(>|t|)`),
`Pr(>|t|)` = NULL)
coef_stats
## Param Estimate Std. Error Lower Upper signif
## 1 (Intercept) 1.029435 0.03622658 0.9564252 1.1024449 ***
## 2 f(conc, 1) -0.633376 0.03594683 -0.7058221 -0.5609299 ***
#-----------------------------------------------------------------------
# Curvas ajustadas com as bandas.
# Valores preditos com intervalo de confiança para a média.
pred <- emmeans(m_fit,
specs = ~conc,
at = list(conc = conc_seq)) %>%
as.data.frame()
# Gráfico.
gg0 +
geom_ribbon(data = pred,
mapping = aes(x = conc,
ymin = lower.CL,
ymax = upper.CL),
size = 0.25, lty = 2,
fill = "gray30", alpha = 0.1, inherit.aes = FALSE) +
geom_line(data = pred,
mapping = aes(x = conc,
y = emmean),
inherit.aes = FALSE)
#