#To determine mycelial EC50 of isolate 2022-065 (being used for the NCSFR Botrytis screenhouse trial 2024)
#because based on the 1 ppm fenhexamid dose that is discriminatory for conidial germination,
#I saw mycelial growth in 2022-065 on 1 ppm...even though it should be sensitive.
#IN SUMMARY: RELATIVE EC50 OF 2022-065 = 1.106921 ppm fenhexamid
# ABSOLUTE EC50 OF 2022-065 = 1.065093 ppm fenhexamid
#Noel et al. (2018) recommend package 'drc'
install.packages("drc")
##
## The downloaded binary packages are in
## /var/folders/5g/_lqy13l94p1gv6wm1sf90hdr0000gn/T//RtmpjqTkY3/downloaded_packages
library(drc)
## Loading required package: MASS
##
## '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
EC50 <- read.csv("2022-065 mycelial growth EC50.csv")
#Make sure the factor variables are being treated as such
EC50$Isolate <- factor(EC50$Isolate)
EC50$Concentration <- factor(EC50$conc)
#Want to only look at 2022-065 isolate, not 2022-860 isolate, which was my resistant check
#Want to exclude rows where conc = "control", because "Eth" is our true control -
#"control" was completely unamended medium.
EC50_065 <- subset(EC50, Isolate == "2022-065", conc != "Control")
#When trying models out with non-transformed "conc", I kept getting fail to converge error messages.
#I need to log-transform to conc + 1 - this is how data is most often displayed in dose response curves.
#Log-logistic is the most commonly used model for dose response curves.
#However, my dataset is small (<15-20 data points per concentration level) so it may not
#be appropriate here per Ritz et al. (2015) outlining this package.
#Trying initial four-parameter log-logistic model
mycelial.m1 <- drm(response ~ log_conc, data = EC50_065, fct = LL.4())
#No error message about failing to converge.
#Selecting appropriate model
# Use mselect() to select the best-fitting model
#Here, I try some of the most commonly suggested in the example for mselect()
mselect(mycelial.m1, list(LL.2(), LL.3(), LL.4(), W1.2(), W1.3(), baro5()))
## Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt, :
## non-finite finite-difference value [1]
## logLik IC Lack of fit Res var
## W1.3 -126.3688 260.7376 9.995635e-01 25.88834
## LL.3 -127.0282 262.0565 7.626356e-01 26.71418
## LL.4 -126.5053 263.0106 8.834683e-01 26.74290
## LL.4 -126.5053 263.0106 8.834683e-01 26.74290
## LL.2 -230.0026 466.0052 4.953004e-38 3510.16111
## W1.2 -230.0026 466.0052 4.953004e-38 3510.16111
## baro5 NA NA NA NA
#For Akaike's information criterion (AIC) and the residual standard error:
#the smaller the better and for lack-of-fit test (against a one-way ANOVA model):
#the larger (the p-value) the better. Note that the residual standard error
#is only available for continuous dose-response data.
#Log likelihood values cannot be used for comparison unless the models are nested.
#Based off of AIC, W1.3 model is best which is a three-parameter Weibull model.
#(I also tried W1.4 and it wasn't better)
#Let's do it
mycelial.model <- drm(response ~ log_conc, data = EC50_065, fct = W1.3())
summary(mycelial.model)
##
## Model fitted: Weibull (type 1) with lower limit at 0 (3 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) 0.7688206 0.0577324 13.317 4.378e-16 ***
## d:(Intercept) 99.9986341 1.4687425 68.085 < 2.2e-16 ***
## e:(Intercept) 0.1015825 0.0077235 13.152 6.545e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 5.088058 (39 degrees of freedom)
#Model fitted: Weibull (type 1) with lower limit at 0 (3 parms)
#note: Weibull is an asymmetric model type, while something like log-logistic is symmetric.
#The three parameters in this model are upper asymptote exp(c)
#(not given here), lower asymptote- exp(d), slope - exp(b),
#and exp(e) - the inflection point, which is the relative EC50!
exp(0.1015825)
## [1] 1.106921
# # # #1.106921 is the relative EC50 of 2022-065. # # # #
#We want both absolute and relative EC50 values.
#Absolute EC50 is where 50% of the maximum response occurs;
#this is in line with how FRAC defines EC50
#and so it is what I will use as my phenotyping dose on media.
#As mentioned above, relative EC50 is the inflection point on a dose-response curve.
#From Noel et al. 2018: "The absolute EC50 was estimated by fitting percent relative
#growth against the log concentration using a four-, three-, and two-parameter
#log-logistic model but specifying type = “absolute” within the “ED” function of
#“drc” (Ritz and Streibig 2005)."
#default in the ED() argument is relative
ED(mycelial.model,respLev = c(50), type = "absolute")
##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 0.0630624 0.0051015
#Now take reverse log of the estimate for EC50
exp(0.0630624)
## [1] 1.065093
# # # # #1.065093 ppm is the absolute EC50 of 2022-065. # # # #
#Plotting the model with absolute and relative EC50 as labelled data points
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(magrittr)
#Use my drm to predict y values for the relative EC50 (we know absolute y is 50)
x_value <- c(0.1015825)
predicted_y <- predict(mycelial.model, newdata = data.frame(LogDose = x_value))
#Make plot based on the drm model, "mycelial.model"
plot(mycelial.model,
xlab = expression(bold("Log[Fenhexamid]")),
ylab = expression(bold("Percent growth relative to control")),
lwd = 2, type = "n", pch = 19, xaxt = "n", yaxt = "n", legend = F, cex.axis = 1.2)
points(0.0630624, 50, pch = 19, cex = 1.5)
points(0.10158215, predicted_y, pch = 15, cex = 1.5)
text(0.0630624, 50, labels=expression(bold("Absolute EC50")), cex= 1, adj = c(-0.1, 1))
text(0.10158215, predicted_y, labels=expression(bold("Relative EC50")), cex= 1, adj = c(-0.1, -0.9))
title("Dose Response Curve for 2022-065")
