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

## ----data---------------------------------------------------------------------
library(oda)

# Cross-classification: rows = attacks (0-7), cols = treatment arm.
#          T1 (0)  T2 (1)
#  0 att:    13       5
#  1 att:     9      13
#  2 att:     4       6
#  3 att:     2       1
#  4 att:     1       2
#  5 att:     1       3
#  6 att:     3       3
#  7 att:     0       1

treatment <- c(
  rep(0L, 13), rep(1L,  5),   # attacks = 0
  rep(0L,  9), rep(1L, 13),   # attacks = 1
  rep(0L,  4), rep(1L,  6),   # attacks = 2
  rep(0L,  2), rep(1L,  1),   # attacks = 3
  rep(0L,  1), rep(1L,  2),   # attacks = 4
  rep(0L,  1), rep(1L,  3),   # attacks = 5
  rep(0L,  3), rep(1L,  3),   # attacks = 6
  rep(0L,  0), rep(1L,  1)    # attacks = 7
)
attacks <- c(
  rep(0L, 18), rep(1L, 22), rep(2L, 10),
  rep(3L,  3), rep(4L,  3), rep(5L,  4),
  rep(6L,  6), rep(7L,  1)
)

table(attacks, treatment,
      dnn = c("Migraine Attacks (0-7)", "Treatment (0=T1, 1=T2)"))

## ----fit-canonical, eval=FALSE------------------------------------------------
# # Canonical reference run (mc_iter = 25000L; not evaluated in CRAN vignette)
# fit <- oda_fit(
#   x         = attacks,
#   y         = treatment,
#   attr_type = "ordered",
#   mc_iter   = 25000L,
#   loo       = "on"
# )

## ----fit----------------------------------------------------------------------
# CRAN-safe run: mc_iter = 500L for vignette rendering speed.
# Training rule, ESS, and confusion matrix are identical to the canonical run.
fit <- oda_fit(
  x         = attacks,
  y         = treatment,
  attr_type = "ordered",
  mc_iter   = 500L,
  mc_seed   = 42L,
  loo       = "on"
)

## ----print-fit----------------------------------------------------------------
print(fit)

## ----confusion----------------------------------------------------------------
# Confusion matrix: actual treatment (rows) x predicted treatment (cols)
conf_mat <- matrix(
  c(fit$confusion$TN, fit$confusion$FP,
    fit$confusion$FN, fit$confusion$TP),
  nrow = 2L, byrow = TRUE,
  dimnames = list(Actual    = c("T1(0)", "T2(1)"),
                  Predicted = c("T1(0)", "T2(1)"))
)
print(conf_mat)

## ----metrics------------------------------------------------------------------
summary(fit)

## ----pv-----------------------------------------------------------------------
# Predictive value: accuracy when the model makes a prediction into each class
pv_t1 <- fit$confusion$TN / (fit$confusion$TN + fit$confusion$FN)
pv_t2 <- fit$confusion$TP / (fit$confusion$TP + fit$confusion$FP)
cat("PV Treatment 1 (0):", round(pv_t1 * 100, 1), "%\n")
cat("PV Treatment 2 (1):", round(pv_t2 * 100, 1), "%\n")