## ----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 = "survival") != "") {
  library(survival)
}
library(dplyr)

## ----eval=TRUE----------------------------------------------------------------
if (system.file(package = "survival") != "") {
  data(reliability, package = "survival")
  capacitor |> setDT()
  df <- copy(capacitor)
} else {
  voltage <- c(200, 200, 200, 200, 250, 250, 250, 250, 300, 300, 300, 300, 350, 350, 350, 350, 200, 200, 200, 200, 250, 250, 250, 250, 300, 300, 300, 300, 350, 350, 350, 350, 200, 200, 200, 200, 250, 250, 250, 250, 300, 300, 300, 300, 350, 350, 350, 350, 200, 200, 200, 200, 250, 250, 250, 250, 300, 300, 300, 300, 350, 350, 350, 350)
  temperature <- c(170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 170, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180, 180)
  time <- c(439, 904, 1092, 1105, 572, 690, 904, 1090, 315, 315, 439, 628, 258, 258, 347, 588, 959, 1065, 1065, 1087, 216, 315, 455, 473, 241, 315, 332, 380, 241, 241, 435, 455, 1105, 1105, 1105, 1105, 1090, 1090, 1090, 1090, 628, 628, 628, 628, 588, 588, 588, 588, 1087, 1087, 1087, 1087, 473, 473, 473, 473, 380, 380, 380, 380, 455, 455, 455, 455)
  status <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
  df <- data.table(
    voltage = voltage,
    temperature = temperature,
    time = time,
    status = status
  )
}

df$voltage <- (df$voltage - 200) / 150
df$temperature <- (df$temperature - 170) / 10
df$time <- (df$time - 216) / (1105 - 216)

control <- list(ncores = 1, maxiter = 100, verbose = 2)

a_n <- c(0.01, 0.01)

e1 <- CoxRun(Cox(time, status) ~ loglinear(temperature, voltage, 0), df,
  a_n = a_n, control = control
)
print(e1, 5)

e2 <- CoxRun(Cox(time, status) ~ loglinear(temperature, 0) + plinear(voltage, 0),
  df,
  a_n = a_n, control = control
)
print(e2, 5)

## ----eval=TRUE----------------------------------------------------------------
names <- c("temperature", "voltage")
tform <- c("loglin", "loglin")
ci_1 <- c(
  e1$beta_0[1] - 1.96 * e1$Standard_Deviation[1],
  e1$beta_0[1] + 1.96 * e1$Standard_Deviation[1]
)
ci_2 <- c(
  e1$beta_0[2] - 1.96 * e1$Standard_Deviation[2],
  e1$beta_0[2] + 1.96 * e1$Standard_Deviation[2]
)


curve_control <- list(
  maxstep = 100,
  alpha = 0.05,
  para_number = 1, manual = TRUE
)
e <- LikelihoodBound(e1, df, curve_control, control = control)
cat("|------------------- Wald Estimate -------------------|")
print(ci_1)
print(e, 5)

curve_control <- list(
  maxstep = 100,
  alpha = 0.05,
  para_number = 2, manual = TRUE
)
e <- LikelihoodBound(e1, df, curve_control, control = control)
cat("|------------------- Likelihood Bound Estimate -------------------|")
print(ci_2)
print(e, 5)

## ----eval=TRUE----------------------------------------------------------------
ci_1 <- c(
  e2$beta_0[1] - 1.96 * e2$Standard_Deviation[1],
  e2$beta_0[1] + 1.96 * e2$Standard_Deviation[1]
)
ci_2 <- c(
  e2$beta_0[2] - 1.96 * e2$Standard_Deviation[2],
  e2$beta_0[2] + 1.96 * e2$Standard_Deviation[2]
)

curve_control <- list(
  maxstep = 100,
  alpha = 0.05,
  para_number = 1, manual = TRUE
)
e <- LikelihoodBound(e2, df, curve_control, control = control)
cat("|------------------- Wald Estimate -------------------|")
print(ci_1)
print(e, 5)

a_n <- c(1.138152, 1.988403)
curve_control <- list(
  maxstep = 100,
  alpha = 0.05,
  para_number = 2, manual = TRUE
)
e <- LikelihoodBound(e2, df, curve_control, control = control)
cat("|------------------- Wald Estimate -------------------|")
print(ci_2)
print(e, 5)

## ----eval=FALSE---------------------------------------------------------------
# fname <- "base_example.csv"
# df <- fread(fname)
# 
# keep_constant <- c(0, 0, 1)
# a_n <- c(-1.493177, 5.020007, 1.438377)
# model <- Cox(entry, exit, event) ~ loglinear(dose0, dose1, 0) + linear(dose0, 1)
# #
# control <- list(
#   ncores = 1, lr = 0.75, maxiters = c(100, 100), halfmax = 5,
#   epsilon = 1e-6, deriv_epsilon = 1e-6, step_max = 1.0,
#   thres_step_max = 100.0, verbose = 2,
#   ties = "breslow"
# )
# 
# v0 <- sort(c((0:50 / 50 - 1.0) * 0.8, 1:50 / 50 * 3, 1.438377, -0.5909))
# for (v in v0) {
#   a_n <- c(-1.493177, 5.020007, v)
#   e <- CoxRun(model, df, a_n = a_n, control = control)
#   ll <- e$LogLik
#   beta <- e$beta_0
#   print(c(ll, beta[3]))
# }

## ----fig.width=7,fig.height=4-------------------------------------------------
x <- c(-0.8, -0.784, -0.768, -0.752, -0.736, -0.72, -0.704, -0.688, -0.672, -0.656, -0.64, -0.624, -0.608, -0.592, -0.5909, -0.576, -0.56, -0.544, -0.528, -0.512, -0.496, -0.48, -0.464, -0.448, -0.432, -0.416, -0.4, -0.384, -0.368, -0.352, -0.336, -0.32, -0.304, -0.288, -0.272, -0.256, -0.24, -0.224, -0.208, -0.192, -0.176, -0.16, -0.144, -0.128, -0.112, -0.096, -0.08, -0.064, -0.048, -0.032, -0.016, 0.0, 0.06, 0.12, 0.18, 0.24, 0.3, 0.36, 0.42, 0.48, 0.54, 0.6, 0.66, 0.72, 0.78, 0.84, 0.9, 0.96, 1.02, 1.08, 1.14, 1.2, 1.26, 1.32, 1.38, 1.438377, 1.44, 1.5, 1.56, 1.62, 1.68, 1.74, 1.8, 1.86, 1.92, 1.98, 2.04, 2.1, 2.16, 2.22, 2.28, 2.34, 2.4, 2.46, 2.52, 2.58, 2.64, 2.7, 2.76, 2.82, 2.88, 2.94, 3.0)
y <- c(-18500.53, -18499.829, -18499.273, -18498.831, -18498.482, -18498.21, -18497.995, -18497.832, -18497.71, -18497.621, -18497.56, -18497.519, -18497.498, -18497.49, -18497.4904, -18497.495, -18497.51, -18497.53, -18497.558, -18497.589, -18497.624, -18497.66, -18497.7, -18497.739, -18497.779, -18497.818, -18497.86, -18497.896, -18497.933, -18497.969, -18498.003, -18498.04, -18498.067, -18498.096, -18498.124, -18498.15, -18498.17, -18498.196, -18498.216, -18498.235, -18498.252, -18498.27, -18498.281, -18498.292, -18498.303, -18498.311, -18498.32, -18498.324, -18498.329, -18498.332, -18498.334, -18498.33, -18498.33, -18498.31, -18498.27, -18498.23, -18498.18, -18498.13, -18498.07, -18498.01, -18497.96, -18497.9, -18497.84, -18497.78, -18497.73, -18497.68, -18497.63, -18497.59, -18497.56, -18497.52, -18497.49, -18497.47, -18497.45, -18497.44, -18497.43, -18497.429487, -18497.43, -18497.43, -18497.44, -18497.45, -18497.47, -18497.5, -18497.52, -18497.56, -18497.6, -18497.64, -18497.69, -18497.74, -18497.8, -18497.86, -18497.93, -18498.0, -18498.07, -18498.15, -18498.23, -18498.32, -18498.41, -18498.5, -18498.6, -18498.7, -18498.81, -18498.92, -18499.03)

df <- data.table(x = x, y = y)
if (system.file(package = "ggplot2") != "") {
  g <- ggplot2::ggplot(df, ggplot2::aes(x = .data$x, y = .data$y)) +
    ggplot2::geom_line(color = "black", alpha = 1, linewidth = 1.5) +
    ggplot2::labs(x = "Linear Parameter Value", y = "Log-Likelihood") +
    ggplot2::ggtitle("Multi-Peak Curve")
} else {
  g <- message("ggplot2 wasn't detected. Please install to see the plot")
}
g

## ----eval=FALSE---------------------------------------------------------------
# fname <- "base_example.csv"
# df <- fread(fname)
# 
# model <- Cox(entry, exit, event) ~ loglinear(dose0, dose1, 0) + linear(dose0, 1)
# keep_constant <- c(0, 0, 0)
# a_n <- c(-1.493177, 5.020007, 1.438377)
# #
# control <- list(
#   ncores = 1, lr = 0.75, maxiter = 100, halfmax = 5,
#   verbose = 2
# )
# coxres <- CoxRun(model, df, a_n = a_n, control = control)
# 
# curve_control <- list(maxstep = 20, alpha = 0.005, para_number = 3, step_size = 0.5, bisect = TRUE)
# e <- LikelihoodBound(coxres, df, curve_control, control = control)