Ajuste de modelos não lineares para eclosão e mortalidade de nematoides
Marlon Henrique Hahn, Henrique Duarte & Walmes Marques Zeviani
2020-12-24
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-21 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
# library(multcomp)
# library(emmeans)
library(nlme)
library(drc)
library(tidyverse)
# source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/pairwise.R")2 Importação e preparo
#-----------------------------------------------------------------------
# Importação.
# Leitura.
tb <- read_csv2("dados_in_vitro.csv", col_types = cols())## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.attr(tb, "spec") <- NULL
str(tb)## tibble [84 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ prod : chr [1:84] "F" "F" "F" "F" ...
## $ exp : chr [1:84] "A" "A" "A" "A" ...
## $ bloco: chr [1:84] "I" "I" "I" "I" ...
## $ xe : num [1:84] 0 0.1 0.2 0.3 0.4 0.5 0.75 0 0.1 0.2 ...
## $ xm : num [1:84] 0 0.5 2.5 5 10 15 20 0 0.5 2.5 ...
## $ ye : num [1:84] 55.9 54.5 40.3 28.1 20.6 ...
## $ ym : num [1:84] 6.06 14.63 32.43 21.74 34.48 ...# Converte variáveis para fator.
tb <- tb %>%
mutate_at(c("prod", "exp", "bloco"), "factor") %>%
mutate(plot = interaction(exp, prod, bloco, drop = TRUE, sep = ":"),
expbloco = interaction(exp, bloco, drop = TRUE, sep = ":"))
# Inspeção das frequências.
tb %>%
xtabs(formula = ~exp + prod + bloco) %>%
ftable()## bloco I II III
## exp prod
## A F 7 7 7
## P 7 7 7
## B F 7 7 7
## P 7 7 73 Análise exploratória
#-----------------------------------------------------------------------
# Gráficos.
# Eclosão.
ggplot(data = tb,
mapping = aes(x = xe, y = ye, color = bloco)) +
facet_grid(facets = exp ~ prod) +
geom_line() +
geom_point() +
labs(x = "Concentration (%)",
y = "Hatched (%)")
# NOTE: O comportamento da eclosão é bastante próximo do linear.
# Mortalidade.
ggplot(data = tb,
mapping = aes(x = xm, y = ym, color = bloco)) +
facet_grid(facets = exp ~ prod) +
geom_line() +
geom_point() +
labs(x = "Concentration (%)",
y = "Mortality (%)")
# NOTE: O comportamento da mortalidade pode ser sigmoidal. Talvez seja
# preciso transformar.4 Análise para a eclosão
4.1 Modelo com fatores qualitativos
#-----------------------------------------------------------------------
# Ajuste do modelo com fatores qualitativos apenas para apreciar a
# interação e efeito de bloco.
# Estrutura de dados agrupados.
tb <- groupedData(formula = ye ~ xe | plot,
data = tb,
order = FALSE)
# Modelo de parcela subdivida ou medidas repetidas.
m0 <- lme(ye ~ expbloco + prod * factor(xe),
random = ~1 | plot,
data = tb)
# Quadro de testes de Wald para os termos de efeito fixo.
anova(m0)## numDF denDF F-value p-value
## (Intercept) 1 60 1058.7137 <.0001
## expbloco 5 5 0.4742 0.7839
## prod 1 5 0.0001 0.9908
## factor(xe) 6 60 25.8974 <.0001
## prod:factor(xe) 6 60 0.0794 0.9980# ATTENTION: não houve efeito de bloco.
# # Teste da aditividade do bloco (rever).
# dae::tukey.1df(aov.obj = aov(ye ~ expbloco + interaction(prod, xe) +
# Error(plot),
# data = tb),
# data = tb)
# IMPORTANT: como não houve efeito de exp:bloco, para simplificar a
# análise, estes termos podem ser abandonados dos modelos. Isso facilita
# o ajuste das curvas. No entanto, alguns autores dizem que os termos de
# blocagem representam restrição à casualização e devem permanecer
# sempre no modelo.4.2 Ajuste de curvas
#-----------------------------------------------------------------------
# Ajuste de vários modelos.
# Para exibir os objetos públicos do pacote.
# ls("package:drc")
sort(unique(tb$xe))## [1] 0.00 0.10 0.20 0.30 0.40 0.50 0.75ye_LL4 <- drm(ye ~ xe, curveid = prod, data = tb, fct = LL.4())
mselect(ye_LL4,
fctList = list(W2.3(), G.4(), LL.3(), LL.5(), LL.4(), W2.4()),
sorted = "IC")## logLik IC Lack of fit Res var
## W2.3 -314.9984 643.9968 0.9906341 113.9994
## LL.3 -315.5494 645.0987 0.9583022 115.5047
## W2.4 -314.9644 647.9289 0.9571503 116.9047
## LL.4 -315.4002 648.8004 0.8920737 118.1240
## LL.4 -315.4002 648.8004 0.8920737 118.1240
## G.4 -316.1868 650.3736 0.7266816 120.3572
## LL.5 -314.8691 651.7382 0.8523797 119.7922# Modelo linear.
logLik(lm(ye ~ prod * xe, data = tb))## 'log Lik.' -318.896 (df=5)#-----------------------------------------------------------------------
# Ajuste do melhor modelo.
# Modelo vencedor.
ye_fit <- drm(ye ~ xe, curveid = prod, data = tb, fct = W2.3())
ye_fit##
## A 'drc' model.
##
## Call:
## drm(formula = ye ~ xe, curveid = prod, data = tb, fct = W2.3())
##
## Coefficients:
## b:F b:P d:F d:P e:F e:P
## -1.1238 -1.0524 61.1072 64.1637 0.3193 0.2814summary(ye_fit)##
## Model fitted: Weibull (type 2) with lower limit at 0 (3 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:F -1.123818 0.232567 -4.8322 6.625e-06 ***
## b:P -1.052408 0.206568 -5.0947 2.373e-06 ***
## d:F 61.107214 3.584040 17.0498 < 2.2e-16 ***
## d:P 64.163712 3.913019 16.3975 < 2.2e-16 ***
## e:F 0.319268 0.044468 7.1797 3.578e-10 ***
## e:P 0.281409 0.042818 6.5723 5.070e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 10.67705 (78 degrees of freedom)# Estimativas dos parâmetros
ye_param <- confint(ye_fit) %>%
as.data.frame() %>%
rename("Lower" = 1, "Upper" = 2) %>%
rownames_to_column(var = "comb") %>%
mutate("Estimate" = coef(ye_fit)) %>%
separate(col = "comb", into = c("param", "prod"), sep = ":")
ye_param## param prod Lower Upper Estimate
## 1 b F -1.5868233 -0.6608127 -1.1238180
## 2 b P -1.4636539 -0.6411624 -1.0524081
## 3 d F 53.9719374 68.2424901 61.1072137
## 4 d P 56.3734907 71.9539337 64.1637122
## 5 e F 0.2307384 0.4077972 0.3192678
## 6 e P 0.1961656 0.3666522 0.2814089ggplot(data = ye_param,
mapping = aes(y = prod, x = Estimate, xmin = Lower, xmax = Upper)) +
facet_wrap(facets = ~param, scale = "free_x", ncol = 1) +
geom_errorbarh(height = 0.2) +
geom_point() +
labs(x = "Estimate", y = "Product")
#-----------------------------------------------------------------------
# Visualização.
# Sequência de valores para a predição.
xe_seq <- with(tb, seq(0, 1.05 * max(xe), length.out = 53))
pred <- with(tb,
crossing(xe = xe_seq,
prod = levels(prod)))
pred## # A tibble: 106 x 2
## xe prod
## <dbl> <chr>
## 1 0 F
## 2 0 P
## 3 0.0151 F
## 4 0.0151 P
## 5 0.0303 F
## 6 0.0303 P
## 7 0.0454 F
## 8 0.0454 P
## 9 0.0606 F
## 10 0.0606 P
## # … with 96 more rows# Predição com bandas de confiança.
pred <- fitted(ye_fit,
newdata = as.data.frame(pred),
interval = "confidence") %>%
as.data.frame() %>%
cbind(pred)
names(pred) <- c("fit", "lwr", "upr", "xe", "prod")
# Gráfico.
gg_fit <- ggplot(data = pred,
mapping = aes(x = xe,
y = fit,
color = prod,
ymin = lwr,
ymax = upr)) +
geom_ribbon(fill = NA, lty = 1, size = 0.25) +
geom_line(size = 1) +
geom_point(data = tb,
mapping = aes(x = xe, y = ye, color = prod),
inherit.aes = FALSE) +
labs(x = "Concentration (%)",
y = "Hatched (%)",
color = "Product") +
theme(legend.position = c(0.975, 0.975),
legend.justification = c(1, 1))
gg_fit
# ATTENTION! Mas afinal qual é a expressão desse modelo?
#-----------------------------------------------------------------------
# Refazendo com o `gnls()` só para conhecimento.
g0 <- gnls(ye ~ d * (1 - exp(-exp(b * (log(xe) - log(e))))),
data = tb,
params = b + d + e ~ 0 + prod,
start = coef(ye_fit))
# Confere as estimativas.
cbind(`drc()` = coef(ye_fit), `gnls()` = coef(g0))## drc() gnls()
## b:F -1.1238180 -1.1237536
## b:P -1.0524081 -1.0525617
## d:F 61.1072137 61.1075544
## d:P 64.1637122 64.1631518
## e:F 0.3192678 0.3192585
## e:P 0.2814089 0.2814140c(`drc()` = logLik(g0), `gnls()` = logLik(ye_fit))## drc() gnls()
## -314.9984 -314.9984#-----------------------------------------------------------------------
# Valores de EC50.
# Determina a EC50 relativa ao máximo.
ye_EC50 <- ED(ye_fit,
respLev = 50,
interval = "delta",
reference = "upper",
type = "relative") %>%
as.data.frame() %>%
mutate(prod = levels(tb$prod),
param = "EC50")##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:F:50 0.442377 0.054456 0.333963 0.550791
## e:P:50 0.398645 0.051249 0.296616 0.500673ye_EC50## Estimate Std. Error Lower Upper prod param
## 1 0.4423769 0.05445628 0.3339627 0.5507910 F EC50
## 2 0.3986445 0.05124871 0.2966162 0.5006729 P EC50ye_EC50 <- cbind(ye_EC50,
value = c(0, 2))
# Adiciona ao gráfico.
gg_fit +
expand_limits(y = c(0, NA)) +
geom_segment(data = ye_EC50,
mapping = aes(x = Lower, xend = Upper,
y = value, yend = value,
color = prod),
arrow = arrow(length = unit(2, "mm"),
ends = "both",
angle = 90),
inherit.aes = FALSE) +
geom_point(data = ye_EC50,
mapping = aes(x = Estimate, y = value, color = prod),
inherit.aes = FALSE)
#-----------------------------------------------------------------------
# Coeficiente de determinação.
ye_R2 <- tb %>%
mutate(fit = fitted(ye_fit)) %>%
group_by(prod) %>%
summarise(R2 = cor(ye, fit)^2)## `summarise()` ungrouping output (override with `.groups` argument)# IC para os parâmetros.
ggplot(data = bind_rows(ye_param, ye_EC50),
mapping = aes(y = prod, x = Estimate, xmin = Lower, xmax = Upper)) +
facet_wrap(facets = ~param, scale = "free_x", ncol = 1) +
geom_errorbarh(height = 0.2) +
geom_point() +
labs(x = "Estimate", y = "Product")
#5 Análise para a mortalidade
5.1 Modelo com fatores qualitativos
#-----------------------------------------------------------------------
# Ajuste do modelo com fatores qualitativos apenas para apreciar a
# interação e efeito de bloco.
# Estrutura de dados agrupados.
tb <- groupedData(formula = ym ~ xm | plot,
data = tb,
order = FALSE)
# Modelo de parcela subdivida ou medidas repetidas.
m0 <- lme(ym ~ expbloco + prod * factor(xm),
random = ~1 | plot,
data = tb)
# Quadro de testes de Wald para os termos de efeito fixo.
anova(m0)## numDF denDF F-value p-value
## (Intercept) 1 60 2535.3156 <.0001
## expbloco 5 5 1.1774 0.4311
## prod 1 5 0.0544 0.8249
## factor(xm) 6 60 233.9859 <.0001
## prod:factor(xm) 6 60 2.0668 0.0705# ATTENTION: não houve efeito de bloco.
# # Teste da aditividade do bloco (rever).
# dae::tukey.1df(aov.obj = aov(ym ~ expbloco + interaction(prod, xm) +
# Error(plot),
# data = tb),
# data = tb)
# IMPORTANT: como não houve efeito de exp:bloco, para simplificar a
# análise, estes termos podem ser abandonados dos modelos. Isso facilita
# o ajuste das curvas. No entanto, alguns autores dizem que os termos de
# blocagem representam restrição à casualização e devem permanecer
# sempre no modelo.5.2 Ajuste de curvas
#-----------------------------------------------------------------------
# Ajuste de vários modelos.
# Para exibir os objetos públicos do pacote.
# ls("package:drc")
sort(unique(tb$xm))## [1] 0.0 0.5 2.5 5.0 10.0 15.0 20.0ym_LL4 <- drm(ym ~ xm, curveid = prod, data = tb, fct = LL.4())
mselect(ym_LL4,
fctList = list(W2.3(), G.4(), LL.3(), LL.5(), LL.4(), W2.4()),
sorted = "IC")## logLik IC Lack of fit Res var
## W2.4 -293.4316 604.8633 1.758766e-02 70.01251
## LL.4 -293.9721 605.9441 1.226914e-02 70.91921
## LL.4 -293.9721 605.9441 1.226914e-02 70.91921
## G.4 -294.3905 606.7810 9.252162e-03 71.62930
## LL.5 -294.4419 610.8838 2.052102e-03 73.65529
## W2.3 -341.0805 696.1611 1.136894e-17 212.12936
## LL.3 -341.2740 696.5479 9.721709e-18 213.10846# Modelo linear.
logLik(lm(ym ~ prod * xm, data = tb))## 'log Lik.' -328.6005 (df=5)#-----------------------------------------------------------------------
# Ajuste do melhor modelo.
# Modelo vencedor.
ym_fit <- drm(ym ~ xm, curveid = prod, data = tb, fct = W2.4())
ym_fit##
## A 'drc' model.
##
## Call:
## drm(formula = ym ~ xm, curveid = prod, data = tb, fct = W2.4())
##
## Coefficients:
## b:F b:P c:F c:P d:F d:P e:F e:P
## 7.976 3.696 21.309 18.684 97.577 100.213 13.503 13.578summary(ym_fit)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:F 7.97617 1.58063 5.0462 2.993e-06 ***
## b:P 3.69648 0.73846 5.0057 3.506e-06 ***
## c:F 21.30900 1.71524 12.4233 < 2.2e-16 ***
## c:P 18.68449 1.88857 9.8934 2.654e-15 ***
## d:F 97.57739 3.41519 28.5716 < 2.2e-16 ***
## d:P 100.21315 4.72982 21.1875 < 2.2e-16 ***
## e:F 13.50307 0.43492 31.0472 < 2.2e-16 ***
## e:P 13.57771 0.60905 22.2931 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 8.367348 (76 degrees of freedom)# Estimativas dos parâmetros
ym_param <- confint(ym_fit) %>%
as.data.frame() %>%
rename("Lower" = 1, "Upper" = 2) %>%
rownames_to_column(var = "comb") %>%
mutate("Estimate" = coef(ym_fit)) %>%
separate(col = "comb", into = c("param", "prod"), sep = ":")
ym_param## param prod Lower Upper Estimate
## 1 b F 4.828075 11.124270 7.976172
## 2 b P 2.225711 5.167244 3.696477
## 3 c F 17.892805 24.725192 21.308999
## 4 c P 14.923068 22.445904 18.684486
## 5 d F 90.775456 104.379333 97.577395
## 6 d P 90.792890 109.633412 100.213151
## 7 e F 12.636850 14.369288 13.503069
## 8 e P 12.364673 14.790748 13.577710ggplot(data = ym_param,
mapping = aes(y = prod, x = Estimate, xmin = Lower, xmax = Upper)) +
facet_wrap(facets = ~param, scale = "free_x", ncol = 1) +
geom_errorbarh(height = 0.2) +
geom_point() +
labs(x = "Estimate", y = "Product")
#-----------------------------------------------------------------------
# Visualização.
# Sequência de valores para a predição.
xm_seq <- with(tb, seq(0, 1.05 * max(xm), length.out = 53))
pred <- with(tb,
crossing(xm = xm_seq,
prod = levels(prod)))
pred## # A tibble: 106 x 2
## xm prod
## <dbl> <chr>
## 1 0 F
## 2 0 P
## 3 0.404 F
## 4 0.404 P
## 5 0.808 F
## 6 0.808 P
## 7 1.21 F
## 8 1.21 P
## 9 1.62 F
## 10 1.62 P
## # … with 96 more rows# Predição com bandas de confiança.
pred <- fitted(ym_fit,
newdata = as.data.frame(pred),
interval = "confidence") %>%
as.data.frame() %>%
cbind(pred)
names(pred) <- c("fit", "lwr", "upr", "xm", "prod")
# Gráfico.
gg_fit <- ggplot(data = pred,
mapping = aes(x = xm,
y = fit,
color = prod,
ymin = lwr,
ymax = upr)) +
geom_ribbon(fill = NA, lty = 1, size = 0.25) +
geom_line(size = 1) +
geom_point(data = tb,
mapping = aes(x = xm, y = ym, color = prod),
inherit.aes = FALSE) +
labs(x = "Concentration (%)",
y = "Mortality (%)",
color = "Product") +
theme(legend.position = c(0.025, 0.975),
legend.justification = c(0, 1))
gg_fit
# ATTENTION! Mas afinal qual é a expressão desse modelo?
#-----------------------------------------------------------------------
# Refazendo com o `gnls()` só para conhecimento.
g0 <- gnls(ym ~ c + (d - c) * (1 - exp(-exp(b * (log(xm) - log(e))))),
data = tb,
params = b + c + d + e ~ 0 + prod,
start = coef(ym_fit))
# Confere as estimativas.
cbind(`drc()` = coef(ym_fit), `gnls()` = coef(g0))## drc() gnls()
## b:F 7.976172 7.976150
## b:P 3.696477 3.696465
## c:F 21.308999 21.309057
## c:P 18.684486 18.684471
## d:F 97.577395 97.577799
## d:P 100.213151 100.213339
## e:F 13.503069 13.503108
## e:P 13.577710 13.577728c(`drc()` = logLik(g0), `gnls()` = logLik(ym_fit))## drc() gnls()
## -293.4316 -293.4316#-----------------------------------------------------------------------
# Valores de EC50.
# Determina a EC50 relativa ao máximo.
ym_EC50 <- ED(ym_fit,
respLev = 50,
interval = "delta",
reference = "upper",
type = "relative") %>%
as.data.frame() %>%
mutate(prod = levels(tb$prod),
param = "EC50")##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:F:50 12.89663 0.46257 11.97534 13.81792
## e:P:50 12.29604 0.49760 11.30499 13.28709ym_EC50## Estimate Std. Error Lower Upper prod param
## 1 12.89663 0.4625726 11.97534 13.81792 F EC50
## 2 12.29604 0.4975964 11.30499 13.28709 P EC50ym_EC50 <- cbind(ym_EC50,
value = c(0, 2))
# Adiciona ao gráfico.
gg_fit +
expand_limits(y = c(0, NA)) +
geom_segment(data = ym_EC50,
mapping = aes(x = Lower, xend = Upper,
y = value, yend = value,
color = prod),
arrow = arrow(length = unit(2, "mm"),
ends = "both",
angle = 90),
inherit.aes = FALSE) +
geom_point(data = ym_EC50,
mapping = aes(x = Estimate, y = value, color = prod),
inherit.aes = FALSE)
#-----------------------------------------------------------------------
# Coeficiente de determinação.
ym_R2 <- tb %>%
mutate(fit = fitted(ym_fit)) %>%
group_by(prod) %>%
summarise(R2 = cor(ym, fit)^2)## `summarise()` ungrouping output (override with `.groups` argument)# IC para os parâmetros.
ggplot(data = bind_rows(ym_param, ym_EC50),
mapping = aes(y = prod, x = Estimate, xmin = Lower, xmax = Upper)) +
facet_wrap(facets = ~param, scale = "free_x", ncol = 1) +
geom_errorbarh(height = 0.2) +
geom_point() +
labs(x = "Estimate", y = "Product")
#