# Primary Packages #
library(tcplfit2)
library(tcpl)
# Data Formatting Packages #
library(data.table)
library(DT)
library(htmlTable)
library(dplyr)
library(stringr)
# Plotting Packages #
library(ggplot2)
library(gridExtra)

# tested concentrations
  conc <- list(.03,.1,.3,1,3,10,30,100)
# observed responses at respective concentrations
  resp <- list(0,.2,.1,.4,.7,.9,.6, 1.2)
# row object with relevant parameters
  row = list(conc = conc, resp = resp, bmed = 0, cutoff = 1, onesd = .5,name="some chemical")
# execute concentration-response modeling through potency estimation
  res <- concRespCore(row,
                      fitmodels = c("cnst", "hill", "gnls",
                                    "poly1", "poly2", "pow", "exp2", "exp3",
                                        "exp4", "exp5"),
                      conthits = T)

# plot the winning curve from example 1, add a title
concRespPlot2(res, log_conc = TRUE) + ggtitle("Example 1: Chemical A")

# read in the data
# Loading in the level 3 example data set from invitrodb
  data("mc3")
  head(mc3)
#>          dtxsid   casrn        name         spid      logc       resp
#> 1 DTXSID7020182 80-05-7 Bisphenol A Tox21_400088  1.477121 36.3584083
#> 2 DTXSID7020182 80-05-7 Bisphenol A Tox21_400088 -2.000000  0.3244336
#> 3 DTXSID7020182 80-05-7 Bisphenol A Tox21_400088  0.000000  9.0799288
#> 4 DTXSID7020182 80-05-7 Bisphenol A Tox21_400088 -3.000000  1.4104484
#> 5 DTXSID7020182 80-05-7 Bisphenol A Tox21_400088 -1.000000  0.5884602
#> 6 DTXSID7020182 80-05-7 Bisphenol A Tox21_400088 -2.000000 -0.5267254
#>                         assay
#> 1 TOX21_ERa_BLA_Agonist_ratio
#> 2 TOX21_ERa_BLA_Agonist_ratio
#> 3 TOX21_ERa_BLA_Agonist_ratio
#> 4 TOX21_ERa_BLA_Agonist_ratio
#> 5 TOX21_ERa_BLA_Agonist_ratio
#> 6 TOX21_ERa_BLA_Agonist_ratio

# determine the background variation
# chosen as logc <= -2 in this example but will be assay/application specific
  temp <- mc3[mc3$logc<= -2,"resp"]
  bmad <- mad(temp)
  onesd <- sd(temp)
  cutoff <- 3*bmad

# select six chemical samples. Note that there may be more than one sample processed for a given chemical
  spid.list <- unique(mc3$spid)
  spid.list <- spid.list[1:6]
  
# create empty objects to store results and plots
  model_fits <- NULL
  result_table <- NULL
  plt_lst <- NULL

# loop over the samples to perform concentration-response modeling & hitcalling
  for(spid in spid.list) {
    # select the data for just this sample
    temp <- mc3[is.element(mc3$spid,spid),]

    # The data file stores concentrations in log10 units, so back-transform
    conc <- 10**temp$logc
    # Save the response values
    resp <- temp$resp

    # pull out all of the chemical identifiers and the assay name
    dtxsid <- temp[1,"dtxsid"]
    casrn <- temp[1,"casrn"]
    name <- temp[1,"name"]
    assay <- temp[1,"assay"]
    
    # Execute curve fitting
    # Input concentrations, responses, cutoff, a list of models to fit, and other model fitting requirements
    # force.fit is set to true so that all models will be fit regardless of cutoff
    # bidirectional = FALSE indicates only fit models in the positive direction.
    # if using bidirectional = TRUE the coff only needs to be specified in the positive direction.
    model_fits[[spid]] <- tcplfit2_core(conc, resp, cutoff, force.fit = TRUE, 
                                        fitmodels = c("cnst", "hill", "gnls", 
                                                      "poly1", "poly2", "pow", 
                                                      "exp2","exp3", "exp4", "exp5"),
                                        bidirectional = FALSE)
    # Get a plot of all curve fits
    plt_lst[[spid]] <- plot_allcurves(model_fits[[spid]], 
                                      conc = conc, resp = resp, log_conc = TRUE)
    
    # Pass the output from 'tcplfit2_core' to 'tcplhit2_core' along with
    # cutoff, onesd, and any identifiers
    out <- tcplhit2_core(model_fits[[spid]], conc, resp, bmed = 0,
                         cutoff = cutoff, onesd = onesd, 
                         identifiers = c(dtxsid = dtxsid, casrn = casrn, 
                                         name = name, assay = assay))
    # store all results in one table
    result_table <- rbind(result_table,out)
  }

# shows the structure of the output object from tcplfit2_core (only top level)
str(model_fits[[1]],max.lev = 1)
#> List of 12
#>  $ cnst      :List of 5
#>  $ hill      :List of 17
#>  $ gnls      :List of 22
#>  $ poly1     :List of 13
#>  $ poly2     :List of 17
#>  $ pow       :List of 15
#>  $ exp2      :List of 15
#>  $ exp3      :List of 17
#>  $ exp4      :List of 15
#>  $ exp5      :List of 17
#>  $ modelnames: chr [1:10] "cnst" "hill" "gnls" "poly1" ...
#>  $ errfun    : chr "dt4"

str(model_fits[[1]][["hill"]])
#> List of 17
#>  $ success: int 1
#>  $ aic    : num 3931
#>  $ cov    : int 1
#>  $ rme    : num 5.4
#>  $ modl   : num [1:675] 2.84e+01 1.96e-03 1.22e+01 1.87e-05 2.04e-01 ...
#>  $ tp     : num 28.5
#>  $ ga     : num 1.15
#>  $ p      : num 2.02
#>  $ er     : num 1.07
#>  $ tp_sd  : num 0.386
#>  $ ga_sd  : num 0.0591
#>  $ p_sd   : num 0.133
#>  $ er_sd  : num 0.0425
#>  $ pars   : chr [1:4] "tp" "ga" "p" "er"
#>  $ sds    : chr [1:4] "tp_sd" "ga_sd" "p_sd" "er_sd"
#>  $ top    : num 28.5
#>  $ ac50   : num 1.15

grid.arrange(grobs=plt_lst,ncol=2)

htmlTable::htmlTable(head(result_table),
        align = 'l',
        align.header = 'l',
        rnames = FALSE  ,
        css.cell =  ' padding-bottom: 5px;  vertical-align:top; padding-right: 10px;min-width: 5em ')

# plot the first row
concRespPlot2(result_table[1,],log_conc = TRUE) + 
  ggtitle(paste(result_table[1,"dtxsid"], result_table[1,"name"]))

# Loading in the Level 0 example data set from invitrodb
data("mc0")
data.table::setDTthreads(2)
dat <- mc0

# Order by the following columns
setkeyv(dat, c('acid', 'srcf', 'apid', 'coli', 'rowi', 'spid', 'conc'))

# Define a temporary replicate ID (rpid) column for test compound wells
# rpid consists of the sample ID, well type (wllt), source file, assay plate ID, and 
# concentration.
nconc <- dat[wllt == "t" , ## denotes test well as the well type (wllt)
             list(n = lu(conc)), #total number of unique concentrations
             by = list(acid, apid, spid)][ , list(nconc = min(n)), by = acid]
dat[wllt == "t" & acid %in% nconc[nconc > 1, acid],
    rpid := paste(acid, spid, wllt, srcf, apid, "rep1", conc, sep = "_")]
dat[wllt == "t" & acid %in% nconc[nconc == 1, acid],
    rpid := paste(acid, spid, wllt, srcf, "rep1", conc, sep = "_")]

# Define rpid column for non-test compound wells
dat[wllt != "t",
    rpid := paste(acid, spid, wllt, srcf, apid, "rep1", conc, sep = "_")]

# set the replicate index (repi) based on rowid 
# increment repi every time a replicate ID is duplicated
dat[, dat_rpid := rowid(rpid)]
dat[, rpid := sub("_rep[0-9]+.*", "",rpid, useBytes = TRUE)]
dat[, rpid := paste0(rpid,"_rep",dat_rpid)]

# For each replicate, define concentration index
# by ranking the unique concentrations
indexfunc <- function(x) as.integer(rank(unique(x))[match(x, unique(x))])
# the := operator is a data.table function to add/update rows
dat[ , cndx := indexfunc(conc), by = list(rpid)]

# If no adjustments are required for the data, the corrected value (cval) should be set as original rval
dat[,cval := rval]

# Poor well quality (wllq) wells should be removed
dat <- dat[!wllq == 0,]

##Fitting generally cannot occur if response values are NA therefore values need to be removed
dat <- dat[!is.na(cval),]

# calculate bval of the median of all the wells that have a type of n
dat[, bval := median(cval[wllt == "n"]), by = list(apid)]
# calculate pval based on the wells that have type of m or o excluding any NA wells
dat[, pval := median(cval[wllt %in% c("m","o")], na.rm = TRUE), by = list(apid, wllt, conc)]
# take pval as the minimum per assay plate (apid)
dat[, pval := min(pval, na.rm = TRUE), by = list(apid)]

# Calculate normalized responses
dat[, resp := ((cval - bval)/(pval - bval) * 100)]

bmad <- mad(dat[cndx %in% c(1, 2) & wllt == "t", resp])
onesd <- sd(dat[cndx %in% c(1, 2) & wllt == "t", resp])

#do tcplfit2 fitting
myfun <- function(y) {
  res <- tcplfit2::tcplfit2_core(y$conc,
                          y$resp,
                          cutoff = 3*bmad,
                          bidirectional = TRUE,
                          verbose = FALSE,
                          force.fit = TRUE,
                          fitmodels = c("cnst", "hill", "gnls", "poly1",
                                        "poly2", "pow", "exp2", "exp3",
                                        "exp4", "exp5")
                          )
  list(list(res)) #use list twice because data.table uses list(.) to look for values to assign to columns
}

# create a subset that contains 6 samples and run curve fitting
subdat <- dat[spid %in% unique(spid)[10:15],]
subdat[wllt == 't',params:= myfun(.SD), by = .(spid)]

myfun2 <- function(y) {
  res <- tcplfit2::tcplhit2_core(params = y$params[[1]],
                                 conc = y$conc,
                                 resp = y$resp,
                                 cutoff = 3*bmad,
                                 onesd = onesd
                                 )
  list(list(res))
}

# continue with hitcalling
res <- subdat[wllt == 't', myfun2(.SD), by = .(spid)]

# pivot wider
res_wide <- rbindlist(Map(cbind, spid = res$spid, res$V1))

# allocate a place-holder object
  plt_list <- NULL
# plot results using `concRespPlot`
  for(i in 1:nrow(res_wide)){
    plt_list[[i]] <- concRespPlot2(res_wide[i,])
  }
# compile and display winning model plots for concentration-response series
  grid.arrange(grobs=plt_list,ncol=2)

# We'll use data from mc3 in this section
data("mc3")

# determine the background variation
# background is defined per the assay.  In this case we use logc <= -2
# However, background should be defined in a way that makes sense for your application
temp <- mc3[mc3$logc<= -2,"resp"]
bmad <- mad(temp)
onesd <- sd(temp)
cutoff <- 3*bmad

# load example data
spid <- unique(mc3$spid)[94]
ex_df <- mc3[is.element(mc3$spid,spid),]

# The data file has stored concentration in log10 form, fix it 
conc <- 10^ex_df$logc # back-transforming concentrations on log10 scale
resp <- ex_df$resp

# modify the data for demonstration purposes 
conc2 <- conc[conc>0.41]
resp2 <- resp[which(conc>0.41)]

# pull out all of the chemical identifiers and the name of the assay
dtxsid <- ex_df[1,"dtxsid"]
casrn <- ex_df[1,"casrn"]
name <- ex_df[1,"name"]
assay <- ex_df[1,"assay"]

# create the row object
row_low <- list(conc = conc2, resp = resp2, bmed = 0, cutoff = cutoff, onesd = onesd,
            assay=assay, dtxsid=dtxsid,casrn=casrn,name=name)

# run the concentration-response modeling for a single sample
res_low <- concRespCore(row_low,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", 
                                          "pow", "exp2", "exp3", "exp4", "exp5"), 
                        bidirectional=F)

concRespPlot2(res_low, log_conc = T) + 
  geom_rect(aes(xmin = log10(res_low[1, "bmdl"]),
                xmax = log10(res_low[1, "bmdu"]),ymin = 0,ymax = 30),
            alpha = 0.05,fill = "skyblue") + 
  geom_segment(aes(x = log10(res_low[, "bmd"]),
                   xend = log10(res_low[, "bmd"]), y = 0, 
                   yend = 30),col = "blue")

# function results
res_low['Min. Conc.'] <- min(conc2)
res_low['Name'] <- name
res_low[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")] <- round(res_low[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")], 3)

# using the argument to set a lower bound for BMD
res_low2 <- concRespCore(row_low,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", 
                                           "pow", "exp2", "exp3", "exp4", "exp5"), 
                         conthits = T, aicc = F, bidirectional=F, bmd_low_bnd = 0.8)
#> Warning in fitpoly2(conc = c(79.9999999999993, 0.600000000000001,
#> 29.9999999999998, : The `bidirectional` argument is ignored when `biphasic =
#> TRUE`.

# print out the new results
# include previous results side by side for comparison 
res_low2['Min. Conc.'] <- min(conc2)
res_low2['Name'] <- paste(name, "after `bounding`", sep = "-")
res_low['Name'] <- paste(name, "before `bounding`", sep = "-")
res_low2[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")] <- round(res_low2[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")], 3)

output_low <- rbind(res_low[1, c('Name', "Min. Conc.", "bmd", "bmdl", "bmdu")], 
                    res_low2[1, c('Name', "Min. Conc.", "bmd", "bmdl", "bmdu")])

# generate some concentration for the fitted curve 
logc_plot <- seq(from=-3,to=2,by=0.05)
conc_plot <- 10**logc_plot

# initiate the plot
plot(conc2,resp2,xlab="conc (uM)",ylab="Response",xlim=c(0.001,100),ylim=c(-5,60),
       log="x",main=paste(name,"\n",assay),cex.main=0.9)

# add vertical lines to mark the minimum concentration in the data and the lower threshold set by bmd_low_bnd
abline(v=min(conc2), lty = 1, col = "brown", lwd = 2)
abline(v=res_low2$bmd, lty = 2, col = "darkviolet", lwd = 2)

# add markers for BMD and its boundaries before `bounding`
lines(c(res_low$bmd,res_low$bmd),c(0,50),col="green",lwd=2)
rect(xleft=res_low$bmdl,ybottom=0,xright=res_low$bmdu,ytop=50,col=rgb(0,1,0, alpha = .5), border = NA)
points(res_low$bmd, 0, pch = "x", col = "green")

# add markers for BMD and its boundaries after `bounding`
lines(c(res_low2$bmd,res_low2$bmd),c(0,50),col="blue",lwd=2)
rect(xleft=res_low2$bmdl,ybottom=0,xright=res_low2$bmdu,ytop=50,col=rgb(0,0,1, alpha = .5), border = NA)
points(res_low2$bmd, 0, pch = "x", col = "blue")

# add the fitted curve
lines(conc_plot, exp4(ps = c(res_low$tp, res_low$ga), conc_plot))
legend(1e-3, 60, legend=c("Lowest Dose Tested", "Boundary", "BMD-before", "BMD-after"),
       col=c("brown", "darkviolet", "green", "blue"), lty=c(1,2,1,1))

# load example data
spid <- unique(mc3$spid)[26]
ex_df <- mc3[is.element(mc3$spid,spid),]

# The data file has stored concentration in log10 form, so fix that
conc <- 10**ex_df$logc # back-transforming concentrations on log10 scale
resp <- ex_df$resp

# pull out all of the chemical identifiers and the name of the assay
dtxsid <- ex_df[1,"dtxsid"]
casrn <- ex_df[1,"casrn"]
name <- ex_df[1,"name"]
assay <- ex_df[1,"assay"]

# create the row object
row_up <- list(conc = conc, resp = resp, bmed = 0, cutoff = cutoff, onesd = onesd,assay=assay,
            dtxsid=dtxsid,casrn=casrn,name=name)

# run the concentration-response modeling for a single sample
res_up <- concRespCore(row_up,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", 
                                         "pow", "exp2", "exp3", "exp4", "exp5"), 
                       conthits = T, aicc = F, bidirectional=F)
#> Warning in fitpoly2(conc = c(0.0199999999999999, 0.3, 0.0599999999999995, : The
#> `bidirectional` argument is ignored when `biphasic = TRUE`.

concRespPlot2(res_up, log_conc = T)

# max conc
res_up['Max Conc.'] <- max(conc)
res_up['Name'] <- name
res_up[1, c("Max Conc.", "bmd", "bmdl", "bmdu")] <- round(res_up[1, c("Max Conc.", "bmd", "bmdl", "bmdu")], 3)
# function results

# using bmd_up_bnd = 2
res_up2 <- concRespCore(row_up,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", 
                                          "pow", "exp2", "exp3", "exp4", "exp5"), 
                        conthits = T, aicc = F, bidirectional=F, bmd_up_bnd = 2)
#> Warning in fitpoly2(conc = c(0.0199999999999999, 0.3, 0.0599999999999995, : The
#> `bidirectional` argument is ignored when `biphasic = TRUE`.

# print out the new results
# include previous results side by side for comparison 
res_up2['Max Conc.'] <- max(conc)
res_up2['Name'] <- paste(name, "after `bounding`", sep = "-")
res_up['Name'] <- paste(name, "before `bounding`", sep = "-")
res_up2[1, c("Max Conc.", "bmd", "bmdl", "bmdu")] <- round(res_up2[1, c("Max Conc.", "bmd", "bmdl", "bmdu")], 3)

output_up <- rbind(res_up[1, c('Name', "Max Conc.", "bmd", "bmdl", "bmdu")], 
                   res_up2[1, c('Name', "Max Conc.", "bmd", "bmdl", "bmdu")])

# generate some concentration for the fitting curve 
logc_plot <- seq(from=-3,to=2,by=0.05)
conc_plot <- 10**logc_plot

# initiate plot
plot(conc,resp,xlab="conc (uM)",ylab="Response",xlim=c(0.001,500),ylim=c(-5,40),
       log="x",main=paste(name,"\n",assay),cex.main=0.9)
# add vertical lines to mark the maximum concentration in the data and the upper boundary set by bmd_up_bnd
abline(v=max(conc), lty = 1, col = "brown", lwd=2)
abline(v=160, lty = 2, col = "darkviolet", lwd=2)

# add marker for BMD and its boundaries before `bounding`
lines(c(res_up$bmd,res_up$bmd),c(0,50),col="green",lwd=2)
rect(xleft=res_up$bmdl,ybottom=0,xright=res_up$bmdu,ytop=50,col=rgb(0,1,0, alpha = .5), border = NA)
points(res_up$bmd, 0, pch = "x", col = "green")

# add marker for BMD and its boundaries after `bounding`
lines(c(res_up2$bmd,res_up2$bmd),c(0,50),col="blue",lwd=2)
rect(xleft=res_up2$bmdl,ybottom=0,xright=res_up2$bmdu,ytop=50,col=rgb(0,0,1, alpha = .5), border = NA)
points(res_up2$bmd, 0, pch = "x", col = "blue")

# add the fitting curve
lines(conc_plot, poly1(ps = c(res_up$a), conc_plot))
legend(1e-3, 40, legend=c("Maximum Dose Tested", "Boundary", "BMD-before", "BMD-after"),
       col=c("brown", "darkviolet", "green", "blue"), lty=c(1,2,1,1))

# using the same data, fit curves 
param <- tcplfit2_core(conc2, resp2, cutoff = cutoff)
hit_res <- tcplhit2_core(param, conc2, resp2, cutoff = cutoff, onesd = onesd, 
                         bmd_low_bnd = 0.8)

# adding the result from tcplhit2_core to the output table for comparison
hit_res["Name"]<-  paste("Chlorothalonil", "tcplhit2_core", sep = "-")
hit_res['Min. Conc.'] <- min(conc2)
hit_res[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")] <- round(hit_res[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")], 3)

output_low <- rbind(output_low, 
                    hit_res[1, c('Name', "Min. Conc.", "bmd", "bmdl", "bmdu")])

res_low3 <- concRespCore(row_low,fitmodels = c("cnst", "hill", "gnls", "poly1", "poly2", 
                                           "pow", "exp2", "exp3", "exp4", "exp5"), 
                         conthits = T, aicc = F, bidirectional=F, bmd_low_bnd = 0.4)
#> Warning in fitpoly2(conc = c(79.9999999999993, 0.600000000000001,
#> 29.9999999999998, : The `bidirectional` argument is ignored when `biphasic =
#> TRUE`.

# print out the new results
# add to previous results for comparison 
res_low3['Min. Conc.'] <- min(conc2)
res_low3['Name'] <- paste("Chlorothalonil", "after `bounding` (two fifths)", sep = "-")
res_low3[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")] <- round(res_low3[1, c("Min. Conc.", "bmd", "bmdl", "bmdu")], 3)

output_low <- rbind(output_low[-3, ], 
                    res_low3[1, c('Name', "Min. Conc.", "bmd", "bmdl", "bmdu")])

# initiate the plot
plot(conc2,resp2,xlab="conc (uM)",ylab="Response",xlim=c(0.001,100),ylim=c(-5,60),
       log="x",main=paste(name,"\n",assay),cex.main=0.9)

# add vertical lines to mark the minimum concentration in the data and the lower boundary set by bmd_low_bnd
abline(v=min(conc2), lty = 1, col = "brown", lwd = 2)
abline(v=0.4*min(conc2), lty = 2, col = "darkviolet", lwd = 2)

# add markers for BMD and its boundaries before `bounding`
lines(c(res_low$bmd,res_low$bmd),c(0,50),col="green",lwd=2)
rect(xleft=res_low$bmdl,ybottom=0,xright=res_low$bmdu,ytop=50,col=rgb(0,1,0, alpha = .5), border = NA)
points(res_low$bmd, 0, pch = "x", col = "green")

# add markers for BMD and its boundaries after `bounding`
lines(c(res_low3$bmd,res_low3$bmd),c(0,50),col="blue",lwd=2)
rect(xleft=res_low3$bmdl,ybottom=0,xright=res_low3$bmdu,ytop=50,col=rgb(0,0,1, alpha = .5), border = NA)
points(res_low3$bmd, 0, pch = "x", col = "blue")

# add the fitted curve
lines(conc_plot, exp4(ps = c(res_low$tp, res_low$ga), conc_plot))
legend(1e-3, 60, legend=c("Lowest Dose Tested", "Boundary Dose", "BMD-before", "BMD-after"),
       col=c("brown", "darkviolet", "green", "blue"), lty=c(1,2,1,1))

# Load the example data set
data("signatures")
# using the first row of signatures data as an example 
signatures[1,]
#>      sample_id        dtxsid                     name
#> 1 TP0001651A05 DTXSID8020337 Clomiphene citrate (1:1)
#>                              signature    cutoff     onesd
#> 1 CMAP mometasone 7.6e-06 100 2494 100 0.2393037 0.1196519
#>                     conc
#> 1 0.03|0.1|0.3|1|10|3|30
#>                                                                                                                                 resp
#> 1 -0.267223789187984|0.0446549190032419|-0.0150009098668537|0.097550989137631|0.420105528343523|0.0353858177275554|0.657069577568875

# using the first row of signature as an example 
conc <- as.numeric(str_split(signatures[1,"conc"],"\\|")[[1]])
resp <- as.numeric(str_split(signatures[1,"resp"],"\\|")[[1]])
cutoff <- signatures[1,"cutoff"]

# run curve fitting
output <- tcplfit2_core(conc, resp, cutoff)
# show the structure of the output 
summary(output)
#>            Length Class  Mode     
#> cnst        5     -none- list     
#> hill       17     -none- list     
#> gnls       22     -none- list     
#> poly1      13     -none- list     
#> poly2      17     -none- list     
#> pow        15     -none- list     
#> exp2       15     -none- list     
#> exp3       17     -none- list     
#> exp4       15     -none- list     
#> exp5       17     -none- list     
#> modelnames 10     -none- character
#> errfun      1     -none- character

# get plots in normal and in log-10 concentration scale
basic <- plot_allcurves(output, conc, resp)
basic_log <- plot_allcurves(output, conc, resp, log_conc = T)
grid.arrange(basic, basic_log)

# prepare the 'row' object for concRespCore
row <- list(conc=conc,
           resp=resp,
           bmed=0,
           cutoff=cutoff,
           onesd=signatures[1,"onesd"],
           name=signatures[1,"name"],
           assay=signatures[1,"signature"])

# run concentration-response modeling 
out <-  concRespCore(row,conthits=F)
# show the output
out
#>                       name                                assay n_gt_cutoff
#> 1 Clomiphene citrate (1:1) CMAP mometasone 7.6e-06 100 2494 100           3
#>      cutoff fit_method top_over_cutoff      rmse  a  b        tp  p  q       ga
#> 1 0.2393037       exp4        3.130382 0.1135895 NA NA 0.7491121 NA NA 9.591142
#>   la      er       bmr     bmdl     bmdu caikwt mll hitcall     ac50 ac50_loss
#> 1 NA -2.6244 0.1614104 2.380694 4.848624     NA  NA       1 9.591142        NA
#>         top       ac5     ac10     ac20      acc    ac1sd      bmd
#> 1 0.7491121 0.7097501 1.457883 3.087658 5.325258 2.408014 3.357835
#>                     conc
#> 1 0.03|0.1|0.3|1|10|3|30
#>                                                                                                                                 resp
#> 1 -0.267223789187984|0.0446549190032419|-0.0150009098668537|0.097550989137631|0.420105528343523|0.0353858177275554|0.657069577568875
#>   errfun
#> 1    dt4

# pass the output to the plotting function
basic_plot <- concRespPlot2(out)
basic_log <- concRespPlot2(out, log_conc = TRUE)
res <- grid.arrange(basic_plot, basic_log)

# Using the fitted result and plot from the example in the last section
# get the cutoff from the output
cutoff <- out[, "cutoff"]

basic_plot + 
  # Cutoff Band - a transparent rectangle
  geom_rect(aes(xmin = 0,xmax = 30,ymin = -cutoff,ymax = cutoff),
            alpha = 0.1,fill = "skyblue") +
  # Titles
  ggtitle(
    label = paste("Best Model Fit",
                  out[, "name"],
                  sep = "\n"),
    subtitle = paste("Assay Endpoint: ",
                     out[, "assay"])) +
  ## Add BMD and BMR labels
  geom_hline(
    aes(yintercept = out[, "bmr"]),
    col = "blue") +
  geom_segment(
    aes(x = out[, "bmd"], xend = out[, "bmd"], y = -0.5, yend = out[, "bmr"]),
    col = "blue"
  ) + geom_point(aes(x = out[, "bmd"], y = out[, "bmr"], fill = "BMD"), shape = 21, cex = 2.5)
#> Warning in geom_segment(aes(x = out[, "bmd"], xend = out[, "bmd"], y = -0.5, : All aesthetics have length 1, but the data has 7 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.
#> Warning in geom_point(aes(x = out[, "bmd"], y = out[, "bmr"], fill = "BMD"), : All aesthetics have length 1, but the data has 7 rows.
#> ℹ Please consider using `annotate()` or provide this layer with data containing
#>   a single row.

# Get all potency estimates and the corresponding y value on the curve
estimate_points <- out %>%
  select(bmd, acc, ac50, ac10, ac5) %>%
  tidyr::pivot_longer(everything(), names_to = "Potency Estimates") %>%
  mutate(`Potency Estimates` = toupper(`Potency Estimates`)) 

y <-  c(out[, "bmr"], out[, "cutoff"], rep(out[, "top"], 3))
y <-  y * c(1, 1, .5, .1, .05)
estimate_points <- cbind(estimate_points, y = y)

# add Potency Estimate Points and set colors
basic_plot + geom_point(
  data = estimate_points,
  aes(x = value, y = y, fill = `Potency Estimates`), shape = 21, cex = 2.5
)

# add Potency Estimate Points and set colors - with plot in log-10 concentration
basic_log + geom_point(
  data = estimate_points,
  aes(x = log10(value), y = y, fill = `Potency Estimates`), shape = 21, cex = 2.5
)

# maybe want to extract and use the same x's in the base plot 
# to calculate predicted responses 
conc_plot <- basic_plot[["layers"]][[2]][["data"]][["conc_plot"]]

basic_plot +
  # fitted parameter values of another curve you want to add
  geom_line(data=data.frame(x=conc_plot, y=tcplfit2::exp5(c(0.5, 10, 1.2), conc_plot)), aes(x,y,color = "exp5"))+
  # add different colors for comparisons 
  scale_colour_manual(values=c("#CC6666", "#9999CC"),
                      labels = c("Curve 1-exp4", "Curve 2-exp5")) +
  labs(title = "Curve 1 v.s. Curve 2")

# some example data
conc <- list(.03, .1, .3, 1, 3, 10, 30, 100)
resp <- list(0, .2, .1, .4, .7, .9, .6, 1.2)
row <- list(conc = conc,
            resp = resp,
            bmed = 0,
            cutoff = 1,
            onesd = .5)

# AUC is included in the output
concRespCore(row, conthits = TRUE, AUC = TRUE)
#>      n_gt_cutoff cutoff fit_method top_over_cutoff      rmse  a  b       tp
#> cnst           1      1       hill        1.225599 0.1750312 NA NA 1.225599
#>              p  q       ga la        er    bmr     bmdl     bmdu       caikwt
#> cnst 0.7752844 NA 2.554272 NA -2.467853 0.6745 2.341811 4.822553 6.064845e-05
#>           mll   hitcall     ac50 ac50_loss      top        ac5      ac10
#> cnst 4.492301 0.9650614 2.554272        NA 1.225599 0.05726215 0.1501198
#>           ac20      acc    ac1sd      bmd                       conc
#> cnst 0.4272681 17.43267 1.580024 3.314775 0.03|0.1|0.3|1|3|10|30|100
#>                               resp errfun      AUC
#> cnst 0|0.2|0.1|0.4|0.7|0.9|0.6|1.2    dt4 106.1212

# This is taken from the example under tcplfit2_core
conc_ex2 <- c(.03, .1, .3, 1, 3, 10, 30, 100)
resp_ex2 <- c(0, .1, 0, .2, .6, .9, 1.1, 1)

# fit all available models in the package
# show all fitted curves 
output_ex2 <- tcplfit2_core(conc_ex2, resp_ex2, .8)
grid.arrange(plot_allcurves(output_ex2, conc_ex2, resp_ex2),
          plot_allcurves(output_ex2, conc_ex2, resp_ex2, log_conc = TRUE), ncol = 2)

fit_method <- "hill"
# extract the parameters 
modpars <- output_ex2[[fit_method]][output_ex2[[fit_method]]$pars]

# plug into get_AUC function 
estimated_auc1 <- get_AUC(fit_method, min(conc_ex2), max(conc_ex2), modpars)
estimated_auc1
#> [1] 97.41475

# extract the predicted responses from the model
pred_resp <- output_ex2[[fit_method]][["modl"]]

# plot to see if the result make sense
# the shaded area is what the function tries to find
plot(conc_ex2, pred_resp)
lines(conc_ex2, pred_resp)
polygon(c(conc_ex2, max(conc_ex2)), c(pred_resp, min(pred_resp)), col=rgb(1, 0, 0,0.5))

# list of models
fitmodels <- c("gnls", "poly1", "poly2", "pow", "exp2", "exp3", "exp4", "exp5")
mylist <- list()
for (model in fitmodels){

  fit_method <- model
  # extract corresponding model parameters
  modpars <- output_ex2[[fit_method]][output_ex2[[fit_method]]$pars]
  
  # get AUC
  mylist[[fit_method]] <- get_AUC(fit_method, min(conc_ex2), max(conc_ex2), modpars)
  
}
# print AUC's for other models 
data.frame(mylist,row.names = "AUC")
#>        gnls    poly1    poly2      pow     exp2     exp3     exp4     exp5
#> AUC 97.4147 58.09263 121.7488 97.43487 55.01408 96.18956 98.80625 98.65734

# Taking the code from example 3 in the vignette 
library(stringr)  # string management package
data("signatures")

# use row 5 in the data
conc <- as.numeric(str_split(signatures[5,"conc"],"\\|")[[1]])
resp <- as.numeric(str_split(signatures[5,"resp"],"\\|")[[1]])
cutoff <- signatures[5,"cutoff"]

# plot all models, this is an example of negative curves 
output_negative <- tcplfit2_core(conc, resp, cutoff)
grid.arrange(plot_allcurves(output_negative, conc, resp),
          plot_allcurves(output_negative, conc, resp, log_conc = TRUE), ncol = 2)

fit_method <- "exp3"

# extract corresponding model parameters and predicted response
modpars <- output_negative[[fit_method]][output_negative[[fit_method]]$pars]
pred_resp <- output_negative[[fit_method]][["modl"]]

estimated_auc2 <- get_AUC(fit_method, min(conc), max(conc), modpars)
estimated_auc2
#> [1] -12.92738

# plot this curve
pred_resp <- pred_resp[order(conc)]
plot(conc[order(conc)], pred_resp)
lines(conc[order(conc)], pred_resp)
polygon(c(conc[order(conc)], max(conc)), c(pred_resp, max(pred_resp)), col=rgb(1, 0, 0,0.5))

get_AUC(fit_method, min(conc), max(conc), modpars, return.abs = TRUE) 
#> [1] 12.92738

# simulate a poly2 curve
conc_sim <- seq(0,3, length.out = 100)
## biphasic poly2 parameters
b1 <- -1.3
b2 <- 0.7
## converted to tcplfit2's poly2 parameters
a <- b1^2/b2
b <- b1/b2

## plot the curve
resp_sim <- poly2(c(a, b, 0.1), conc_sim)
plot(conc_sim, resp_sim, type = "l")
abline(h = 0)

# get AUC for the simulated Polynomial 2 curve 
get_AUC("poly2", min(conc_sim), max(conc_sim), ps = c(a, b))
#> [1] 0.45

out <- tcplhit2_core(output_ex2, conc_ex2, resp_ex2, 0.8, onesd = 0.4)
out
#>      n_gt_cutoff cutoff fit_method top_over_cutoff       rmse  a  b       tp
#> cnst           3    0.8       hill        1.279049 0.05022924 NA NA 1.023239
#>             p  q       ga la        er    bmr    bmdl     bmdu       caikwt
#> cnst 1.592714 NA 2.453014 NA -3.295307 0.5396 2.33877 2.979982 1.446965e-08
#>           mll   hitcall     ac50 ac50_loss      top       ac5      ac10
#> cnst 12.71495 0.9999999 2.453014        NA 1.023239 0.3862094 0.6174049
#>          ac20      acc    ac1sd      bmd                       conc
#> cnst 1.027285 5.466819 1.856853 2.627574 0.03|0.1|0.3|1|3|10|30|100
#>                           resp errfun
#> cnst 0|0.1|0|0.2|0.6|0.9|1.1|1    dt4
post_hit_AUC(out)
#> [1] 97.41475

poly1_plot <- ggplot(mapping=aes(ex_conc)) +  
  geom_line(aes(y = 55*ex_conc, color = "a=55")) +
  geom_line(aes(y = 10*ex_conc, color = "a=10")) +
  geom_line(aes(y = 0.05*ex_conc, color = "a=0.05")) +
  geom_line(aes(y = -5*ex_conc, color = "a=(-5)")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.1,0.8)) +
  scale_color_manual(name='a values',
                     breaks=c('a=(-5)', 'a=0.05', 'a=10', 'a=55'),
                     values=c('a=(-5)'='black', 'a=0.05' = 'red', 'a=10'='blue', 'a=55'='darkviolet'))

poly1_plot

fits_poly <- data.frame(
  # change a 
  y1 = poly2(ps = c(a = 40, b = 2),x = ex_conc),
  y2 = poly2(ps = c(a = 6, b = 2),x = ex_conc),
  y3 = poly2(ps = c(a = 0.1, b = 2),x = ex_conc),
  y4 = poly2(ps = c(a = -2, b = 2),x = ex_conc),
  y5 = poly2(ps = c(a = -20, b = 2),x = ex_conc),
  
  # change b 
  y6 = poly2(ps = c(a = 4,b = 1.8),x = ex_conc),
  y7 = poly2(ps = c(a = 4,b = 7),x = ex_conc),
  y8 = poly2(ps = c(a = 4,b = 16),x = ex_conc)
)

# shows how changes in parameter 'a' affect the shape of the curve 
poly2_plot1 <- ggplot(fits_poly, aes(ex_conc)) +
  geom_line(aes(y = y1, color = "a=40")) +
  geom_line(aes(y = y2, color = "a=6")) +
  geom_line(aes(y = y3, color = "a=0.1")) +
  geom_line(aes(y = y4, color = "a=(-2)")) +
  geom_line(aes(y = y5, color = "a=(-20)")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='a values',
                     breaks=c('a=(-20)', 'a=(-2)', 'a=0.1', 'a=6', 'a=40'),
                     values=c('a=(-20)'='black', 'a=(-2)'='red', 'a=0.1'='blue', 'a=6'='darkviolet', 'a=40'='darkgoldenrod1'))

# shows how changes in parameter 'b' affect the shape of the curve 
poly2_plot2 <- ggplot(fits_poly, aes(ex_conc)) +  
  geom_line(aes(y = y6, color = "b=1.8")) +
  geom_line(aes(y = y7, color = "b=7")) +
  geom_line(aes(y = y8, color = "b=16")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='b values',
                     breaks=c('b=1.8', 'b=7', 'b=16'),
                     values=c('b=1.8'='black', 'b=7'='red', 'b=16'='blue'))

grid.arrange(poly2_plot1, poly2_plot2, ncol = 2)

fits_pow <- data.frame(
  # change a
  y1 = pow(ps = c(a = 0.48,p = 1.45),x = ex2_conc),
  y2 = pow(ps = c(a = 7.2,p = 1.45),x = ex2_conc),
  y3 = pow(ps = c(a = -3.2,p = 1.45),x = ex2_conc),
  
  # change p
  y4 = pow(ps = c(a = 1.2,p = 0.3),x = ex2_conc),
  y5 = pow(ps = c(a = 1.2,p = 1.6),x = ex2_conc),
  y6 = pow(ps = c(a = 1.2,p = 3.2),x = ex2_conc)
)

# shows how changes in parameter 'a' affect the shape of the curve
pow_plot1 <- ggplot(fits_pow, aes(ex2_conc)) +  
  geom_line(aes(y = y1, color = "a=0.48")) +
  geom_line(aes(y = y2, color = "a=7.2")) +
  geom_line(aes(y = y3, color = "a=(-3.2)")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='a values',
                     breaks=c('a=(-3.2)', 'a=0.48', 'a=7.2'),
                     values=c('a=(-3.2)'='black', 'a=0.48'='red', 'a=7.2'='blue'))

# shows how changes in parameter 'p' affect the shape of the curve
pow_plot2 <- ggplot(fits_pow, aes(ex2_conc)) +  
  geom_line(aes(y = y4, color = "p=0.3")) +
  geom_line(aes(y = y5, color = "p=1.6")) +
  geom_line(aes(y = y6, color = "p=3.2")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='p values',
                     breaks=c('p=0.3', 'p=1.6', 'p=3.2'),
                     values=c('p=0.3'='black', 'p=1.6'='red', 'p=3.2'='blue'))



grid.arrange(pow_plot1, pow_plot2, ncol = 2)

fits_hill <- data.frame(
  # change tp
  y1 = hillfn(ps = c(tp = -200,ga = 5,p = 1.76), x = ex_conc),
  y2 = hillfn(ps = c(tp = 200,ga = 5,p = 1.76), x = ex_conc),
  y3 = hillfn(ps = c(tp = 850,ga = 5,p = 1.76), x = ex_conc),

  # change ga
  y4 = hillfn(ps = c(tp = 120,ga = 4,p = 1.76), x = ex_conc),
  y5 = hillfn(ps = c(tp = 120,ga = 12,p = 1.76), x = ex_conc),
  y6 = hillfn(ps = c(tp = 120,ga = 20,p = 1.76), x = ex_conc),
  
  # change p
  y7 = hillfn(ps = c(tp = 120,ga = 5,p = 0.5), x = ex_conc),
  y8 = hillfn(ps = c(tp = 120,ga = 5,p = 2), x = ex_conc),
  y9 = hillfn(ps = c(tp = 120,ga = 5,p = 5), x = ex_conc)
  
)

# shows how changes in parameter 'tp' affect the shape of the curve
hill_plot1 <- ggplot(fits_hill, aes(log10(ex_conc))) +  
  geom_line(aes(y = y1, color = "tp=(-200)")) +
  geom_line(aes(y = y2, color = "tp=200")) +
  geom_line(aes(y = y3, color = "tp=850")) +
  labs(x = "Concentration in Log-10 Scale", y = "Response") +
  theme(legend.position = c(0.2,0.8),
        legend.key.size = unit(0.5, 'cm')) +
  scale_color_manual(name='tp values',
                     breaks=c('tp=(-200)', 'tp=200', 'tp=850'),
                     values=c('tp=(-200)'='black', 'tp=200'='red', 'tp=850'='blue'))

# shows how changes in parameter 'ga' affect the shape of the curve
hill_plot2 <- ggplot(fits_hill, aes(log10(ex_conc))) + 
  geom_line(aes(y = y4, color = "ga=4")) +
  geom_line(aes(y = y5, color = "ga=12")) +
  geom_line(aes(y = y6, color = "ga=20")) +
  labs(x = "Concentration in Log-10 Scale", y = "Response") +
  theme(legend.position = c(0.8,0.25),
        legend.key.size = unit(0.4, 'cm')) +
  scale_color_manual(name='ga values',
                     breaks=c('ga=4', 'ga=12', 'ga=20'),
                     values=c('ga=4'='black', 'ga=12'='red', 'ga=20'='blue'))

# shows how changes in parameter 'p' affect the shape of the curve
hill_plot3 <- ggplot(fits_hill, aes(log10(ex_conc))) +  
  geom_line(aes(y = y7, color = "p=0.5")) +
  geom_line(aes(y = y8, color = "p=2")) +
  geom_line(aes(y = y9, color = "p=5")) +
  labs(x = "Concentration in Log-10 Scale", y = "Response") +
  theme(legend.position = c(0.8,0.2),
        legend.key.size = unit(0.4, 'cm')) +
  scale_color_manual(name='p values',
                     breaks=c('p=0.5', 'p=2', 'p=5'),
                     values=c('p=0.5'='black', 'p=2'='red', 'p=5'='blue'))


grid.arrange(hill_plot1, hill_plot2, hill_plot3, ncol = 2, nrow = 2)

fits_gnls <- data.frame(
  # change la
  y1 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 17,q = 1.34), x = ex_conc),
  y2 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 50,q = 1.34), x = ex_conc),
  y3 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 100,q = 1.34), x = ex_conc),

  # change q
  y4 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 20,q = 0.3), x = ex_conc),
  y5 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 20,q = 1.2), x = ex_conc),
  y6 = gnls(ps = c(tp = 750,ga = 15,p = 1.45,la = 20,q = 8), x = ex_conc)
  
)

# shows how changes in parameter 'la' affect the shape of the curve
gnls_plot1 <- ggplot(fits_gnls, aes(log10(ex_conc))) +  
  geom_line(aes(y = y1, color = "la=17")) +
  geom_line(aes(y = y2, color = "la=50")) +
  geom_line(aes(y = y3, color = "la=100")) +
  labs(x = "Concentration in Log-10 Scale", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='la values',
                     breaks=c('la=17', 'la=50', 'la=100'),
                     values=c('la=17'='black', 'la=50'='red', 'la=100'='blue'))

# shows how changes in parameter 'q' affect the shape of the curve
gnls_plot2 <- ggplot(fits_gnls, aes(log10(ex_conc))) +  
  geom_line(aes(y = y4, color = "q=0.3")) +
  geom_line(aes(y = y5, color = "q=1.2")) +
  geom_line(aes(y = y6, color = "q=8")) +
  labs(x = "Concentration in Log-10 Scale", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='q values',
                     breaks=c('q=0.3', 'q=1.2', 'q=8'),
                     values=c('q=0.3'='black', 'q=1.2'='red', 'q=8'='blue'))
  
  
grid.arrange(gnls_plot1, gnls_plot2, ncol = 2)

fits_exp2 <- data.frame(
  # change a
  y1 = exp2(ps = c(a = 20,b = 12), x = ex2_conc),
  y2 = exp2(ps = c(a = 9,b = 12), x = ex2_conc),
  y3 = exp2(ps = c(a = 0.1,b = 12), x = ex2_conc),
  y4 = exp2(ps = c(a = -3,b = 12), x = ex2_conc),
  
  # change b
  y5 = exp2(ps = c(a = 0.45,b = 4), x = ex2_conc),
  y6 = exp2(ps = c(a = 0.45,b = 9), x = ex2_conc),
  y7 = exp2(ps = c(a = 0.45,b = 20), x = ex2_conc)
  
)

# shows how changes in parameter 'a' affect the shape of the curve 
exp2_plot1 <- ggplot(fits_exp2, aes(ex2_conc)) +  
  geom_line(aes(y = y1, color = "a=20")) +
  geom_line(aes(y = y2, color = "a=9")) +
  geom_line(aes(y = y3, color = "a=0.1")) +
  geom_line(aes(y = y4, color = "a=(-3)")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='a values',
                     breaks=c('a=(-3)', 'a=0.1', 'a=9', 'a=20'),
                     values=c('a=(-3)'='black', 'a=0.1'='red', 'a=9'='blue', 'a=20'='darkviolet'))

# shows how changes in parameter 'b' affect the shape of the curve 
exp2_plot2 <- ggplot(fits_exp2, aes(ex2_conc)) +  
  geom_line(aes(y = y5, color = "b=4")) +
  geom_line(aes(y = y6, color = "b=9")) +
  geom_line(aes(y = y7, color = "b=20")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='b values',
                     breaks=c('b=4', 'b=9', 'b=20'),
                     values=c('b=4'='black', 'b=9'='red', 'b=20'='blue'))

grid.arrange(exp2_plot1, exp2_plot2, ncol = 2)

fits_exp3 <- data.frame(
  # change p
  y1 = exp3(ps = c(a = 1.67,b = 12.5,p = 0.3), x = ex2_conc),
  y2 = exp3(ps = c(a = 1.67,b = 12.5,p = 0.9), x = ex2_conc),
  y3 = exp3(ps = c(a = 1.67,b = 12.5,p = 1.2), x = ex2_conc)
  
)

# shows how changes in parameter 'p' affect the shape of the curve 
exp3_plot <- ggplot(fits_exp3, aes(ex2_conc)) +  
  geom_line(aes(y = y1, color = "p=0.3")) +
  geom_line(aes(y = y2, color = "p=0.9")) +
  geom_line(aes(y = y3, color = "p=1.2")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.15,0.8)) +
  scale_color_manual(name='p values',
                     breaks=c('p=0.3', 'p=0.9', 'p=1.2'),
                     values=c('p=0.3'='black', 'p=0.9'='red', 'p=1.2'='blue'))


exp3_plot

fits_exp4 <- data.frame(
  # change tp  
  y1 = exp4(ps = c(tp = 895,ga = 15),x = ex_conc),
  y2 = exp4(ps = c(tp = 200,ga = 15),x = ex_conc),
  y3 = exp4(ps = c(tp = -500,ga = 15),x = ex_conc),
  
  # change ga
  y4 = exp4(ps = c(tp = 500,ga = 0.4),x = ex_conc),
  y5 = exp4(ps = c(tp = 500,ga = 10),x = ex_conc),
  y6 = exp4(ps = c(tp = 500,ga = 20),x = ex_conc)
  
)

# shows how changes in parameter 'tp' affect the shape of the curve 
exp4_plot1 <- ggplot(fits_exp4, aes(ex_conc)) +  
  geom_line(aes(y = y1, color = "tp=895")) +
  geom_line(aes(y = y2, color = "tp=200")) +
  geom_line(aes(y = y3, color = "tp=(-500)")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.1,0.8)) +
  scale_color_manual(name='tp values',
                     breaks=c('tp=(-500)', 'tp=200', 'tp=895'),
                     values=c('tp=(-500)'='black', 'tp=200'='red', 'tp=895'='blue'))


# shows how changes in parameter 'ga' affect the shape of the curve 
exp4_plot2 <- ggplot(fits_exp4, aes(ex_conc)) +  
  geom_line(aes(y = y4, color = "ga=0.4")) +
  geom_line(aes(y = y5, color = "ga=10")) +
  geom_line(aes(y = y6, color = "ga=20")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.8,0.2)) +
  scale_color_manual(name='ga values',
                     breaks=c('ga=0.4', 'ga=10', 'ga=20'),
                     values=c('ga=0.4'='black', 'ga=10'='red', 'ga=20'='blue'))


grid.arrange(exp4_plot1, exp4_plot2, ncol = 2)

fits_exp5 <- data.frame(
  # change p
  y1 = exp5(ps = c(tp = 793,ga = 6.25,p = 0.3), x = ex_conc),
  y2 = exp5(ps = c(tp = 793,ga = 6.25,p = 3.4), x = ex_conc),
  y3 = exp5(ps = c(tp = 793,ga = 6.25,p = 8), x = ex_conc)
  
)

# shows how changes in parameter 'p' affect the shape of the curve 
exp5_plot <- ggplot(fits_exp5, aes(ex_conc)) +  
  geom_line(aes(y = y1, color = "p=0.3")) +
  geom_line(aes(y = y2, color = "p=3.4")) +
  geom_line(aes(y = y3, color = "p=8")) +
  labs(x = "Concentration", y = "Response") +
  theme(legend.position = c(0.8,0.2)) +
  scale_color_manual(name='p values',
                     breaks=c('p=0.3', 'p=3.4', 'p=8'),
                     values=c('p=0.3'='black', 'p=3.4'='red', 'p=8'='blue'))


exp5_plot