EC50

library(gsheet)
library(tidyverse)

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)) 
  

FGT152 <- dat2 |> 
  filter(Code == "FGT152")

FGT152 |> 
  ggplot(aes(factor(Dose), GC_mean))+
  geom_point()+
  geom_line()

dat2 |> 
   ggplot(aes(factor(Dose), GC_mean))+
  geom_point()+
  geom_line()+
  facet_wrap(~Code)

EC50 com pacote DRM

library(drc)
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)

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