install.packages("drda")

# install.packages("remotes")
remotes::install_github("albertopessia/drda")

library(drda)

?voropm2

head(voropm2)

# by default `drda` uses a 4-parameter logistic function for model fitting

# common R API for fitting models
fit <- drda(response ~ log_dose, data = voropm2)

# get a general overview of the results
summary(fit)

# get parameter estimates by using generic functions...
coef(fit)
sigma(fit)

# ... or accessing the variables directly
fit$coefficients
fit$sigma

# compare the estimated model against a flat horizontal line, or the full
# 5-parameter logistic model, using AIC, BIC, and the Likelihood Ratio Test
# (LRT)
#
# note that the LRT is testing the null hypothesis of a flat horizontal line
# being as a good fit as the chosen model, therefore we expect the test to be
# significant
#
# if the test is not significant, a horizontal line is probably a better model
anova(fit)

# use the `mean_function` argument to select a different model
fit_l2 <- drda(response ~ log_dose, data = voropm2, mean_function = "logistic2")
fit_l4 <- drda(response ~ log_dose, data = voropm2, mean_function = "logistic4")
fit_l5 <- drda(response ~ log_dose, data = voropm2, mean_function = "logistic5")
fit_gz <- drda(response ~ log_dose, data = voropm2, mean_function = "gompertz")

# which model should be chosen?
anova(fit_l2, fit_l4, fit_l5, fit_gz)

# 5-parameter logistic function provides the best fit (AIC and BIC are minimum)

# it is possible to give each observation its own weight
fit_weighted <- drda(response ~ log_dose, data = voropm2, weights = weight)

# all the commands shown so far are available for a weighted fit as well

# it is possible to fix parameter values by setting the `lower_bound` and
# `upper_bound` appropriately
#
# unconstrained parameters have a lower bound of `-Inf` and an upper bound of
# `Inf`
#
# Important: be careful when deciding the constraints, because the optimization
#            problem might become very difficult to solve within a reasonable
#            number of iterations.
#
# In this particular example we are:
#   - fixing the alpha parameter to 1
#   - fixing the delta parameter to -1
#   - constraining the growth rate to be between 1 and 5
#   - not constraining the `phi` parameter, i.e. the `log(EC50)`
lb <- c(1, -1, 1, -Inf)
ub <- c(1, -1, 5,  Inf)

fit <- drda(
  response ~ log_dose, data = voropm2, lower_bound = lb, upper_bound = ub,
  max_iter = 260
)

summary(fit)

# if the algorithm does not converge, we can try to increase the maximum number
# of iterations or provide our own starting point
fit <- drda(
  response ~ log_dose, data = voropm2, lower_bound = lb, upper_bound = ub,
  start = c(1, -1, 2.6, 5), max_iter = 10000
)

summary(fit)

fit_l5 <- drda(response ~ log_dose, data = voropm2, mean_function = "logistic5")

# plot the data used for fitting, the maximum likelihood curve, and
# *approximate* confidence intervals for the curve
plot(fit_l5)

# combine all curves in the same plot
fit_l2 <- drda(response ~ log_dose, data = voropm2, mean_function = "logistic2")
fit_l4 <- drda(response ~ log_dose, data = voropm2, mean_function = "logistic4")
plot(fit_l2, fit_l4, fit_l5)

# modify default plotting options
# use `legend_show = FALSE` to remove the legend altogether
plot(
  fit_l5, base = "10", col = "magenta", xlab = "x", ylab = "y", level = 0.9,
  midpoint = FALSE, main = "Example plot", legend_location = "topright",
  legend = "5-parameter logistic function"
)