presentation.Rmd

---
title: "Workshop: Limits of Agreement"
subtitle: "Bland-Altman methods for assessing agreement of clinical measurements"
author: "Sam Gardiner"
institute: "Cell & Molecular Therapies, Royal Prince Alfred Hospital"
date: "14 April 2021"
output:
  xaringan::moon_reader:
    lib_dir: libs
    css: [default, metropolis, metropolis-fonts, tweaks.css]
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
---

```{r setup, include=FALSE}
# Libraries
options(htmltools.dir.version = FALSE)
library(tidyverse)
library(RefManageR)
library(bibtex)
library(ggtext)
library(patchwork)

# Bibliography
BibOptions(check.entries = FALSE,
           bib.style = "authoryear",
           cite.style = "authoryear",
           style = "markdown",
           hyperlink = FALSE,
           dashed = FALSE)
bib <- ReadBib("./references/references.bib", check = FALSE)

# Knitting
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE,
                      dev = "svglite",
                      fig.align = "center")

# Plotting
theme_ba <- theme_classic() +
  theme(plot.subtitle = element_markdown(),
        plot.background = element_rect(fill = "transparent"),
        plot.margin = margin(5,5,5,5, "mm"),
        legend.position = "none")
theme_set(theme_ba)

# Data
pefr_wide <- read_csv("data/pefr_wide.csv")
fix <- read_csv("data/fix_points.csv")
```

```{r functions}
blandize <- function(data, x, y) {
  data %>%
    transmute(magnitude = ({{x}} + {{y}}) / 2,
              difference = {{x}} - {{y}})
}

blandstats <- function(bland_data,
                       magnitude = magnitude,
                       difference = difference,
                       alpha = 0.05) {
  with(bland_data,
       lst(
         bias = mean({{difference}}),
         n = nrow(bland_data),
         sd = sd({{difference}}),
         se = sqrt(var({{difference}}) / n),
         loa.se = sqrt(3 * var({{difference}}) / n),
         stat = qt(alpha / 2, df = n - 1, lower.tail = FALSE),
         bias.upper = bias + se * stat,
         bias.lower = bias - se * stat,
         limit.upper = bias + 1.96 * sd,
         limit.lower = bias - 1.96 * sd,
         limit.upper.ci.upper = limit.upper + loa.se * stat,
         limit.upper.ci.lower = limit.upper - loa.se * stat,
         limit.lower.ci.upper = limit.lower + loa.se * stat,
         limit.lower.ci.lower = limit.lower - loa.se * stat
       ))
}

gg_ba <- function(data, ba_stats) {
  fix_ba_abs <- ggplot(data, aes(magnitude, difference)) +
    geom_point(alpha = 1/2) +
    geom_hline(yintercept = 0, linetype = "dotted") +
    geom_hline(yintercept = ba_stats$bias, linetype = "dashed", colour = "firebrick") +
    geom_hline(yintercept = c(ba_stats$limit.lower, ba_stats$limit.upper),
               linetype = "dashed",
               colour = "dodgerblue")
}
```


# The _Limits of Agreement_ method

```{r results = "asis"}
NoCite(bib, "Bland1986")
print(bib[key = "Bland1986"])
```

- The 29th most-cited paper of all time! `r Citep(bib, "Noorden2014")`
- Still the gold standard for measuring agreement between continuous clinical measurements.
- Simple enough to do by hand (in Excel) if needed, but also available in almost all statistical software: R, GraphPad Prism, SAS etc.

---
class: middle, center, inverse

# Agreement

---

# Agreement

- It is often useful to compare two methods of measuring some clinical parameter. For example:
  - One-stage vs. chromogenic FIX activity 
  - Axillary vs tympanic temperature 
  - NucleoCounter vs CELL-DYN cell counts
- If the two methods "agree" (within clinically meaningful limits), you might be able to retire the more expensive, more laborious or otherwise less convenient method.

---

# Not agreement

## Correlation

- What about $r$, the standard (Pearson product-moment) correlation coefficient? 
--

  - $r$ measures linear correlation between two variables, not agreement.
  - Two measurement methods can be perfectly linearly correlated, but not agree.
  - Being correlated just means that two variables tend to go up or down together.
  - Correlation $r$ is a function of the variability of the data: two variables that cover a wide range will have larger $r$ than similar variables which cover a small range, even if the degree of agreement is the same.

---
class: middle

# Not agreement

## Perfectly correlated, but not in agreement

```{r fig.height = 5}
set.seed(1987)

eg1_data <- tibble(
  x = runif(50),
  y1 = x * 1.5,
  y2 = x + 0.25
)

eg1_1_stats <- with(eg1_data, cor.test(x, y1))
eg1_2_stats <- with(eg1_data, cor.test(x, y2))

eg1_1 <- ggplot(eg1_data, aes(x, y1)) +
  geom_point(alpha = 1/2) +
  coord_equal() +
  scale_x_continuous(breaks = c(0, 0.5, 1)) +
  labs(subtitle = str_glue("_r_ = {eg1_1_stats$estimate}; _p_ = {format.pval(eg1_1_stats$p.value)}"),
       x = "Method 1",
       y = "Method 2")

eg1_2 <- ggplot(eg1_data, aes(x, y2)) +
  geom_point(alpha = 1/2) +
  coord_equal() +
  scale_x_continuous(breaks = c(0, 0.5, 1)) +
  scale_y_continuous(breaks = c(0, 0.5, 1, 1.5), limits = c(0, 1.5)) +
  labs(subtitle = str_glue("_r_ = {eg1_2_stats$estimate}; _p_ = {format.pval(eg1_2_stats$p.value)}"),
       x = "Method 1",
       y = "Method 2")

eg1_1 + eg1_2
```

---

# Not agreement

## Perfectly correlated, but not in agreement

```{r fig.height=5}
line_equal <- geom_abline(slope = 1, linetype = "dashed", colour = "firebrick")
eg1_1 + line_equal + eg1_2 + line_equal
```

---

# Not agreement

## Even worse:

```{r fig.height = 4, width = 4}
anscombe_long <- anscombe %>%
  pivot_longer(everything(),
               names_to = c("dimension", "set"),
               names_pattern = "([xy])([1234])") %>%
  pivot_wider(names_from = dimension, values_from = value) %>%
  unnest(cols = c(x, y))

ggplot(anscombe_long, aes(x, y)) +
  geom_point() +
  facet_wrap(vars(set), nrow = 2, ncol = 2) +
  geom_smooth(method = "lm", alpha = 1/2, se = FALSE) +
  annotate("text", x = 5, y = 12.5, label = "r = 0.816", hjust = 0) +
  theme(strip.background = element_blank(), strip.text = element_blank()) +
  scale_x_continuous(breaks = c(5, 10, 15, 20), limits = c(2.5, 20)) +
  scale_y_continuous(breaks = c(0, 5, 10, 15), limits = c(2.5, 15)) +
  coord_equal()
```

Data: `r Citet(bib, "Anscombe1973")`

---

# Not agreement

## Calibration

- Is measuring agreement the same as calibration?
--

  - Generally, **no**.
  - Calibration compares a single method against a ground truth. 
  - Agreement compares two imperfect methods (which are assumed to have measurement error) with each other.
  - If the "ground truth" isn't particularly precise, agreement and calibration may be the same concept.

---

# Not quite agreement

## Repeatability

- Repeatability is a closely-related concept: if a measurement method agrees with itself over repeated measurements, it is _repeatable_.
- The Bland-Altman _Limits of Agreement_ methods work well for assessing repeatability, as well.

---
class: middle, center, inverse

# Assessing agreement

---

# Eyeball the data

.pull-left[
- Plot:
  - each method against the other
  - the line of equality (the line with slope 1, passing through the origin)
- Do the observations lie approximately along the line of equality?
- Are there any obvious systematic differences?
]

.pull-right[ 
```{r, fig.width = 4, fig.height = 4}
ggplot(pefr_wide, aes(Wright1, Mini1)) +
  geom_point() +
  scale_x_continuous(limits = c(0, 800)) +
  scale_y_continuous(limits = c(0, 800)) +
  line_equal +
  coord_equal() +
  labs(x = "PEFR by large meter (L/min)",
       y = "PEFR by mini meter (L/min)")
```

PEFR: Peak expiratory flow rate, a measure of lung function.
]

---

# The _Limits of Agreement_ method


0. Decide on a clinically acceptable threshold of agreement. Use clinical reasoning or published evidence. For example, you might consider methods in agreement if they are within
  - 5mmHg for blood pressure
  - 0.1 for blood pH
  - 5% clotting activity for a FIX assay
0. Visualise the magnitude of the measurements against the difference of the two methods.
  - magnitude: estimate with the mean of the two methods
  - difference: subtract one method from the other
0. Find the bias and its standard deviation. The bias is the average differences between methods.
0. Find the limits of agreement:
  - $\text{Limits} = \text{Bias} \pm \text{SD(Bias)} \times 1.96$
0. Critically appraise: 
  - are there systematic differences between the methods?
  - is the scale of the difference the same over the range of the measurements?
  - are the 95% limits of agreement within the predefined clinically meaningful threshold?

---

# Why 1.96?

```{r fig.height = 5}
ggplot() +
  scale_x_continuous(limits = c(-3, 3)) +
  scale_y_continuous(limits = c(0, 0.5), expand = c(0, 0)) +
  stat_function(geom = "line", fun = dnorm) +
  stat_function(geom = "area", fun = dnorm, xlim = c(-1.96, 1.96), fill = "firebrick", alpha = 1/2 ) +
  geom_vline(xintercept = c(-1.96, 1.96), linetype = "dashed", colour = "firebrick") +
  labs(x = "Standardised difference", y = "Density")
```


---

# Anatomy of a Bland-Altman plot

.pull-left[
```{r}
pefr_wide %>%
  select(Subject, "Large meter" = Wright1, "Mini meter" = Mini1) %>%
  head(10) %>%
  knitr::kable()
```
]

.pull-right[
## Example dataset:
Comparison of **p**eak **e**xpiratory **f**low **r**ate (PEFR in L/minute) by a large Wright peak flow meter and a mini Wright meter, measure in the same subject. `r Citet(bib, "Bland1986")`.
]

---

# Anatomy of a Bland-Altman plot

```{r fig.height = 5}
pefr_bland <- blandize(pefr_wide, Wright1, Mini1)
pefr_stats <- blandstats(pefr_bland)

anatomy0 <- ggplot(pefr_bland, aes(magnitude, difference)) +
  scale_y_continuous(limits = c(-120, 120)) +
  labs(x = "Magnitude: Mean of Large and Mini (L/min)",
       y = "Difference: Large - Mini (L/min)")

anatomy0
```

---
# Anatomy of a Bland-Altman plot

```{r fig.height = 5}
anatomy1 <- anatomy0 + 
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dotted")

anatomy1
```

---

# Anatomy of a Bland-Altman plot

```{r fig.height = 5}
anatomy2 <- anatomy1 +
  geom_hline(yintercept = pefr_stats$bias, linetype = "solid", colour = "firebrick")

anatomy2
```

---

# Anatomy of a Bland-Altman plot

```{r fig.height = 5}
anatomy3 <- anatomy2 +
  geom_hline(yintercept = c(pefr_stats$limit.upper, pefr_stats$limit.lower),
             linetype = "dashed",
             colour = "dodgerblue")
  
anatomy3
```

---
# Anatomy of a Bland-Altman plot

```{r fig.height = 5}
anatomy4 <- anatomy3 +
  annotate("ribbon",
           x = c(-Inf, Inf),
           ymin = pefr_stats$bias.lower,
           ymax = pefr_stats$bias.upper,
           fill = "firebrick",
           alpha = 1/10)
anatomy4
```

---
# Anatomy of a Bland-Altman plot

```{r fig.height = 5}
anatomy5 <- anatomy4 + 
    annotate("ribbon",
           x = c(-Inf, Inf),
           ymin = pefr_stats$limit.upper.ci.lower,
           ymax = pefr_stats$limit.upper.ci.upper,
           fill = "dodgerblue",
           alpha = 1/10) +
      annotate("ribbon",
           x = c(-Inf, Inf),
           ymin = pefr_stats$limit.lower.ci.lower,
           ymax = pefr_stats$limit.lower.ci.upper,
           fill = "dodgerblue",
           alpha = 1/10)

anatomy5
```
---

# Assessing agreement

## Assumptions

The PEFR example relies on some assumptions about the data:
- That there is no systematic change to the degree of agreement over the range of the measurements.
- That the measurement error (the difference between the two measurements) is normally distributed.
  - The limits of agreement and confidence intervals rely on this assumption to be accurate, but should be OK with other distributions as long as the sample size isn't tiny.

---

# Assessing agreement

## What if there _is_ a systematic difference?

Bland and Altman suggest two remedies:

- Working with percentage difference instead of absolute difference
- Log-transforming your data.

---

# Systematic difference

.pull-left[
- Paired plasma samples from the SPK-9001-101 participants were measured by a chromogenic FIX activity assay at the trial central laboratory, and by a one-stage assay at each site's local laboratory.
- $n=15$ participants with $147$ measurements.
]

.pull-right[

```{r fig.height = 5, fig.width = 4}
participants <- unique(fix$fill)
scale_part <- scale_colour_manual(values = set_names(participants))

fix_plot <- ggplot(fix, aes(x, y, colour = fill)) +
coord_equal() +
  scale_x_continuous(limits = c(0, NA), breaks = scales::breaks_width(10)) +
  scale_y_continuous(limits = c(0, 100), breaks = scales::breaks_width(10)) +
  theme(legend.position = "none") +
  scale_part +
  geom_point(alpha = 1/2) +
  labs(x = "Central Laboratory FIX:C Value (%)",
       y = "Local Laboratory FIX:C Value (%)")
  
fix_plot
```
.center[`r Citet(bib, "Robinson2021", .opts = list(max.names = 3, longnamesfirst = FALSE))`]]

---

# Systematic difference

```{r fig.height=6}
fix_plot +
  line_equal +
  geom_smooth(method = "lm",
              se = FALSE,
              aes(group = NA))
```

---

# Systematic difference

## Plotting the absolute difference would be a mistake

```{r fig.height = 5}

fix_bland <- blandize(fix, x, y) %>%
  cbind(fix)

fix_stats <- blandstats(fix_bland)

fix_ba_abs <- gg_ba(fix_bland, fix_stats) +
  labs(x = "Magnitude (average FIX activity by both methods)",
       y = "Difference (central lab - local lab result)") +
  geom_point(aes(colour = fill)) +
  scale_part

fix_ba_abs
```

---

# Systematic difference

## Plotting the absolute difference would be a mistake

```{r fig.height = 5}
fix_ba_abs +
  geom_smooth(method = "lm")
```

---

# Systematic difference

## Plot the percentage or ratio difference, instead

```{r fig.height = 5}
fix_pct_bland <- blandize(fix, x, y) %>%
  mutate(difference = abs(difference) / magnitude) %>%
  cbind(fix)

fix_pct_stats <- blandstats(fix_pct_bland)

fix_log_ba <- gg_ba(fix_pct_bland, fix_pct_stats) +
  geom_point(aes(colour = fill)) +
  labs(x = "Magnitude (average FIX %)",
       y = "Ratio of difference to magnitude") +
  scale_part

fix_log_ba 
```

---

# More to learn...

- Repeated measures versions where each method is used to measure a sample or individual multiple times.
  - Paired or unpaired? 
  - Constant underlying true value, or time-dependent? `r Citep(bib, "Bland2007")`
- Applications to transcriptomics (gene expression) data: the MA plot `r Citep(bib, "Dudoit2002")`.

---

# References and further reading

```{r results="asis"}
PrintBibliography(bib, end = 5)
```
---

# References and further reading

```{r results="asis"}
PrintBibliography(bib, start = 6)
```