## ----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(-2, 0.01, 0.01)

poismodel <- get_form(Pois(time, status) ~ loglinear(CONST, temperature, voltage), df)$model

poisres <- PoisRun(poismodel, df,
  a_n = a_n, control = control
)

print(poisres)

## ----eval=TRUE----------------------------------------------------------------
res_p <- Residual(poisres, df, pearson = TRUE, deviance = TRUE)
print(res_p$Risk[1:5])
print(res_p$Raw_Residual[1:5])
print(res_p$Pearson_Residual[1:5])
print(res_p$Deviance_Residual[1:5])

a_n <- poisres$beta_0
res_d <- Residual(poismodel, df, a_n = a_n, pearson = TRUE, deviance = TRUE)
print(res_p$Risk[1:5])
print(res_p$Raw_Residual[1:5])
print(res_p$Pearson_Residual[1:5])
print(res_p$Deviance_Residual[1:5])