## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
Sys.setenv(OMP_THREAD_LIMIT = 1) # Reducing core use, to avoid accidental use of too many cores
library(Colossus)
library(data.table)
if (system.file(package = "ggplot2") != "") {
  library(ggplot2)
}

## ----eval=TRUE----------------------------------------------------------------
# Suppose we have a linear-quadratic model
# an initial slope of 0.2 and a threshold of 5
para_0 <- c(0.2, 5)
# We might also have a linear-quadratic model
# same slope and threshold
# also an exponential slope of 0.1
para_1 <- c(0.2, 5, 0.1)

paras <- list(cov_0 = para_0, cov_1 = para_1)

# We pass a list of subterm formula options,
# either 'quad' or 'exp' to choose subterm options.

# The names of the tform and para lists should match.
tforms <- list(cov_0 = "quad", cov_1 = "exp")

res <- Linked_Dose_Formula(tforms, paras)
res_quad <- round(res$cov_0, 3)
res_exp <- round(res$cov_1, 3)

print("Linear-Quadratic Model")
print(paste0("Threshold:", res_quad[1]))
print(paste0("linear slope:", res_quad[2]))
print(paste0("quadratic slope:", res_quad[3]))
print(paste0("piecewise intercept:", res_quad[4]))

print("Linear-Exponential Model")
print(paste0("Threshold:", res_exp[1]))
print(paste0("linear slope:", res_exp[2]))
print(paste0("piece-wise intercept:", res_exp[3]))
print(paste0("exponential slope:", res_exp[4]))
print(paste0("exponential intercept:", res_exp[5]))

x <- (0:50) / 50 * 15
y0 <- ifelse(x < res_quad[1], res_quad[2] * x, res_quad[3] * x^2 + res_quad[4])
y1 <- ifelse(x < res_exp[1], res_exp[2] * x, res_exp[3] + exp(res_exp[4] * x + res_exp[5]))
y <- c(y0, y1)
c <- c(rep("LQ", length(x)), rep("LEXP", length(x)))
df <- data.table("x" = x, "y" = y, "model" = c)
if (system.file(package = "ggplot2") != "") {
  g <- ggplot2::ggplot(df, ggplot2::aes(x = .data$x, y = .data$y, group = .data$model, color = .data$model)) +
    ggplot2::geom_line(linewidth = 1.2) +
    labs(x = "Dose", y = "Response")
} else {
  g <- message("ggplot2 wasn't detected. Please install to see the plot")
}
g

## ----eval=TRUE----------------------------------------------------------------
# Suppose we used the same slope and intercept,
# our asymptote can be any value greater than 5*0.2
# Suppose we want the asymptote to be 3
threshold <- 5
slope <- 0.2

res <- Linked_Lin_Exp_Para(threshold, slope, 3)
print(paste0("Exponential slope: ", round(res, 3)))

x <- c()
y <- c()
c <- c()
i <- 0
for (slope in c(-0.2, 0.2)) {
  for (max in c(2, 4)) {
    i <- i + 1
    res <- Linked_Lin_Exp_Para(threshold, slope, sign(slope) * max)
    res <- Linked_Dose_Formula(list(temp = "exp"), list(temp = c(slope, threshold, res)))$temp
    xt <- 0:50
    yt <- ifelse(xt < res[1], res[2] * xt, res[3] - sign(slope) * exp(res[4] * xt + res[5]))
    x <- c(x, xt)
    y <- c(y, yt)
    c <- c(c, rep(i, length(xt)))
  }
}
df <- data.table("x" = x, "y" = y, "model" = c)
if (system.file(package = "ggplot2") != "") {
  g <- ggplot2::ggplot(df, ggplot2::aes(x = .data$x, y = .data$y, group = .data$model, color = .data$model)) +
    ggplot2::geom_line(linewidth = 1.2) +
    labs(x = "Dose", y = "Response")
} else {
  g <- message("ggplot2 wasn't detected. Please install to see the plot")
}
g