#Analysis of Dose-Response Curves
rm(list = ls())
#install.packages(“drc”, dependencies=TRUE)
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
library(ggplot2)
#IC50 is the “half-maximal inhibitory concentration” for an entity (typically a drug) against a biological process or function (eg, enzyme activity, cell number, metabolic activity). 
#If you have an inhibitor of a given enzyme and an in vitro assay for that enzyme, the IC50 is the concentration of inhibitor at which that enzyme runs at 50% of its normal activity in that assay.
#EC50 is, comparably, the “half-maximal effective concentration”. The effective is in there because not all entities will completely wipe out the biological process that you are assaying. 
#Fitting data to the EC50 equation
## Displaying the data set
head(ryegrass)
##      rootl conc
## 1 7.580000    0
## 2 8.000000    0
## 3 8.328571    0
## 4 7.250000    0
## 5 7.375000    0
## 6 7.962500    0
table(ryegrass$conc)
## 
##    0 0.94 1.88 3.75  7.5   15   30 
##    6    3    3    3    3    3    3
dim(ryegrass)
## [1] 24  2
## Fitting a four-parameter log-logistic model
## with user-defined parameter names
ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, 
                   fct = LL.4(names = c("Slope", "Lower Limit", "Upper Limit", "ED50")))

LL.4’ provide the four-parameter log-logistic function

#b = the Hill coefficient;c = min_value;d = max_value;e = EC50.
summary(ryegrass.m1)
## 
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
## 
## Parameter estimates:
## 
##                         Estimate Std. Error t-value   p-value    
## Slope:(Intercept)        2.98222    0.46506  6.4125 2.960e-06 ***
## Lower Limit:(Intercept)  0.48141    0.21219  2.2688   0.03451 *  
## Upper Limit:(Intercept)  7.79296    0.18857 41.3272 < 2.2e-16 ***
## ED50:(Intercept)         3.05795    0.18573 16.4644 4.268e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  0.5196256 (20 degrees of freedom)
confint(ryegrass.m1, level = 0.90)
##                              5 %      95 %
## Slope:(Intercept)       2.180118 3.7843205
## Lower Limit:(Intercept) 0.115441 0.8473854
## Upper Limit:(Intercept) 7.467733 8.1181836
## ED50:(Intercept)        2.737621 3.3782892
#We can convert the inferred values to an IC50.
coefs <- setNames(
  ryegrass.m1$coefficients,
  c("hill", "min_value", "max_value", "ec_50")
)
coefs 
##      hill min_value max_value     ec_50 
## 2.9822191 0.4814132 7.7929583 3.0579550
ic_50 <- with(
  as.list(coefs),
  exp(
    log(ec_50) + (1 / hill) * log(max_value / (max_value - 2 * min_value))
  )
)
ic_50 
## [1] 3.196215
##############################plot
#plot_1
plot(ryegrass.m1, broken=F, xlab="Dose (mM)", ylab="Root length (cm)", lwd=1.2, 
     cex=1.2, cex.axis=1.2, cex.lab=1.2)

#plot_2
p1 <- ggplot(data = ryegrass,
       aes(x = conc, y = rootl)) + 
  geom_point(col = "yellowgreen") + 
  geom_smooth(method = drm,col = "skyblue", method.args = list(fct = L.4()), se = F) +
  scale_x_log10() +
  labs(title= "", x = "Dose (mM)",  y = "Root length (cm)") + 
  theme(legend.position = "",
        panel.background = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=0),
        panel.grid = element_blank(),
        axis.text.x   = element_text(size= 12, color = "darkblue",family = "sans",hjust = 0.5),
        axis.text.y   = element_text(size= 12, color = "darkblue",family = "sans",vjust = 0.5,hjust = 0.5),
        axis.title  = element_text(size=12, color = "darkred",family = "sans"),
        axis.ticks =  element_line(size= 0.5),
        axis.ticks.length = unit(3, "pt"))
p1
## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 6 rows containing non-finite values (stat_smooth).

#output p1
ggsave(paste0(Sys.Date(),"-plot2.tiff"), 
       plot = p1, device = NULL, 
       scale = 1, width = 10, height = 10, units ="cm",dpi = 300, 
       limitsize = TRUE,compression = "lzw")
## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 6 rows containing non-finite values (stat_smooth).
#plot_2 mean + sd
head(ryegrass)
##      rootl conc
## 1 7.580000    0
## 2 8.000000    0
## 3 8.328571    0
## 4 7.250000    0
## 5 7.375000    0
## 6 7.962500    0
ryegrass_1 <- aggregate(ryegrass[1],list(conc = ryegrass$conc),FUN = mean)
ryegrass_1$sd <- aggregate(ryegrass[1],list(conc = ryegrass$conc),FUN = sd)$rootl
ryegrass_1
##    conc     rootl        sd
## 1  0.00 7.7493452 0.4151924
## 2  0.94 7.6732804 0.7236913
## 3  1.88 6.4145503 0.4755951
## 4  3.75 3.0146825 1.1595582
## 5  7.50 1.0339286 0.1663975
## 6 15.00 0.6791667 0.1501735
## 7 30.00 0.3033333 0.1193035
sd(ryegrass$rootl[ryegrass$conc == 0])
## [1] 0.4151924
p2 <- ggplot(data = ryegrass_1,
             aes(x = conc, y = rootl)) + 
  geom_point(col = "red", size=3, shape=21, fill="green",position = "identity") + 
  geom_errorbar(aes(ymin=rootl-sd, ymax=rootl+sd), width=.05,position = "identity")+
  geom_smooth(method = drm,col = "skyblue", method.args = list(fct = L.4()), se = F) +
  labs(title= "", x = "Dose (mM)",  y = "Root length (cm)") + 
  theme(legend.position = "",
        panel.background = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=0),
        panel.grid = element_blank(),
        axis.text.x   = element_text(size= 12, color = "darkblue",family = "sans",hjust = 0.5),
        axis.text.y   = element_text(size= 12, color = "darkblue",family = "sans",vjust = 0.5,hjust = 0.5),
        axis.title  = element_text(size=12, color = "darkred",family = "sans"),
        axis.ticks =  element_line(size= 0.5),
        axis.ticks.length = unit(3, "pt")) +
        scale_x_log10() 
p2
## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 1 rows containing non-finite values (stat_smooth).

#output p3
ggsave(paste0(Sys.Date(),"-plot2.tiff"), 
       plot = p2, device = NULL, 
       scale = 1, width = 10, height = 10, units ="cm",dpi = 300, 
       limitsize = TRUE,compression = "lzw")
## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous x-axis
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
#
#We can predict values using predict on the newly defined curved_fit variable.
additional_conc <- data.frame(
  # This gives a sequence of 10 evenly spaced values between 0 and 10
  conc = seq(0, 10, 1)
)
predicted_data <- predict(ryegrass.m1, newdata = additional_conc)
predicted_data
##  [1] 7.7929583 7.5411183 6.1851478 4.2414604 2.7468200 1.8523337 1.3452816
##  [8] 1.0517508 0.8744747 0.7625279 0.6888809
data.frame(additional_conc,predicted_data)
##    conc predicted_data
## 1     0      7.7929583
## 2     1      7.5411183
## 3     2      6.1851478
## 4     3      4.2414604
## 5     4      2.7468200
## 6     5      1.8523337
## 7     6      1.3452816
## 8     7      1.0517508
## 9     8      0.8744747
## 10    9      0.7625279
## 11   10      0.6888809
###############ref and acknowledgement
#https://rpubs.com/russH/378543