1 Dose resposta

Essa figura representa uma típica curva de dose resposta:

Nesse caso a figura representa a resposta de doses de herbicidas em uma planta.

2 Importando os pacotes

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)

3 Dados

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
df_path <- url("https://raw.githubusercontent.com/maxwelco/workshop-esalq/main/data/dose_resposta.csv")

# Lendo via read_csv
dose <- read_csv(df_path)

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

4 Visualização de dados

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.

4.1 Gráfico de pontos

ggplot(dose, aes(x = rate, y = control, color = biotype)) + 
  geom_point() +
  ylim(0, 100)

4.2 Gráfico de dispersão

ggplot(dose, aes(x = factor(rate), y = control, color = biotype)) + 
  geom_boxplot()

5 Modelo

library(drc)
model <- drm(control ~ # variável resposta
               rate, # variável explanatória
             biotype, # fatores
             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.

5.1 Parâmetros

A função tidy mostra os parâmetros da equação usada: l4 (4 parâmetros log-logistic).

broom.mixed::tidy(model)
## # 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?

5.2 Visualizando o modelo

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.

5.3 Gerando valor a sua análise

5.3.1 EDs

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!

5.3.2 Comparando os EDs

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.

6 Gráfico

# gerar uma novo dataset
nd <- expand.grid(rate=exp(seq(log(0.5), log(1680), length=1680)),
                       biotype = c("S", "R"))

# predict the model
pm <- predict(model, newdata=nd, interval="confidence")

# predict intervalo de confiança
nd$p <- pm[,1] 
nd$pmin <- pm[,2] 
nd$pmax <- pm[,3] 

# ajustando o 0
dose$rate0 <- dose$rate
dose$rate0[dose$rate==0] <- 0.5
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))

7 Leitura complementar