plots.Rmd

---
title: "Plots"
date: "`r format(Sys.time(), '%B %d, %Y')`"
author: "Richard Wen <rrwen.dev@gmail.com>"
---

```{r, echo=FALSE}
knitr::opts_chunk$set(
  message=FALSE,
  warning=FALSE
)
options(readr.show_col_types = FALSE)
```

```{r, echo=FALSE, eval=FALSE}
#install.packages("openVA")
#install.packages("tidyverse")
#install.packages("dplyr")
#install.packages("ggrepel")
#install.packages("patchwork")
#install.packages("glue")
#install.packages("kableExtra)
#install.packages("boot")
#install.pacakges("scales")
```

# Libraries

```{r}
library(tidyverse)
library(dplyr)
library(ggrepel)
library(ggpattern)
library(patchwork)
library(glue)
library(kableExtra)
library(boot)
library(scales)
```

```{r, echo=FALSE}
cat("R Version\n---")
R.version

cat("\nPackage Versions\n---")
cat(paste0("\ntidyverse ", packageVersion("tidyverse")))
```

# Functions

## Data Functions

```{r}

# Order textual age ranges
order_arange <- function(ages) {
    out <- ages %>%
        data.frame %>%
        rename("age_range" = ".") %>%
        separate( # sep into range and unit
            age_range,
            into = c("age_range", "unit"),
            sep = " "
        ) %>%
        separate( # sep into min and max age
            age_range,
            into = c("age_min", "age_max"),
            sep = "-"
        ) %>%
        mutate( # format into correct data types
            age_min = as.integer(age_min),
            age_max = as.integer(age_max),
            unit = str_to_lower(unit) %>%
                factor(., levels = c(
                    "hours",
                    "days",
                    "weeks",
                    "months",
                    "years"
                ))
        ) %>%
        arrange(unit, age_min, age_max) %>% # sort ages
        mutate( # create ordered ages
            age_range = paste0(age_min, "-", age_max, " ", unit)
        ) %>%
        pull(age_range)
    return(out)
}

```

## Metric Functions

Functions for PCCC and CSMF Accuracy.

```{r}

# Calc PCCC
calc_pccc <- function(actual, pred, k = 1, N = NULL, assume_zero = T) {
    
    # Calc N num of causes if not known
    N <- if (is.null(N)) length(unique(na.omit(c(actual, pred)))) else N
    
    # Calc frac of deaths in top k causes
    TP <- actual == pred
    TP[is.na(TP) | is.null(TP)] <- FALSE # for no preds
    C <- sum(TP) / length(actual)
    
    # Calc PCCC
    if (!all(TP)) {
        out <- (C - (k/N)) / (1 - (k/N))
    } else { # avoid division by 0 when max pccc
        out <- 1
    }
    
    # Scale negative numbers to 0 for easier plotting
    out <- if (assume_zero & out < 0) 0 else out
    return(out)
}

# Calc CSMF accuracy
calc_csmf_acc <- function(actual, pred) {
    
    # Get all unique causes
    causes <- unique(c(actual, pred))
    
    # Get csmfs
    cases <- length(actual)
    csmf_true <- table(actual) / cases
    csmf_pred <- table(pred) / cases
    
    # Correct for missing causes in either actual or pred
    csmf_true <- vapply(
        causes,
        function(x) if (x %in% names(csmf_true)) csmf_true[x] else 0,
        FUN.VALUE = numeric(1)
    )
    csmf_pred <- vapply(
        causes,
        function(x) if (x %in% names(csmf_pred)) csmf_pred[x] else 0,
        FUN.VALUE = numeric(1)
    )
    
    # Calc csmf max error
    csmf_max_error <- 2 * (1 - min(csmf_true))
    
    # Calc csmf acc
    out <- 1 - (sum(abs(csmf_true - csmf_pred)) / csmf_max_error)
    return(out)
}

# Calc bootstrapped cis for a metric
calc_ci <- function(
        x,
        f = calc_pccc,
        n_boot = 1000,
        conf_level = 0.95,
        col_actual = "cod",
        col_pred = "model_cghr10",
        ...
    ) {
    
    # Function to compute metric on a bootstrap sample
    f_boot <- function(df, i) {
        f(df[[col_actual]][i], df[[col_pred]][i], ...)
    }
    
    # Run bootstrap using the 'boot' package
    boot_res <- boot(data = x, statistic = f_boot, R = n_boot)
    
    # Compute percentile-based confidence interval
    ci <- boot.ci(boot_res, type = "perc", conf = conf_level)$percent[4:5]
    
    # Return named vector with metric and CI
    out <- c(ci_mean = mean(boot_res$t), ci_low = ci[1], ci_hi = ci[2])
    return(out)
}

```

## Plot Functions

Functions for generating plots.

```{r}

# Plot performance of models using lines and points
plot_perf <- function(
        df_long,
        f = calc_pccc,
        group = "cod_name",
        color = "model_name",
        text_scale = 1,
        title = NULL,
        title_size = 10,
        x_title = "Partial Chance Corrected Concordance (PCCC)",
        y_title = NULL,
        text_mid_size = 2.25 * text_scale,
        text_mid_vjust = 2.5,
        text_start_size = 2.75 * text_scale,
        text_end_size = 3.5 * text_scale,
        text_hjust = 0.2,
        add_min_points = F,
        add_max_points = F,
        mid_points_size = 4,
        min_points_size = 4,
        max_points_size = 4,
        x_limits = c(0, 1.15),
        x_breaks = seq(0, 1, by = 0.1),
        x_expand = expansion(add = 0.1),
        x_color_labels = c(
            "G3" = "GPT-3.5 (G3)",
            "G4" = "GPT-4 (G4)",
            "G5" = "GPT-5 (G5)",
            "ISV" = "InSilicoVA (ISV)",
            "I5" = "InterVA-5 (I5)"
        ),
        x_shapes = 21:26,
        x_title_size = 9 * text_scale,
        y_title_size = 9 * text_scale,
        x_text_size = 8 * text_scale,
        y_text_size = 9 * text_scale,
        sort_group = T,
        legend_position = "bottom",
        legend_title = element_blank(),
        legend_margin = margin(c(-8, -2, -2, -50)),
        legend_add = c("R = Range"),
        ...
    ) {
    
    #  Group data
    df_plot <- df_long %>%
        group_by(.data[[group]], model_name, .drop = F)
    
    # Calc metric
    if (identical(f, calc_pccc)) {
        df_plot <- df_plot %>%
            summarize(
                metric = f(cod, model_cghr10, N = unique(ncod)[1], ...)
            )
    } else {
        df_plot <- df_plot %>%
            summarize(
                metric = f(cod, model_cghr10, ...)
            )
    }
        
    # Add metric value to model name
    df_plot <- df_plot %>%
        mutate( 
            model_label = paste0(model_name, ": ", round(metric, 2))
        )
    
    # Add min and max values for each group
    df_ends <- df_plot %>%
        group_by(.data[[group]]) %>%
        summarize(
            min_metric = min(metric),
            max_metric = max(metric),
            model_min = model_name[which.min(metric)],
            model_max = model_name[which.max(metric)],
            group_range = diff(range(metric))
        ) %>%
        arrange(max_metric)
    
    # Sort group by max model metric
    if (sort_group) {
        df_plot[[group]] <- factor(df_plot[[group]], levels = df_ends[[group]])
        df_ends[[group]] <- factor(df_ends[[group]], levels = df_ends[[group]])
    }
    
    # Filter out the mid models
    df_mids <- df_plot %>%
        left_join(df_ends, by = group) %>%
        filter(model_name != model_max & model_name != model_min)
    
    # Add blank legend items
    if (!is.null(legend_add)) {
        
        # Add to manual color labels
        for (i in legend_add) {
            x_color_labels[i] <- i
        }
        
        # Add to factor levels
        df_plot[[color]] <- factor(
            df_plot[[color]],
            levels = c(levels(df_plot[[color]]), legend_add)
        )
    }
    
    # Get custom colors for models
    default_colors <- scales::hue_pal()(length(unique(df_plot[[color]])))
    default_colors <- c(default_colors, rep(NA, length(legend_add)))
    names(default_colors) <- c(unique(as.character(df_plot[[color]])), legend_add)
    
    # Plot metric for each model
    out <- ggplot(df_plot, aes(x = metric, y = .data[[group]])) +
        geom_segment( # min max lines
            data = df_ends,
            aes(
                x = min_metric,
                xend = max_metric,
                y = .data[[group]],
                yend = .data[[group]]
            ),
            linewidth = 0.75
        ) +
        geom_point( # points
            color = "black",
            aes(
                shape = .data[[color]],
                fill = .data[[color]]
            ),
            size = mid_points_size
        ) +
        geom_text( # min model text
            data = df_ends,
            aes(
                x = min_metric,
                y = .data[[group]],
                label = paste0(model_min, ": ", round(min_metric, 2))
            ),
            hjust = 1 + text_hjust,
            vjust = 0.5,
            size = text_start_size
        ) +
        geom_text( # max model text
            data = df_ends,
            aes(
                x = max_metric,
                y = .data[[group]],
                label = paste0(
                    model_max, ": ", round(max_metric, 2),
                    " (R=", round(group_range, 2), ")"
                )
            ),
            hjust = 0 - text_hjust,
            vjust = 0.5,
            size = text_end_size,
            fontface = "bold"
        ) +
        geom_text_repel( # mid model text
            data = df_mids,
            aes(label = model_label),
            color = "gray60",
            vjust = text_mid_vjust,
            max.overlaps = Inf,
            size = text_mid_size
        ) +
        scale_fill_manual(
            values = default_colors,
            labels = x_color_labels,
            drop = ifelse(is.null(legend_add), T, F)
        ) +
        scale_shape_manual(
            values = x_shapes,
            labels = x_color_labels,
            drop = ifelse(is.null(legend_add), T, F)
        ) +
        scale_x_continuous(
            limits = x_limits,
            breaks = x_breaks,
            expand = x_expand
        ) +
        labs(
            title = title,
            x = x_title,
            y = y_title
        ) +
        theme_minimal() +
        theme(
            legend.position = legend_position,
            legend.title = legend_title,
            legend.margin = legend_margin,
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            plot.title = element_text(size = title_size),
            axis.title.x = element_text(size = x_title_size),
            axis.title.y = element_text(size = y_title_size),
            axis.text.x = element_text(size = x_text_size),
            axis.text.y = element_text(size = y_text_size)
        )
    
    # Add min and max points
    if (add_min_points) {
        out <- out + geom_point( # min points
            data = df_ends,
            aes(x = min_metric),
            size = 2
        )
    }
    if (add_max_points) {
        out <- out + geom_point( # max points
            data = df_ends,
            aes(x = max_metric),
            size = 2
        )
    }
    return(out)
}

# Plot performance of models using grouped bars
plot_perf_bar <- function(
        df_long,
        f = calc_pccc,
        group = "cod_name",
        color = "model_name",
        text_scale = 1,
        title = NULL,
        title_size = 10,
        x_title = "Partial Chance Corrected Concordance (PCCC)",
        y_title = NULL,
        x_hpos = -0.47 * text_scale,
        x_limits = c(0, 1),
        x_breaks = seq(0, 1, by = 0.1),
        x_expand = expansion(add = ((x_hpos * -1) + 0.03)),
        x_shapes = 21:26,
        x_title_size = 8 * text_scale,
        y_title_size = 8 * text_scale,
        x_text_size = 7 * text_scale,
        x_color_size = 2.5 * text_scale,
        y_text_size = 8 * text_scale,
        legend_position = "bottom",
        legend_title = element_blank(),
        legend_margin = margin(c(-8, -2, -2, -50)),
        legend_key_size = 0.5,
        legend_text_size = 7 * text_scale,
        legend_spacing_x = unit(58, "pt"),
        n_boot = 1000,
        error_bar = T,
        error_bar_width = 0.05,
        error_bar_size = 0.3,
        error_bar_label = "95% Confidence Interval (Bootstrap, 1000 Resamples)"
    ) {
    
    #  Group data
    df_plot <- df_long %>%
        group_by(.data[[group]], model_name, .drop = F)
    
    # Calc metric
    if (identical(f, calc_pccc)) {
        df_plot <- df_plot %>%
            group_map(~{
                
                # Calc metric
                out <- c(
                    metric = f(.x$cod, .x$model_cghr10, N = unique(.x$ncod)[1])
                )
                
                # Include ci if error bar
                if (error_bar) {
                    out <- c(.y, out, calc_ci(.x, n_boot = n_boot))
                }
                return(out)
                
            }, .keep = T) %>%
            bind_rows
        
    } else {
        df_plot <- df_plot %>%
            group_map(~{
                
                # Calc metric
                out <- c(
                    metric = f(.x$cod, .x$model_cghr10)
                )
                
                # Include ci if error bar
                if (error_bar) {
                    out <- c(.y, out, calc_ci(.x, n_boot = n_boot))
                }
                return(out)
                
            }, .keep = T) %>%
            bind_rows
    }
    
    # Get custom colors for models
    default_colors <- scales::hue_pal()(length(unique(df_plot[[color]])))
    names(default_colors) <- c(unique(as.character(df_plot[[color]])))
    
    # Format plot data
    df_plot <- df_plot %>%
        mutate(
            x = metric,
            y = .data[[group]],
            color = .data[[color]]
        )

    # Sort bars within groups (ascending)
    df_pos <- df_plot %>%
        group_by(y) %>%
        arrange(x, .by_group = TRUE) %>%  
        mutate(
            rank_in_group = row_number(),
            n_colors = n()
        ) %>%
        ungroup()
    
    # Reverse sort groups by highest-performing bar
    group_levels_order <- df_pos %>%
        group_by(y) %>%
        summarize(max_x = max(x)) %>%
        arrange(max_x) %>%
        pull(y)
    df_pos <- df_pos %>%
        mutate(
            y = factor(y, levels = group_levels_order),
            group_center = as.numeric(y)
        )
    
    # Bar width and offset
    group_width <- 0.8
    max_items <- max(df_pos$n_colors)
    bar_width <- group_width / max_items
    df_pos <- df_pos %>%
        mutate(
            offset = (rank_in_group - (n_colors + 1)/2) * bar_width,
            x_pos = group_center + offset,
            label_text = paste0(color, ": ", format(round(x, 2), nsmall = 2)),
            label_x = x_pos
        )
    
    # Add ci to model label
    if (error_bar) {
        df_pos <- df_pos %>%
            mutate(
                label_text = paste0(label_text, " (", round(ci_low, 2), ",", round(ci_hi, 2), ")")
            )
    }
    
    # Identify top-performing model in each group
    df_pos <- df_pos %>%
        group_by(y) %>%
        mutate(
            top_model = x == max(x),
            fontface = ifelse(top_model, "bold", "plain"),
            alpha_val = ifelse(top_model, 1, 0.35),       # transparency for non-top bars
            outline_color = ifelse(top_model, "black", NA) # black outline for top bars
        ) %>%
        ungroup()
    
    # Compute range and modify group labels using square brackets
    group_labels <- df_pos %>%
        group_by(y) %>%
        summarize(range_x = max(x) - min(x)) %>%
        mutate(label = paste0(as.character(y), "\n[Range=", format(round(range_x, 2), nsmall = 2), "]")) %>%
        arrange(match(y, levels(df_pos$y)))
    
    # Replace factor levels with new labels including range
    levels(df_pos$y) <- group_labels$label

    # Plot metric for each model
    df_pos$ci_label <- error_bar_label # dummy label for legend
    out <- ggplot(df_pos, aes(x = x_pos, y = x, fill = color)) +
        geom_col(
            width = bar_width * 0.9,
            position = "identity",
            color = df_pos$outline_color,
            alpha = df_pos$alpha_val
        ) +
        geom_text(
            aes(
                x = label_x,
                label = label_text,
                fontface = fontface
            ),
            y = x_hpos,
            hjust = 0,
            vjust = 0.5,
            size = x_color_size
        ) +
        coord_flip() +
        scale_fill_manual(
            values = default_colors
        ) +
        scale_x_continuous(
            breaks = seq_along(levels(df_pos$y)),
            labels = levels(df_pos$y),
            expand = c(0, 0)
        ) +
        scale_y_continuous(
            expand = x_expand,
            breaks = x_breaks,
            limits = x_limits
        ) +
        labs(
            title = title,
            x = y_title,
            y = x_title
        ) +
        theme_minimal() +
        theme(
            legend.position = legend_position,
            legend.title = legend_title,
            legend.margin = legend_margin,
            legend.key.size = unit(legend_key_size, "lines"),
            legend.text = element_text(size = legend_text_size),
            legend.spacing.x = legend_spacing_x,
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            plot.title = element_text(size = title_size),
            axis.title.x = element_text(size = x_title_size),
            axis.title.y = element_text(size = y_title_size),
            axis.text.x = element_text(size = x_text_size),
            axis.text.y = element_text(size = y_text_size)
        )
    
    # Show error bars
    if (error_bar) {
        out <- out +
            geom_errorbar(
                aes(ymin = ci_low, ymax = ci_hi, color = ci_label),
                width = error_bar_width,
                size = error_bar_size
            ) +
            scale_color_manual(name = "", values = "black") +
            guides(
                fill = guide_legend(order = 1), # bar legend first
                color = guide_legend(order = 2) # ci legend last
            )
    }
    return(out)
}


# Plot perf of models using trend lines
plot_perf_trend <- function(
        df_long,
        f = calc_pccc,
        group = "age_range",
        text_scale = 1,
        title = NULL,
        x_title = NULL,
        y_title = NULL,
        title_size = 10,
        x_title_size = 9 * text_scale,
        y_title_size = 9 * text_scale,
        x_text_size = 8 * text_scale,
        y_text_size = 9 * text_scale,
        legend_text_size = 9 * text_scale,
        legend_position = "bottom",
        legend_direction = "horizontal",
        legend_spacing = unit(0, "pt"),
        y_limits = c(0, 1.15),
        y_breaks = seq(0, 1, by = 0.1),
        x_color_labels = c(
            "G3" = "GPT-3.5",
            "G4" = "GPT-4",
            "G5" = "GPT-5",
            "ISV" = "InSilicoVA",
            "I5" = "InterVA-5"
        ),
        x_linetypes = c(
            "solid",
            "solid",
            "dashed",
            "dotted",
            "dotted"
        ),
        x_linealpha = c(
            1,
            0.5,
            0.9,
            0.9,
            0.7
        ),
        x_linewidths = c(
            2,
            1,
            1.15,
            2,
            1.15
        ),
        plot_mean = F
    ) {
    
    # Calc metric
    df_plot <- df_long %>%
        group_by(.data[[group]], model_name, .drop = F) %>%
        summarize( # calc pccc
            age = paste0(unique(age), collapse = ","),
            metric = f(cod, model_cghr10, N = unique(ncod)[1])
        ) %>%
        mutate( # add pccc to model name
            model_label = paste0(model_name, " (", round(metric, 2), ")")
        )

    # Get custom colors for models
    default_colors <- scales::hue_pal()(length(unique(df_plot$model_name)))
    names(default_colors) <- unique(df_plot$model_name)
    default_colors["Mean Across Models"] <- "black"
    
    # Plot trend lines for each model
    out <- ggplot(
        df_plot,
        aes(
            x = .data[[group]],
            y = metric,
            color = model_name,
            group = model_name
        )
    ) +
        geom_line(
            aes(
                linetype = model_name,
                alpha = model_name,
                linewidth = model_name
            )
        ) + 
        scale_y_continuous(
            limits = y_limits,
            breaks = y_breaks
        ) +
        scale_color_manual(
            values = default_colors,
            labels = x_color_labels
        ) +
        scale_linetype_manual(
            values = x_linetypes,
            labels = x_color_labels
        ) +
        scale_linewidth_manual(
           values = x_linewidths,
           labels = x_color_labels
        ) +
        scale_alpha_manual(
            values = x_linealpha,
            labels = x_color_labels
        ) +
        labs(
            title = title,
            x = x_title,
            y = y_title
        ) +
        theme_minimal() +
        theme(
            legend.key.width = unit(32, "pt"),
            legend.title = element_blank(),
            legend.text = element_text(size = legend_text_size),
            legend.position = legend_position, 
            legend.direction = legend_direction,
            legend.box.spacing = legend_spacing,
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            plot.title = element_text(size = title_size),
            axis.title.x = element_text(size = x_title_size),
            axis.title.y = element_text(size = y_title_size),
            axis.text.x = element_text(size = x_text_size),
            axis.text.y = element_text(size = y_text_size)
        )
    
    # Add mean line
    if (plot_mean) {
        out <- out +
            geom_line(
                data = df_plot %>%
                    group_by(.data[[group]]) %>%
                    summarize(
                        metric = mean(metric),
                        model_name = "Mean Across Models"
                    ),
                aes(color = model_name),
                se = F,
                linewidth = 1.5
            )
    }
    return(out)
    
}

```

# Data

## Raw Data

Load data from `data` folder.

```{r}

# Read preds
df_raw <- read_csv("../data/healsl_rd1to2_cod_v2.csv")

# Display preds
cat("Cases: ", nrow(df_raw), "\n")
print(head(df_raw))

```

## Physician Cases

Filter for cases that were coded by physicians.

```{r}

# Filter for phy cases
df <- df_raw %>%
    filter(!is.na(physician_cghr10) & !is.null(physician_cghr10))

# Preview phy cases
cat("Non-coded Cases: ", nrow(df_raw) - nrow(df))
cat("\nPhysician Coded Cases: ", nrow(df))
df

```

## Inspect Data

Display case counts for raw data and prepared data after combining models and filtering for physician coded cases.

```{r}

cat("\nRaw Data\n--------\n\n")

# Physician info
raw_physicians <- df_raw$physician_cghr10
cat(paste0(
    "Cases (",
    length(unique(na.omit(raw_physicians))), " CODs): ",
    length(raw_physicians),
    "\n"
))

# Physician agreed cases
raw_physicians_agree <- df_raw %>%
    filter(is_agreed == TRUE) %>%
    pull(physician_cghr10)
cat(paste0(
    "Agreed Cases (",
    length(unique(na.omit(raw_physicians_agree))), " CODs): ",
    length(raw_physicians_agree),
    "\n"
))

# Models cases
raw_models <- df_raw %>% select(ends_with("_cghr10"), -physician_cghr10)
for (mcol in colnames(raw_models)) {
    m <- raw_models[[mcol]]
    cat(paste0(
        mcol, " Predicted Cases (", 
        length(unique(na.omit(m))), " CODs): ",
        length(m) - sum(is.na(m) | is.null(m)),
        "\n"
    ))
}

cat("\nPrepared Data\n--------------\n")

# Physician cases
physicians <- df$physician_cghr10
cat(paste0(
    "Physician Coded Cases (",
    length(unique(na.omit(physicians))), " CODs): ",
    length(physicians) - sum(is.na(physicians) | is.null(physicians)),
    "\n"
))

# Physician agreed cases
physicians_agree <- df %>%
    filter(is_agreed == TRUE) %>%
    pull(physician_cghr10)
cat(paste0(
    "Physician Agreed Coded Cases (",
    length(unique(na.omit(physicians_agree))), " CODs): ",
    length(physicians_agree) - sum(is.na(physicians_agree) | is.null(physicians_agree)),
    "\n"
))

# Models cases
models <- df %>% select(ends_with("_cghr10"), -physician_cghr10)
for (mcol in colnames(models)) {
    m <- models[[mcol]]
    cat(paste0(
        mcol, " Predicted Cases (", 
        length(unique(na.omit(m))), " CODs): ",
        length(m) - sum(is.na(m) | is.null(m)),
        "\n"
    ))
}

# Models agreed cases
models_agree <- df %>%
    filter(is_agreed == TRUE) %>%
    select(ends_with("_cghr10"), -physician_cghr10)
for (mcol in colnames(models_agree)) {
    m <- models_agree[[mcol]]
    cat(paste0(
        mcol, " Agreed Predicted Cases (", 
        length(unique(na.omit(m))), " CODs): ",
        length(m) - sum(is.na(m) | is.null(m)),
        "\n"
    ))
}
```

# Preparation

Prepare data for plotting.

## Agreed Cases

Filter for agreed cases only.

```{r}

# Agreed cases only
df <- df %>% filter(is_agreed == TRUE)

```

## Rename COD Names

Rename COD names into simpler names.

```{r}

# Rename COD names
df <- df %>%
    mutate(across(
        c(
            physician_cghr10,
            gpt3_cghr10,
            gpt4_cghr10,
            gpt5_cghr10,
            insilicova_cghr10,
            interva5_cghr10
        ),
        ~ case_when(
            .x == "Maternal conditions" ~ "Maternal",
            .x == "Acute respiratory infections" ~ "Pneumonia",
            .x == "Neonatal infections" ~ "Infections",
            .x == "Road and transport injuries" ~ "Road/transport injuries",
            .x == "Liver and alcohol related diseases" ~ "Liver/alcohol",
            .x == "Birth asphyxia and trauma" ~ "Birth asphyxia/trauma",
            .default = .x
        ) %>%
            str_replace_all(" diseases| disease", "")
    ))

# Rename raw COD names
df_raw <- df_raw %>%
    mutate(across(
        c(
            physician_cghr10,
            gpt3_cghr10,
            gpt4_cghr10,
            gpt5_cghr10,
            insilicova_cghr10,
            interva5_cghr10
        ),
        ~ case_when(
            .x == "Maternal conditions" ~ "Maternal",
            .x == "Acute respiratory infections" ~ "Pneumonia",
            .x == "Neonatal infections" ~ "Infections",
            .x == "Road and transport injuries" ~ "Road/transport injuries",
            .x == "Liver and alcohol related diseases" ~ "Liver/alcohol",
            .x == "Birth asphyxia and trauma" ~ "Birth asphyxia/trauma",
            .default = .x
        ) %>%
            str_replace_all(" diseases| disease", "")
    ))

```

## Shorten COD Names

Shorten long names into new lines.

```{r}

# Shorten COD names
df <- df %>%
    mutate(across(
        c(
            physician_cghr10,
            gpt3_cghr10,
            gpt4_cghr10,
            gpt5_cghr10,
            insilicova_cghr10,
            interva5_cghr10
        ),
        ~ str_wrap(.x, width = 18)
    ))

# Shorten raw COD names
df_raw <- df_raw %>%
    mutate(across(
        c(
            physician_cghr10,
            gpt3_cghr10,
            gpt4_cghr10,
            gpt5_cghr10,
            insilicova_cghr10,
            interva5_cghr10
        ),
        ~ str_wrap(.x, width = 18)
    ))

```

## Standardize Age Ranges

Standardize age ranges to follow age range and unit format in lowercase.

```{r}

# Standardize raw ages to age range + unit
df_raw <- df_raw %>%
    mutate(
        age_range_orig = age_range,
        age_range = str_to_lower(age_range),
        age_range = str_replace( # neo ages
            age_range,
            " \\(early neonatal\\)| \\(late neonatal\\)",
            ""
        ),
        age_range = case_when( # adult ages
            age == "adult" ~ paste0(age_range, " years"),
            .default = age_range
        )
    )

# Standardize ages to age range + unit
df <- df %>%
    mutate(
        age_range = str_to_lower(age_range),
        age_range = str_replace( # neo ages
            age_range,
            " \\(early neonatal\\)| \\(late neonatal\\)",
            ""
        ),
        age_range = case_when( # adult ages
            age == "adult" ~ paste0(age_range, " years"),
            .default = age_range
        )
    )
df$age_range %>% unique

```

## Order Age Ranges

Order age ranges in sequential order.

```{r}

# Get age ranges for adult
ages_adult <- df %>%
    filter(age == "adult") %>%
    pull(age_range) %>%
    unique %>%
    str_replace(., "10-14 years", "12-14 years") %>%
    order_arange

# Get age ranges for child
ages_child <- df %>%
    filter(age == "child") %>%
    pull(age_range) %>%
    unique %>%
    order_arange

# Get age ranges for neo
ages_neo <- df %>%
    filter(age == "neo") %>%
    pull(age_range) %>%
    unique %>%
    order_arange

# Order raw age ranges as factor
df_raw <- df_raw %>%
    mutate(
        age_range = factor(
            age_range %>%
                str_replace(., "10-14 years", "12-14 years"),
            levels = c(ages_neo, ages_child, ages_adult)
        )
    )

# Order age ranges as factor
df <- df %>%
    mutate(
        age_range = factor(
            age_range %>%
                str_replace(., "10-14 years", "12-14 years"),
            levels = c(ages_neo, ages_child, ages_adult)
        )
    )
df$age_range %>% unique

```

## COD Case Count

Count the number and percentage of cases for each CGHR-10 COD.

```{r}

# Count cases per cod
df <- df %>%
    left_join( # count cases and add as column
        df %>%
            group_by(age, physician_cghr10) %>%
            count(name = "physician_cghr10_n") %>%
            group_by(age) %>%
            mutate(
                physician_cghr10_perc = (physician_cghr10_n / sum(physician_cghr10_n)) * 100
            ),
        by = c("age", "physician_cghr10")
    )
df %>% select(rowid, age, physician_cghr10, physician_cghr10_n)

```

## Format Data

Format data into long form with each record consisting of one model.

Also calculates supplementary information for plots and tables.

```{r}

# Number of causes for each age group
ncod_adult <- 19
ncod_child <- 10
ncod_neo <- 7
ncod_all <- ncod_adult + ncod_child + ncod_neo

# Model full name
gpt3_name <- "GPT-3.5"
gpt4_name <- "GPT-4"
gpt5_name <- "GPT-5"
isva_name <- "InSilicoVA"
iva5_name <- "InterVA-5"

# Pivot data into long format for calc
df_long <- df %>%
    filter( # remove low count cod
        physician_cghr10_n >= 10
    ) %>%
    pivot_longer( # long format
        ends_with("_cghr10") & !starts_with("physician_"),
        names_to = "model",
        values_to = "model_cghr10"
    ) %>%
    mutate(
        model = str_replace( # remove cghr10
            model,
            "_cghr10",
            ""
        ),
        model_name = str_replace_all(model, "_", " or ") %>% # use full model names
            str_replace(., "gpt3", gpt3_name) %>%
            str_replace(., "gpt4", gpt4_name) %>%
            str_replace(., "gpt5", gpt5_name) %>%
            str_replace(., "insilicova", isva_name) %>%
            str_replace(., "interva5", iva5_name) %>%
            factor(., levels = c(
                gpt5_name,
                gpt4_name,
                gpt3_name,
                isva_name,
                iva5_name
            )),
        age_name = case_when( # use full age name
            age == "adult" ~ "Adult",
            age == "child" ~ "Child",
            age == "neo" ~ "Neonate",
            .default = NA
        ),
        cod = physician_cghr10,
        cod_name = paste0(
            cod,
            "\n(n=", physician_cghr10_n, ", ", round(physician_cghr10_perc), "%)"
        ),
        ncod = case_when( # num of cods for each age group
            age == "adult" ~ ncod_adult,
            age == "child" ~ ncod_child,
            age == "neo" ~ ncod_neo,
            .default = NA
        )
    ) %>%
    filter( # remove combine models
        !str_detect(model_name, " or ")
    )

# Calc num of records for each age
n_all <- df %>% group_by(age) %>% count
n_adult <- n_all %>% filter(age == "adult") %>% pull(n)
n_child <- n_all %>% filter(age == "child") %>% pull(n)
n_neo <- n_all %>% filter(age == "neo") %>% pull(n)

# Calc perc of records for each age
perc_adult <- (n_adult / nrow(df)) * 100
perc_child <- (n_child / nrow(df)) * 100
perc_neo <- (n_neo / nrow(df)) * 100

# Calc num of cods for each age
cod_all <- df_long %>% group_by(age) %>% summarize(n = n_distinct(cod))
cod_adult <- cod_all %>% filter(age == "adult") %>% pull(n)
cod_child <- cod_all %>% filter(age == "child") %>% pull(n)
cod_neo <- cod_all %>% filter(age == "neo") %>% pull(n)

# Preview
df_long

```

# Plots

Generate plots from the predictions.

## PCCC and CSMF

```{r}

# Prep data
df_plot <- df_long %>%
    group_by(model_name) %>%
    summarize(
        pccc = calc_pccc(cod, model_cghr10, N = ncod_all),
        csmf = calc_csmf_acc(cod, model_cghr10)
    ) %>%
    arrange(pccc) %>%
    mutate(
        model_name = model_name %>%
            case_when(
                . == "G3" ~ "GPT-3.5",
                . == "G4" ~ "GPT-4",
                . == "G5" ~ "GPT-5",
                . == "ISV" ~ "InSilicoVA",
                . == "I5" ~ "InterVA-5",
                .default = .
            ) %>%
            factor(., levels = .)
    ) %>%
    pivot_longer(
        c(pccc, csmf),
        names_to = "metric",
        values_to = "metric_value"
    ) %>%
    mutate(
        metric = str_to_upper(metric) %>%
            str_replace(., "CSMF", "Cause Specific Mortality Fraction (CSMF) Accuracy") %>%
            str_replace(., "PCCC", "Partial Chance Corrected Concordance (PCCC)") %>%
            factor(., levels = rev(c("Partial Chance Corrected Concordance (PCCC)", "Cause Specific Mortality Fraction (CSMF) Accuracy"))),
        metric_label = paste0(metric, "=", round(metric_value, 2)) %>%
            str_replace(., " Accuracy", ""),
        metric_pattern = case_when(
            metric == "Partial Chance Corrected Concordance (PCCC)" ~ "stripe",
            metric == "Cause Specific Mortality Fraction (CSMF) Accuracy" ~ "none"
        )
    )

# Plot data
fig_pccc_csmf <- ggplot(
    df_plot,
    aes(
        x = metric_value,
        y = model_name,
        fill = metric
    )
) +
    geom_bar_pattern(
        aes(pattern = metric_pattern),
        stat = "identity",
        position = position_dodge(),
        color = "black",
        width = 0.7,
        pattern_fill = "black",
        pattern_angle = 45,
        pattern_spacing = 0.05,
        pattern_colour = "black"
    ) +
    geom_text(
        data = df_plot,
        aes(
            label = round(metric_value, 2)
        ),
        position = position_dodge(width = 0.7),
        color = "black",
        fill = "white",
        size = 4.5,
        label.r = unit(0, "pt"),
        vjust = 0.5,
        hjust = -0.25
    ) +
    scale_x_continuous(
        limits = c(0, 1),
        breaks = seq(0, 1, 0.1)
    ) +
    scale_pattern_manual(values = c(
        "stripe" = "stripe",
        "none" = "none"
    )) +
    labs(
        title = NULL,
        x = NULL,
        y = NULL
    ) +
    guides(
        pattern = "none",
        fill = guide_legend(
            reverse = TRUE,
            override.aes = list(
                pattern = c("stripe", "none")
            )
        )
    ) +
    theme_minimal() +
    theme(
        legend.position = "bottom",
        legend.margin = margin(c(-8, -2, -2, -60)),
        legend.title = element_blank(),
        legend.text = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(size = 10),
        axis.title.x = element_text(size = 9),
        axis.title.y = element_text(size = 9),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12)
    )

# Save plot
ggsave(
    "../manuscript/figures/fig-pccc-csmf.pdf",
    fig_pccc_csmf,
    width = 170,
    height = 100,
    units = "mm",
    dpi = 300
)
ggsave(
    "../manuscript/figures/png/fig-pccc-csmf.png",
    fig_pccc_csmf,
    width = 170,
    height = 100,
    units = "mm",
    dpi = 300
)

# Display plot
fig_pccc_csmf

```

## PCCC By Age

```{r}

# Plot perf for pccc by age groups
fig_pccc <- plot_perf_bar(
    df_long %>%
        mutate(
            age_label = case_when(
                age_name == "Adult"
                    ~ paste0(age_name, "\n12+ Years\n(n=", n_adult, ")"),
                age_name == "Child"
                    ~ paste0(age_name, "\n28 Days to 11 Years\n(n=", n_child, ")"),
                age_name == "Neonate"
                    ~ paste0(age_name, "\n<28 Days\n(n=", n_neo, ")"),
                .default = age_name
            ) %>%
                factor(., levels = unique(.) %>% sort %>% rev)
        ),
    calc_pccc,
    group = "age_label",
    text_scale = 0.95
)

# Save plot
ggsave(
    "../manuscript/figures/fig-pccc.pdf",
    fig_pccc,
    width = 170,
    height = 75,
    units = "mm",
    dpi = 300
)
ggsave(
    "../manuscript/figures/png/fig-pccc.png",
    fig_pccc,
    width = 170,
    height = 75,
    units = "mm",
    dpi = 300
)

# Display plot
fig_pccc

```

## PCCC by Age and COD

### Adult

```{r}

# Plot adult perf for pccc by cod
fig_pccc_cod_adult <- plot_perf_bar(
    df_long %>% filter(age == "adult"),
    calc_pccc,
    group = "cod_name",
    title = paste0("Adult Causes of Death, 12+ Years (n=", n_adult, ")"),
    text_scale = 0.7,
    x_hpos = -0.27
)
    
# Save plot
ggsave(
    "../manuscript/figures/fig-pccc-cod-adult.pdf",
    fig_pccc_cod_adult,
    width = 170, # min 170 mm full / 85 mm half
    height = 225, # max 225 mm
    units = "mm",
    dpi = 300
)
ggsave(
    "../manuscript/figures/png/fig-pccc-cod-adult.png",
    fig_pccc_cod_adult,
    width = 170, # min 170 mm full / 85 mm half
    height = 225, # max 225 mm
    units = "mm",
    dpi = 300
)

# Display plot
fig_pccc_cod_adult

```

### Child

```{r}

# Plot child perf for pccc by cod
fig_pccc_cod_child <- plot_perf_bar(
    df_long %>% filter(age == "child"),
    calc_pccc,
    group = "cod_name",
    title = paste0("Child Causes of Death, 28 Days to 11 Years (n=", n_child, ")"),
    text_scale = 0.92,
    x_hpos = -0.43
)

# Save plot
ggsave(
    "../manuscript/figures/fig-pccc-cod-child.pdf",
    fig_pccc_cod_child,
    width = 170,
    height = 150,
    units = "mm",
    dpi = 300
)
ggsave(
    "../manuscript/figures/png/fig-pccc-cod-child.png",
    fig_pccc_cod_child,
    width = 170,
    height = 150,
    units = "mm",
    dpi = 300
)

# Display plot
fig_pccc_cod_child

```

### Neonate

```{r}

# Plot neo perf for pccc by cod
fig_pccc_cod_neo <- plot_perf_bar(
    df_long %>% filter(age == "neo"),
    calc_pccc,
    group = "cod_name",
    title = paste0("Neonatal Causes of Death, <28 Days (n=", n_neo, ")"),
    x_title = "Partial Chance Corrected Concordance (PCCC)",
    text_scale = 0.92,
    x_hpos = -0.43
)

# Save plot
ggsave(
    "../manuscript/figures/fig-pccc-cod-neo.pdf",
    fig_pccc_cod_neo,
    width = 170,
    height = 100,
    units = "mm",
    dpi = 300
)
ggsave(
    "../manuscript/figures/png/fig-pccc-cod-neo.png",
    fig_pccc_cod_neo,
    width = 170,
    height = 100,
    units = "mm",
    dpi = 300
)

# Display plot
fig_pccc_cod_neo

```

## PCCC by Age Range

```{r}

# Format levels of age range with breaks
age_levels <- levels(df_long$age_range) %>%
    str_replace(., " ", "\n")

# Plot perf trend over age ranges
fig_pccc_age <- plot_perf_trend(
    df_long %>%
        mutate(
            age_range = str_replace(age_range, " ", "\n") %>%
                factor(., levels = age_levels)
        ),
    y_title = "Partial Chance Corrected Concordance (PCCC)",
    text_scale = 0.9,
    y_limits = c(0.2, 0.9),
    y_breaks = seq(0, 1, by = 0.05),
)

# Save plot
ggsave(
    "../manuscript/figures/fig-pccc-age.pdf",
    fig_pccc_age,
    width = 170,
    height = 100,
    units = "mm",
    dpi = 300
)
ggsave(
    "../manuscript/figures/png/fig-pccc-age.png",
    fig_pccc_age,
    width = 170,
    height = 100,
    units = "mm",
    dpi = 300
)

# Display plot
fig_pccc_age

```

# Tables

```{r}

# Get total records
tab_all <- df_raw %>%
    filter(!is.na(physician_cghr10) & !is.null(physician_cghr10)) # remove uncoded cases
tab_all_n <- tab_all %>% nrow

# Get n records by age group
tab_agegroup <- df %>%
    group_by(age) %>%
    summarize(
        n = n(),
        n_cod_uniq = n_distinct(physician_cghr10)
    ) %>%
    ungroup %>%
    left_join(
        tab_all %>%
            group_by(age) %>%
            summarize(n_all = n()) %>%
            ungroup,
        by = "age"
    ) %>%
    mutate(
        n_agree_perc = round(n / n_all * 100, 1) 
    )

# Get n records by age group and sex
tab_agegroup_sex <- df %>%
    group_by(age, sex) %>%
    summarize(n_sex = n()) %>%
    pivot_wider(
        names_from = sex,
        values_from = n_sex,
        names_glue = "n_{str_to_lower(sex)}"
    )

```

## Records by COD

```{r}

# Get n records by cod
tab_cod_n <- df %>%
    group_by(age, physician_cghr10) %>%
    summarize(n_cod = n()) %>%
    left_join(
        tab_all %>%
        group_by(age, physician_cghr10) %>%
        summarize(n_cod_all = n()) %>%
        ungroup,
        by = c("age", "physician_cghr10")
    ) %>%
    mutate(
        n_cod_agree_perc = round(n_cod / n_cod_all * 100 , 1)
    )

# Get n records by cod and sex
tab_cod_sex <- df %>%
    group_by(age, physician_cghr10, sex) %>%
    summarize(n_cod_sex = n()) %>%
    pivot_wider(
        names_from = sex,
        values_from = n_cod_sex,
        names_glue = "n_cod_{str_to_lower(sex)}"
    )

# Combine all
tab_cod <- tab_cod_n %>%
    left_join(tab_agegroup, by = "age") %>%
    left_join(tab_agegroup_sex, by = "age") %>%
    left_join(tab_cod_sex, by = c("age", "physician_cghr10"))

# Calc perc
tab_cod <- tab_cod %>%
    mutate(
        n_perc = round(n / nrow(df) * 100, 1),
        n_male_perc = round(n_male / n * 100, 1),
        n_female_perc = round(n_female / n * 100, 1),
        n_cod_perc = round(n_cod / n * 100, 1),
        n_cod_male_perc = round(n_cod_male / n_cod * 100, 1),
        n_cod_female_perc = round(n_cod_female / n_cod * 100, 1)
    ) %>%
    ungroup

# Add labels to combined
tab_cod <- tab_cod %>%
    mutate(
        age = if_else(age == "neo", "Neonate", str_to_title(age)),
        `Age Group` = glue(
            "{age} (n={n}, {n_perc}%, {n_cod_uniq} CODs)\n- Female (n={n_female}, {n_female_perc}%)\n- Male (n={n_male}, {n_male_perc}%)\n- Agree ({n_agree_perc}% of {n_all})"
        ),
        `Cause of Death (COD)` = glue(
            "{str_to_title(physician_cghr10)}"
        ),
        `Agree` = glue("{n_cod_agree_perc}% of {n_cod_all}"),
        `Female` = glue("{n_cod_female} ({n_cod_female_perc}%)"),
        `Male` = glue("{n_cod_male} ({n_cod_male_perc}%)"),
        `Total` = glue("{n_cod} ({n_cod_perc}%)")
    ) %>%
    mutate(across(
        c(`Female`, `Male`, `Total`),
        ~ if_else(.x == "NA (NA%)", NA, .x)
    ))

# Save table
write_csv(tab_cod, "../data/healsl_rd1to2_cod_profile_v1.csv", na = "")

# Display table
tab_cod
```

```{r}

# Prepare and format the data
tab_cod_format <- tab_cod %>%
    ungroup %>%
    select(`Age Group`, `Cause of Death (COD)`, Agree, Female, Male, Total) %>%
    mutate(
        `Age Group` = str_replace_all(`Age Group`, "\\\\n|\\n", "<br>"),
        `Cause of Death (COD)` = str_replace_all(`Cause of Death (COD)`, "\\\\n|\\n", "<br>"),
        Agree =  str_replace_all(Agree, "\\\\n|\\n", "<br>")
    ) %>%
    group_by(`Age Group`) %>%
    mutate(`Age Group` = ifelse(row_number() == 1, `Age Group`, "")) %>%
    ungroup()

# Build the table
tab_cod_html <- tab_cod_format %>%
    kable(format = "html", escape = FALSE, align = "l") %>%
    kable_styling(full_width = FALSE, position = "left", html_font = "Arial") %>%
    row_spec(0, extra_css = "border-top: 2px solid; border-bottom: 1px solid;") %>%     # thick top border
    row_spec(nrow(tab_cod_format), extra_css = "border-bottom: 2px solid;")
tab_cod_html

```

## Records by Age

```{r}

# Get n records by age
tab_age_n <- df %>%
    group_by(age, age_range) %>%
    summarize(n_age = n()) %>%
    left_join(
        tab_all %>%
        group_by(age, age_range) %>%
        summarize(n_age_all = n()) %>%
        ungroup,
        by = c("age", "age_range")
    ) %>%
    mutate(
        n_age_agree_perc = round(n_age / n_age_all * 100 , 1)
    )

# Get n records by age and sex
tab_age_sex <- df %>%
    group_by(age, age_range, sex) %>%
    summarize(n_age_sex = n()) %>%
    pivot_wider(
        names_from = sex,
        values_from = n_age_sex,
        names_glue = "n_age_{str_to_lower(sex)}"
    )

# Combine all
tab_age <- tab_age_n %>%
    left_join(tab_agegroup, by = "age") %>%
    left_join(tab_agegroup_sex, by = "age") %>%
    left_join(tab_age_sex, by = c("age", "age_range"))

# Calc perc
tab_age <- tab_age %>%
    mutate(
        n_perc = round(n / nrow(df) * 100, 1),
        n_male_perc = round(n_male / n * 100, 1),
        n_female_perc = round(n_female / n * 100, 1),
        n_age_perc = round(n_age / n * 100, 1),
        n_age_male_perc = round(n_age_male / n_age * 100, 1),
        n_age_female_perc = round(n_age_female / n_age * 100, 1)
    ) %>%
    ungroup

# Add labels to combined
tab_age <- tab_age %>%
    mutate(
        age = if_else(age == "neo", "Neonate", str_to_title(age)),
        `Age Group` = glue(
            "{age} (n={n}, {n_perc}%, {n_cod_uniq} CODs)\n- Female (n={n_female}, {n_female_perc}%)\n- Male (n={n_male}, {n_male_perc}%)\n- Agree ({n_agree_perc}% of {n_all})"
        ),
        `Age Range` = glue(
            "{str_to_title(age_range)}"
        ),
        `Agree` = glue("{n_age_agree_perc}% of {n_age_all}"),
        `Female` = glue("{n_age_female} ({n_age_female_perc}%)"),
        `Male` = glue("{n_age_male} ({n_age_male_perc}%)"),
        `Total` = glue("{n_age} ({n_age_perc}%)")
    ) %>%
    mutate(across(
        c(`Female`, `Male`, `Total`),
        ~ if_else(.x == "NA (NA%)", NA, .x)
    ))

# Save table
write_csv(tab_age, "../data/healsl_rd1to2_age_profile_v1.csv", na = "")

# Display table
tab_age

```

```{r}

# Prepare and format the data
tab_age_format <- tab_age %>%
    ungroup %>%
    select(`Age Group`, `Age Range`, Agree, Female, Male, Total) %>%
    mutate(
        `Age Group` = str_replace_all(`Age Group`, "\\\\n|\\n", "<br>"),
        `Age Range` = str_replace_all(`Age Range`, "\\\\n|\\n", "<br>"),
        `Agree` = str_replace_all(`Agree`, "\\\\n|\\n", "<br>")
    ) %>%
    group_by(`Age Group`) %>%
    mutate(`Age Group` = ifelse(row_number() == 1, `Age Group`, "")) %>%
    ungroup()

# Build the table
tab_age_html <- tab_age_format %>%
    kable(format = "html", escape = FALSE, align = "l") %>%
    kable_styling(full_width = FALSE, position = "left", html_font = "Arial") %>%
    row_spec(0, extra_css = "border-top: 2px solid; border-bottom: 1px solid;") %>%     # thick top border
    row_spec(nrow(tab_age_format), extra_css = "border-bottom: 2px solid;")
tab_age_html

```