# =============================================================================
# MRMC AGREEMENT ANALYSIS
# Kappa, Gwet's AC1, prevalence sensitivity, and repeated-reader agreement
#
# Developed by: Planimeter Kft.
# Version:      1.0.0
# Date:         2026-09-22
# Copyright:    (c) 2026 Planimeter Kft.
#
# License:
#   This program is free software: you may redistribute it and/or modify it
#   under the terms of the GNU General Public License as published by the
#   Free Software Foundation, either version 3 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful, but
#   WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#
# SPDX-License-Identifier: GPL-3.0-or-later
#
# AI use statement:
#   Generative AI tools were used to support code development, review, and
#   validation. The analytical design, final code, interpretation, and
#   responsibility for the released version remain with Planimeter Kft.
#
# Purpose:
#   1) Read the fixed MRMC simulation database.
#   2) Reproduce baseline inter-reader agreement results.
#   3) Generate publication-ready figures.
#   4) Evaluate agreement statistics across prevalence levels.
#   5) Summarise intra-reader agreement.
#   6) Save derived tables and figures.
#
# Expected project structure:
#
#   project/
#   |-- mrmc_agreement_analysis_planimeter_v1.0.0.R
#   |-- data/
#   |   `-- mrmc_full_simulation.csv
#   |-- figures/        # created automatically
#   `-- output/         # created automatically
#
# Required input:
#   data/mrmc_full_simulation.csv
#
# Main outputs:
#   figures/01_cover_mrmc_agreement.png
#   figures/02_mrmc_design.png
#   figures/03_chance_agreement.png
#   figures/04_prevalence_sensitivity.png
#   figures/05_intra_reader_agreement.png
#   output/baseline_summary.csv
#   output/inter_reader_agreement_results.csv
#   output/prevalence_sensitivity_results.csv
#   output/reader_intra_agreement.csv
# =============================================================================


# ============================================================
# 0. PACKAGES
# ============================================================

required_packages <- c(
  "tidyverse",
  "scales"
)

missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]

if (length(missing_packages) > 0) {
  install.packages(missing_packages)
}

# irrCAC may not be available from CRAN in all R installations.
# If it is missing, install the maintained GitHub version.

if (!requireNamespace("irrCAC", quietly = TRUE)) {

  if (!requireNamespace("remotes", quietly = TRUE)) {
    install.packages("remotes")
  }

  remotes::install_github("kgwet/irrCAC")
}

library(tidyverse)
library(scales)
library(irrCAC)


# ============================================================
# 1. PATHS AND GLOBAL SETTINGS
# ============================================================

data_path <- file.path("data", "mrmc_full_simulation.csv")

figure_dir <- "figures"
output_dir <- "output"

dir.create(figure_dir, showWarnings = FALSE, recursive = TRUE)
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)

# Visual palette used for the generated figures.

COL_DARK   <- "#17212B"
COL_ACCENT <- "#2E6372"
COL_LIGHT  <- "#D8E2E6"
COL_MID    <- "#7B8790"
COL_PALE   <- "#F2F5F6"


theme_planimeter <- function(base_size = 14) {

  theme_minimal(
    base_size = base_size,
    base_family = "sans"
  ) +

    theme(
      plot.title = element_text(
        face = "bold",
        size = base_size * 1.45,
        colour = COL_DARK,
        margin = margin(b = 8)
      ),

      plot.subtitle = element_text(
        size = base_size,
        colour = COL_MID,
        lineheight = 1.15,
        margin = margin(b = 16)
      ),

      plot.caption = element_text(
        size = base_size * 0.75,
        colour = COL_MID,
        hjust = 0,
        margin = margin(t = 14)
      ),

      axis.title = element_text(
        face = "bold",
        colour = COL_DARK
      ),

      axis.text = element_text(
        colour = COL_DARK
      ),

      panel.grid.minor = element_blank(),

      panel.grid.major.x = element_blank(),

      legend.position = "top",

      legend.title = element_blank(),

      plot.margin = margin(
        20, 24, 20, 24
      )
    )
}


save_blog_plot <- function(
    p,
    filename,
    width = 12,
    height = 7
) {

  ggsave(
    filename = file.path(
      figure_dir,
      filename
    ),
    plot = p,
    width = width,
    height = height,
    dpi = 320,
    bg = "white"
  )
}


# ============================================================
# 2. READ FIXED SIMULATION DATABASE
# ============================================================

if (!file.exists(data_path)) {
  stop(
    paste0(
      "Input file not found: ",
      data_path,
      "\nPlace mrmc_full_simulation.csv in the relative data/ directory."
    )
  )
}

mrmc <- read_csv(
  data_path,
  show_col_types = FALSE
)


required_columns <- c(
  "case_id",
  "reader_id",
  "reading",
  "truth",
  "rating"
)

if (!all(required_columns %in% names(mrmc))) {
  stop(
    "The input database does not contain all required columns."
  )
}


# ============================================================
# 3. BASIC MRMC METRICS
# ============================================================

pairwise_inter_agreement <- function(dat) {

  dat %>%

    filter(
      reading == 1
    ) %>%

    group_by(
      case_id
    ) %>%

    summarise(

      n0 = sum(rating == 0),

      n1 = sum(rating == 1),

      n_readers = n(),

      pair_agreement =
        (
          choose(n0, 2) +
          choose(n1, 2)
        ) /
        choose(n_readers, 2),

      .groups = "drop"
    ) %>%

    summarise(
      agreement =
        mean(pair_agreement)
    ) %>%

    pull(agreement)
}


intra_agreement_overall <- function(dat) {

  dat %>%

    select(
      case_id,
      reader_id,
      reading,
      rating
    ) %>%

    pivot_wider(
      names_from = reading,
      values_from = rating,
      names_prefix = "reading_"
    ) %>%

    summarise(
      agreement =
        mean(
          reading_1 == reading_2
        )
    ) %>%

    pull(agreement)
}


baseline_summary <- tibble(

  metric = c(
    "Sensitivity",
    "Specificity",
    "Inter-reader raw agreement",
    "Intra-reader raw agreement"
  ),

  value = c(

    mean(
      mrmc$rating[
        mrmc$truth == 1
      ] == 1
    ),

    mean(
      mrmc$rating[
        mrmc$truth == 0
      ] == 0
    ),

    pairwise_inter_agreement(
      mrmc
    ),

    intra_agreement_overall(
      mrmc
    )
  )
)


write_csv(
  baseline_summary,
  file.path(
    output_dir,
    "baseline_summary.csv"
  )
)


# ============================================================
# 4. INTER-READER AGREEMENT:
#    READING 1 ONLY
# ============================================================

dat_inter <- mrmc %>%

  filter(
    reading == 1
  ) %>%

  select(
    case_id,
    reader_id,
    rating
  )


ratings_wide <- dat_inter %>%

  pivot_wider(
    names_from = reader_id,
    values_from = rating
  ) %>%

  arrange(
    case_id
  )


ratings <- ratings_wide %>%

  select(
    -case_id
  ) %>%

  as.data.frame()


stopifnot(
  !anyNA(ratings)
)


pa_result <- irrCAC::pa.coeff.raw(
  ratings
)$est


fleiss_result <- irrCAC::fleiss.kappa.raw(

  ratings,

  weights = "unweighted",

  categ.labels = c(0, 1)

)$est


conger_result <- irrCAC::conger.kappa.raw(

  ratings,

  weights = "unweighted",

  categ.labels = c(0, 1)

)$est


gwet_result <- irrCAC::gwet.ac1.raw(

  ratings,

  weights = "unweighted",

  categ.labels = c(0, 1)

)$est


kripp_result <- irrCAC::krippen.alpha.raw(

  ratings,

  weights = "unweighted",

  categ.labels = c(0, 1)

)$est


bp_result <- irrCAC::bp.coeff.raw(

  ratings,

  weights = "unweighted",

  categ.labels = c(0, 1)

)$est


extract_agreement_result <- function(
    x,
    method_name
) {

  tibble(

    method =
      method_name,

    observed_agreement =
      x$pa,

    chance_agreement =
      x$pe,

    estimate =
      x$coeff.val,

    SE =
      x$coeff.se,

    CI_95 =
      x$conf.int,

    p_value =
      x$p.value
  )
}


agreement_results <- bind_rows(

  extract_agreement_result(
    pa_result,
    "Percent agreement"
  ),

  extract_agreement_result(
    fleiss_result,
    "Fleiss' kappa"
  ),

  extract_agreement_result(
    conger_result,
    "Conger's kappa"
  ),

  extract_agreement_result(
    gwet_result,
    "Gwet's AC1"
  ),

  extract_agreement_result(
    kripp_result,
    "Krippendorff's alpha"
  ),

  extract_agreement_result(
    bp_result,
    "Brennan-Prediger"
  )
)


write_csv(
  agreement_results,
  file.path(
    output_dir,
    "inter_reader_agreement_results.csv"
  )
)


# ============================================================
# 5. VALUES USED IN THE BLOG COVER
# ============================================================

observed_agreement <- agreement_results %>%

  filter(
    method == "Percent agreement"
  ) %>%

  pull(estimate)


fleiss_kappa <- agreement_results %>%

  filter(
    method == "Fleiss' kappa"
  ) %>%

  pull(estimate)


gwet_ac1 <- agreement_results %>%

  filter(
    method == "Gwet's AC1"
  ) %>%

  pull(estimate)


# ============================================================
# FIGURE 1 / COVER IMAGE
# ============================================================

cover_plot <- ggplot() +

  annotate(
    "rect",
    xmin = 0,
    xmax = 10,
    ymin = 0,
    ymax = 10,
    fill = "white",
    colour = NA
  ) +

  annotate(
    "text",
    x = 5,
    y = 8.9,
    label = "HIGH AGREEMENT. SURPRISINGLY LOW KAPPA.",
    size = 6.5,
    fontface = "bold",
    colour = COL_DARK
  ) +

  annotate(
    "text",
    x = 2.0,
    y = 5.7,
    label = percent(
      observed_agreement,
      accuracy = 0.1
    ),
    size = 15,
    fontface = "bold",
    colour = COL_ACCENT
  ) +

  annotate(
    "text",
    x = 2.0,
    y = 4.5,
    label = "Observed agreement",
    size = 4.8,
    colour = COL_MID
  ) +

  annotate(
    "text",
    x = 5.0,
    y = 5.7,
    label = sprintf(
      "%.2f",
      fleiss_kappa
    ),
    size = 15,
    fontface = "bold",
    colour = COL_DARK
  ) +

  annotate(
    "text",
    x = 5.0,
    y = 4.5,
    label = "Fleiss' kappa",
    size = 4.8,
    colour = COL_MID
  ) +

  annotate(
    "text",
    x = 8.0,
    y = 5.7,
    label = sprintf(
      "%.2f",
      gwet_ac1
    ),
    size = 15,
    fontface = "bold",
    colour = COL_ACCENT
  ) +

  annotate(
    "text",
    x = 8.0,
    y = 4.5,
    label = "Gwet's AC1",
    size = 4.8,
    colour = COL_MID
  ) +

  annotate(
    "segment",
    x = 2.9,
    xend = 4.0,
    y = 5.4,
    yend = 5.4,
    linewidth = 1.0,
    colour = COL_LIGHT
  ) +

  annotate(
    "segment",
    x = 6.0,
    xend = 7.0,
    y = 5.4,
    yend = 5.4,
    linewidth = 1.0,
    colour = COL_LIGHT
  ) +

  annotate(
    "text",
    x = 5,
    y = 2.1,
    label =
      "Same cases. Same readers. Same observed agreement.\nDifferent definitions of chance agreement.",
    size = 5.1,
    lineheight = 1.15,
    colour = COL_DARK
  ) +

  coord_cartesian(
    xlim = c(0, 10),
    ylim = c(0, 10),
    clip = "off"
  ) +

  theme_void() +

  theme(
    plot.margin =
      margin(
        25, 25, 25, 25
      )
  )


save_blog_plot(
  cover_plot,
  "01_cover_mrmc_agreement.png",
  width = 13.5,
  height = 7.6
)


# ============================================================
# FIGURE 2 / MRMC DESIGN
# ============================================================

design_plot <- ggplot() +

  annotate(
    "rect",
    xmin = 0.5,
    xmax = 3.0,
    ymin = 4.0,
    ymax = 6.4,
    fill = COL_PALE,
    colour = COL_ACCENT,
    linewidth = 1
  ) +

  annotate(
    "text",
    x = 1.75,
    y = 5.6,
    label = "500 CASES",
    size = 7,
    fontface = "bold",
    colour = COL_DARK
  ) +

  annotate(
    "text",
    x = 1.75,
    y = 4.8,
    label = "Binary truth\n10% prevalence",
    size = 4.6,
    lineheight = 1.1,
    colour = COL_MID
  ) +

  annotate(
    "rect",
    xmin = 4.0,
    xmax = 6.5,
    ymin = 4.0,
    ymax = 6.4,
    fill = COL_PALE,
    colour = COL_ACCENT,
    linewidth = 1
  ) +

  annotate(
    "text",
    x = 5.25,
    y = 5.6,
    label = "5 READERS",
    size = 7,
    fontface = "bold",
    colour = COL_DARK
  ) +

  annotate(
    "text",
    x = 5.25,
    y = 4.8,
    label = "≈85% sensitivity\n≈95% specificity",
    size = 4.6,
    lineheight = 1.1,
    colour = COL_MID
  ) +

  annotate(
    "rect",
    xmin = 7.5,
    xmax = 10.0,
    ymin = 4.0,
    ymax = 6.4,
    fill = COL_PALE,
    colour = COL_ACCENT,
    linewidth = 1
  ) +

  annotate(
    "text",
    x = 8.75,
    y = 5.6,
    label = "2 SESSIONS",
    size = 7,
    fontface = "bold",
    colour = COL_DARK
  ) +

  annotate(
    "text",
    x = 8.75,
    y = 4.8,
    label = "Repeated reads\nfor each reader",
    size = 4.6,
    lineheight = 1.1,
    colour = COL_MID
  ) +

  annotate(
    "segment",
    x = 3.1,
    xend = 3.9,
    y = 5.2,
    yend = 5.2,
    arrow = grid::arrow(
      length = grid::unit(
        0.18,
        "inches"
      )
    ),
    linewidth = 0.9,
    colour = COL_MID
  ) +

  annotate(
    "segment",
    x = 6.6,
    xend = 7.4,
    y = 5.2,
    yend = 5.2,
    arrow = grid::arrow(
      length = grid::unit(
        0.18,
        "inches"
      )
    ),
    linewidth = 0.9,
    colour = COL_MID
  ) +

  annotate(
    "text",
    x = 5.25,
    y = 2.6,
    label =
      "Case difficulty + reader heterogeneity + stable reader × case effects\n+ small session-specific variation",
    size = 5.0,
    lineheight = 1.2,
    colour = COL_DARK
  ) +

  annotate(
    "text",
    x = 5.25,
    y = 1.4,
    label =
      "Fully crossed design: every reader evaluates every case in both sessions",
    size = 4.5,
    fontface = "italic",
    colour = COL_MID
  ) +

  coord_cartesian(
    xlim = c(0, 10.5),
    ylim = c(0.5, 7),
    clip = "off"
  ) +

  labs(
    title =
      "A reproducible multi-reader, multi-case simulation",
    subtitle =
      "The data-generating mechanism separates case, reader and repeat-reading variation."
  ) +

  theme_void() +

  theme(
    plot.title =
      element_text(
        face = "bold",
        size = 21,
        colour = COL_DARK
      ),

    plot.subtitle =
      element_text(
        size = 13.5,
        colour = COL_MID,
        margin = margin(
          b = 20
        )
      ),

    plot.margin =
      margin(
        25, 25, 25, 25
      )
  )


save_blog_plot(
  design_plot,
  "02_mrmc_design.png",
  width = 13,
  height = 7.3
)


# ============================================================
# FIGURE 3 / CHANCE-AGREEMENT COMPARISON
# ============================================================

chance_plot_data <- agreement_results %>%

  filter(
    method %in%
      c(
        "Fleiss' kappa",
        "Gwet's AC1"
      )
  ) %>%

  select(
    method,
    observed_agreement,
    chance_agreement,
    estimate
  ) %>%

  pivot_longer(
    cols =
      c(
        observed_agreement,
        chance_agreement
      ),
    names_to = "agreement_type",
    values_to = "value"
  ) %>%

  mutate(

    agreement_type =
      recode(
        agreement_type,
        observed_agreement =
          "Observed agreement",
        chance_agreement =
          "Expected by chance"
      ),

    method =
      factor(
        method,
        levels =
          c(
            "Fleiss' kappa",
            "Gwet's AC1"
          )
      )
  )


coefficient_labels <- agreement_results %>%

  filter(
    method %in%
      c(
        "Fleiss' kappa",
        "Gwet's AC1"
      )
  ) %>%

  transmute(
    method =
      factor(
        method,
        levels =
          c(
            "Fleiss' kappa",
            "Gwet's AC1"
          )
      ),
    label =
      paste0(
        if_else(
          method == "Fleiss' kappa",
          "κ = ",
          "AC1 = "
        ),
        sprintf(
          "%.2f",
          estimate
        )
      )
  )


chance_plot <- ggplot(
  chance_plot_data,
  aes(
    x = method,
    y = value,
    fill = agreement_type
  )
) +

  geom_col(
    position =
      position_dodge(
        width = 0.74
      ),
    width = 0.62
  ) +

  geom_text(
    aes(
      label =
        percent(
          value,
          accuracy = 0.1
        )
    ),
    position =
      position_dodge(
        width = 0.74
      ),
    vjust = -0.45,
    size = 4.5,
    fontface = "bold",
    colour = COL_DARK
  ) +

  geom_text(
    data =
      coefficient_labels,
    aes(
      x = method,
      y = 1.02,
      label = label
    ),
    inherit.aes = FALSE,
    size = 5.4,
    fontface = "bold",
    colour = COL_ACCENT
  ) +

  scale_fill_manual(
    values =
      c(
        "Observed agreement" =
          COL_ACCENT,
        "Expected by chance" =
          COL_LIGHT
      )
  ) +

  scale_y_continuous(
    labels =
      percent_format(
        accuracy = 1
      ),
    limits =
      c(
        0,
        1.08
      ),
    expand =
      expansion(
        mult =
          c(
            0,
            0
          )
      )
  ) +

  labs(
    title =
      "The data do not change. The definition of chance agreement does.",
    subtitle =
      "Both coefficients start from the same 88% observed agreement, but use very different chance corrections.",
    x = NULL,
    y = "Agreement",
    caption =
      "Reading 1 only; 500 cases and 5 readers."
  ) +

  theme_planimeter() +

  theme(
    legend.position = "top"
  )


save_blog_plot(
  chance_plot,
  "03_chance_agreement.png",
  width = 11,
  height = 7
)


# ============================================================
# 6. PREVALENCE SENSITIVITY SIMULATION
#
# Same reader model, different prevalence.
#
# The same seed is used across prevalence values within each
# Monte Carlo replicate to keep the latent case and reader
# structure aligned across prevalence scenarios.
# ============================================================


simulate_mrmc <- function(
    n_cases = 500,
    n_readers = 5,
    n_readings = 2,

    prevalence = 0.10,

    mean_sensitivity = 0.87,
    mean_specificity = 0.96,

    reader_sd = 0.25,
    case_difficulty_sd = 0.50,
    reader_case_sd = 0.45,

    reader_occasion_sd = 0.03,
    occasion_sd = 0.10,

    rho_intra = 0.65,

    seed = 20260812
) {

  set.seed(seed)


  cases <- tibble(

    case_id =
      sprintf(
        "C%04d",
        seq_len(n_cases)
      ),

    truth_u =
      runif(
        n_cases
      ),

    case_difficulty =
      rnorm(
        n_cases,
        mean = 0,
        sd = case_difficulty_sd
      )
  ) %>%

    mutate(
      truth =
        as.integer(
          truth_u < prevalence
        )
    )


  readers <- tibble(

    reader_id =
      sprintf(
        "R%02d",
        seq_len(n_readers)
      ),

    reader_sensitivity =
      plogis(
        qlogis(
          mean_sensitivity
        ) +
          rnorm(
            n_readers,
            0,
            reader_sd
          )
      ),

    reader_specificity =
      plogis(
        qlogis(
          mean_specificity
        ) +
          rnorm(
            n_readers,
            0,
            reader_sd
          )
      )
  )


  reader_case <- crossing(

    case_id =
      cases$case_id,

    reader_id =
      readers$reader_id

  ) %>%

    mutate(

      reader_case_effect =
        rnorm(
          n(),
          0,
          reader_case_sd
        ),

      z_reading1 =
        rnorm(
          n()
        ),

      z_reading2 =
        rho_intra *
        z_reading1 +
        sqrt(
          1 -
            rho_intra^2
        ) *
        rnorm(
          n()
        ),

      threshold_reading1 =
        pnorm(
          z_reading1
        ),

      threshold_reading2 =
        pnorm(
          z_reading2
        )
    )


  reader_occasion <- crossing(

    reader_id =
      readers$reader_id,

    reading =
      seq_len(
        n_readings
      )

  ) %>%

    mutate(
      reader_occasion_effect =
        rnorm(
          n(),
          0,
          reader_occasion_sd
        )
    )


  crossing(

    case_id =
      cases$case_id,

    reader_id =
      readers$reader_id,

    reading =
      seq_len(
        n_readings
      )

  ) %>%

    left_join(
      cases,
      by = "case_id"
    ) %>%

    left_join(
      readers,
      by = "reader_id"
    ) %>%

    left_join(
      reader_case,
      by =
        c(
          "case_id",
          "reader_id"
        )
    ) %>%

    left_join(
      reader_occasion,
      by =
        c(
          "reader_id",
          "reading"
        )
    ) %>%

    mutate(

      occasion_noise =
        rnorm(
          n(),
          0,
          occasion_sd
        ),

      linear_predictor =
        if_else(

          truth == 1,

          qlogis(
            reader_sensitivity
          ) -
            case_difficulty +
            reader_case_effect +
            reader_occasion_effect +
            occasion_noise,

          qlogis(
            1 -
              reader_specificity
          ) +
            case_difficulty +
            reader_case_effect +
            reader_occasion_effect +
            occasion_noise
        ),

      prob_positive =
        plogis(
          linear_predictor
        ),

      decision_threshold =
        if_else(
          reading == 1,
          threshold_reading1,
          threshold_reading2
        ),

      rating =
        as.integer(
          decision_threshold <
            prob_positive
        )
    )
}


# Closed-form calculations for a complete binary multiple-rater design.
# These implement the chance-correction definitions used for
# Fleiss' kappa and Gwet's AC1 in the prevalence analysis.

fast_binary_agreement <- function(dat) {

  wide <- dat %>%

    filter(
      reading == 1
    ) %>%

    select(
      case_id,
      reader_id,
      rating
    ) %>%

    pivot_wider(
      names_from = reader_id,
      values_from = rating
    ) %>%

    arrange(
      case_id
    )


  mat <- as.matrix(
    wide %>%
      select(
        -case_id
      )
  )


  n_readers <-
    ncol(
      mat
    )


  n1 <-
    rowSums(
      mat == 1
    )


  n0 <-
    n_readers -
    n1


  pa <-
    mean(
      (
        choose(
          n1,
          2
        ) +
          choose(
            n0,
            2
          )
      ) /
        choose(
          n_readers,
          2
        )
    )


  p_positive <-
    mean(
      mat == 1
    )


  p_negative <-
    1 -
    p_positive


  pe_fleiss <-
    p_positive^2 +
    p_negative^2


  pe_gwet <-
    2 *
    p_positive *
    p_negative


  kappa <-
    (
      pa -
      pe_fleiss
    ) /
    (
      1 -
      pe_fleiss
    )


  ac1 <-
    (
      pa -
      pe_gwet
    ) /
    (
      1 -
      pe_gwet
    )


  sensitivity <-
    mean(
      dat$rating[
        dat$truth == 1
      ] == 1
    )


  specificity <-
    mean(
      dat$rating[
        dat$truth == 0
      ] == 0
    )


  tibble(
    actual_prevalence =
      mean(
        dat$truth
      ),
    sensitivity =
      sensitivity,
    specificity =
      specificity,
    observed_agreement =
      pa,
    fleiss_kappa =
      kappa,
    gwet_ac1 =
      ac1,
    positive_rating_rate =
      p_positive
  )
}


prevalence_grid <- c(
  0.01,
  0.02,
  0.05,
  0.10,
  0.15,
  0.20,
  0.30,
  0.40,
  0.50
)


B_prevalence <- 30


prevalence_results <- map_dfr(

  seq_len(
    B_prevalence
  ),

  function(b) {

    map_dfr(

      prevalence_grid,

      function(p) {

        # The seed is held constant across prevalence levels within
        # each replicate to align latent case and reader components.

        dat <-
          simulate_mrmc(

            prevalence = p,

            seed =
              200000 +
              b
          )


        fast_binary_agreement(
          dat
        ) %>%

          mutate(
            target_prevalence = p,
            simulation = b,
            .before = 1
          )
      }
    )
  }
)


prevalence_summary <- prevalence_results %>%

  group_by(
    target_prevalence
  ) %>%

  summarise(

    actual_prevalence =
      mean(
        actual_prevalence
      ),

    sensitivity =
      mean(
        sensitivity,
        na.rm = TRUE
      ),

    specificity =
      mean(
        specificity,
        na.rm = TRUE
      ),

    observed_agreement =
      mean(
        observed_agreement
      ),

    fleiss_kappa =
      mean(
        fleiss_kappa
      ),

    gwet_ac1 =
      mean(
        gwet_ac1
      ),

    .groups = "drop"
  )


write_csv(
  prevalence_results,
  file.path(
    output_dir,
    "prevalence_sensitivity_results_all_simulations.csv"
  )
)


write_csv(
  prevalence_summary,
  file.path(
    output_dir,
    "prevalence_sensitivity_results.csv"
  )
)


# ============================================================
# FIGURE 4 / PREVALENCE SENSITIVITY
# ============================================================

prevalence_plot_data <- prevalence_summary %>%

  select(
    target_prevalence,
    observed_agreement,
    fleiss_kappa,
    gwet_ac1
  ) %>%

  pivot_longer(

    cols =
      c(
        observed_agreement,
        fleiss_kappa,
        gwet_ac1
      ),

    names_to =
      "method",

    values_to =
      "estimate"
  ) %>%

  mutate(

    method =
      recode(
        method,
        observed_agreement =
          "Observed agreement",
        fleiss_kappa =
          "Fleiss' kappa",
        gwet_ac1 =
          "Gwet's AC1"
      )
  )


prevalence_plot <- ggplot(

  prevalence_plot_data,

  aes(
    x = target_prevalence,
    y = estimate,
    group = method,
    linetype = method,
    colour = method
  )

) +

  geom_line(
    linewidth = 1.25
  ) +

  geom_point(
    size = 2.7
  ) +

  scale_colour_manual(
    values =
      c(
        "Observed agreement" =
          COL_MID,
        "Fleiss' kappa" =
          COL_DARK,
        "Gwet's AC1" =
          COL_ACCENT
      )
  ) +

  scale_linetype_manual(
    values =
      c(
        "Observed agreement" =
          "dotted",
        "Fleiss' kappa" =
          "solid",
        "Gwet's AC1" =
          "longdash"
      )
  ) +

  scale_x_continuous(
    labels =
      percent_format(
        accuracy = 1
      ),
    breaks =
      prevalence_grid
  ) +

  scale_y_continuous(
    limits =
      c(
        0,
        1
      ),
    breaks =
      seq(
        0,
        1,
        0.2
      )
  ) +

  labs(
    title =
      "The same reader model can produce very different kappa values",
    subtitle =
      "Only prevalence changes; the underlying reader-performance model is held constant.",
    x =
      "True prevalence",
    y =
      "Agreement coefficient",
    caption =
      paste0(
        "Monte Carlo mean across ",
        B_prevalence,
        " simulations per prevalence level."
      )
  ) +

  theme_planimeter() +

  theme(
    legend.position = "top"
  )


save_blog_plot(
  prevalence_plot,
  "04_prevalence_sensitivity.png",
  width = 12,
  height = 7
)


# ============================================================
# 7. INTRA-READER AGREEMENT
# ============================================================

intra_by_reader <- mrmc %>%

  select(
    case_id,
    reader_id,
    reading,
    rating
  ) %>%

  pivot_wider(
    names_from = reading,
    values_from = rating,
    names_prefix = "reading_"
  ) %>%

  group_by(
    reader_id
  ) %>%

  summarise(

    intra_agreement =
      mean(
        reading_1 ==
          reading_2
      ),

    .groups = "drop"
  )


overall_inter <-
  pairwise_inter_agreement(
    mrmc
  )


overall_intra <-
  mean(
    intra_by_reader$intra_agreement
  )


write_csv(
  intra_by_reader,
  file.path(
    output_dir,
    "reader_intra_agreement.csv"
  )
)


# ============================================================
# FIGURE 5 / INTRA-READER REPRODUCIBILITY
# ============================================================

intra_plot <- ggplot(

  intra_by_reader,

  aes(
    x = reader_id,
    y = intra_agreement
  )

) +

  geom_col(
    width = 0.62,
    fill = COL_ACCENT
  ) +

  geom_text(
    aes(
      label =
        percent(
          intra_agreement,
          accuracy = 0.1
        )
    ),
    vjust = -0.45,
    size = 4.1,
    fontface = "bold",
    colour = COL_DARK
  ) +

  geom_hline(
    yintercept =
      overall_inter,
    linewidth = 0.9,
    linetype = "dashed",
    colour = COL_DARK
  ) +

  annotate(
    "text",
    x = Inf,
    y =
      overall_inter +
      0.008,
    label =
      paste0(
        "Overall inter-reader agreement = ",
        percent(
          overall_inter,
          accuracy = 0.1
        )
      ),
    hjust = 1.05,
    vjust = 0,
    size = 3.9,
    colour = COL_DARK
  ) +

  scale_y_continuous(
    labels =
      percent_format(
        accuracy = 1
      ),
    limits =
      c(
        0.75,
        1.00
      ),
    breaks =
      seq(
        0.75,
        1.00,
        0.05
      )
  ) +

  labs(
    title =
      "Readers can be more consistent with themselves than with each other",
    subtitle =
      paste0(
        "Overall intra-reader agreement = ",
        percent(
          overall_intra,
          accuracy = 0.1
        ),
        ". The dashed line shows overall inter-reader agreement."
      ),
    x =
      NULL,
    y =
      "Intra-reader agreement",
    caption =
      "Agreement between reading 1 and reading 2 for the same reader and case."
  ) +

  theme_planimeter()


save_blog_plot(
  intra_plot,
  "05_intra_reader_agreement.png",
  width = 11,
  height = 7
)


# ============================================================
# 8. CONSOLE SUMMARY
# ============================================================

cat("\n")
cat("============================================================\n")
cat("MRMC BLOG ANALYSIS COMPLETE\n")
cat("============================================================\n\n")

cat(
  "Observed inter-reader agreement: ",
  sprintf(
    "%.3f",
    observed_agreement
  ),
  "\n",
  sep = ""
)

cat(
  "Fleiss' kappa: ",
  sprintf(
    "%.3f",
    fleiss_kappa
  ),
  "\n",
  sep = ""
)

cat(
  "Gwet's AC1: ",
  sprintf(
    "%.3f",
    gwet_ac1
  ),
  "\n",
  sep = ""
)

cat(
  "Overall intra-reader agreement: ",
  sprintf(
    "%.3f",
    overall_intra
  ),
  "\n\n",
  sep = ""
)

cat(
  "Figures saved to: ",
  figure_dir,
  "/\n",
  sep = ""
)

cat(
  "Derived tables saved to: ",
  output_dir,
  "/\n",
  sep = ""
)

cat("\n")
cat("============================================================\n")