Aula17

Author

ARLAM 

library(gsheet)
library(tidyverse)
Warning: package 'ggplot2' was built under R version 4.2.3
Warning: package 'tibble' was built under R version 4.2.3
Warning: package 'dplyr' was built under R version 4.2.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
dat <- gsheet2tbl("https://docs.google.com/spreadsheets/d/15pCj0zljvd-TGECe67OMt6sa21xO8BqUgv4d-kU8qi8/edit#gid=0")

options(scipen=999)
dat2 <- dat |> 
  select(-Isolate, -Population) |> 
  group_by(Code, Year, Dose) |> 
  summarise(GC_mean = mean(GC)) 
`summarise()` has grouped output by 'Code', 'Year'. You can override using the
`.groups` argument.
FGT152 <- dat2 |> 
  filter(Code == "FGT152")

FGT152 |> 
  ggplot(aes(factor(Dose), GC_mean))+
  geom_point()+
  geom_line()
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

dat2 |> 
   ggplot(aes(factor(Dose), GC_mean))+
  geom_point()+
  geom_line()+
  facet_wrap(~Code)
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

EC50 com pacote DRM

library(drc)
Warning: package 'drc' was built under R version 4.2.3
Carregando pacotes exigidos: MASS
Warning: package 'MASS' was built under R version 4.2.3

Attaching package: 'MASS'
The following object is masked from 'package:dplyr':

    select

'drc' has been loaded.
Please cite R and 'drc' if used for a publication,
for references type 'citation()' and 'citation('drc')'.

Attaching package: 'drc'
The following objects are masked from 'package:stats':

    gaussian, getInitial
drc1 <- drm(GC_mean ~ Dose, data = FGT152,
    fct = LL.3())
AIC(drc1)
[1] 33.60846
summary(drc1)

Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (3 parms)

Parameter estimates:

               Estimate Std. Error t-value     p-value    
b:(Intercept)  0.401905   0.053427  7.5225    0.001672 ** 
d:(Intercept) 47.540342   1.459890 32.5643 0.000005302 ***
e:(Intercept)  7.220130   2.340119  3.0854    0.036739 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 1.993805 (4 degrees of freedom)
plot(drc1)

ED(drc1, 50, interval = "delta")

Estimated effective doses

       Estimate Std. Error    Lower    Upper
e:1:50  7.22013    2.34012  0.72292 13.71734
library(ec50estimator)
Warning: package 'ec50estimator' was built under R version 4.2.3
df_ec50 <- estimate_EC50(GC_mean ~ Dose,
                         data = dat2,
                         isolate_col = "Code",
                         interval = "delta",
                         fct = drc::LL.3())


## Função criada pelo Cha

df_ED50 <- function(formula, isolate_col, fct = LL.3(), ed_perc = 50, interval = "delta", data) {
  # get the unique values of the isolate column
  isolates <- unique(data %>% pull(isolate_col))

  # apply drm() and ED() to each isolate
  quiet_map <- quietly(map_df)
  results <- quiet_map(
    isolates, 
    ~{
      # subset the data for the current isolate
      isolate_data <- data %>% filter(!!sym(isolate_col) == .x)
      
      # apply drm()
      drm_result <- drm(formula, data = isolate_data, fct = fct)
      
      # apply ED() and extract the EC50 value and its confidence interval
      ed_result <- as.data.frame(ED(drm_result, ed_perc, interval = interval))
      
      # return as a data frame
      tibble(
        isolate = .x,
        EC50 = ed_result["e:1:50", "Estimate"],
        lower = ed_result["e:1:50", "Lower"],
        upper = ed_result["e:1:50", "Upper"]
      )
    }
  )$result

  return(results)
}

results <- df_ED50(GC_mean ~ Dose, 
                   isolate_col = "Code", 
                   data = dat2)


# print the results
results
# A tibble: 12 × 4
   isolate    EC50    lower  upper
   <chr>     <dbl>    <dbl>  <dbl>
 1 FGT152   7.22     0.723  13.7  
 2 FGT169   0.577    0.0549  1.10 
 3 FGT170   2.28     1.41    3.15 
 4 FGT186   3.29    -1.15    7.72 
 5 FGT187   0.0985  -0.0226  0.220
 6 FGT189   9.30     2.29   16.3  
 7 FGT201   4.93     0.172   9.68 
 8 FGT202  13.8    -17.4    45.0  
 9 FGT215   0.505    0.213   0.797
10 FGT229   0.568   -0.241   1.38 
11 FGT231   1.83     0.153   3.50 
12 FGT232   0.340    0.105   0.575