3.7 Agreement and reliability

Two methods that measure the same thing can be highly correlated and still disagree badly. Agreement asks how close the values are, not whether they rise together, and reliability asks how much of the variation in a measurement is real difference between subjects rather than noise.

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes

from statly import style

style.use()
rng = np.random.default_rng(1)

# Simulated: blood pressure of 80 patients measured by a cuff (the old method) and a new wrist device.
N = 80
true_bp = rng.normal(135, 18, N)
cuff = true_bp + rng.normal(0, 5, N)
wrist = (
    0.75 * true_bp + 38 + rng.normal(0, 5, N)
)  # reads high for low pressures, low for high ones

3.7.1 Correlation is not agreement

A new device is compared with the standard one on the same patients. A correlation of 0.9 sounds reassuring, but correlation only says that the two rise and fall together. A device that always reads 10 too high has a correlation of 1 with the standard and is useless without correction.

fig, axes = style.subplots(2, height=3.4, width=10)
for ax, other, title in (
    (axes[0], wrist, "Wrist device"),
    (axes[1], cuff + 10, "Always 10 too high"),
):
    ax.scatter(cuff, other, s=12, color=style.BLUE, alpha=0.7)
    line = np.array([90, 185])
    ax.plot(line, line, color=style.INK, ls="--", lw=1)
    style.note(ax, f"r = {np.corrcoef(cuff, other)[0, 1]:.2f}", x=0.03, y=0.95, color=style.INK)
    ax.set(title=title, xlabel="cuff", ylabel="other method", aspect="equal")
fig.tight_layout()
plt.show()

Figure 3.7.1: Blood pressure of 80 simulated patients, measured with a cuff and with a new wrist device. Left: the two correlate strongly, but the points do not lie on the line of perfect agreement (dashed): the device reads too high for low pressures and too low for high ones. Right: an extreme case, a device that is always 10 too high, with a correlation of 1 and no agreement at all.

3.7.2 The Bland-Altman plot

The standard picture of agreement plots, for each subject, the difference between the two methods against their mean. The average difference is the bias, and the limits of agreement, the bias ± 1.96 standard deviations of the differences, give the range in which 95% of differences fall:

$$\text{limits of agreement} = \bar d \pm 1.96\, s_d$$

Whether those limits are acceptable is a clinical question, not a statistical one: would a difference of that size change a decision?

difference = wrist - cuff
mean_of_two = (wrist + cuff) / 2
bias, sd = difference.mean(), difference.std(ddof=1)

fig, ax = plt.subplots(figsize=(8, 3.4))
ax.scatter(mean_of_two, difference, s=14, color=style.BLUE, alpha=0.7)
for value, text in (
    (bias, f"bias {bias:+.1f}"),
    (bias + 1.96 * sd, f"upper limit {bias + 1.96 * sd:+.1f}"),
    (bias - 1.96 * sd, f"lower limit {bias - 1.96 * sd:+.1f}"),
):
    ax.axhline(
        value,
        color=style.CORAL if value == bias else style.INK,
        ls="-" if value == bias else "--",
        lw=1.3,
    )
    ax.text(mean_of_two.max() + 1, value, text, va="center", fontsize=8.5)
trend = np.polyfit(mean_of_two, difference, 1)
xs = np.array([mean_of_two.min(), mean_of_two.max()])
ax.plot(xs, np.polyval(trend, xs), color=style.GREY, lw=1)
ax.axhline(0, color=style.GREY, lw=0.6)
ax.set(
    xlabel="mean of the two methods", ylabel="wrist − cuff", xlim=(None, mean_of_two.max() + 30)
)
plt.show()

Figure 3.7.2: Bland-Altman plot for the wrist device against the cuff. Each dot is a patient: the difference (wrist − cuff) against the mean of the two. The middle line is the bias, the outer lines the 95% limits of agreement. The differences also slope downwards: the device’s error depends on the level, overestimating low pressures and underestimating high ones, which a single bias hides.

3.7.3 Agreement on categories: Cohen’s kappa

When two raters sort the same cases into categories (two pathologists grading tumours, two coders labelling texts), the share of cases they agree on is flattering, because some agreement happens by chance. Cohen’s kappa compares the observed agreement $p_o$ with the agreement $p_e$ expected if the two rated independently with their own frequencies:

$$\kappa = \frac{p_o - p_e}{1 - p_e}$$

κ = 1 is perfect agreement, 0 is no better than chance. Rough labels: below 0.4 poor, 0.4 to 0.6 moderate, 0.6 to 0.8 good, above 0.8 very good.

table = np.array([[156, 8], [12, 24]])  # rows: rater A benign/suspicious; columns: rater B
total = table.sum()
p_observed = np.trace(table) / total
p_expected = (table.sum(axis=1) / total) @ (table.sum(axis=0) / total)
kappa = (p_observed - p_expected) / (1 - p_expected)

fig, axes = style.subplots(2, height=3.2, width=10)
axes[0].imshow(table, cmap="Blues")
for i in range(2):
    for j in range(2):
        axes[0].text(j, i, table[i, j], ha="center", va="center", fontsize=12)
axes[0].set_xticks([0, 1], ["benign", "suspicious"])
axes[0].set_yticks([0, 1], ["benign", "suspicious"])
axes[0].set(xlabel="rater B", ylabel="rater A")
axes[0].grid(False)
axes[1].barh(
    ["expected by chance", "observed"],
    [p_expected, p_observed],
    color=[style.GREY, style.BLUE],
    alpha=0.7,
)
axes[1].axvline(1, color=style.INK, lw=0.8)
for i, value in enumerate((p_expected, p_observed)):
    axes[1].text(value + 0.01, i, f"{value:.0%}", va="center", fontsize=9)
style.note(
    axes[1],
    f"κ = ({p_observed:.2f} − {p_expected:.2f}) / (1 − {p_expected:.2f}) = {kappa:.2f}",
    x=0.03,
    y=0.95,
    color=style.INK,
)
axes[1].xaxis.set_major_formatter(style.PERCENT)
axes[1].set(xlim=(0, 1.15), xlabel="agreement")
fig.tight_layout()
plt.show()

Figure 3.7.3: Two raters classify 200 skin images as benign or suspicious (simulated). Left: the table of their ratings. Right: the observed agreement is 90%, which sounds excellent; but because most images are benign, two raters guessing with the same frequencies would already agree 72% of the time. Kappa measures how much of the room above chance they achieve: 0.65, good but far from the 90% suggests.

For ordered categories, weighted kappa gives partial credit to near misses. With more than two raters, Fleiss’ kappa generalises the idea.

3.7.4 Reliability: the intraclass correlation

How reliable is a measurement? If the same subjects are measured several times (by different raters, on different days), part of the variation is real difference between subjects and part is noise. The intraclass correlation is the share that is real:

$$\text{ICC} = \frac{\sigma^2_{\text{subjects}}}{\sigma^2_{\text{subjects}} + \sigma^2_{\text{error}}}$$

It is the same quantity as the ICC of mixed models (Chapter 4.7), seen as a property of a measurement. An ICC of 0.9 means 90% of the spread in single measurements reflects true differences.

def icc_one_way(scores: np.ndarray) -> float:
    """One-way random-effects ICC from a subjects × raters table."""
    n, k = scores.shape
    between = k * np.sum((scores.mean(axis=1) - scores.mean()) ** 2) / (n - 1)
    within = np.sum((scores - scores.mean(axis=1, keepdims=True)) ** 2) / (n * (k - 1))
    return (between - within) / (between + (k - 1) * within)


subjects = rng.normal(50, 10, 15)
order = np.argsort(subjects)
fig, axes = style.subplots(2, height=3.2, width=10)
for ax, noise, title in ((axes[0], 3, "Reliable"), (axes[1], 10, "Unreliable")):
    scores = subjects[:, None] + rng.normal(0, noise, (15, 3))
    for r, color in enumerate((style.BLUE, style.CORAL, style.GREEN)):
        ax.scatter(np.arange(15), scores[order, r], s=22, color=color, alpha=0.8)
    style.note(ax, f"ICC = {icc_one_way(scores):.2f}", x=0.03, y=0.95, color=style.INK)
    ax.set(
        title=title, xlabel="subject (sorted by true value)", ylabel="measurement", xticks=[]
    )
fig.tight_layout()
plt.show()

Figure 3.7.4: Fifteen subjects, each measured by three raters (dots of one colour per rater). Left: raters agree closely relative to how much subjects differ; the ICC is high. Right: the same subjects with noisier raters; the measurements of one subject spread as much as the subjects themselves, and the ICC is low.

Low reliability has a quiet cost: noise in a predictor weakens every correlation and regression slope involving it (attenuation). The observed correlation between two measures is at most the square root of the product of their reliabilities.

3.7.5 Questionnaires: Cronbach’s alpha

A questionnaire scale adds up several items meant to measure the same thing. Cronbach’s alpha estimates the reliability of the sum from how strongly the items correlate:

$$\alpha = \frac{k}{k - 1}\Big(1 - \frac{\sum_{i} s_i^2}{s_{\text{total}}^2}\Big)$$

with k items, item variances $s_i^2$ and the variance of the total score $s^2_{\text{total}}$. More items, or more strongly related items, give a higher alpha:

def cronbach_alpha(items: np.ndarray) -> float:
    """Cronbach's alpha from a respondents × items table."""
    k = items.shape[1]
    return k / (k - 1) * (1 - items.var(axis=0, ddof=1).sum() / items.sum(axis=1).var(ddof=1))


def scale(k: int, r: float, n: int = 500) -> np.ndarray:
    """k items that each correlate r with one another, through one common trait."""
    trait = rng.normal(size=(n, 1))
    return np.sqrt(r) * trait + np.sqrt(1 - r) * rng.normal(size=(n, k))


ITEMS = np.arange(2, 21)
fig, ax = plt.subplots(figsize=(8, 3.2))
for r, color in ((0.2, style.CORAL), (0.4, style.GOLD), (0.6, style.BLUE)):
    ax.plot(
        ITEMS,
        [cronbach_alpha(scale(k, r)) for k in ITEMS],
        "o-",
        color=color,
        ms=4,
        lw=2,
        label=f"items correlate {r}",
    )
ax.axhline(0.7, color=style.INK, ls="--", lw=1)
ax.legend(fontsize=8, loc="lower right")
ax.set(xlabel="number of items", ylabel="Cronbach's alpha", ylim=(0, 1))
plt.show()

Figure 3.7.5: Cronbach’s alpha of simulated scales, by number of items, for items that correlate 0.2, 0.4 and 0.6 with each other (each point from 500 respondents). Even weakly related items add up to a reliable scale if there are enough of them. A very high alpha can also mean the items are near-copies of each other, which adds length without adding information.

3.7.6 In one picture

• Correlation is not agreement: compare methods with a Bland-Altman plot; bias and limits of agreement $\bar d \pm 1.96 s_d$; look for a trend in the differences.
• For categories, Cohen’s kappa $(p_o - p_e)/(1 - p_e)$ corrects agreement for chance.
• The ICC is the share of variation that is real difference between subjects.
• Unreliable measures attenuate correlations and slopes.
• Cronbach’s alpha grows with the number of items and their correlation.

Next: checking the assumptions behind the tests.