B-Reader Accuracy Analysis for Pneumoconiosis Classification

Introduction

This document presents an optimized analysis of B-reader accuracy for pneumoconiosis classification. The analysis evaluates reader performance across multiple outcomes including parenchymal abnormalities, pleural abnormalities, opacity classifications, and small opacity detection.

Objectives

• Analyze B-reader accuracy across four key outcomes
• Calculate sensitivity, specificity, and overall accuracy metrics
• Generate comprehensive summary statistics
• Provide reproducible analysis workflow

Setup and Configuration

Load Required Libraries

# Load required libraries
suppressPackageStartupMessages({
  library(readxl)
  library(dplyr)
  library(tidyr)
  library(purrr)
})

Configuration

# Set working directory (use relative paths for portability)
setwd("D:\\OneDrive - Michigan State University\\images_bl\\readings\\reading_accuracy\\")

# File paths
DATA_FILES <- list(
  detailed_readings = "detailed_B-readings_8_2_2024.xlsx",
  detailed_readings_2nd = "detailed_readings_NIOSH_2nd_batch.csv",
  test_consensus = "test_consensus_4_7_2025.csv"
)

Helper Functions

Accuracy Metrics Calculation

#' Calculate accuracy metrics from confusion matrix
#' @param confusion_matrix A 2x2 confusion matrix
#' @return List with accuracy, sensitivity, and specificity
calculate_metrics <- function(confusion_matrix) {
  if (is.null(confusion_matrix) || length(confusion_matrix) == 0) {
    return(list(accuracy = NA, sensitivity = NA, specificity = NA))
  }
  
  # Handle different matrix dimensions
  n_rows <- nrow(confusion_matrix)
  n_cols <- ncol(confusion_matrix)
  
  if (n_rows == 2 && n_cols == 2) {
    # Standard 2x2 matrix
    accuracy <- (confusion_matrix[1,1] + confusion_matrix[2,2]) / sum(confusion_matrix)
    sensitivity <- confusion_matrix[2,2] / sum(confusion_matrix[2,])
    specificity <- confusion_matrix[1,1] / sum(confusion_matrix[1,])
  } else if (n_rows == 1 && n_cols == 2) {
    # Only one class present
    if (rownames(confusion_matrix)[1] == "1") {
      accuracy <- sensitivity <- confusion_matrix[1,2] / sum(confusion_matrix)
      specificity <- NA
    } else {
      accuracy <- specificity <- confusion_matrix[1,1] / sum(confusion_matrix)
      sensitivity <- NA
    }
  } else if (n_rows == 2 && n_cols == 1) {
    # Only one prediction class
    if (colnames(confusion_matrix)[1] == "1") {
      accuracy <- confusion_matrix[2,1] / sum(confusion_matrix)
      sensitivity <- 1
      specificity <- 0
    } else {
      accuracy <- confusion_matrix[1,1] / sum(confusion_matrix)
      sensitivity <- 0
      specificity <- 1
    }
  } else if (n_rows == 1 && n_cols == 1) {
    # Perfect prediction
    accuracy <- sensitivity <- specificity <- 1
  } else {
    # Handle other cases
    accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
    sensitivity <- specificity <- NA
  }
  
  list(accuracy = accuracy, sensitivity = sensitivity, specificity = specificity)
}

#' Calculate multi-class accuracy
#' @param confusion_matrix A confusion matrix
#' @return Accuracy value
calculate_multiclass_accuracy <- function(confusion_matrix) {
  if (is.null(confusion_matrix) || length(confusion_matrix) == 0) {
    return(NA)
  }
  
  n_rows <- nrow(confusion_matrix)
  n_cols <- ncol(confusion_matrix)
  
  if (n_rows == n_cols) {
    # Square matrix
    sum(diag(confusion_matrix)) / sum(confusion_matrix)
  } else if (n_rows < n_cols) {
    # More predicted classes than actual
    sum(diag(confusion_matrix[, rownames(confusion_matrix)])) / sum(confusion_matrix)
  } else if (n_rows > n_cols) {
    # More actual classes than predicted
    sum(diag(confusion_matrix[colnames(confusion_matrix), ])) / sum(confusion_matrix)
  } else {
    # Single class
    1
  }
}

Data Loading and Preprocessing

Load Detailed Readings Data

#' Load and preprocess detailed readings data
load_detailed_readings <- function() {
  # Load first batch
  detailed_readings <- read_excel(DATA_FILES$detailed_readings)
  
  # Load second batch
  detailed_readings_2nd <- read.csv(DATA_FILES$detailed_readings_2nd)
  
  # Process second batch
  detailed_readings_2nd_processed <- detailed_readings_2nd %>%
    select(1:5, 26, 27, 108) %>%
    mutate(
      SmallOpacities_Profusion_major = as.numeric(substr(SMALL_OPAC_CATEGORY, 1, 1)),
      largeopacity_n = case_when(
        LARGE_OPAC_STAGE == "O" ~ 0,
        LARGE_OPAC_STAGE == "A" ~ 1,
        LARGE_OPAC_STAGE == "B" ~ 2,
        LARGE_OPAC_STAGE == "C" ~ 3,
        TRUE ~ NA_real_
      )
    ) %>%
    select(1, 5, 9, 10, 8)
  
  # Process first batch
  detailed_readings_processed <- detailed_readings %>%
    select(1, 2, 14, 17, 13)
  
  # Ensure column names match
  names(detailed_readings_2nd_processed) <- names(detailed_readings_processed)
  
  # Combine datasets
  rbind(detailed_readings_processed, detailed_readings_2nd_processed) %>%
    rename(ID = 1)
}

# Load the data
detailed_reading <- load_detailed_readings()
test <- read.csv(DATA_FILES$test_consensus)
test_detailed_reading <- merge(test, detailed_reading, by = "ID", all.x = TRUE)

cat("Data loaded successfully!\n")

## Data loaded successfully!

cat("Number of records:", nrow(test_detailed_reading), "\n")

## Number of records: 1078

cat("Number of unique readers:", length(unique(test_detailed_reading$READER_ID)), "\n")

## Number of unique readers: 46

Create Reading Classifications

#' Create reading classifications
create_reading_classifications <- function(test_detailed_reading) {
  test_detailed_reading %>%
    mutate(
      parenchymal_reading = as.numeric(SmallOpacities_Profusion_major > 0 | largeopacity_n > 0),
      PARENCHYMAL_PLEURAL_reading = as.numeric(parenchymal_reading > 0 | PleuralAbnormal > 0),
      Opacity_reading = case_when(
        SmallOpacities_Profusion_major > 0 & largeopacity_n == 0 ~ 1,
        SmallOpacities_Profusion_major > 0 & largeopacity_n > 0 ~ 2,
        TRUE ~ 0
      )
    )
}

#' Create small opacity classifications
create_small_opacity_classifications <- function(test_detailed_reading) {
  test_detailed_reading %>%
    filter(small_opacity <= 1) %>%
    mutate(small_opacity_reading = as.numeric(SmallOpacities_Profusion_major == 1))
}

# Apply classifications
test_detailed_reading <- create_reading_classifications(test_detailed_reading)
test_detailed_reading2 <- create_small_opacity_classifications(test_detailed_reading)

cat("Classifications created successfully!\n")

## Classifications created successfully!

Main Analysis

Reader Accuracy Analysis Function

#' Analyze reader accuracy for all outcomes
analyze_reader_accuracy <- function(test_detailed_reading, test_detailed_reading2) {
  reader_ids <- unique(test_detailed_reading$READER_ID)
  
  # Pre-allocate results
  results <- list()
  
  for (i in seq_along(reader_ids)) {
    reader_id <- reader_ids[i]
    cat("Processing reader", i, "of", length(reader_ids), ":", reader_id, "\n")
    
    # Get reader data
    reader_data <- test_detailed_reading[test_detailed_reading$READER_ID == reader_id, ]
    reader_data2 <- test_detailed_reading2[test_detailed_reading2$READER_ID == reader_id, ]
    
    # Calculate metrics for each outcome
    results[[i]] <- list(
      reader_id = reader_id,
      
      # Outcome 1: Parenchymal abnormal
      outcome1 = calculate_metrics(
        table(reader_data$Parenchymal_abnormal_two_class, reader_data$parenchymal_reading)
      ),
      
      # Outcome 2: Parenchymal pleural abnormal
      outcome2 = calculate_metrics(
        table(reader_data$PARENCHYMAL_PLEURAL_abnormal_two_class, reader_data$PARENCHYMAL_PLEURAL_reading)
      ),
      
      # Outcome 3: Opacity three class
      outcome3 = list(
        accuracy = calculate_multiclass_accuracy(
          table(reader_data$Opacity_three_class, reader_data$Opacity_reading)
        ),
        sensitivity = NA,
        specificity = NA
      ),
      
      # Outcome 4: Small opacity
      outcome4 = calculate_metrics(
        table(reader_data2$small_opacity, reader_data2$small_opacity_reading)
      )
    )
  }
  
  results
}

Perform Analysis

cat("Starting reader accuracy analysis...\n")

## Starting reader accuracy analysis...

# Perform analysis
results <- analyze_reader_accuracy(test_detailed_reading, test_detailed_reading2)

## Processing reader 1 of 46 : 2573 
## Processing reader 2 of 46 : 1848 
## Processing reader 3 of 46 : 176 
## Processing reader 4 of 46 : 1194 
## Processing reader 5 of 46 : 1532 
## Processing reader 6 of 46 : 2112 
## Processing reader 7 of 46 : 3385 
## Processing reader 8 of 46 : 3115 
## Processing reader 9 of 46 : 3352 
## Processing reader 10 of 46 : 2800 
## Processing reader 11 of 46 : 1679 
## Processing reader 12 of 46 : 153 
## Processing reader 13 of 46 : 1196 
## Processing reader 14 of 46 : 3114 
## Processing reader 15 of 46 : 2967 
## Processing reader 16 of 46 : BT 
## Processing reader 17 of 46 : ReaderC 
## Processing reader 18 of 46 : ReaderA 
## Processing reader 19 of 46 : ReaderB 
## Processing reader 20 of 46 : 6 
## Processing reader 21 of 46 : 7 
## Processing reader 22 of 46 : PANEL 
## Processing reader 23 of 46 : 12 
## Processing reader 24 of 46 : 2 
## Processing reader 25 of 46 : 10 
## Processing reader 26 of 46 : 13 
## Processing reader 27 of 46 : 8 
## Processing reader 28 of 46 : 14 
## Processing reader 29 of 46 : 3 
## Processing reader 30 of 46 : 11 
## Processing reader 31 of 46 : 9 
## Processing reader 32 of 46 : 21 
## Processing reader 33 of 46 : 15 
## Processing reader 34 of 46 : 19 
## Processing reader 35 of 46 : 18 
## Processing reader 36 of 46 : 5 
## Processing reader 37 of 46 : 16 
## Processing reader 38 of 46 : 20 
## Processing reader 39 of 46 : 24 
## Processing reader 40 of 46 : 28 
## Processing reader 41 of 46 : 1 
## Processing reader 42 of 46 : 27 
## Processing reader 43 of 46 : 4 
## Processing reader 44 of 46 : 25 
## Processing reader 45 of 46 : 26 
## Processing reader 46 of 46 : 17

cat("Analysis completed successfully!\n")

## Analysis completed successfully!

cat("Number of readers analyzed:", length(results), "\n")

## Number of readers analyzed: 46

Results Processing

Create Summary Table

#' Convert results to summary table
create_summary_table <- function(results) {
  # Extract metrics
  metrics <- map_dfr(results, function(result) {
    data.frame(
      reader_id = result$reader_id,
      accuracy1 = result$outcome1$accuracy,
      sensitivity1 = result$outcome1$sensitivity,
      specificity1 = result$outcome1$specificity,
      accuracy2 = result$outcome2$accuracy,
      sensitivity2 = result$outcome2$sensitivity,
      specificity2 = result$outcome2$specificity,
      accuracy3 = result$outcome3$accuracy,
      accuracy4 = result$outcome4$accuracy,
      sensitivity4 = result$outcome4$sensitivity,
      specificity4 = result$outcome4$specificity
    )
  })
  
  # Calculate summary statistics
  summary_table <- matrix(0, 4, 7)
  
  # Count non-NA values
  summary_table[, 1] <- c(
    sum(!is.na(metrics$accuracy1)),
    sum(!is.na(metrics$accuracy2)),
    sum(!is.na(metrics$accuracy3)),
    sum(!is.na(metrics$accuracy4))
  )
  
  # Mean accuracy
  summary_table[, 2] <- c(
    mean(metrics$accuracy1, na.rm = TRUE),
    mean(metrics$accuracy2, na.rm = TRUE),
    mean(metrics$accuracy3, na.rm = TRUE),
    mean(metrics$accuracy4, na.rm = TRUE)
  )
  
  # SD accuracy
  summary_table[, 3] <- c(
    sd(metrics$accuracy1, na.rm = TRUE),
    sd(metrics$accuracy2, na.rm = TRUE),
    sd(metrics$accuracy3, na.rm = TRUE),
    sd(metrics$accuracy4, na.rm = TRUE)
  )
  
  # Mean sensitivity
  summary_table[, 4] <- c(
    mean(metrics$sensitivity1, na.rm = TRUE),
    mean(metrics$sensitivity2, na.rm = TRUE),
    NA,
    mean(metrics$sensitivity4, na.rm = TRUE)
  )
  
  # SD sensitivity
  summary_table[, 5] <- c(
    sd(metrics$sensitivity1, na.rm = TRUE),
    sd(metrics$sensitivity2, na.rm = TRUE),
    NA,
    sd(metrics$sensitivity4, na.rm = TRUE)
  )
  
  # Mean specificity
  summary_table[, 6] <- c(
    mean(metrics$specificity1, na.rm = TRUE),
    mean(metrics$specificity2, na.rm = TRUE),
    NA,
    mean(metrics$specificity4, na.rm = TRUE)
  )
  
  # SD specificity
  summary_table[, 7] <- c(
    sd(metrics$specificity1, na.rm = TRUE),
    sd(metrics$specificity2, na.rm = TRUE),
    NA,
    sd(metrics$specificity4, na.rm = TRUE)
  )
  
  summary_table
}

# Create summary table
summary_table <- create_summary_table(results)

Display Results

# Create a formatted results table
outcome_names <- c("Parenchymal Abnormal", "Parenchymal-Pleural Abnormal", 
                   "Opacity Three-Class", "Small Opacity")

results_df <- data.frame(
  Outcome = outcome_names,
  N_Readers = summary_table[, 1],
  Mean_Accuracy = round(summary_table[, 2], 3),
  SD_Accuracy = round(summary_table[, 3], 3),
  Mean_Sensitivity = round(summary_table[, 4], 3),
  SD_Sensitivity = round(summary_table[, 5], 3),
  Mean_Specificity = round(summary_table[, 6], 3),
  SD_Specificity = round(summary_table[, 7], 3)
)

# Display results
knitr::kable(results_df, 
             caption = "B-Reader Accuracy Summary Statistics",
             digits = 3)

B-Reader Accuracy Summary Statistics

Outcome	N_Readers	Mean_Accuracy	SD_Accuracy	Mean_Sensitivity	SD_Sensitivity	Mean_Specificity	SD_Specificity
Parenchymal Abnormal	46	0.930	0.092	0.912	0.174	0.910	0.222
Parenchymal-Pleural Abnormal	46	0.931	0.091	0.916	0.171	0.908	0.224
Opacity Three-Class	46	0.858	0.117	NA	NA	NA	NA
Small Opacity	46	0.808	0.208	0.693	0.292	0.902	0.237

Save Results

# Save results to CSV
write.csv(summary_table, "B-reader_test_6_12_2025.csv", row.names = FALSE, na = "")

cat("Results saved to B-reader_test_6_12_2025.csv\n")

## Results saved to B-reader_test_6_12_2025.csv

Summary

This analysis evaluated B-reader accuracy across four key outcomes:

1. Parenchymal Abnormal: Detection of parenchymal abnormalities
2. Parenchymal-Pleural Abnormal: Detection of combined parenchymal and pleural abnormalities
   
3. Opacity Three-Class: Multi-class classification of opacity types
4. Small Opacity: Detection of small opacities

The analysis provides comprehensive metrics including accuracy, sensitivity, and specificity for each outcome, enabling assessment of reader performance and identification of areas for improvement in pneumoconiosis classification.

Key Findings

• The analysis processed 46 readers
• Results are saved in the file B-reader_test_6_12_2025.csv
• Summary statistics show mean performance across all readers for each outcome

For detailed individual reader results, refer to the results object which contains metrics for each reader across all outcomes.