Essa figura representa uma típica curva de dose resposta:
Nesse caso a figura representa a resposta de doses de herbicidas em uma planta.
Vamos importar alguns pacotes até o final dessa aula, porém vamos trabalhar na base tidyverse. O tidyverse é um pacote que contém vários pacotes que auxiliam na manipulação de dados e construção de gráficos. Um exemplo que vamos usar é o pipe (%>%) do pacote magrittr, que faz parte do tidyverse. Assim como o pacote dplyr
para manipulaccção de dados.
library(tidyverse)
options(scipen=TRUE)
Parte dos resultados de Inheritance of mesotrione resistance in an Amaranthus tuberculatus (var. rudis) population from Nebraska, USA de Oliveira et al. (2018).
Vamos usar o pacote RCurl
para baixar o pacote via GitHub.
library(RCurl)
# Baixando via GitHub
<- url("https://raw.githubusercontent.com/maxwelco/workshop-esalq/main/data/dose_resposta.csv")
df_path
# Lendo via read_csv
<- read_csv(df_path)
dose
dose
## # A tibble: 144 x 6
## biotype run rep rate control biomass
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 R I 1 0 0 0.580
## 2 R I 2 0 0 0.54
## 3 R I 3 0 0 0.41
## 4 R I 4 0 0 0.5
## 5 R I 1 13.1 5 0.5
## 6 R I 2 13.1 5 0.64
## 7 R I 3 13.1 10 0.47
## 8 R I 4 13.1 10 0.72
## 9 R I 1 26.2 15 0.52
## 10 R I 2 26.2 10 0.54
## # … with 134 more rows
Colunas | Significado |
---|---|
biotype | biótipo |
run | repetição dos experimentos |
rep | replicação das unidades experimentais |
rate | dose do herbicida |
control | % de controle de planta daninha |
biomass | biomassa seca 21 dias após tratamento |
Sempre visualize seus dados antes de proceder com a análise estatística. Aqui vamos usar o pacote ggplot
do grupo tidyverse
. Você já importou o pacote ggplot
via tidyverse
.
ggplot(dose, aes(x = rate, y = control, color = biotype)) +
geom_point() +
ylim(0, 100)
ggplot(dose, aes(x = factor(rate), y = control, color = biotype)) +
geom_boxplot()
library(drc)
<- drm(control ~ # variável resposta
model # variável explanatória
rate, # fatores
biotype, fct = l4(), # curva usada (existem várias)
data = dose) # onde estão os dados
Veja em detalhes os modelos que podem ser usados aqui DRC CRAN https://cran.r-project.org/web/packages/drc/drc.pdf). Existem várias equações que podem ser usadas em dose-resposta. Qual a melhor? Depende de seus dados!
No DRC CRAN você encontra todas as informações do pacote drc
. Vale a pena explorá-lo antes de começar a usar o pacote.
A função tidy
mostra os parâmetros da equação usada: l4 (4 parâmetros log-logistic).
::tidy(model) broom.mixed
## # A tibble: 8 x 6
## term curve estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 b R -0.826 0.149 -5.55 1.44e- 7
## 2 b S -2.52 0.275 -9.15 8.25e- 16
## 3 c R 1.78 2.59 0.686 4.94e- 1
## 4 c S -0.265 2.68 -0.0987 9.21e- 1
## 5 d R 162. 37.5 4.33 2.88e- 5
## 6 d S 99.4 1.25 79.7 2.66e-116
## 7 e R 1062. 600. 1.77 7.88e- 2
## 8 e S 17.3 0.790 21.9 8.40e- 47
#tidy é uma função do pacote broom.mixed
Parâmetro | Significado |
---|---|
b | inclinação |
c | limite inferior |
d | limite superior |
e | ponto de infecção da curva |
Existe um problema com o resultado de um dos parâmetros!
Visualize o limite superior! É possível isso ocorrer biologicamente?
O pacote drc
possui a função plot
que pode ser usada junto com o modelo usado. Veja abaixo.
plot(model)
Você pode adicionar vários atributos na função plot
como cores, legendas, tamanho etc. Veja como no CRAN do drc
. Mas a função plot
é obseleta no sentido de reportar a figura. Porém, muito usada para rápida visualização de adequação do modelo usado.
A função ED
do pacote drc gera a dose na qual causa 50 ou 90% de controle da planta daninha. Veja que aqui a variável resposta é % de controle de plantas daninhas.
ED(model, c(50, 90), interval = "delta")
##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:R:50 1062.47359 599.87212 -123.80998 2248.75717
## e:R:90 15203.89625 15580.07840 -15606.65569 46014.44819
## e:S:50 17.33522 0.79036 15.77224 18.89820
## e:S:90 41.48941 4.25617 33.07257 49.90625
ED(model, c(50, 90), interval = "delta",
type = "absolute")
##
## Estimated effective doses
##
## Estimate Std. Error Lower Upper
## e:R:50 382.32000 152.96053 79.83127 684.80874
## e:R:90 1355.52014 820.73737 -267.53788 2978.57817
## e:S:50 17.45215 0.79482 15.88036 19.02395
## e:S:90 42.53316 4.46712 33.69916 51.36715
Observe que os resultados mudam (e muito!) quando usamos type = "absolute"
. Existem duas tipos aqui type = "absolute"
ou type = "relative"
. Tipo relativo é o padrão usado, então quando usamos a função ED
sem o atributo type, o resultado será relativo. Relativo ao limite superior da curva. Por outro lado, quando especificamos type = "absolute"
, o resultado será o valor absoluto especificado, como 50 ou 90% de controle (valores do eixo-y).
Se quiser ver todos os atributos para a função ED
, escreva no console ?ED(). Vai aparecer em Help todas as informações. Isso é válido para qualquer função no R.
É muito importante saber diferenciar os tipos de EDs. Cuidado!
função EDcomp
faz comparações entre os EDs. Por exemplo, a código abaixo compara a dose efetiva para causar 50% de injúrias (controle) do herbicida nos biótipos de planta daninhas.
EDcomp(model, c(50, 50), type = "absolute")
##
## Estimated ratios of effect doses
##
## Estimate Std. Error t-value p-value
## R/S:50/50 21.906751 8.821166 2.370067 0.019189
Se você se lembrar, estamos com dois biótipos de plantas daninhas: R e S. Nesse caso uma população estudada é suspectível (S) ao herbicida e a outra resistênte (R). Os resultados nos dão confiança que o efeito do herbicida nos dois biótipos são diferentes. Por exemplo, o fator de resistêcia é 22 (p-value = 0.019189)
Você também pode comparar os quatros parâmetros entre os biótipos. Veja como no CRAN do drc.
# gerar uma novo dataset
<- expand.grid(rate=exp(seq(log(0.5), log(1680), length=1680)),
nd biotype = c("S", "R"))
# predict the model
<- predict(model, newdata=nd, interval="confidence")
pm
# predict intervalo de confiança
$p <- pm[,1]
nd$pmin <- pm[,2]
nd$pmax <- pm[,3]
nd
# ajustando o 0
$rate0 <- dose$rate
dose$rate0[dose$rate==0] <- 0.5 dose
ggplot(dose, aes(x = rate, y = control,
color = biotype,
fill = biotype)) +
geom_line(data=nd, aes(x=rate, y=p,
linetype = biotype),
size=1.3) +
geom_ribbon(data=nd, aes(x=rate, y=p, ymin=pmin, ymax=pmax,
color = NULL),
alpha=0.2) +
theme_bw() +
coord_trans(x="log") +
scale_x_continuous(breaks=c(0, 1, 5, 12, 26, 53, 105,
210, 420, 840, 1680))
Keshtkar et al. (2021) Perspective: common errors in dose–response analysis and how to avoid them
Knezevic et al. (2007) Utilizing R Software Package for Dose-Response Studies: The Concept and Data Analysis
Ritz et al. (2015) Dose-Response Analysis Using R