---
title: "Agreement Tutorial"
author: "Dr BASHIR SSUNA"
date: "2025-05-20"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, warning=FALSE, message=FALSE}
# Install if needed:
# install.packages(c("irr", "DescTools", "ggplot2", "BlandAltmanLeh"))

library(irr)            # for kappa2()
library(DescTools)      # for CCC()
library(ggplot2)        # for plotting
library(BlandAltmanLeh) # for bland.altman.plot()
library(tidyverse)
library(readxl)
library(kableExtra)
library(broom)
library(gt)         # for nice looking tables

```

Genearte sample datasets for the tutorial

1. For categorical data

```{r}
set.seed(123)
n_cat <- 1000

# Simulate two raters with three categories
rater1 <- sample(c("A","B","C"), n_cat, replace = TRUE, prob = c(0.4,0.3,0.3))


# Introduce some disagreement
rater2 <- ifelse(runif(n_cat) < 0.8,
                 rater1,
                 sample(c("A","B","C"), n_cat, replace = TRUE))

cat_data <- data.frame(
  rater1 = factor(rater1, levels = c("A","B","C")),
  rater2 = factor(rater2, levels = c("A","B","C"))
)

head(cat_data)


```

2. For continuous data

Here you might want to test is a test gives similar results to a gold standard test or gives similar
results each time you measure the same thing.

```{r}
set.seed(456)
n_cont <- 1000

# Method 1: gold standard (or True measurement or can be 1st measurement)
method1 <- rnorm(n_cont, mean = 10, sd = 2)


# Method 2: noisy measurement with slight bias (or 2nd measurement)
method2 <- method1 + rnorm(n_cont, mean = 0.5, sd = 1)

cont_data <- data.frame(method1, method2)
head(cont_data)


```

# Categorical data
An unweighted Cohen’s kappa of 0.817 indicates “almost perfect” agreement beyond chance between the two raters 
(Landis & Koch scale: 0.81–1.00) .


```{r}


# Unweighted Cohen's kappa for two raters
kappa_res <- kappa2(cat_data, weight = "unweighted")
kappa_res

```

# Continuous data

```{r}


# Lets first reshape to long form to have box plots separate for each method
cont_long <- cont_data %>%
  pivot_longer(
    cols      = c(method1, method2),
    names_to  = "Method",
    values_to = "Value"
  )

# Then do a side-by-side boxplots for easy comparison
ggplot(cont_long, aes(x = Method, y = Value, fill = Method)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    x     = "Method",
    y     = "Values",
    title = "Comparison of Method 1 and Method 2"
  ) +
  theme_minimal() +
  theme(legend.position = "none")


```


```{r, warning=FALSE}
# Tests have a table of summary statistics 

# Compute summary stats
boxplot_summary <- data.frame(
  Pair           = "Method1 vs Method2",
  Median_Method1 = median(cont_data$method1),
  Median_Method2 = median(cont_data$method2),
  IQR_Method1    = IQR(cont_data$method1),
  IQR_Method2    = IQR(cont_data$method2)
)

# Render with gt for nice formatting

boxplot_table <- gt(boxplot_summary) %>%
  tab_header(
    title    = "Box Plot Summary Statistics",
    subtitle = "Median and Interquartile Range (IQR) for Each Pair"
  ) %>%
  cols_label(
    Pair            = "Pair Name",
    Median_Method1  = "Median (Method 1)",
    Median_Method2  = "Median (Method 2)",
    IQR_Method1     = "IQR (Method 1)",
    IQR_Method2     = "IQR (Method 2)"
  ) %>%
  fmt_number(
    columns = vars(Median_Method1, Median_Method2, IQR_Method1, IQR_Method2),
    decimals = 3
  ) %>%
  opt_table_font(font = "Arial", size = "medium") %>%
  opt_align_table_header(align = "center") %>%
  tab_options(
    table.width     = pct(80),
    table.font.size = px(12)
  )

# Print the table
boxplot_table


```

```{r, warning=FALSE}
# Mean absolute error
# Compute Mean Absolute Error (MAE)
mae_val <- mean(abs(cont_data$method1 - cont_data$method2))

mae_table <- tibble(
  Metric = "Mean Absolute Error",
  Value  = mae_val
)

mae_gt <- mae_table %>%
  gt() %>%
  tab_header(
    title    = "MAE Summary",
    subtitle = "Mean Absolute Error Between Method 1 and Method 2"
  ) %>%
  cols_label(
    Metric = "Metric",
    Value  = "MAE"
  ) %>%
  fmt_number(
    columns = vars(Value),
    decimals = 3
  ) %>%
  opt_table_font(font = "Arial", size = "medium") %>%
  opt_align_table_header(align = "center") %>%
  tab_options(
    table.width     = pct(50),
    table.font.size = px(12)
  )

# Print MAE table
mae_gt

```


```{r, warning=FALSE}
# Paired t-test
ttest_res <- t.test(cont_data$method1, cont_data$method2, paired = TRUE)
ttest_tbl  <- tidy(ttest_res) %>%
  select(
    `Mean Difference` = estimate,
    `t Statistic`    = statistic,
    DF               = parameter,
    `p-value`        = p.value,
    `CI Lower`       = conf.low,
    `CI Upper`       = conf.high
  )

ttest_gt <- ttest_tbl %>%
  gt() %>%
  tab_header(
    title    = "Paired t-Test Results",
    subtitle = "Method 1 vs Method 2"
  ) %>%
  fmt_number(
    columns = vars(`Mean Difference`, `t Statistic`, DF, `p-value`, `CI Lower`, `CI Upper`),
    decimals = 3
  ) %>%
  opt_table_font(font = "Arial", size = "medium") %>%
  opt_align_table_header(align = "center") %>%
  tab_options(
    table.width     = pct(80),
    table.font.size = px(12)
  )

# Print Paired t-test table
ttest_gt

```


```{r}
# Mean absolute error plots
# Compute MAE
mae_val <- mean(abs(cont_data$method1 - cont_data$method2))

# Scatter + vertical “error” lines + identity line
ggplot(cont_data, aes(x = method1)) +
  geom_segment(aes(xend = method1, y = method1, yend = method2),
               color = "grey70") +
  geom_point(aes(y = method2), color = "steelblue", size = 2) +
  # 45° reference line
  geom_abline(slope = 1, intercept = 0,
              linetype = "dashed", color = "red") +
  
  labs(
    title = paste0("Method1 vs. Method2 (MAE = ",
                   round(mae_val, 2), ")"),
    x     = "Method1",
    y     = "Method2"
  ) +
  theme_minimal()
```


```{r}
# Bland-Altman plot and limits of agreement

# Compute means and differences
bland_data <- cont_data %>%
  mutate(
    mean_val = (method1 + method2) / 2,
    diff     = method2 - method1
  )

# Calculate bias and limits of agreement
bias      <- mean(bland_data$diff)
sd_diff   <- sd(bland_data$diff)
loa_lower <- bias - 1.96 * sd_diff
loa_upper <- bias + 1.96 * sd_diff

# Print summary
bland_summary <- tibble(
  Bias           = bias,
  SD_of_Differences = sd_diff,
  LoA_Lower      = loa_lower,
  LoA_Upper      = loa_upper
)
print(bland_summary)



```

```{r}
# Bland-Altman plot

ggplot(bland_data, aes(x = mean_val, y = diff)) +
  geom_point(color = "steelblue", size = 2, alpha = 0.7) +
  geom_hline(yintercept = bias,      linetype = "solid",  color = "red") +
  geom_hline(yintercept = loa_lower, linetype = "dashed", color = "darkgrey") +
  geom_hline(yintercept = loa_upper, linetype = "dashed", color = "darkgrey") +
  labs(
    title = "Bland–Altman Plot: Method 1 vs Method 2",
    x     = "Mean of Method 1 and 2",
    y     = "Difference (Method 2 – Method 1)"
  ) +
  annotate(
    "text",
    x = min(bland_data$mean_val),
    y = bias,
    label = paste0("Bias = ", round(bias, 3)),
    hjust = 0, vjust = -1, color = "red"
  ) +
  annotate(
    "text",
    x = min(bland_data$mean_val),
    y = loa_upper,
    label = paste0("LoA upper = ", round(loa_upper, 3)),
    hjust = 0, vjust = -1, color = "darkgrey"
  ) +
  annotate(
    "text",
    x = min(bland_data$mean_val),
    y = loa_lower,
    label = paste0("LoA lower = ", round(loa_lower, 3)),
    hjust = 0, vjust = 1, color = "darkgrey"
  ) +
  theme_minimal()

```


```{r, warning=FALSE}
# Compute Lin’s CCC 
ccc_res <- CCC(cont_data$method1, cont_data$method2)
ccc_val <- ccc_res$rho.c[1]

# Plot 
ggplot(cont_data, aes(x = method1, y = method2)) +
  geom_point(color = "steelblue", size = 2, alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0,
              linetype = "dashed", color = "red") +
  geom_label(
    aes(x = Inf, y = Inf,
        label = paste0("CCC = ", sprintf("%.2f", ccc_val))),
    hjust = 1.1,    
    vjust = 1.1,    
    size  = 4,
    fill  = "white",
    label.size = 0  
  ) +
  labs(
    title = "Concordance Plot: Method 1 vs Method 2",
    x     = "Method 1",
    y     = "Method 2"
  ) +
  theme_minimal()

```

```{r}


```


```{r}


```

```{r}


```


```{r}


```


```{r}


```


```{r}


```

```{r}


```


```{r}


```

```{r}


```


```{r}


```

```{r}


```


```{r}


```


```{r}


```


```{r}


```

```{r}


```


```{r}


```

```{r}


```


```{r}


```

```{r}


```


```{r}


```


```{r}


```


```{r}


```

```{r}


```


```{r}


```

```{r}


```


```{r}


```

```{r}


```


```{r}


```