library(readxl)
library(ggplot2)
library(tidyverse)
library(scales)
library(forcats)
library(stringr)
library(plotly)
library(readr)
library(emmeans)
library(multcomp)
library(multcompView)
library(DHARMa)
library(car)
library(glmmTMB)
library(lme4)
library(performance)
library(MASS)
library(dunn.test)
Análise dos dados do experimento de eficiência de inseticidas no controle de lagartas da traça-do-tomateiro
Objetivos e hipóteses
Por meio da análise dos dados, objetivou-se determinar a eficiência de inseticidas no controle de lagartas de Phthorimaea absoluta.
Hipóteses testadas
Hipótese nula (H0): Não há diferença significativa entre os tratamentos. A média da eficiência dos tratamentos no controle é igual para todos os inseticidas.
Hipótese alternativa (Ha): Pelo menos um dos inseticidas apresenta um desempenho significativamente melhor do que o controle, em comparação com os demais.
Análise das variáveis resposta
Pré-análise dos dados
Pacotes: os seguintes pacotes R foram utilizados para as análises.
Carregamento dos dados: como o dataframe contendo os dados de eficiência dos inseticidas estava localizado em uma planilha Excel no desktop, para o carregamento dos dados, foi necessário utilizar a função read_excel
do pacote readxl
. O dataframe foi atribuído ao objeto denominado dados
.
<- read_excel("teste_lagarta_traça.xlsx")
dados
dados
# A tibble: 96 × 6
tratamento repetição tempo vivos mortos total
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 abamex 1 48 0 10 10
2 abamex 2 48 0 10 10
3 abamex 3 48 0 10 10
4 abamex 4 48 0 10 10
5 abamex 5 48 0 10 10
6 abamex 6 48 0 10 10
7 actara 1 48 6 4 10
8 actara 2 48 4 6 10
9 actara 3 48 6 4 10
10 actara 4 48 5 5 10
# ℹ 86 more rows
Análise visual dos dados
Usando o pacote ggplot2
, foi explorado visualmente, por meio de diferentes gráficos, os dados de eficiência dos inseticidas.
Gráfico de barras: primeiramente, foi visualizada a proporção de mortalidade de insetos em cada tratamento. A proporção de insetos mortos foi calculada pela fórmula:
\(PM = TM / (TV + TM)\)
PM = proporção de mortos
TM = total de insetos mortos no tratamento
TV = total de insetos vivos no tratamento
<- dados |>
proporcao_mortos group_by(tratamento) |>
summarise(
TM = sum(mortos, na.rm = TRUE),
TV = sum(vivos, na.rm = TRUE)
|>
) mutate(
tratamento = str_to_title(tratamento),
PM = TM / (TM + TV)
|>
) ungroup() |>
arrange(PM) |>
mutate(tratamento = factor(tratamento, levels = tratamento))
print(proporcao_mortos)
# A tibble: 16 × 4
tratamento TM TV PM
<fct> <dbl> <dbl> <dbl>
1 Controle 1 59 0.0167
2 Actara 30 30 0.5
3 Benevia 43 17 0.717
4 Premio 44 16 0.733
5 Nomolt 23 7 0.767
6 Hayate 51 9 0.85
7 Tracer 53 7 0.883
8 Avatar 56 4 0.933
9 Belt 57 3 0.95
10 Vertimec 58 2 0.967
11 Talstar 59 1 0.983
12 Abamex 60 0 1
13 Cartap 60 0 1
14 Delegate 60 0 1
15 Joiner 60 0 1
16 Pirate 60 0 1
<- ggplot(proporcao_mortos, aes(x = tratamento, y = PM)) +
grafico_proporcao geom_col(fill = "darkgreen") +
labs(
title = "Proporção de insetos mortos por tratamento",
x = "Tratamentos",
y = "Proporção de mortos"
+
) ylim(0, 1) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
<- plotly::ggplotly(grafico_proporcao) |>
grafico_interativo ::config(displayModeBar = FALSE)
plotly
grafico_interativo
Adição de facetas: a visualização de cada tratamento individualmente foi possibilitada pela adição de facetas, permitindo analisar o comportamento de cada repetição (número de insetos mortos). Isso foi feito com a função facet_wrap
.
|>
dados ggplot(aes(repetição, mortos))+
geom_col(fill = "darkgreen")+
labs(x = "Repetições",
y = "Número de insetos mortos")+
theme_bw()+
theme(panel.grid.minor = element_blank())+
scale_y_continuous(limits = c(0, 10), n.breaks = 3)+
facet_wrap(~ tratamento, labeller = labeller(tratamento = tools::toTitleCase))
Mortalidade
A mortalidade é descrita como:
\(MT = M / T\)
MT = mortalidade
M = mortos
T = total
<- dados |>
dados mutate(
mortalidade = mortos / total)
<- dados |>
media_mortalidade group_by(tratamento) |>
summarise(
mortalidade_media = mean(mortalidade, na.rm = TRUE),
.groups = "drop")
<- media_mortalidade |>
mortalidade_controle filter(grepl("controle", tolower(tratamento))) |>
pull(mortalidade_media)
print(mortalidade_controle)
[1] 0.01666667
print(media_mortalidade)
# A tibble: 16 × 2
tratamento mortalidade_media
<chr> <dbl>
1 abamex 1
2 actara 0.5
3 avatar 0.933
4 belt 0.95
5 benevia 0.717
6 cartap 1
7 controle 0.0167
8 delegate 1
9 hayate 0.85
10 joiner 1
11 nomolt 0.767
12 pirate 1
13 premio 0.733
14 talstar 0.983
15 tracer 0.883
16 vertimec 0.967
Análises estatísticas
Análise das pressuposições da ANOVA
Para que a ANOVA seja realizada é necessário que ela se enquandre nos seguintes pressupostos:
Normalidade dos resíduos
Homogeneidade de variâncias (homocedasticidade)
Aditividade do modelo
<- aov(mortalidade ~ tratamento, data = dados)
modelo_anova
summary(modelo_anova)
Df Sum Sq Mean Sq F value Pr(>F)
tratamento 15 6.070 0.4046 65.18 <2e-16 ***
Residuals 80 0.497 0.0062
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Histograma dos resíduos: foi criado o histograma dos resíduos da ANOVA, criado com a função hist
a partir do objeto modelo_anova
, serve para visualizar a distribuição desses resíduos. Isso permite verificar se eles são aproximadamente normais, uma suposição fundamental para a validade dos resultados da ANOVA.
hist(residuals(modelo_anova))
plot(simulateResiduals(modelo_anova))
Normalidade dos resíduos: foi verificada pelo Shapiro teste, executado a partir da função shapiro.test
a partir do objeto modelo_anova
.
shapiro.test(residuals(modelo_anova))
Shapiro-Wilk normality test
data: residuals(modelo_anova)
W = 0.90715, p-value = 4.649e-06
Como p-valor < 0.05, os resíduos não possuem normalidade.
Homogeneidade das variâncias: foi testada pelo teste de Levene, por meio da função leveneTest
aplicada a mortalidade por tratamento.
::leveneTest(mortalidade ~ tratamento, data = dados) car
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 15 3.598 9.517e-05 ***
80
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Como o p-valor < 0.05, as variâncias são heterogêneas.
Assim, observamos que os dados não atendem as pressuposições da ANOVA. Assim, se torna necessária a transformação dos dados ou a utilização de outros testes e modelos para a análise de dados. Com isso, apresentamos as duas formas de análise dos dados.
Transformação dos dados
Os dados foram transformados pela função anscombe = asin(sqrt(x))
(transformação angular (ou de Anscombe)). É uma transformação usada para dados de proporção
<- dados %>%
dados mutate(
anscombe = asin(sqrt((mortos + 3/8) / (total + 3/4)))
)
print(dados)
# A tibble: 96 × 8
tratamento repetição tempo vivos mortos total mortalidade anscombe
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 abamex 1 48 0 10 10 1 1.38
2 abamex 2 48 0 10 10 1 1.38
3 abamex 3 48 0 10 10 1 1.38
4 abamex 4 48 0 10 10 1 1.38
5 abamex 5 48 0 10 10 1 1.38
6 abamex 6 48 0 10 10 1 1.38
7 actara 1 48 6 4 10 0.4 0.692
8 actara 2 48 4 6 10 0.6 0.879
9 actara 3 48 6 4 10 0.4 0.692
10 actara 4 48 5 5 10 0.5 0.785
# ℹ 86 more rows
Aplicando a ANOVA novamente
<- aov(anscombe ~ tratamento, data = dados) modelo_anova2
Histograma dos resíduos
hist(residuals(modelo_anova2))
Normalidade dos resíduos
shapiro.test(residuals(modelo_anova2))
Shapiro-Wilk normality test
data: residuals(modelo_anova2)
W = 0.93441, p-value = 0.0001243
Como o p-valor < 0.05, os resíduos não são normais.
Homogeneidade das variâncias
::leveneTest(anscombe ~ tratamento, data = dados) car
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 15 3.0498 0.0006496 ***
80
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Como o p-valor < 0.05, as variâncias são heterogêneas.
Assim, observamos que a transformação dos dados não foi suficiente para que eles atendessem as pressuposições da ANOVA.
Teste de Kruskal-Wallis
O teste de Kruskal-Wallis é um teste não-paramétrico que verifica se há diferenças significativas entre os tratamentos analisados. Assim, ele foi aplicado aos dados de mortalidade.
kruskal.test(mortalidade ~ tratamento, data = dados)
Kruskal-Wallis rank sum test
data: mortalidade by tratamento
Kruskal-Wallis chi-squared = 77.032, df = 15, p-value = 2.427e-10
Como p-valor < 0.05, há diferenças significativas na mortalidade entre os tratamentos.
Comparações múltiplas
Para realizar as comparações dos tratamentos foi realizado o teste de Dunn.
<- dunn.test(dados$mortalidade, dados$tratamento, method = "holm") teste_dunn
Kruskal-Wallis rank sum test
data: x and group
Kruskal-Wallis chi-squared = 77.0316, df = 15, p-value = 0
Comparison of x by group
(Holm)
Col Mean-|
Row Mean | abamex actara avatar belt benevia cartap
---------+------------------------------------------------------------------
actara | 4.145820
| 0.0019*
|
avatar | 1.463230 -2.682589
| 1.0000 0.3141
|
belt | 1.097423 -3.048397 -0.365807
| 1.0000 0.1081 1.0000
|
benevia | 3.297811 -0.848008 1.834581 2.200388
| 0.0512 1.0000 1.0000 1.0000
|
cartap | 0.000000 -4.145820 -1.463230 -1.097423 -3.297811
| 1.0000 0.0019* 1.0000 1.0000 0.0507
|
controle | 4.622479 0.476658 3.159248 3.525056 1.324667 4.622479
| 0.0002* 1.0000 0.0751 0.0227* 1.0000 0.0002*
|
delegate | 0.000000 -4.145820 -1.463230 -1.097423 -3.297811 0.000000
| 1.0000 0.0019* 1.0000 1.0000 0.0502 1.0000
|
hayate | 2.023027 -2.122793 0.559796 0.925604 -1.274784 2.023027
| 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
|
joiner | 0.000000 -4.145820 -1.463230 -1.097423 -3.297811 0.000000
| 1.0000 0.0019* 1.0000 1.0000 0.0497 1.0000
|
nomolt | 2.810068 -1.335752 1.346837 1.712645 -0.487743 2.810068
| 0.2254 1.0000 1.0000 1.0000 1.0000 0.2229
|
pirate | 0.000000 -4.145820 -1.463230 -1.097423 -3.297811 0.000000
| 1.0000 0.0019* 1.0000 1.0000 0.0492 1.0000
|
premio | 3.181418 -0.964402 1.718187 2.083995 -0.116393 3.181418
| 0.0733 1.0000 1.0000 1.0000 1.0000 0.0725
|
talstar | 0.365807 -3.780012 -1.097423 -0.731615 -2.932004 0.365807
| 1.0000 0.0085* 1.0000 1.0000 0.1566 1.0000
|
tracer | 2.333409 -1.812410 0.870178 1.235986 -0.964402 2.333409
| 0.7949 1.0000 1.0000 1.0000 1.0000 0.7851
|
vertimec | 0.731615 -3.414205 -0.731615 -0.365807 -2.566196 0.731615
| 1.0000 0.0339 1.0000 1.0000 0.4318 1.0000
Col Mean-|
Row Mean | controle delegate hayate joiner nomolt pirate
---------+------------------------------------------------------------------
delegate | -4.622479
| 0.0002*
|
hayate | -2.599451 2.023027
| 0.3968 1.0000
|
joiner | -4.622479 0.000000 -2.023027
| 0.0002* 1.0000 1.0000
|
nomolt | -1.812410 2.810068 0.787040 2.810068
| 1.0000 0.2204 1.0000 0.2179
|
pirate | -4.622479 0.000000 -2.023027 0.000000 -2.810068
| 0.0002* 1.0000 1.0000 0.5000 0.2155
|
premio | -1.441060 3.181418 1.158391 3.181418 0.371350 3.181418
| 1.0000 0.0718 1.0000 0.0711 1.0000 0.0703
|
talstar | -4.256671 0.365807 -1.657219 0.365807 -2.444260 0.365807
| 0.0012* 1.0000 1.0000 1.0000 0.5951 1.0000
|
tracer | -2.289069 2.333409 0.310382 2.333409 -0.476658 2.333409
| 0.8389 0.7753 1.0000 0.7654 1.0000 0.7556
|
vertimec | -3.890863 0.731615 -1.291412 0.731615 -2.078452 0.731615
| 0.0054* 1.0000 1.0000 1.0000 1.0000 1.0000
Col Mean-|
Row Mean | premio talstar tracer
---------+---------------------------------
talstar | -2.815610
| 0.2239
|
tracer | -0.848008 1.967602
| 1.0000 1.0000
|
vertimec | -2.449803 0.365807 -1.601794
| 0.5932 1.0000 1.0000
alpha = 0.05
Reject Ho if p <= alpha/2
<- teste_dunn$comparisons
comparacoes <- teste_dunn$P.adjusted
pvals <- unique(as.character(dados$tratamento))
grupos <- matrix(1, nrow = length(grupos), ncol = length(grupos),
matriz_p dimnames = list(grupos, grupos))
for (i in seq_along(comparacoes)) {
<- unlist(strsplit(comparacoes[i], " - "))
par 1], par[2]] <- pvals[i]
matriz_p[par[2], par[1]] <- pvals[i]
matriz_p[par[
}
<- multcompLetters(matriz_p)$Letters
letras <- data.frame(tratamento = names(letras), letra = letras)
letras_df
<- dados %>%
summary_ic group_by(tratamento) %>%
summarise(
media = mean(mortalidade, na.rm = TRUE),
n = n(),
sd = sd(mortalidade, na.rm = TRUE),
se = sd / sqrt(n),
ic = qt(0.975, df = n - 1) * se
)
<- summary_ic %>%
plot_df left_join(letras_df, by = "tratamento") %>%
mutate(tratamento = str_to_title(tratamento))
ggplot(plot_df, aes(x = tratamento, y = media)) +
geom_errorbar(aes(ymin = media - ic, ymax = media + ic), width = 0.2) +
geom_point(size = 3) +
geom_text(aes(label = letra, y = media + ic + 0.03), size = 5) +
labs(
title = "Mortalidade por tratamento (Teste de Dunn)",
x = "Tratamento",
y = "Mortalidade"
+
) theme_minimal() +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)
)
Modelos alternativos
Uma forma de analisar os dados sem transformação é a utilização de modelos que não exigem características semelhantes a da ANOVA. Dentre eles pode ser utilizado o GLM e GLMM.
Ajuste dos dados
Como possuimos dados com valores de 0 ou 100% (ou seja, tudo morto ou tudo vivo), o modelo encontra problemas. Assim, os dados foram ajustado adicionando 0.5
ao número de mortos e vivos. A função cbind
foi utilizada, pois estamos tratando de contagem de eventos (morto x vivo).
$tratamento <- as.factor(dados$tratamento)
dados
<- dados |>
dados_ajustados ::filter(mortos + vivos > 0) |>
dplyr::mutate(
dplyrmortos_ajustados = mortos + 0.5,
vivos_ajustados = vivos + 0.5,
resposta_binomial = cbind(mortos_ajustados, vivos_ajustados))
GLM - Modelo linear generalizado
Inicialmente, foi testado o modelo com distribuição binomial. Esse modelo é indicado quando os dados são proporcionais (insetos mortos / vivos). Com isso, é possível verificar se há diferença significativa entre os inseticidas testados quanto à mortalidade dos insetos. O argumento family = binomial
foi utilizada por ter apenas dois resultados possíveis (sucesso ou fracasso). O link logit
foi utilizado por se tratar de dados de proporção.
<- glm(resposta_binomial ~ tratamento,
modelo_glm family = binomial(link = "logit"),
data = dados_ajustados)
summary(modelo_glm)
Call:
glm(formula = resposta_binomial ~ tratamento, family = binomial(link = "logit"),
data = dados_ajustados)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.045e+00 5.909e-01 5.152 2.58e-07 ***
tratamentoactara -3.045e+00 6.402e-01 -4.756 1.98e-06 ***
tratamentoavatar -9.129e-01 7.134e-01 -1.280 0.200698
tratamentobelt -7.419e-01 7.297e-01 -1.017 0.309291
tratamentobenevia -2.212e+00 6.488e-01 -3.409 0.000652 ***
tratamentocartap -8.319e-15 8.357e-01 0.000 1.000000
tratamentocontrole -5.785e+00 7.844e-01 -7.375 1.64e-13 ***
tratamentodelegate -5.482e-15 8.357e-01 0.000 1.000000
tratamentohayate -1.540e+00 6.716e-01 -2.294 0.021808 *
tratamentojoiner -6.571e-15 8.357e-01 0.000 1.000000
tratamentonomolt -2.089e+00 6.983e-01 -2.991 0.002776 **
tratamentopirate -6.974e-15 8.357e-01 0.000 1.000000
tratamentopremio -2.139e+00 6.505e-01 -3.288 0.001008 **
tratamentotalstar -3.037e-01 7.844e-01 -0.387 0.698651
tratamentotracer -1.322e+00 6.834e-01 -1.934 0.053106 .
tratamentovertimec -5.431e-01 7.521e-01 -0.722 0.470209
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 331.240 on 95 degrees of freedom
Residual deviance: 26.177 on 80 degrees of freedom
AIC: 274.57
Number of Fisher Scoring iterations: 4
<- anova(modelo_glm, data = dados_ajustados)
anova_glm print(anova_glm)
Analysis of Deviance Table
Model: binomial, link: logit
Response: resposta_binomial
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 95 331.24
tratamento 15 305.06 80 26.18 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(modelo_glm))
qqnorm(residuals(modelo_glm))
qqline(residuals(modelo_glm))
<- sum(residuals(modelo_glm, type = "deviance")^2)
res_dev_binomial <- df.residual(modelo_glm)
df_binomial <- res_dev_binomial / df_binomial
dispersion_param_binomial print(dispersion_param_binomial)
[1] 0.3272128
A análise do Modelo Linear Generalizado (GLM) com família binomial revelou um parâmetro de dispersão residual de 0.3272128. Este valor, que é menor que 1, indicando a presença de subdispersão nos dados. Embora o teste de dispersão do DHARMa não tenha sido estatisticamente significativo (p = 0.184), a magnitude do parâmetro estimado (0.327) sugere que a variabilidade observada é consideravelmente menor do que a esperada sob a distribuição binomial padrão. Assim, a family foi ajustada para quasibinomial
.
Distribuição quasibinomial: o modelo GLM foi ajustado usando o argumento family = quasibinomial
. Ele permite que o parâmetro de dispersão dos dados seja estimado, corrigindo os erros padrão e os testes de significância.
<- glm(resposta_binomial ~ tratamento,
modelo_glm_qb family = quasibinomial(link = "logit"),
data = dados_ajustados)
summary(modelo_glm_qb)
Call:
glm(formula = resposta_binomial ~ tratamento, family = quasibinomial(link = "logit"),
data = dados_ajustados)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.045e+00 3.334e-01 9.132 4.84e-14 ***
tratamentoactara -3.045e+00 3.612e-01 -8.430 1.16e-12 ***
tratamentoavatar -9.129e-01 4.025e-01 -2.268 0.026029 *
tratamentobelt -7.419e-01 4.117e-01 -1.802 0.075304 .
tratamentobenevia -2.212e+00 3.660e-01 -6.042 4.57e-08 ***
tratamentocartap -8.319e-15 4.715e-01 0.000 1.000000
tratamentocontrole -5.785e+00 4.426e-01 -13.072 < 2e-16 ***
tratamentodelegate -5.482e-15 4.715e-01 0.000 1.000000
tratamentohayate -1.540e+00 3.789e-01 -4.065 0.000111 ***
tratamentojoiner -6.571e-15 4.715e-01 0.000 1.000000
tratamentonomolt -2.089e+00 3.940e-01 -5.302 9.91e-07 ***
tratamentopirate -6.974e-15 4.715e-01 0.000 1.000000
tratamentopremio -2.139e+00 3.670e-01 -5.828 1.13e-07 ***
tratamentotalstar -3.037e-01 4.426e-01 -0.686 0.494580
tratamentotracer -1.322e+00 3.856e-01 -3.428 0.000963 ***
tratamentovertimec -5.431e-01 4.243e-01 -1.280 0.204264
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasibinomial family taken to be 0.3183159)
Null deviance: 331.240 on 95 degrees of freedom
Residual deviance: 26.177 on 80 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4
<- emmeans(modelo_glm_qb, ~ tratamento, type = "response")
m print(m)
tratamento prob SE df asymp.LCL asymp.UCL
abamex 0.9545 0.0145 Inf 0.9161 0.976
actara 0.5000 0.0347 Inf 0.4324 0.568
avatar 0.8939 0.0214 Inf 0.8442 0.929
belt 0.9091 0.0200 Inf 0.8617 0.941
benevia 0.6970 0.0319 Inf 0.6311 0.756
cartap 0.9545 0.0145 Inf 0.9161 0.976
controle 0.0606 0.0166 Inf 0.0352 0.102
delegate 0.9545 0.0145 Inf 0.9161 0.976
hayate 0.8182 0.0268 Inf 0.7597 0.865
joiner 0.9545 0.0145 Inf 0.9161 0.976
nomolt 0.7222 0.0421 Inf 0.6328 0.797
pirate 0.9545 0.0145 Inf 0.9161 0.976
premio 0.7121 0.0314 Inf 0.6468 0.770
talstar 0.9394 0.0166 Inf 0.8976 0.965
tracer 0.8485 0.0249 Inf 0.7930 0.891
vertimec 0.9242 0.0184 Inf 0.8794 0.953
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
<- cld(m, adjust = "Tukey", Letters = LETTERS)
cld_results print(cld_results)
tratamento prob SE df asymp.LCL asymp.UCL .group
controle 0.0606 0.0166 Inf 0.0266 0.132 A
actara 0.5000 0.0347 Inf 0.3990 0.601 B
benevia 0.6970 0.0319 Inf 0.5957 0.782 C
premio 0.7121 0.0314 Inf 0.6115 0.795 CD
nomolt 0.7222 0.0421 Inf 0.5834 0.828 CD
hayate 0.8182 0.0268 Inf 0.7258 0.884 CDE
tracer 0.8485 0.0249 Inf 0.7598 0.908 DEF
avatar 0.8939 0.0214 Inf 0.8126 0.942 EFG
belt 0.9091 0.0200 Inf 0.8307 0.953 EFG
vertimec 0.9242 0.0184 Inf 0.8491 0.964 EFG
talstar 0.9394 0.0166 Inf 0.8679 0.973 FG
cartap 0.9545 0.0145 Inf 0.8871 0.982 G
pirate 0.9545 0.0145 Inf 0.8871 0.982 G
delegate 0.9545 0.0145 Inf 0.8871 0.982 G
joiner 0.9545 0.0145 Inf 0.8871 0.982 G
abamex 0.9545 0.0145 Inf 0.8871 0.982 G
Confidence level used: 0.95
Conf-level adjustment: sidak method for 16 estimates
Intervals are back-transformed from the logit scale
P value adjustment: tukey method for comparing a family of 16 estimates
Tests are performed on the log odds ratio scale
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.
<- as.data.frame(cld_results) |>
cld_df ::rename(group_sig = .group,
dplyremmean = prob,
lower.CL = asymp.LCL,
upper.CL = asymp.UCL)
<- cld_df |>
cld_df mutate(tratamento = str_to_sentence(tratamento))
ggplot(cld_df, aes(x = reorder(tratamento, emmean), y = emmean, fill = group_sig)) +
geom_bar(stat = "identity", fill = "darkgreen", color = "black") +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2, position = position_dodge(0.9)) +
geom_text(aes(y = upper.CL, label = group_sig),
vjust = -0.5, size = 4, color = "black")+
labs(title = "Probabilidade estimada de resposta por tratamento (GLM Quasibinomial)",
x = "Tratamentos",
y = "Probabilidade estimada (EMM)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1)) +
ylim(0, max(cld_df$upper.CL) * 1)
<- summary(modelo_glm_qb)$dispersion
dispersion_param_quasibinomial print(dispersion_param_quasibinomial)
[1] 0.3183159
O valor 0.3183159 indica que os dados têm menos variabilidade (subdispersão) do que o observado pelo modelo binomial simples. O modelo quasibinomial
utiliza esse valor para ajustar os p-valores e erros padrão, tornando a análise mais precisa.
GLMM - Modelo linear generalizado misto
O GLMM é apropriado para os dados, porque a variável é binomial (mortos / vivos), tem um efeito fixo a ser testado (tratamentos) e possui efeitos aleatórios (repetições).
$repetição <- as.factor(dados$repetição)
dados
$resposta_binomial <- cbind(dados$mortos, dados$vivos)
dados
<- glmer(resposta_binomial ~ tratamento + (1 | repetição),
modelo_glmm family = binomial(link = "logit"),
data = dados_ajustados)
summary(modelo_glmm)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: resposta_binomial ~ tratamento + (1 | repetição)
Data: dados_ajustados
AIC BIC logLik -2*log(L) df.resid
261.4 305.0 -113.7 227.4 79
Scaled residuals:
Min 1Q Median 3Q Max
-1.2917 -0.1526 0.5160 0.5299 2.3622
Random effects:
Groups Name Variance Std.Dev.
repetição (Intercept) 6.904e-05 0.008309
Number of obs: 96, groups: repetição, 6
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.6733 0.5009 5.337 9.43e-08 ***
tratamentoactara -2.7346 0.5581 -4.900 9.61e-07 ***
tratamentoavatar -0.3887 0.6569 -0.592 0.554011
tratamentobelt -0.3751 0.6584 -0.570 0.568923
tratamentobenevia -1.9130 0.5663 -3.378 0.000729 ***
tratamentocartap -0.3704 0.6590 -0.562 0.574026
tratamentocontrole -6.3354 0.9336 -6.786 1.15e-11 ***
tratamentodelegate -0.3598 0.6602 -0.545 0.585745
tratamentohayate -1.3913 0.5831 -2.386 0.017021 *
tratamentojoiner -0.3946 0.6562 -0.601 0.547622
tratamentonomolt -1.7774 0.6211 -2.862 0.004216 **
tratamentopirate -0.3672 0.6593 -0.557 0.577625
tratamentopremio -1.8524 0.5677 -3.263 0.001102 **
tratamentotalstar -0.4441 0.6508 -0.682 0.494997
tratamentotracer -0.6936 0.6269 -1.106 0.268522
tratamentovertimec -0.3731 0.6587 -0.566 0.571119
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.573179 (tol = 0.002, component 1)
<- emmeans(modelo_glmm, ~ tratamento, type = "response")
emm <- cld(emm, adjust = "tukey", Letters = LETTERS)
letras
print(letras)
tratamento prob SE df asymp.LCL asymp.UCL .group
controle 0.025 0.0192 Inf 0.00251 0.208 A
actara 0.485 0.0615 Inf 0.31272 0.660 B
benevia 0.681 0.0574 Inf 0.49538 0.823 BC
premio 0.694 0.0567 Inf 0.50829 0.833 BC
nomolt 0.710 0.0756 Inf 0.45340 0.879 BC
hayate 0.783 0.0508 Inf 0.59916 0.897 C
tracer 0.879 0.0402 Inf 0.70443 0.957 C
talstar 0.903 0.0365 Inf 0.73186 0.969 C
joiner 0.907 0.0357 Inf 0.73667 0.971 C
avatar 0.908 0.0356 Inf 0.73723 0.972 C
belt 0.909 0.0354 Inf 0.73852 0.972 C
vertimec 0.909 0.0354 Inf 0.73871 0.972 C
cartap 0.909 0.0354 Inf 0.73896 0.972 C
pirate 0.909 0.0353 Inf 0.73926 0.973 C
delegate 0.910 0.0352 Inf 0.73994 0.973 C
abamex 0.935 0.0303 Inf 0.76797 0.984 C
Confidence level used: 0.95
Conf-level adjustment: sidak method for 16 estimates
Intervals are back-transformed from the logit scale
P value adjustment: tukey method for comparing a family of 16 estimates
Tests are performed on the log odds ratio scale
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.
<- letras |>
letras mutate(tratamento = str_to_sentence(tratamento))
ggplot(as.data.frame(letras), aes(x = reorder (tratamento, prob), y = prob)) +
geom_col(fill = "darkgreen", color = "black") +
geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2) +
geom_text(aes(y = asymp.UCL, label = .group),
vjust = -0.5, size = 4, color = "black") +
labs(title = "Proporção de insetos mortos por tratamento (GLMM)",
y = "Proporção estimada de mortos", x = "Tratamentos") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1)) +
ylim(0, max(cld_df$upper.CL) * 1.05)
<- sum(residuals(modelo_glmm, type = "pearson")^2) / df.residual(modelo_glmm)
dispersion_param_glmm print(dispersion_param_glmm)
[1] 0.4918094
O modelo apresentou um parâmetro de dispersão de 0.4918094. Este valor, é inferior a 1, indicando claramente a presença de subdispersão nos dados, ou seja, a variabilidade observada é menor do que a esperada pela distribuição binomial padrão. Assim é necessário o ajuste do argumento.
Dispformula: o modelo GLMM foi ajustado utilizando o argumento dispformula = ~1
, que permite estimar um um valor de dispersão diferente do padrão. Assim, ele permite inferências mais precisas.
$mortos_ajustados <- ceiling(dados_ajustados$mortos_ajustados)
dados_ajustados$vivos_ajustados <- ceiling(dados_ajustados$vivos_ajustados)
dados_ajustados
<- glmmTMB(cbind(mortos_ajustados, vivos_ajustados) ~ tratamento + (1 | repetição),
modelo_glmm_disp family = binomial(link = "logit"),
data = dados_ajustados,
dispformula = ~1)
summary(modelo_glmm_disp)
Family: binomial ( logit )
Formula: cbind(mortos_ajustados, vivos_ajustados) ~ tratamento + (1 |
repetição)
Data: dados_ajustados
AIC BIC logLik -2*log(L) df.resid
269.5 313.1 -117.7 235.5 79
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
repetição (Intercept) 1.973e-10 1.405e-05
Number of obs: 96, groups: repetição, 6
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.398e+00 4.264e-01 5.624 1.87e-08 ***
tratamentoactara -2.398e+00 4.872e-01 -4.922 8.58e-07 ***
tratamentoavatar -5.733e-01 5.458e-01 -1.050 0.293543
tratamentobelt -4.520e-01 5.557e-01 -0.813 0.416010
tratamentobenevia -1.642e+00 4.957e-01 -3.312 0.000927 ***
tratamentocartap -1.989e-05 6.030e-01 0.000 0.999974
tratamentocontrole -4.626e+00 5.831e-01 -7.933 2.13e-15 ***
tratamentodelegate 1.715e-05 6.030e-01 0.000 0.999977
tratamentohayate -1.063e+00 5.158e-01 -2.061 0.039327 *
tratamentojoiner -1.189e-06 6.030e-01 0.000 0.999998
tratamentonomolt -1.596e+00 5.415e-01 -2.947 0.003214 **
tratamentopirate 3.850e-06 6.030e-01 0.000 0.999995
tratamentopremio -1.577e+00 4.973e-01 -3.171 0.001518 **
tratamentotalstar -1.694e-01 5.831e-01 -0.291 0.771415
tratamentotracer -8.853e-01 5.251e-01 -1.686 0.091777 .
tratamentovertimec -3.185e-01 5.678e-01 -0.561 0.574924
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
<- performance::check_overdispersion(modelo_glmm_disp)
overdisp_check print(overdisp_check)
# Overdispersion test
dispersion ratio = 0.863
p-value = 0.088
<- emmeans(modelo_glmm_disp, ~ tratamento, type = "response")
emm print(emm)
tratamento prob SE df asymp.LCL asymp.UCL
abamex 0.9167 0.0326 Inf 0.8267 0.962
actara 0.5000 0.0589 Inf 0.3865 0.613
avatar 0.8611 0.0408 Inf 0.7607 0.924
belt 0.8750 0.0390 Inf 0.7769 0.934
benevia 0.6806 0.0549 Inf 0.5649 0.778
cartap 0.9167 0.0326 Inf 0.8267 0.962
controle 0.0972 0.0349 Inf 0.0471 0.190
delegate 0.9167 0.0326 Inf 0.8267 0.962
hayate 0.7917 0.0479 Inf 0.6827 0.870
joiner 0.9167 0.0326 Inf 0.8267 0.962
nomolt 0.6905 0.0713 Inf 0.5370 0.811
pirate 0.9167 0.0326 Inf 0.8267 0.962
premio 0.6944 0.0543 Inf 0.5792 0.790
talstar 0.9028 0.0349 Inf 0.8098 0.953
tracer 0.8194 0.0453 Inf 0.7134 0.892
vertimec 0.8889 0.0370 Inf 0.7932 0.943
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
<- cld(emm, adjust = "tukey", Letters = LETTERS)
letras print(letras)
tratamento prob SE df asymp.LCL asymp.UCL .group
controle 0.0972 0.0349 Inf 0.0323 0.258 A
actara 0.5000 0.0589 Inf 0.3330 0.667 B
benevia 0.6806 0.0549 Inf 0.5028 0.818 BC
nomolt 0.6905 0.0713 Inf 0.4547 0.856 BC
premio 0.6944 0.0543 Inf 0.5167 0.829 BC
hayate 0.7917 0.0479 Inf 0.6176 0.899 C
tracer 0.8194 0.0453 Inf 0.6478 0.918 C
avatar 0.8611 0.0408 Inf 0.6942 0.944 C
belt 0.8750 0.0390 Inf 0.7100 0.952 C
vertimec 0.8889 0.0370 Inf 0.7259 0.960 C
talstar 0.9028 0.0349 Inf 0.7419 0.968 C
cartap 0.9167 0.0326 Inf 0.7579 0.975 C
joiner 0.9167 0.0326 Inf 0.7579 0.975 C
abamex 0.9167 0.0326 Inf 0.7579 0.975 C
pirate 0.9167 0.0326 Inf 0.7579 0.975 C
delegate 0.9167 0.0326 Inf 0.7579 0.975 C
Confidence level used: 0.95
Conf-level adjustment: sidak method for 16 estimates
Intervals are back-transformed from the logit scale
P value adjustment: tukey method for comparing a family of 16 estimates
Tests are performed on the log odds ratio scale
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.
<- as.data.frame(letras) |>
letras_df ::rename(group_sig = .group,
dplyremmean = prob,
lower.CL = asymp.LCL,
upper.CL = asymp.UCL)
<- letras_df |>
letras_df mutate(
grupo = case_when(
== "Controle" ~ "Controle",
tratamento < 0.3 ~ "Grupo Baixa Mortalidade",
emmean >= 0.3 & emmean < 0.6 ~ "Grupo Média Mortalidade",
emmean >= 0.6 ~ "Grupo Alta Mortalidade",
emmean TRUE ~ "Outros"
)
)
<- letras_df |>
letras_df mutate(tratamento = str_to_sentence(tratamento))
<- ggplot(letras_df, aes(x = reorder(tratamento, emmean), y = emmean, fill = grupo)) +
graf_glmm geom_bar(stat = "identity", color = "black") +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2, position = position_dodge(0.9)) +
geom_text(aes(y = upper.CL, label = group_sig),
vjust = -0.5, size = 4, color = "black", fontface = "bold") +
labs(title = "Proporção de insetos mortos por tratamento (GLMM Binomial com dispersão)",
x = "Tratamentos",
y = "Proporção estimada de mortos",
fill = "Grupos de tratamento") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1, size = 9)) +
ylim(0, max(letras_df$upper.CL) * 1.05) +
coord_flip() +
scale_fill_manual(values = c(
"Controle" = "grey",
"Grupo Baixa Mortalidade" = "#99d594",
"Grupo Média Mortalidade" = "#67a9cf",
"Grupo Alta Mortalidade" = "#af8dc3",
"Outros" = "lightgrey"
))
print(graf_glmm)
<- sigma(modelo_glmm_disp)
dispersion_param_modelado print(paste("Parâmetro de dispersão (sigma) para Binomial com dispformula:", dispersion_param_modelado))
[1] "Parâmetro de dispersão (sigma) para Binomial com dispformula: 1"
No modelo GLMM ajustado com glmmTMB
e família binomial possui parâmetro de dispersão de 1. Isso indica que, ao considerar os efeitos fixos e aleatórios, a variabilidade em seus dados está perfeitamente alinhada com o que a distribuição binomial prevê. Portanto, não há superdispersão nem subdispersão residual, e as inferências (p-valores e erros padrão) do seu modelo são consideradas válidas e precisas.
Comparação dos modelos
Comparação do índice de dispersão
Foi realizada a comparação do modelo GLM com distribuição quasibinomial e o modelo GLMM binomial com dispformula
. O índice de dispersão do GLM foi de 0.318, indicando subdispersão, ou seja, variação abaixo da esperada, o que pode comprometer a validade estatística dos testes. Já o GLMM apresentou índice de dispersão de 0.863, próximo de 1, o que indica um ajuste mais adequado à variância dos dados. Portanto, o GLMM com dispersão é mais adequado por fornecer estimativas mais robustas e realistas para a inferência dos efeitos dos tratamentos.
<- sum(residuals(modelo_glm_qb, type = "pearson")^2) / df.residual(modelo_glm_qb)
disp_glm
<- check_overdispersion(modelo_glmm_disp)
disp_glmm
cat("GLM quasibinomial - Índice de dispersão:", round(disp_glm, 3), "\n")
GLM quasibinomial - Índice de dispersão: 0.318
cat("GLMM binomial - Índice de dispersão:", round(disp_glmm$dispersion_ratio, 3), "\n")
GLMM binomial - Índice de dispersão: 0.863
Modelo selecionado
print(graf_glmm)
Curva de Probit
A Curva de Probit foi feita para os inseticidas abamex, avatar, benevia, delegate e joiner. A Curva de Probit é um modelo de regressão onde a reposta binária (morto / vivo) é transformada pela função probit
para ajustar uma relação linear com a dose transformada em log. Comum para cálculo da CL50.
Abamex
<- data.frame(
dados_ab conc = c(25, 50, 5, 40, 10, 1),
total = c(60, 60, 60, 60, 60, 60),
mortos = c(36, 56, 17, 41, 27, 8)
)
$mort_prop <- dados_ab$mortos / dados_ab$total
dados_ab$lconc <- log10(dados_ab$conc)
dados_ab
print("Dados Iniciais:")
[1] "Dados Iniciais:"
print(dados_ab)
conc total mortos mort_prop lconc
1 25 60 36 0.6000000 1.39794
2 50 60 56 0.9333333 1.69897
3 5 60 17 0.2833333 0.69897
4 40 60 41 0.6833333 1.60206
5 10 60 27 0.4500000 1.00000
6 1 60 8 0.1333333 0.00000
<- lm(mort_prop ~ lconc, data = dados_ab)
modelo_reg_ab
print(anova(modelo_reg_ab))
Analysis of Variance Table
Response: mort_prop
Df Sum Sq Mean Sq F value Pr(>F)
lconc 1 0.37296 0.37296 36.247 0.003834 **
Residuals 4 0.04116 0.01029
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(modelo_reg_ab))
Call:
lm(formula = mort_prop ~ lconc, data = dados_ab)
Residuals:
1 2 3 4 5 6
-0.05453 0.15114 -0.07476 -0.05776 -0.03576 0.07167
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.06166 0.08577 0.719 0.51197
lconc 0.42410 0.07044 6.021 0.00383 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1014 on 4 degrees of freedom
Multiple R-squared: 0.9006, Adjusted R-squared: 0.8758
F-statistic: 36.25 on 1 and 4 DF, p-value: 0.003834
<- glm(
modelo_probit_ab cbind(mortos, total - mortos) ~ lconc,
data = dados_ab,
family = binomial(link = "probit")
)
print(summary(modelo_probit_ab))
Call:
glm(formula = cbind(mortos, total - mortos) ~ lconc, family = binomial(link = "probit"),
data = dados_ab)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3277 0.1727 -7.69 1.47e-14 ***
lconc 1.2559 0.1391 9.03 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 112.489 on 5 degrees of freedom
Residual deviance: 13.834 on 4 degrees of freedom
AIC: 42.645
Number of Fisher Scoring iterations: 4
<- sum(residuals(modelo_probit_ab, type = "pearson")^2)
pearson_chisq_ab <- nrow(dados_ab) - length(coef(modelo_probit_ab))
df_chisq_ab <- pchisq(pearson_chisq_ab, df = df_chisq_ab, lower.tail = FALSE)
p_valor_chisq_ab
cat(paste("Estatística Qui-quadrado de Pearson:", round(pearson_chisq_ab, 4), "\n"))
Estatística Qui-quadrado de Pearson: 12.0059
cat(paste("Graus de Liberdade (DF):", df_chisq_ab, "\n"))
Graus de Liberdade (DF): 4
cat(paste("Pr > ChiSq (p-valor):", round(p_valor_chisq_ab, 4), "\n\n"))
Pr > ChiSq (p-valor): 0.0173
<- c(
probabilidades_ab seq(from = 0.01, to = 0.10, by = 0.01),
seq(from = 0.15, to = 0.90, by = 0.05),
seq(from = 0.91, to = 0.99, by = 0.01)
)
<- dose.p(modelo_probit_ab, p = probabilidades_ab)
cl_log_ab <- attr(cl_log_ab, "SE")
erros_padrao_log_ab <- qnorm(0.975)
z_valor_ab <- cl_log_ab - z_valor_ab * erros_padrao_log_ab
limite_inf_log_ab <- cl_log_ab + z_valor_ab * erros_padrao_log_ab
limite_sup_log_ab
<- data.frame(
tabela_escala_log10_ab Probabilidade = probabilidades_ab,
Log10_conc = cl_log_ab,
Limite_Inferior_95 = limite_inf_log_ab,
Limite_Superior_95 = limite_sup_log_ab
)
<- data.frame(
tabela_escala_conc_ab Probabilidade = probabilidades_ab,
Concentracao = 10^cl_log_ab,
Limite_Inferior_95 = 10^limite_inf_log_ab,
Limite_Superior_95 = 10^limite_sup_log_ab
)
options(digits = 6)
print(tabela_escala_log10_ab)
Probabilidade Log10_conc Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 -0.7951538 -1.22739826 -0.362909326
p = 0.02: 0.02 -0.5780975 -0.96506165 -0.191133407
p = 0.03: 0.03 -0.4403823 -0.79886222 -0.081902361
p = 0.04: 0.04 -0.3367845 -0.67399734 0.000428382
p = 0.05: 0.05 -0.2525157 -0.57255226 0.067520960
p = 0.06: 0.06 -0.1807896 -0.48630846 0.124729171
p = 0.07: 0.07 -0.1179000 -0.41077828 0.174978340
p = 0.08: 0.08 -0.0615898 -0.34322997 0.220050443
p = 0.09: 0.09 -0.0103779 -0.28187131 0.261115540
p = 0.10: 0.10 0.0367627 -0.22545999 0.298985453
p = 0.15: 0.15 0.2319375 0.00720882 0.456666155
p = 0.20: 0.20 0.3870562 0.19070952 0.583402870
p = 0.25: 0.25 0.5201342 0.34662370 0.693644762
p = 0.30: 0.30 0.6396424 0.48490696 0.794377881
p = 0.35: 0.35 0.7503846 0.61098060 0.889788545
p = 0.40: 0.40 0.8554680 0.72809962 0.982836463
p = 0.45: 0.45 0.9571375 0.83836664 1.075908390
p = 0.50: 0.50 1.0571951 0.94328837 1.171101749
p = 0.55: 0.55 1.1572526 1.04417649 1.270328710
p = 0.60: 0.60 1.2589221 1.14247116 1.375372992
p = 0.65: 0.65 1.3640055 1.23996952 1.488041569
p = 0.70: 0.70 1.4747477 1.33897868 1.610516713
p = 0.75: 0.75 1.5942559 1.44253545 1.745976327
p = 0.80: 0.80 1.7273339 1.55497045 1.899697401
p = 0.85: 0.85 1.8824526 1.68342186 2.081483404
p = 0.90: 0.90 2.0776274 1.84247724 2.312777530
p = 0.91: 0.91 2.1247680 1.88059292 2.368943083
p = 0.92: 0.92 2.1759799 1.92189677 2.430062992
p = 0.93: 0.93 2.2322901 1.96720203 2.497378147
p = 0.94: 0.94 2.2951798 2.01768047 2.572679060
p = 0.95: 0.95 2.3669058 2.07511632 2.658695220
p = 0.96: 0.96 2.4511746 2.14243816 2.759911038
p = 0.97: 0.97 2.5547724 2.22500506 2.884539759
p = 0.98: 0.98 2.6924876 2.33448953 3.050485764
p = 0.99: 0.99 2.9095439 2.50656436 3.312523467
print(tabela_escala_conc_ab)
Probabilidade Concentracao Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 0.160268 0.0592382 0.433601
p = 0.02: 0.02 0.264182 0.1083773 0.643971
p = 0.03: 0.03 0.362759 0.1589051 0.828128
p = 0.04: 0.04 0.460485 0.2118374 1.000987
p = 0.05: 0.05 0.559093 0.2675764 1.168210
p = 0.06: 0.06 0.659493 0.3263560 1.332690
p = 0.07: 0.07 0.762255 0.3883486 1.496161
p = 0.08: 0.08 0.867781 0.4537013 1.659780
p = 0.09: 0.09 0.976387 0.5225510 1.824381
p = 0.10: 0.10 1.088335 0.5950316 1.990607
p = 0.15: 0.15 1.705837 1.0167374 2.861977
p = 0.20: 0.20 2.438126 1.5513490 3.831800
p = 0.25: 0.25 3.312335 2.2213843 4.939065
p = 0.30: 0.30 4.361566 3.0542668 6.228420
p = 0.35: 0.35 5.628395 4.0830115 7.758693
p = 0.40: 0.40 7.169156 5.3468700 9.612502
p = 0.45: 0.45 9.060194 6.8923392 11.909908
p = 0.50: 0.50 11.407620 8.7758334 14.828655
p = 0.55: 0.55 14.363246 11.0707360 18.634971
p = 0.60: 0.60 18.151899 13.8826112 23.734112
p = 0.65: 0.65 23.120943 17.3767889 30.763913
p = 0.70: 0.70 29.836488 21.8262276 40.786526
p = 0.75: 0.75 39.287635 27.7035515 55.715538
p = 0.80: 0.80 53.374513 35.8897514 79.377497
p = 0.85: 0.85 76.287368 48.2416175 120.637799
p = 0.90: 0.90 119.571420 69.5788491 205.483772
p = 0.91: 0.91 133.280927 75.9613937 233.853074
p = 0.92: 0.92 149.961537 83.5404427 269.192522
p = 0.93: 0.93 170.722235 92.7261068 314.324438
p = 0.94: 0.94 197.323933 104.1550825 373.834226
p = 0.95: 0.95 232.758619 118.8820601 455.716989
p = 0.96: 0.96 282.601588 138.8155627 575.322075
p = 0.97: 0.97 358.733892 167.8823573 766.548713
p = 0.98: 0.98 492.592334 216.0177976 1123.274149
p = 0.99: 0.99 811.977346 321.0438527 2053.635989
<- data.frame(lconc = seq(0, max(dados_ab$lconc) + 0.1, length.out = 200))
lconc_grafico_ab
<- predict(modelo_probit_ab, newdata = lconc_grafico_ab, type = "link", se.fit = TRUE)
predicoes_ab
<- qnorm(0.975)
z_valor_ab $prob_predita_ab <- pnorm(predicoes_ab$fit)
lconc_grafico_ab$prob_superior_ab <- pnorm(predicoes_ab$fit + z_valor_ab * predicoes_ab$se.fit)
lconc_grafico_ab$prob_inferior_ab <- pnorm(predicoes_ab$fit - z_valor_ab * predicoes_ab$se.fit)
lconc_grafico_ab
<- as.numeric(dose.p(modelo_probit_ab, p = 0.50))
lc50_log_ab
<- 10^lc50_log_ab
lc50_conc_ab
<- -1.5
x_lim_inferior_ab <- 1
x_lim_superior_ab
<- data.frame(lconc = seq(x_lim_inferior_ab, x_lim_superior_ab, length.out = 200))
lconc_grafico_ab
<- predict(modelo_probit_ab, newdata = lconc_grafico_ab, type = "link", se.fit = TRUE)
predicoes_ab <- qnorm(0.975)
z_valor_ab $prob_predita_ab <- pnorm(predicoes_ab$fit)
lconc_grafico_ab$prob_superior_ab <- pnorm(predicoes_ab$fit + z_valor_ab * predicoes_ab$se.fit)
lconc_grafico_ab$prob_inferior <- pnorm(predicoes_ab$fit - z_valor_ab * predicoes_ab$se.fit)
lconc_grafico_ab
<- ggplot() +
grafico_probit_final_ab geom_ribbon(data = lconc_grafico_ab,
aes(x = lconc, ymin = prob_inferior * 100, ymax = prob_superior_ab * 100),
fill = "skyblue", alpha = 0.5) +
geom_line(data = lconc_grafico_ab,
aes(x = lconc, y = prob_predita_ab * 100),
color = "blue", size = 1) +
geom_point(data = dados_ab,
aes(x = lconc, y = mort_prop * 100),
color = "red", size = 4) +
geom_hline(yintercept = 50, linetype = "dashed", color = "black", size = 0.7) +
geom_vline(xintercept = lc50_log_ab, linetype = "dashed", color = "black", size = 0.7) +
labs(
title = "Curva dose-resposta (Modelo Probit) - Abamex",
x = "Log10 (Concentração)",
y = "Mortalidade observada e prevista (%)"
+
) annotate("text",
x = lc50_log_ab,
y = 45,
label = paste("LC50 =", format(lc50_conc_ab, digits = 4)),
hjust = -0.1,
vjust = 1,
fontface = "bold",
color = "black") +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
scale_x_continuous(limits = c(x_lim_inferior_ab, x_lim_superior_ab), breaks = seq(-4, 1, 0.5)) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
print(grafico_probit_final_ab)
Avatar
<- data.frame(
dados_av conc = c(0.1, 1, 5, 10, 25, 40, 0.3, 0.5, 0.7),
total = c(60, 60, 60, 60, 60, 60, 60, 60, 60),
mortos = c(0, 4, 16, 1, 15, 14, 3, 8, 6)
)
$mort_prop_av <- dados_av$mortos / dados_av$total
dados_av$lconc_av <- log10(dados_av$conc)
dados_av
print("Dados Iniciais:")
[1] "Dados Iniciais:"
print(dados_av)
conc total mortos mort_prop_av lconc_av
1 0.1 60 0 0.0000000 -1.000000
2 1.0 60 4 0.0666667 0.000000
3 5.0 60 16 0.2666667 0.698970
4 10.0 60 1 0.0166667 1.000000
5 25.0 60 15 0.2500000 1.397940
6 40.0 60 14 0.2333333 1.602060
7 0.3 60 3 0.0500000 -0.522879
8 0.5 60 8 0.1333333 -0.301030
9 0.7 60 6 0.1000000 -0.154902
<- lm(mort_prop_av ~ lconc_av, data = dados_av)
modelo_reg_av
print(anova(modelo_reg_av))
Analysis of Variance Table
Response: mort_prop_av
Df Sum Sq Mean Sq F value Pr(>F)
lconc_av 1 0.03911 0.03911 6.03 0.0438 *
Residuals 7 0.04540 0.00649
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(modelo_reg_av))
Call:
lm(formula = mort_prop_av ~ lconc_av, data = dados_av)
Residuals:
Min 1Q Median 3Q Max
-0.16122 -0.02365 0.00902 0.04143 0.11200
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1008 0.0285 3.54 0.0095 **
lconc_av 0.0771 0.0314 2.46 0.0438 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0805 on 7 degrees of freedom
Multiple R-squared: 0.463, Adjusted R-squared: 0.386
F-statistic: 6.03 on 1 and 7 DF, p-value: 0.0438
<- glm(
modelo_probit_av cbind(mortos, total - mortos) ~ lconc_av,
data = dados_av,
family = binomial(link = "probit")
)
print(summary(modelo_probit_av))
Call:
glm(formula = cbind(mortos, total - mortos) ~ lconc_av, family = binomial(link = "probit"),
data = dados_av)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3420 0.0864 -15.54 < 2e-16 ***
lconc_av 0.3982 0.0858 4.64 3.4e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 53.183 on 8 degrees of freedom
Residual deviance: 30.874 on 7 degrees of freedom
AIC: 63.15
Number of Fisher Scoring iterations: 4
<- sum(residuals(modelo_probit_av, type = "pearson")^2)
pearson_chisq_av <- nrow(dados_av) - length(coef(modelo_probit_av))
df_chisq_av <- pchisq(pearson_chisq_av, df = df_chisq_av, lower.tail = FALSE)
p_valor_chisq_av
cat(paste("Estatística Qui-quadrado de Pearson:", round(pearson_chisq_av, 4), "\n"))
Estatística Qui-quadrado de Pearson: 24.7946
cat(paste("Graus de Liberdade (DF):", df_chisq_av, "\n"))
Graus de Liberdade (DF): 7
cat(paste("Pr > ChiSq (p-valor):", round(p_valor_chisq_av, 4), "\n\n"))
Pr > ChiSq (p-valor): 8e-04
<- c(
probabilidades_av seq(from = 0.01, to = 0.10, by = 0.01),
seq(from = 0.15, to = 0.90, by = 0.05),
seq(from = 0.91, to = 0.99, by = 0.01)
)
<- dose.p(modelo_probit_av, p = probabilidades_av)
cl_log_av <- attr(cl_log_av, "SE")
erros_padrao_log_av <- qnorm(0.975)
z_valor_av <- cl_log_av - z_valor_av * erros_padrao_log_av
limite_inf_log_av <- cl_log_av + z_valor_av * erros_padrao_log_av
limite_sup_log_av
<- data.frame(
tabela_escala_log10_av Probabilidade = probabilidades_av,
Log10_conc = cl_log_av,
Limite_Inferior_95 = limite_inf_log_av,
Limite_Superior_95 = limite_sup_log_av
)
<- data.frame(
tabela_escala_conc_av Probabilidade = probabilidades_av,
Concentracao = 10^cl_log_av,
Limite_Inferior_95 = 10^limite_inf_log_av,
Limite_Superior_95 = 10^limite_sup_log_av
)
options(digits = 6)
print(tabela_escala_log10_av)
Probabilidade Log10_conc Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 -2.47191568 -3.801554 -1.1422769
p = 0.02: 0.02 -1.78732722 -2.841054 -0.7336000
p = 0.03: 0.03 -1.35297778 -2.235993 -0.4699627
p = 0.04: 0.04 -1.02623362 -1.784726 -0.2677412
p = 0.05: 0.05 -0.76045246 -1.421591 -0.0993141
p = 0.06: 0.06 -0.53423094 -1.116726 0.0482641
p = 0.07: 0.07 -0.33587892 -0.854099 0.1823410
p = 0.08: 0.08 -0.15827835 -0.624221 0.3076640
p = 0.09: 0.09 0.00324226 -0.421105 0.4275900
p = 0.10: 0.10 0.15192226 -0.240763 0.5446074
p = 0.15: 0.15 0.76749714 0.405663 1.1293311
p = 0.20: 0.20 1.25673655 0.799861 1.7136125
p = 0.25: 0.25 1.67646040 1.090396 2.2625246
p = 0.30: 0.30 2.05338538 1.333682 2.7730886
p = 0.35: 0.35 2.40266256 1.551269 3.2540557
p = 0.40: 0.40 2.73409244 1.753585 3.7146003
p = 0.45: 0.45 3.05475469 1.946831 4.1626785
p = 0.50: 0.50 3.37033296 2.135362 4.6053037
p = 0.55: 0.55 3.68591122 2.322716 5.0491063
p = 0.60: 0.60 4.00657347 2.512194 5.5009526
p = 0.65: 0.65 4.33800335 2.707318 5.9686887
p = 0.70: 0.70 4.68728053 2.912343 6.4622180
p = 0.75: 0.75 5.06420551 3.133058 6.9953534
p = 0.80: 0.80 5.48392936 3.378323 7.5895356
p = 0.85: 0.85 5.97316877 3.663688 8.2826493
p = 0.90: 0.90 6.58874366 4.022142 9.1553451
p = 0.91: 0.91 6.73742365 4.108639 9.3662080
p = 0.92: 0.92 6.89894426 4.202576 9.5953122
p = 0.93: 0.93 7.07654484 4.305831 9.8472583
p = 0.94: 0.94 7.27489685 4.421112 10.1286812
p = 0.95: 0.95 7.50111837 4.552546 10.4496909
p = 0.96: 0.96 7.76689954 4.706907 10.8268919
p = 0.97: 0.97 8.09364370 4.896601 11.2906866
p = 0.98: 0.98 8.52799313 5.148655 11.9073310
p = 0.99: 0.99 9.21258160 5.545711 12.8794519
print(tabela_escala_conc_av)
Probabilidade Concentracao Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 3.37353e-03 1.57923e-04 7.20648e-02
p = 0.02: 0.02 1.63182e-02 1.44193e-03 1.84672e-01
p = 0.03: 0.03 4.43631e-02 5.80774e-03 3.38873e-01
p = 0.04: 0.04 9.41383e-02 1.64163e-02 5.39832e-01
p = 0.05: 0.05 1.73599e-01 3.78799e-02 7.95584e-01
p = 0.06: 0.06 2.92260e-01 7.64318e-02 1.11754e+00
p = 0.07: 0.07 4.61446e-01 1.39927e-01 1.52174e+00
p = 0.08: 0.08 6.94579e-01 2.37563e-01 2.03079e+00
p = 0.09: 0.09 1.00749e+00 3.79223e-01 2.67664e+00
p = 0.10: 0.10 1.41880e+00 5.74430e-01 3.50435e+00
p = 0.15: 0.15 5.85460e+00 2.54486e+00 1.34689e+01
p = 0.20: 0.20 1.80608e+01 6.30755e+00 5.17145e+01
p = 0.25: 0.25 4.74745e+01 1.23139e+01 1.83031e+02
p = 0.30: 0.30 1.13080e+02 2.15617e+01 5.93046e+02
p = 0.35: 0.35 2.52733e+02 3.55852e+01 1.79496e+03
p = 0.40: 0.40 5.42116e+02 5.67002e+01 5.18323e+03
p = 0.45: 0.45 1.13437e+03 8.84771e+01 1.45438e+04
p = 0.50: 0.50 2.34603e+03 1.36572e+02 4.02999e+04
p = 0.55: 0.55 4.85189e+03 2.10240e+02 1.11971e+05
p = 0.60: 0.60 1.01525e+04 3.25233e+02 3.16922e+05
p = 0.65: 0.65 2.17773e+04 5.09704e+02 9.30441e+05
p = 0.70: 0.70 4.86721e+04 8.17228e+02 2.89880e+06
p = 0.75: 0.75 1.15933e+05 1.35849e+03 9.89358e+06
p = 0.80: 0.80 3.04740e+05 2.38959e+03 3.88629e+07
p = 0.85: 0.85 9.40089e+05 4.60987e+03 1.91712e+08
p = 0.90: 0.90 3.87921e+06 1.05231e+04 1.43003e+09
p = 0.91: 0.91 5.46290e+06 1.28422e+04 2.32385e+09
p = 0.92: 0.92 7.92400e+06 1.59432e+04 3.93833e+09
p = 0.93: 0.93 1.19274e+07 2.02223e+04 7.03491e+09
p = 0.94: 0.94 1.88320e+07 2.63701e+04 1.34487e+10
p = 0.95: 0.95 3.17043e+07 3.56899e+04 2.81638e+10
p = 0.96: 0.96 5.84655e+07 5.09222e+04 6.71262e+10
p = 0.97: 0.97 1.24063e+08 7.88135e+04 1.95293e+11
p = 0.98: 0.98 3.37282e+08 1.40817e+05 8.07850e+11
p = 0.99: 0.99 1.63148e+09 3.51327e+05 7.57621e+12
<- data.frame(lconc_av = seq(0, max(dados_av$lconc_av) + 0.1, length.out = 200))
lconc_grafico_av
<- predict(modelo_probit_av, newdata = lconc_grafico_av, type = "link", se.fit = TRUE)
predicoes_av
<- qnorm(0.975)
z_valor_av $prob_predita_av <- pnorm(predicoes_av$fit)
lconc_grafico_av$prob_superior_av <- pnorm(predicoes_av$fit + z_valor_av * predicoes_av$se.fit)
lconc_grafico_av$prob_inferior_ab <- pnorm(predicoes_av$fit - z_valor_av * predicoes_av$se.fit)
lconc_grafico_av
<- as.numeric(dose.p(modelo_probit_av, p = 0.50))
lc50_log_av
<- 10^lc50_log_av
lc50_conc_av
<- -1.5
x_lim_inferior_av <- 1
x_lim_superior_av
<- data.frame(lconc_av = seq(x_lim_inferior_av, x_lim_superior_av, length.out = 200))
lconc_grafico_av
<- predict(modelo_probit_av, newdata = lconc_grafico_av, type = "link", se.fit = TRUE)
predicoes_av <- qnorm(0.975)
z_valor_av $prob_predita_av <- pnorm(predicoes_av$fit)
lconc_grafico_av$prob_superior_av <- pnorm(predicoes_av$fit + z_valor_av * predicoes_av$se.fit)
lconc_grafico_av$prob_inferior_av <- pnorm(predicoes_av$fit - z_valor_av * predicoes_av$se.fit)
lconc_grafico_av
<- ggplot() +
grafico_probit_final_av geom_ribbon(data = lconc_grafico_av,
aes(x = lconc_av, ymin = prob_inferior_av * 100, ymax = prob_superior_av * 100),
fill = "skyblue", alpha = 0.5) +
geom_line(data = lconc_grafico_av,
aes(x = lconc_av, y = prob_predita_av * 100),
color = "blue", size = 1) +
geom_point(data = dados_av,
aes(x = lconc_av, y = mort_prop_av * 100),
color = "red", size = 4) +
geom_hline(yintercept = 50, linetype = "dashed", color = "black", size = 0.7) +
geom_vline(xintercept = lc50_log_av, linetype = "dashed", color = "black", size = 0.7) +
labs(
title = "Curva dose-resposta (Modelo Probit) - Avatar",
x = "Log10 (Concentração)",
y = "Mortalidade observada e prevista (%)"
+
) annotate("text",
x = lc50_log_av,
y = 45,
label = paste("LC50 =", format(lc50_conc_av, digits = 4)),
hjust = -0.1,
vjust = 1,
fontface = "bold",
color = "black") +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
scale_x_continuous(limits = c(x_lim_inferior_av, x_lim_superior_av), breaks = seq(-4, 1, 0.5)) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
print(grafico_probit_final_av)
Benevia
<- data.frame(
dados_be conc = c(0.1, 1, 5, 10, 25, 40, 60),
total = c(60, 60, 60, 60, 60, 50, 50),
mortos = c(1, 2, 8, 5, 44, 6, 8)
)
$mort_prop_be <- dados_be$mortos / dados_be$total
dados_be$lconc_be <- log10(dados_be$conc)
dados_be
print("Dados Iniciais:")
[1] "Dados Iniciais:"
print(dados_be)
conc total mortos mort_prop_be lconc_be
1 0.1 60 1 0.0166667 -1.00000
2 1.0 60 2 0.0333333 0.00000
3 5.0 60 8 0.1333333 0.69897
4 10.0 60 5 0.0833333 1.00000
5 25.0 60 44 0.7333333 1.39794
6 40.0 50 6 0.1200000 1.60206
7 60.0 50 8 0.1600000 1.77815
<- lm(mort_prop_be ~ lconc_be, data = dados_be)
modelo_reg_be
print(anova(modelo_reg_be))
Analysis of Variance Table
Response: mort_prop_be
Df Sum Sq Mean Sq F value Pr(>F)
lconc_be 1 0.07406 0.07406 1.252 0.314
Residuals 5 0.29577 0.05915
print(summary(modelo_reg_be))
Call:
lm(formula = mort_prop_be ~ lconc_be, data = dados_be)
Residuals:
1 2 3 4 5 6 7
0.0338 -0.0618 -0.0402 -0.1239 0.4814 -0.1548 -0.1346
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0951 0.1208 0.79 0.47
lconc_be 0.1122 0.1003 1.12 0.31
Residual standard error: 0.243 on 5 degrees of freedom
Multiple R-squared: 0.2, Adjusted R-squared: 0.0403
F-statistic: 1.25 on 1 and 5 DF, p-value: 0.314
<- glm(
modelo_probit_be cbind(mortos, total - mortos) ~ lconc_be,
data = dados_be,
family = binomial(link = "probit")
)
print(summary(modelo_probit_be))
Call:
glm(formula = cbind(mortos, total - mortos) ~ lconc_be, family = binomial(link = "probit"),
data = dados_be)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.535 0.153 -10.03 <2e-16 ***
lconc_be 0.653 0.118 5.54 3e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 123.61 on 6 degrees of freedom
Residual deviance: 82.54 on 5 degrees of freedom
AIC: 109.9
Number of Fisher Scoring iterations: 4
<- sum(residuals(modelo_probit_be, type = "pearson")^2)
pearson_chisq_be <- nrow(dados_be) - length(coef(modelo_probit_be))
df_chisq_be <- pchisq(pearson_chisq_be, df = df_chisq_be, lower.tail = FALSE)
p_valor_chisq_be
cat(paste("Estatística Qui-quadrado de Pearson:", round(pearson_chisq_be, 4), "\n"))
Estatística Qui-quadrado de Pearson: 88.801
cat(paste("Graus de Liberdade (DF):", df_chisq_be, "\n"))
Graus de Liberdade (DF): 5
cat(paste("Pr > ChiSq (p-valor):", round(p_valor_chisq_be, 4), "\n\n"))
Pr > ChiSq (p-valor): 0
<- c(
probabilidades_be seq(from = 0.01, to = 0.10, by = 0.01),
seq(from = 0.15, to = 0.90, by = 0.05),
seq(from = 0.91, to = 0.99, by = 0.01)
)
<- dose.p(modelo_probit_be, p = probabilidades_be)
cl_log_be <- attr(cl_log_be, "SE")
erros_padrao_log_be <- qnorm(0.975)
z_valor_be <- cl_log_be - z_valor_be * erros_padrao_log_be
limite_inf_log_be <- cl_log_be + z_valor_be * erros_padrao_log_be
limite_sup_log_be
<- data.frame(
tabela_escala_log10_be Probabilidade = probabilidades_be,
Log10_conc = cl_log_be,
Limite_Inferior_95 = limite_inf_log_be,
Limite_Superior_95 = limite_sup_log_be
)
<- data.frame(
tabela_escala_conc_be Probabilidade = probabilidades_be,
Concentracao = 10^cl_log_be,
Limite_Inferior_95 = 10^limite_inf_log_be,
Limite_Superior_95 = 10^limite_sup_log_be
)
options(digits = 6)
print(tabela_escala_log10_be)
Probabilidade Log10_conc Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 -1.2111644 -2.0676481 -0.3546806
p = 0.02: 0.02 -0.7940017 -1.5096138 -0.0783896
p = 0.03: 0.03 -0.5293253 -1.1571410 0.0984903
p = 0.04: 0.04 -0.3302196 -0.8932204 0.2327811
p = 0.05: 0.05 -0.1682625 -0.6796323 0.3431072
p = 0.06: 0.06 -0.0304116 -0.4988685 0.4380453
p = 0.07: 0.07 0.0904568 -0.3413933 0.5223069
p = 0.08: 0.08 0.1986799 -0.2014262 0.5987861
p = 0.09: 0.09 0.2971046 -0.0751984 0.6694075
p = 0.10: 0.10 0.3877046 0.0398793 0.7355299
p = 0.15: 0.15 0.7628130 0.4985345 1.0270915
p = 0.20: 0.20 1.0609372 0.8279282 1.2939463
p = 0.25: 0.25 1.3167013 1.0744930 1.5589095
p = 0.30: 0.30 1.5463853 1.2697424 1.8230282
p = 0.35: 0.35 1.7592218 1.4353964 2.0830472
p = 0.40: 0.40 1.9611828 1.5840611 2.3383045
p = 0.45: 0.45 2.1565824 1.7229089 2.5902559
p = 0.50: 0.50 2.3488840 1.8564217 2.8413463
p = 0.55: 0.55 2.5411856 1.9878208 3.0945505
p = 0.60: 0.60 2.7365852 2.1198163 3.3533542
p = 0.65: 0.65 2.9385462 2.2550843 3.6220082
p = 0.70: 0.70 3.1513827 2.3967012 3.9060642
p = 0.75: 0.75 3.3810668 2.5487306 4.2134029
p = 0.80: 0.80 3.6368308 2.7172984 4.5563632
p = 0.85: 0.85 3.9349550 2.9130737 4.9568364
p = 0.90: 0.90 4.3100634 3.1586175 5.4615093
p = 0.91: 0.91 4.4006635 3.2178217 5.5835052
p = 0.92: 0.92 4.4990881 3.2821013 5.7160749
p = 0.93: 0.93 4.6073113 3.3527384 5.8618841
p = 0.94: 0.94 4.7281796 3.4315819 6.0247773
p = 0.95: 0.95 4.8660305 3.5214484 6.2106126
p = 0.96: 0.96 5.0279877 3.6269631 6.4290122
p = 0.97: 0.97 5.2270934 3.7565926 6.6975941
p = 0.98: 0.98 5.4917697 3.9287841 7.0547553
p = 0.99: 0.99 5.9089324 4.1999355 7.6179293
print(tabela_escala_conc_be)
Probabilidade Concentracao Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 6.14944e-02 8.55760e-03 4.41895e-01
p = 0.02: 0.02 1.60694e-01 3.09304e-02 8.34854e-01
p = 0.03: 0.03 2.95580e-01 6.96400e-02 1.25456e+00
p = 0.04: 0.04 4.67499e-01 1.27873e-01 1.70915e+00
p = 0.05: 0.05 6.78793e-01 2.09107e-01 2.20347e+00
p = 0.06: 0.06 9.32370e-01 3.17053e-01 2.74186e+00
p = 0.07: 0.07 1.23156e+00 4.55624e-01 3.32895e+00
p = 0.08: 0.08 1.58008e+00 6.28889e-01 3.96996e+00
p = 0.09: 0.09 1.98200e+00 8.41011e-01 4.67097e+00
p = 0.10: 0.10 2.44177e+00 1.09617e+00 5.43914e+00
p = 0.15: 0.15 5.79179e+00 3.15162e+00 1.06437e+01
p = 0.20: 0.20 1.15063e+01 6.72865e+00 1.96764e+01
p = 0.25: 0.25 2.07349e+01 1.18712e+01 3.62168e+01
p = 0.30: 0.30 3.51872e+01 1.86098e+01 6.65316e+01
p = 0.35: 0.35 5.74410e+01 2.72519e+01 1.21073e+02
p = 0.40: 0.40 9.14498e+01 3.83761e+01 2.17924e+02
p = 0.45: 0.45 1.43411e+02 5.28334e+01 3.89274e+02
p = 0.50: 0.50 2.23298e+02 7.18492e+01 6.93979e+02
p = 0.55: 0.55 3.47685e+02 9.72346e+01 1.24323e+03
p = 0.60: 0.60 5.45237e+02 1.31770e+02 2.25608e+03
p = 0.65: 0.65 8.68053e+02 1.79922e+02 4.18801e+03
p = 0.70: 0.70 1.41704e+03 2.49288e+02 8.05498e+03
p = 0.75: 0.75 2.40473e+03 3.53778e+02 1.63457e+04
p = 0.80: 0.80 4.33342e+03 5.21553e+02 3.60050e+04
p = 0.85: 0.85 8.60905e+03 8.18604e+02 9.05391e+04
p = 0.90: 0.90 2.04204e+04 1.44085e+03 2.89407e+05
p = 0.91: 0.91 2.51573e+04 1.65128e+03 3.83270e+05
p = 0.92: 0.92 3.15564e+04 1.91470e+03 5.20086e+05
p = 0.93: 0.93 4.04866e+04 2.25288e+03 7.27586e+05
p = 0.94: 0.94 5.34785e+04 2.70136e+03 1.05871e+06
p = 0.95: 0.95 7.34566e+04 3.32237e+03 1.62410e+06
p = 0.96: 0.96 1.06657e+05 4.23607e+03 2.68542e+06
p = 0.97: 0.97 1.68692e+05 5.70943e+03 4.98418e+06
p = 0.98: 0.98 3.10291e+05 8.48758e+03 1.13437e+07
p = 0.99: 0.99 8.10835e+05 1.58466e+04 4.14886e+07
<- data.frame(lconc_be = seq(0, max(dados_be$lconc_be) + 0.1, length.out = 200))
lconc_grafico_be
<- predict(modelo_probit_be, newdata = lconc_grafico_be, type = "link", se.fit = TRUE)
predicoes_be
<- qnorm(0.975)
z_valor_be $prob_predita_be <- pnorm(predicoes_be$fit)
lconc_grafico_be$prob_superior_be <- pnorm(predicoes_be$fit + z_valor_be * predicoes_be$se.fit)
lconc_grafico_be$prob_inferior_be <- pnorm(predicoes_be$fit - z_valor_be * predicoes_be$se.fit)
lconc_grafico_be
<- as.numeric(dose.p(modelo_probit_be, p = 0.50))
lc50_log_be
<- 10^lc50_log_be
lc50_conc_be
<- -1.5
x_lim_inferior_be <- 1
x_lim_superior_be
<- data.frame(lconc_be = seq(x_lim_inferior_be, x_lim_superior_be, length.out = 200))
lconc_grafico_be
<- predict(modelo_probit_be, newdata = lconc_grafico_be, type = "link", se.fit = TRUE)
predicoes_be <- qnorm(0.975)
z_valor_be $prob_predita_be <- pnorm(predicoes_be$fit)
lconc_grafico_be$prob_superior_be <- pnorm(predicoes_be$fit + z_valor_be * predicoes_be$se.fit)
lconc_grafico_be$prob_inferior_be <- pnorm(predicoes_be$fit - z_valor_be * predicoes_be$se.fit)
lconc_grafico_be
<- ggplot() +
grafico_probit_final_be geom_ribbon(data = lconc_grafico_be,
aes(x = lconc_be, ymin = prob_inferior_be * 100, ymax = prob_superior_be * 100),
fill = "skyblue", alpha = 0.5) +
geom_line(data = lconc_grafico_be,
aes(x = lconc_be, y = prob_predita_be * 100),
color = "blue", size = 1) +
geom_point(data = dados_be,
aes(x = lconc_be, y = mort_prop_be * 100),
color = "red", size = 4) +
geom_hline(yintercept = 50, linetype = "dashed", color = "black", size = 0.7) +
geom_vline(xintercept = lc50_log_be, linetype = "dashed", color = "black", size = 0.7) +
labs(
title = "Curva dose-resposta (Modelo Probit) - Benevia",
x = "Log10 (Concentração)",
y = "Mortalidade observada e prevista (%)"
+
) annotate("text",
x = lc50_log_be,
y = 45,
label = paste("LC50 =", format(lc50_conc_be, digits = 4)),
hjust = -0.1,
vjust = 1,
fontface = "bold",
color = "black") +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
scale_x_continuous(limits = c(x_lim_inferior_be, x_lim_superior_be), breaks = seq(-4, 1, 0.5)) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
print(grafico_probit_final_be)
Delegate
<- data.frame(
dados_de conc = c(0.1, 1),
total = c(60, 60),
mortos = c(0,0)
)
$mort_prop_de <- dados_de$mortos / dados_de$total
dados_de$lconc_de <- log10(dados_de$conc)
dados_de
print("Dados Iniciais:")
[1] "Dados Iniciais:"
print(dados_de)
conc total mortos mort_prop_de lconc_de
1 0.1 60 0 0 -1
2 1.0 60 0 0 0
<- lm(mort_prop_de ~ lconc_de, data = dados_de)
modelo_reg_de
print(anova(modelo_reg_de))
Analysis of Variance Table
Response: mort_prop_de
Df Sum Sq Mean Sq F value Pr(>F)
lconc_de 1 0 0 NaN NaN
Residuals 0 0 NaN
print(summary(modelo_reg_de))
Call:
lm(formula = mort_prop_de ~ lconc_de, data = dados_de)
Residuals:
ALL 2 residuals are 0: no residual degrees of freedom!
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0 NaN NaN NaN
lconc_de 0 NaN NaN NaN
Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared: NaN, Adjusted R-squared: NaN
F-statistic: NaN on 1 and 0 DF, p-value: NA
<- glm(
modelo_probit_de cbind(mortos, total - mortos) ~ lconc_de,
data = dados_de,
family = binomial(link = "probit")
)
print(summary(modelo_probit_de))
Call:
glm(formula = cbind(mortos, total - mortos) ~ lconc_de, family = binomial(link = "probit"),
data = dados_de)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.99e+00 9.44e+03 0 1
lconc_de 1.45e-15 1.33e+04 0 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 0.0000e+00 on 1 degrees of freedom
Residual deviance: 3.3515e-10 on 0 degrees of freedom
AIC: 4
Number of Fisher Scoring iterations: 22
<- sum(residuals(modelo_probit_de, type = "pearson")^2)
pearson_chisq_de <- nrow(dados_de) - length(coef(modelo_probit_de))
df_chisq_de <- pchisq(pearson_chisq_de, df = df_chisq_de, lower.tail = FALSE)
p_valor_chisq_de
cat(paste("Estatística Qui-quadrado de Pearson:", round(pearson_chisq_de, 4), "\n"))
Estatística Qui-quadrado de Pearson: 0
cat(paste("Graus de Liberdade (DF):", df_chisq_de, "\n"))
Graus de Liberdade (DF): 0
cat(paste("Pr > ChiSq (p-valor):", round(p_valor_chisq_de, 4), "\n\n"))
Pr > ChiSq (p-valor): 0
<- c(
probabilidades_de seq(from = 0.01, to = 0.10, by = 0.01),
seq(from = 0.15, to = 0.90, by = 0.05),
seq(from = 0.91, to = 0.99, by = 0.01)
)
<- dose.p(modelo_probit_de, p = probabilidades_de)
cl_log_de <- attr(cl_log_de, "SE")
erros_padrao_log_de <- qnorm(0.975)
z_valor_de <- cl_log_de - z_valor_de * erros_padrao_log_de
limite_inf_log_de <- cl_log_de + z_valor_de * erros_padrao_log_de
limite_sup_log_de
<- data.frame(
tabela_escala_log10_de Probabilidade = probabilidades_de,
Log10_conc = cl_log_de,
Limite_Inferior_95 = limite_inf_log_de,
Limite_Superior_95 = limite_sup_log_de
)
<- data.frame(
tabela_escala_conc_de Probabilidade = probabilidades_de,
Concentracao = 10^cl_log_de,
Limite_Inferior_95 = 10^limite_inf_log_de,
Limite_Superior_95 = 10^limite_sup_log_de
)
options(digits = 6)
print(tabela_escala_log10_de)
Probabilidade Log10_conc Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 3.22082e+15 -5.82243e+34 5.82243e+34
p = 0.02: 0.02 3.40917e+15 -6.16293e+34 6.16293e+34
p = 0.03: 0.03 3.52868e+15 -6.37896e+34 6.37896e+34
p = 0.04: 0.04 3.61858e+15 -6.54147e+34 6.54147e+34
p = 0.05: 0.05 3.69170e+15 -6.67366e+34 6.67366e+34
p = 0.06: 0.06 3.75394e+15 -6.78618e+34 6.78618e+34
p = 0.07: 0.07 3.80851e+15 -6.88483e+34 6.88483e+34
p = 0.08: 0.08 3.85738e+15 -6.97317e+34 6.97317e+34
p = 0.09: 0.09 3.90182e+15 -7.05350e+34 7.05350e+34
p = 0.10: 0.10 3.94272e+15 -7.12745e+34 7.12745e+34
p = 0.15: 0.15 4.11209e+15 -7.43362e+34 7.43362e+34
p = 0.20: 0.20 4.24669e+15 -7.67695e+34 7.67695e+34
p = 0.25: 0.25 4.36217e+15 -7.88571e+34 7.88571e+34
p = 0.30: 0.30 4.46588e+15 -8.07318e+34 8.07318e+34
p = 0.35: 0.35 4.56198e+15 -8.24690e+34 8.24690e+34
p = 0.40: 0.40 4.65316e+15 -8.41174e+34 8.41174e+34
p = 0.45: 0.45 4.74139e+15 -8.57123e+34 8.57123e+34
p = 0.50: 0.50 4.82821e+15 -8.72819e+34 8.72819e+34
p = 0.55: 0.55 4.91504e+15 -8.88515e+34 8.88515e+34
p = 0.60: 0.60 5.00326e+15 -9.04464e+34 9.04464e+34
p = 0.65: 0.65 5.09445e+15 -9.20948e+34 9.20948e+34
p = 0.70: 0.70 5.19055e+15 -9.38320e+34 9.38320e+34
p = 0.75: 0.75 5.29425e+15 -9.57067e+34 9.57067e+34
p = 0.80: 0.80 5.40973e+15 -9.77943e+34 9.77943e+34
p = 0.85: 0.85 5.54434e+15 -1.00228e+35 1.00228e+35
p = 0.90: 0.90 5.71370e+15 -1.03289e+35 1.03289e+35
p = 0.91: 0.91 5.75461e+15 -1.04029e+35 1.04029e+35
p = 0.92: 0.92 5.79905e+15 -1.04832e+35 1.04832e+35
p = 0.93: 0.93 5.84791e+15 -1.05715e+35 1.05715e+35
p = 0.94: 0.94 5.90249e+15 -1.06702e+35 1.06702e+35
p = 0.95: 0.95 5.96473e+15 -1.07827e+35 1.07827e+35
p = 0.96: 0.96 6.03785e+15 -1.09149e+35 1.09149e+35
p = 0.97: 0.97 6.12775e+15 -1.10774e+35 1.10774e+35
p = 0.98: 0.98 6.24725e+15 -1.12935e+35 1.12935e+35
p = 0.99: 0.99 6.43561e+15 -1.16340e+35 1.16340e+35
print(tabela_escala_conc_de)
Probabilidade Concentracao Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 Inf 0 Inf
p = 0.02: 0.02 Inf 0 Inf
p = 0.03: 0.03 Inf 0 Inf
p = 0.04: 0.04 Inf 0 Inf
p = 0.05: 0.05 Inf 0 Inf
p = 0.06: 0.06 Inf 0 Inf
p = 0.07: 0.07 Inf 0 Inf
p = 0.08: 0.08 Inf 0 Inf
p = 0.09: 0.09 Inf 0 Inf
p = 0.10: 0.10 Inf 0 Inf
p = 0.15: 0.15 Inf 0 Inf
p = 0.20: 0.20 Inf 0 Inf
p = 0.25: 0.25 Inf 0 Inf
p = 0.30: 0.30 Inf 0 Inf
p = 0.35: 0.35 Inf 0 Inf
p = 0.40: 0.40 Inf 0 Inf
p = 0.45: 0.45 Inf 0 Inf
p = 0.50: 0.50 Inf 0 Inf
p = 0.55: 0.55 Inf 0 Inf
p = 0.60: 0.60 Inf 0 Inf
p = 0.65: 0.65 Inf 0 Inf
p = 0.70: 0.70 Inf 0 Inf
p = 0.75: 0.75 Inf 0 Inf
p = 0.80: 0.80 Inf 0 Inf
p = 0.85: 0.85 Inf 0 Inf
p = 0.90: 0.90 Inf 0 Inf
p = 0.91: 0.91 Inf 0 Inf
p = 0.92: 0.92 Inf 0 Inf
p = 0.93: 0.93 Inf 0 Inf
p = 0.94: 0.94 Inf 0 Inf
p = 0.95: 0.95 Inf 0 Inf
p = 0.96: 0.96 Inf 0 Inf
p = 0.97: 0.97 Inf 0 Inf
p = 0.98: 0.98 Inf 0 Inf
p = 0.99: 0.99 Inf 0 Inf
<- data.frame(lconc_de = seq(0, max(dados_de$lconc_de) + 0.1, length.out = 200))
lconc_grafico_de
<- predict(modelo_probit_de, newdata = lconc_grafico_de, type = "link", se.fit = TRUE)
predicoes_de
<- qnorm(0.975)
z_valor_de $prob_predita_de <- pnorm(predicoes_de$fit)
lconc_grafico_de$prob_superior_de <- pnorm(predicoes_de$fit + z_valor_de * predicoes_de$se.fit)
lconc_grafico_de$prob_inferior_de <- pnorm(predicoes_de$fit - z_valor_de * predicoes_de$se.fit)
lconc_grafico_de
<- as.numeric(dose.p(modelo_probit_de, p = 0.50))
lc50_log_de
<- 10^lc50_log_de
lc50_conc_de
<- -1.5
x_lim_inferior_de <- 1
x_lim_superior_de
<- data.frame(lconc_de = seq(x_lim_inferior_de, x_lim_superior_de, length.out = 200))
lconc_grafico_de
<- predict(modelo_probit_de, newdata = lconc_grafico_de, type = "link", se.fit = TRUE)
predicoes_de <- qnorm(0.975)
z_valor_de $prob_predita_de <- pnorm(predicoes_de$fit)
lconc_grafico_de$prob_superior_de <- pnorm(predicoes_de$fit + z_valor_de * predicoes_de$se.fit)
lconc_grafico_de$prob_inferior_de <- pnorm(predicoes_de$fit - z_valor_de * predicoes_de$se.fit)
lconc_grafico_de
<- ggplot() +
grafico_probit_final_de geom_ribbon(data = lconc_grafico_de,
aes(x = lconc_de, ymin = prob_inferior_de * 100, ymax = prob_superior_de * 100),
fill = "skyblue", alpha = 0.5) +
geom_line(data = lconc_grafico_de,
aes(x = lconc_de, y = prob_predita_de * 100),
color = "blue", size = 1) +
geom_point(data = dados_de,
aes(x = lconc_de, y = mort_prop_de * 100),
color = "red", size = 4) +
geom_hline(yintercept = 50, linetype = "dashed", color = "black", size = 0.7) +
geom_vline(xintercept = lc50_log_de, linetype = "dashed", color = "black", size = 0.7) +
labs(
title = "Curva dose-resposta (Modelo Probit) - Delegate",
x = "Log10 (Concentração)",
y = "Mortalidade observada e prevista (%)"
+
) annotate("text",
x = lc50_log_de,
y = 45,
label = paste("LC50 =", format(lc50_conc_de, digits = 4)),
hjust = -0.1,
vjust = 1,
fontface = "bold",
color = "black") +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
scale_x_continuous(limits = c(x_lim_inferior_de, x_lim_superior_de), breaks = seq(-4, 1, 0.5)) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
print(grafico_probit_final_de)
Joiner
<- data.frame(
dados_jo conc = c(0.00375, 0.0075, 0.01, 0.015, 0.03, 0.07, 0.1, 1, 3, 5, 10),
total = c(60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60),
mortos = c(43, 6, 14, 2, 12, 15, 26, 36, 45, 49, 46)
)
$mort_prop_jo <- dados_jo$mortos / dados_jo$total
dados_jo$lconc_jo <- log10(dados_jo$conc)
dados_jo
print("Dados Iniciais:")
[1] "Dados Iniciais:"
print(dados_jo)
conc total mortos mort_prop_jo lconc_jo
1 0.00375 60 43 0.7166667 -2.425969
2 0.00750 60 6 0.1000000 -2.124939
3 0.01000 60 14 0.2333333 -2.000000
4 0.01500 60 2 0.0333333 -1.823909
5 0.03000 60 12 0.2000000 -1.522879
6 0.07000 60 15 0.2500000 -1.154902
7 0.10000 60 26 0.4333333 -1.000000
8 1.00000 60 36 0.6000000 0.000000
9 3.00000 60 45 0.7500000 0.477121
10 5.00000 60 49 0.8166667 0.698970
11 10.00000 60 46 0.7666667 1.000000
<- lm(mort_prop_jo ~ lconc_jo, data = dados_jo)
modelo_reg_jo
print(anova(modelo_reg_jo))
Analysis of Variance Table
Response: mort_prop_jo
Df Sum Sq Mean Sq F value Pr(>F)
lconc_jo 1 0.4404 0.4404 9.357 0.0136 *
Residuals 9 0.4236 0.0471
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(modelo_reg_jo))
Call:
lm(formula = mort_prop_jo ~ lconc_jo, data = dados_jo)
Residuals:
Min 1Q Median 3Q Max
-0.2546 -0.1380 -0.0015 0.0380 0.5311
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5981 0.0823 7.27 4.7e-05 ***
lconc_jo 0.1701 0.0556 3.06 0.014 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.217 on 9 degrees of freedom
Multiple R-squared: 0.51, Adjusted R-squared: 0.455
F-statistic: 9.36 on 1 and 9 DF, p-value: 0.0136
<- glm(
modelo_probit_jo cbind(mortos, total - mortos) ~ lconc_jo,
data = dados_jo,
family = binomial(link = "probit")
)
print(summary(modelo_probit_jo))
Call:
glm(formula = cbind(mortos, total - mortos) ~ lconc_jo, family = binomial(link = "probit"),
data = dados_jo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.2635 0.0647 4.07 4.7e-05 ***
lconc_jo 0.4501 0.0447 10.06 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 233.58 on 10 degrees of freedom
Residual deviance: 124.18 on 9 degrees of freedom
AIC: 172.9
Number of Fisher Scoring iterations: 4
<- sum(residuals(modelo_probit_jo, type = "pearson")^2)
pearson_chisq_jo <- nrow(dados_jo) - length(coef(modelo_probit_jo))
df_chisq_jo <- pchisq(pearson_chisq_jo, df = df_chisq_jo, lower.tail = FALSE)
p_valor_chisq_jo
cat(paste("Estatística Qui-quadrado de Pearson:", round(pearson_chisq_jo, 4), "\n"))
Estatística Qui-quadrado de Pearson: 138.0422
cat(paste("Graus de Liberdade (DF):", df_chisq_jo, "\n"))
Graus de Liberdade (DF): 9
cat(paste("Pr > ChiSq (p-valor):", round(p_valor_chisq_jo, 4), "\n\n"))
Pr > ChiSq (p-valor): 0
<- c(
probabilidades_jo seq(from = 0.01, to = 0.10, by = 0.01),
seq(from = 0.15, to = 0.90, by = 0.05),
seq(from = 0.91, to = 0.99, by = 0.01)
)
<- dose.p(modelo_probit_jo, p = probabilidades_jo)
cl_log_jo <- attr(cl_log_jo, "SE")
erros_padrao_log_jo <- qnorm(0.975)
z_valor_jo <- cl_log_jo - z_valor_jo * erros_padrao_log_jo
limite_inf_log_jo <- cl_log_jo + z_valor_jo * erros_padrao_log_jo
limite_sup_log_jo
<- data.frame(
tabela_escala_log10_jo Probabilidade = probabilidades_jo,
Log10_conc = cl_log_jo,
Limite_Inferior_95 = limite_inf_log_jo,
Limite_Superior_95 = limite_sup_log_jo
)
<- data.frame(
tabela_escala_conc_jo Probabilidade = probabilidades_jo,
Concentracao = 10^cl_log_jo,
Limite_Inferior_95 = 10^limite_inf_log_jo,
Limite_Superior_95 = 10^limite_sup_log_jo
)
options(digits = 6)
print(tabela_escala_log10_jo)
Probabilidade Log10_conc Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 -5.753995 -6.7297193 -4.7782698
p = 0.02: 0.02 -5.148341 -6.0096998 -4.2869818
p = 0.03: 0.03 -4.764073 -5.5534174 -3.9747285
p = 0.04: 0.04 -4.475003 -5.2105449 -3.7394615
p = 0.05: 0.05 -4.239867 -4.9319371 -3.5477974
p = 0.06: 0.06 -4.039730 -4.6950470 -3.3844122
p = 0.07: 0.07 -3.864248 -4.4875630 -3.2409331
p = 0.08: 0.08 -3.707125 -4.3019909 -3.1122596
p = 0.09: 0.09 -3.564228 -4.1334141 -2.9950426
p = 0.10: 0.10 -3.432692 -3.9784252 -2.8869579
p = 0.15: 0.15 -2.888094 -3.3392408 -2.4369471
p = 0.20: 0.20 -2.455265 -2.8355782 -2.0749519
p = 0.25: 0.25 -2.083936 -2.4086452 -1.7592275
p = 0.30: 0.30 -1.750472 -2.0317978 -1.4691457
p = 0.35: 0.35 -1.441467 -1.6910879 -1.1918463
p = 0.40: 0.40 -1.148252 -1.3785425 -0.9179613
p = 0.45: 0.45 -0.864563 -1.0887619 -0.6403636
p = 0.50: 0.50 -0.585371 -0.8166736 -0.3540693
p = 0.55: 0.55 -0.306180 -0.5564770 -0.0558833
p = 0.60: 0.60 -0.022491 -0.3017753 0.2567933
p = 0.65: 0.65 0.270724 -0.0459488 0.5873971
p = 0.70: 0.70 0.579729 0.2180368 0.9414209
p = 0.75: 0.75 0.913193 0.4985814 1.3278055
p = 0.80: 0.80 1.284522 0.8074736 1.7615706
p = 0.85: 0.85 1.717351 1.1644782 2.2702238
p = 0.90: 0.90 2.261949 1.6107074 2.9131898
p = 0.91: 0.91 2.393485 1.7181363 3.0688346
p = 0.92: 0.92 2.536382 1.8347220 3.2380427
p = 0.93: 0.93 2.693505 1.9627845 3.4242258
p = 0.94: 0.94 2.868987 2.1056680 3.6323053
p = 0.95: 0.95 3.069124 2.2684671 3.8697815
p = 0.96: 0.96 3.304260 2.4595460 4.1489745
p = 0.97: 0.97 3.593330 2.6942157 4.4924443
p = 0.98: 0.98 3.977598 3.0058340 4.9493618
p = 0.99: 0.99 4.583252 3.4963811 5.6701222
print(tabela_escala_conc_jo)
Probabilidade Concentracao Limite_Inferior_95 Limite_Superior_95
p = 0.01: 0.01 1.76200e-06 1.86329e-07 1.66621e-05
p = 0.02: 0.02 7.10656e-06 9.77913e-07 5.16438e-05
p = 0.03: 0.03 1.72158e-05 2.79629e-06 1.05992e-04
p = 0.04: 0.04 3.34963e-05 6.15822e-06 1.82196e-04
p = 0.05: 0.05 5.75616e-05 1.16967e-05 2.83271e-04
p = 0.06: 0.06 9.12579e-05 2.01815e-05 4.12656e-04
p = 0.07: 0.07 1.36695e-04 3.25415e-05 5.74205e-04
p = 0.08: 0.08 1.96279e-04 4.98895e-05 7.72219e-04
p = 0.09: 0.09 2.72754e-04 7.35505e-05 1.01148e-03
p = 0.10: 0.10 3.69240e-04 1.05093e-04 1.29730e-03
p = 0.15: 0.15 1.29392e-03 4.57888e-04 3.65639e-03
p = 0.20: 0.20 3.50538e-03 1.46023e-03 8.41488e-03
p = 0.25: 0.25 8.24259e-03 3.90261e-03 1.74089e-02
p = 0.30: 0.30 1.77635e-02 9.29399e-03 3.39511e-02
p = 0.35: 0.35 3.61854e-02 2.03663e-02 6.42915e-02
p = 0.40: 0.40 7.10801e-02 4.18271e-02 1.20792e-01
p = 0.45: 0.45 1.36596e-01 8.15151e-02 2.28895e-01
p = 0.50: 0.50 2.59794e-01 1.52520e-01 4.42518e-01
p = 0.55: 0.55 4.94106e-01 2.77666e-01 8.79259e-01
p = 0.60: 0.60 9.49531e-01 4.99143e-01 1.80631e+00
p = 0.65: 0.65 1.86519e+00 8.99604e-01 3.86720e+00
p = 0.70: 0.70 3.79952e+00 1.65210e+00 8.73818e+00
p = 0.75: 0.75 8.18829e+00 3.15197e+00 2.12719e+01
p = 0.80: 0.80 1.92541e+01 6.41909e+00 5.77525e+01
p = 0.85: 0.85 5.21616e+01 1.46042e+01 1.86305e+02
p = 0.90: 0.90 1.82788e+02 4.08044e+01 8.18823e+02
p = 0.91: 0.91 2.47449e+02 5.22560e+01 1.17175e+03
p = 0.92: 0.92 3.43861e+02 6.83474e+01 1.72999e+03
p = 0.93: 0.93 4.93748e+02 9.17877e+01 2.65599e+03
p = 0.94: 0.94 7.39583e+02 1.27546e+02 4.28850e+03
p = 0.95: 0.95 1.17253e+03 1.85553e+02 7.40937e+03
p = 0.96: 0.96 2.01493e+03 2.88102e+02 1.40921e+04
p = 0.97: 0.97 3.92040e+03 4.94556e+02 3.10774e+04
p = 0.98: 0.98 9.49725e+03 1.01352e+03 8.89942e+04
p = 0.99: 0.99 3.83047e+04 3.13604e+03 4.67867e+05
<- data.frame(lconc_jo = seq(0, max(dados_jo$lconc_jo) + 0.1, length.out = 200))
lconc_grafico_jo
<- predict(modelo_probit_jo, newdata = lconc_grafico_jo, type = "link", se.fit = TRUE)
predicoes_jo
<- qnorm(0.975)
z_valor_jo $prob_predita_jo <- pnorm(predicoes_jo$fit)
lconc_grafico_jo$prob_superior_jo <- pnorm(predicoes_jo$fit + z_valor_jo * predicoes_jo$se.fit)
lconc_grafico_jo$prob_inferior_jo <- pnorm(predicoes_jo$fit - z_valor_jo * predicoes_jo$se.fit)
lconc_grafico_jo
<- as.numeric(dose.p(modelo_probit_jo, p = 0.50))
lc50_log_jo
<- 10^lc50_log_jo
lc50_conc_jo
<- -1.5
x_lim_inferior_jo <- 1
x_lim_superior_jo
<- data.frame(lconc_jo = seq(x_lim_inferior_jo, x_lim_superior_jo, length.out = 200))
lconc_grafico_jo
<- predict(modelo_probit_jo, newdata = lconc_grafico_jo, type = "link", se.fit = TRUE)
predicoes_jo <- qnorm(0.975)
z_valor_jo $prob_predita_jo <- pnorm(predicoes_jo$fit)
lconc_grafico_jo$prob_superior_jo <- pnorm(predicoes_jo$fit + z_valor_jo * predicoes_jo$se.fit)
lconc_grafico_jo$prob_inferior_jo <- pnorm(predicoes_jo$fit - z_valor_jo * predicoes_jo$se.fit)
lconc_grafico_jo
<- ggplot() +
grafico_probit_final_jo geom_ribbon(data = lconc_grafico_jo,
aes(x = lconc_jo, ymin = prob_inferior_jo * 100, ymax = prob_superior_jo * 100),
fill = "skyblue", alpha = 0.5) +
geom_line(data = lconc_grafico_jo,
aes(x = lconc_jo, y = prob_predita_jo * 100),
color = "blue", size = 1) +
geom_point(data = dados_jo,
aes(x = lconc_jo, y = mort_prop_jo * 100),
color = "red", size = 4) +
geom_hline(yintercept = 50, linetype = "dashed", color = "black", size = 0.7) +
geom_vline(xintercept = lc50_log_jo, linetype = "dashed", color = "black", size = 0.7) +
labs(
title = "Curva dose-resposta (Modelo Probit) - Joiner",
x = "Log10 (Concentração)",
y = "Mortalidade observada e prevista (%)"
+
) annotate("text",
x = lc50_log_jo,
y = 45,
label = paste("LC50 =", format(lc50_conc_jo, digits = 4)),
hjust = -0.1,
vjust = 1,
fontface = "bold",
color = "black") +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
scale_x_continuous(limits = c(x_lim_inferior_jo, x_lim_superior_jo), breaks = seq(-4, 1, 0.5)) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
print(grafico_probit_final_jo)
Interpretação da curva dose-resposta: a CL50 corresponde no gráfico do ponto de interseção entre a linha horizontal de referência a 50% da mortalidade dos insetos, e a curva gerada pelo modelo. Foi possível estimar a CL50 apenas para o joiner (CL50 = 0.2598). Isso se deve ao fato de possuir poucos pontos para a estimação do parâmetro.
Conclusão
Considerando os agrupamentos dos tratamentos, de modo geral, a maioria dos inseticidas testados apresentam alta mortalidade contra as lagartas da traça do tomateiro. Diante dos resultados obtidos, observamos a necessidade de complementação dos dados.