Submission. Knit your Rmd file into HTML and submit it via Brightspace. Include a link to your git repository (GitHub or git.cs.dal.ca) as a comment on your submission.

0 Setup

This lab uses a few different/new libraries to handle signal processing and to cover all 3 of our analysis paradigms for signal data: a time-domain model, a frequency-domain model, and a state-space model.

install.packages(c(
  "signal",      # digital filters (Butterworth band-pass / notch)
  "forecast",    # ARIMA time-series models (auto.arima)
  "depmixS4"     # Gaussian (continuous) hidden Markov models (state-space)
))

If the block below runs without error you are ready to go. Instead of attaching like normal (using library()) we are just checking if signal and forecast are both installed (using requiredNamespace). This is these libraries contain functions (like filter()) that collide with each other and dplyr. Instead we will call both explicitly in the code with the signal:: / forecast:: prefix where needed.

library(tidyverse)
library(here)
library(fs)
library(depmixS4)
library(yardstick)
library(scales)
requireNamespace('signal')
requireNamespace('forecast')

1 Introduction

The electrocardiogram (ECG) records the small voltages produced as the heart’s chambers depolarise to contract (systole) and repolarise to relax (diastole).

A single normal beat has a very distinctive/standard shape - the P wave (atrial contracting), the QRS complex (ventricular contracting), and the T wave (ventricular relaxing). The time between the main R spikes in successive QRS peaks is the R-R interval.

When a heart doesn’t beat in a nice regular rhythm we call this a cardiac arrhythmia. Atrial fibrillation (AF) is the most common sustained arrhythmia and a major cause of stroke and heart failure. As AF is often asymptomatic, a lot of work has gone into identifying it accurately from ECGs. Typically, AF is recognised in an ECG by the disappearance of the P wave (as the atria are contracting chaotically rather than in an organised wave) and irregular spacing and timing of the QRS complex (as the ventricles contract at irregular intervals).

In this lab we will build a different AF detector for each signal data analysis paradigm and then compare them against one another:

  • a time-domain model that asks how predictable the next R-R interval is (using ARIMA via forecast::auto.arima, the autoregressive / moving-average / integrated family from the lecture),
  • a frequency-domain model that asks how the signal’s energy is spread across frequencies (using the Fourier transform), and
  • a state-space model that treats the rhythm as a sequence of hidden states (a hidden Markov model).

2 The Data

We use a small prepared subset of the MIT-BIH Atrial Fibrillation Database (the same public PhysioNet database Babaeizadeh et al. validated on): two-channel Holter ECGs sampled at 250 Hz with expert beat and rhythm annotations. Reading the raw PhysioNet files needs some slightly awkward platform-specific libraries, so I’ve simplified things a little (including extract R-R intervals) and converted to a more usable format for you. We’ve also simplified away the “temporal dependence” complexity of how you cut-up time-series data to avoid misleading yourself.

As in the earlier labs, we download a prepared archive once and cache it.

data_dir  <- here("data")
cache_dir <- fs::path(data_dir, "cache")
fs::dir_create(cache_dir)

data_url <- "https://maguire-lab.github.io/health_data_science_research_2026/static_files/practicals/ecg_afib.zip"
zip_path <- fs::path(cache_dir, "ecg_afib.zip")
ecg_root <- fs::path(data_dir, "ecg_afib")
if (!fs::dir_exists(ecg_root)) {
  if (!fs::file_exists(zip_path)) {
    download.file(data_url, zip_path, mode = "wb", quiet = TRUE)
  }
  unzip(zip_path, exdir = data_dir)
}

fs::dir_ls(ecg_root)
## /home/fin/Documents/teaching/2025_2026/health_data_science_research_2026/data/ecg_afib/rr_intervals.csv
## /home/fin/Documents/teaching/2025_2026/health_data_science_research_2026/data/ecg_afib/waveforms.csv

The archive contains 2 tables:

  • waveforms.csv - a handful of 10-second raw ECG snippets (2500 samples each at 250 Hz), labelled AF or Normal where record indicates an individual patient.
  • rr_intervals.csv - sequences of consecutive R-R intervals split into fixed-length windows (128 beats each, roughly two minutes) and labelled by rhythm.
waveforms <- read_csv(fs::path(ecg_root, "waveforms.csv")) |>
  mutate(rhythm = factor(rhythm, levels = c("AF", "Normal")))

rr <- read_csv(fs::path(ecg_root, "rr_intervals.csv")) |>
  mutate(rhythm = factor(rhythm, levels = c("AF", "Normal")))

glimpse(waveforms)
## Rows: 30,000
## Columns: 5
## $ record  <chr> "04015", "04015", "04015", "04015", "04015", "04015", "04015",…
## $ segment <chr> "04015_Normal", "04015_Normal", "04015_Normal", "04015_Normal"…
## $ rhythm  <fct> Normal, Normal, Normal, Normal, Normal, Normal, Normal, Normal…
## $ time    <dbl> 0.000, 0.004, 0.008, 0.012, 0.016, 0.020, 0.024, 0.028, 0.032,…
## $ ecg     <dbl> -0.115, -0.090, -0.100, -0.130, -0.160, -0.215, -0.230, -0.285…
glimpse(rr)
## Rows: 55,560
## Columns: 5
## $ record <chr> "04015", "04015", "04015", "04015", "04015", "04015", "04015", …
## $ window <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ rhythm <fct> Normal, Normal, Normal, Normal, Normal, Normal, Normal, Normal,…
## $ beat   <dbl> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ rr     <dbl> 0.856, 0.840, 0.836, 0.820, 0.828, 0.824, 0.780, 0.844, 0.824, …

Q1 (2 marks): Write a short data dictionary for waveforms and rr with your best guess explanation of the values and format of these datasets.

3 Train / Test Split

Before any cleaning or exploration we split the data, so that nothing we decide later can be influenced by the test set.

We split on whole records and keep the train and test data in completely separate objects (rr_train / rr_test and waveforms_train / waveforms_test).

set.seed(GLOBAL_SEED)

all_records  <- rr |> distinct(record) |> pull(record)
test_records <- sample(all_records, size = round(0.3 * length(all_records)))

rr_train <- rr |> filter(!record %in% test_records)
rr_test  <- rr |> filter(record %in% test_records)

waveforms_train <- waveforms |> filter(!record %in% test_records)
waveforms_test  <- waveforms |> filter(record %in% test_records)

bind_rows(
  rr_train |> distinct(record, window, rhythm) |> mutate(split = "train"),
  rr_test  |> distinct(record, window, rhythm) |> mutate(split = "test")
) |>
  count(split, rhythm)

Q2 (2 marks): Why do we split on the record instead of the segment? What could the models end up learning to recognise instead of AF if we split by segment?

4 Cleaning - Frequency Filtering

Raw ECG is rarely analysis-ready. Two classic problems are baseline wander (slow drift from breathing and electrode movement, below the heartbeat frequency) and mains interference(a sharp hum at the power-line frequency - 60 Hz as these are US recordings - table of national frequencies).

The lecture covered two standard fixes: a band-pass filter keeps only frequencies likely to be connected to heart bioelectical signals, and a notch (band-stop) filter removes narrow bands of interference like the mains frequency. The the actual way these filters are implemented varies a lot (with different nuances on how they flatten/amplify things). In this lab we implement the two filters using something called a Butterworth filter from the signal and apply them using filtfilt() (which filters forwards and backwards so features are not shifted in time).

fs <- 250

# keep things in conceivable "normal heart rate" range (0.5 Hz to 40 Hz) to smooth out high and low frequency noise 
bandpass <- signal::butter(2, c(0.5, 40) / (fs / 2), type = "pass")

# cut out 59-61Hz signal from the electrical grid
notch    <- signal::butter(2, c(59, 61) / (fs / 2), type = "stop")

example_wave <- waveforms_train |>
  filter(rhythm == "Normal") |>
  filter(segment == first(segment))

cleaned <- example_wave |>
  mutate(
    bp_filtered = signal::filtfilt(bandpass, ecg),
    filtered    = signal::filtfilt(notch, bp_filtered)
  )

cleaned |>
  filter(time <= 4) |>
  dplyr::select(time, Raw = ecg, Filtered = filtered) |>   # dplyr::select - MASS also exports select()
  pivot_longer(c(Raw, Filtered), names_to = "stage", values_to = "ecg") |>
  mutate(stage = factor(stage, levels = c("Raw", "Filtered"))) |>
  ggplot(aes(time, ecg)) +
  geom_line(linewidth = 0.3) +
  facet_wrap(~stage, ncol = 1) +
  labs(x = "Time (s)", y = "ECG (mV)",
       title = "Before and after band-pass + notch filtering")

As we can see, this cleans things up a little (although isn’t perfect) so now let’s apply the same two filters to every waveform in both the training and test sets.

We keep the original trace as ecg_raw and overwrite ecg with the cleaned signal.

clean_waveforms <- function(df) {
  df |>
    group_by(record, segment) |>
    arrange(time, .by_group = TRUE) |>
    mutate(
      ecg_raw = ecg,
      ecg     = signal::filtfilt(notch, signal::filtfilt(bandpass, ecg))
    ) |>
    ungroup()
}

waveforms_train <- clean_waveforms(waveforms_train)
waveforms_test  <- clean_waveforms(waveforms_test)

glimpse(waveforms_train)
## Rows: 20,000
## Columns: 6
## $ record  <chr> "04746", "04746", "04746", "04746", "04746", "04746", "04746",…
## $ segment <chr> "04746_AF", "04746_AF", "04746_AF", "04746_AF", "04746_AF", "0…
## $ rhythm  <fct> AF, AF, AF, AF, AF, AF, AF, AF, AF, AF, AF, AF, AF, AF, AF, AF…
## $ time    <dbl> 0.000, 0.004, 0.008, 0.012, 0.016, 0.020, 0.024, 0.028, 0.032,…
## $ ecg     <dbl> 0.07331563, 0.12585480, 0.15472892, 0.16593037, 0.17014741, 0.…
## $ ecg_raw <dbl> 0.190, 0.195, 0.200, 0.220, 0.225, 0.225, 0.240, 0.280, 0.280,…

Q3 (1 marks): As part of these filters we are dividing the sampling rate in half. What does fs/2 represent in terms of the signal we’ve captured?

5 Exploratory Data Analysis

Remember, we only explore the training set!

5.1 Class balance and a first look

rr_train |>
  distinct(record, window, rhythm) |>
  count(rhythm) |>
  ggplot(aes(rhythm, n, fill = rhythm)) +
  geom_col(show.legend = FALSE) +
  scale_fill_manual(values = c("AF" = "firebrick", "Normal" = "grey60")) +
  labs(x = NULL, y = "Number of windows (train)", title = "Window class balance")

A normal and an AF waveform snippet side by side - note the slightly more regular spacing in the Normal vs AF signal; the p-wave feature is less obviously different in this case:

waveforms_train |>
  filter(time <= 4) |>
  group_by(rhythm) |>
  filter(segment == first(segment)) |>
  ungroup() |>
  ggplot(aes(time, ecg)) +
  geom_line(linewidth = 0.3) +
  facet_wrap(~rhythm, ncol = 1) +
  labs(x = "Time (s)", y = "ECG (mV)", title = "Normal vs AF (training examples)")

5.2 Heart-rate-variability features

A few standard summaries of each window’s R-R series:

window_features <- rr_train |>
  group_by(record, window, rhythm) |>
  summarise(
    mean_rr = mean(rr),                       # average interval (1 / mean_rr = mean rate)
    sdnn    = sd(rr),                          # overall variability
    .groups = "drop"
  )

window_features |>
  pivot_longer(c(mean_rr, sdnn),
               names_to = "feature", values_to = "value") |>
  ggplot(aes(rhythm, value, fill = rhythm)) +
  geom_boxplot(show.legend = FALSE) +
  facet_wrap(~feature, scales = "free_y") +
  scale_fill_manual(values = c("AF" = "firebrick", "Normal" = "grey60")) +
  labs(x = NULL, y = NULL, title = "HRV features by rhythm (train)")

Q4 (2 marks): Explain what this plot tells you about the R-R intervals in Normal vs AF samples.

6 Time Domain Model (ARIMA)

To detect AF with a time-domain approach we are going to work out how predictable is the next interval given the preceding one(s) across each window?. A regular sinus rhythm is highly predictable, so our model should be able to forecast the next intervals well; whereas an AF window is close to noise, so the forecast error is large. This means we can use our model’s forecast error as the AF score.

We’ll use an ARIMA model to do this. ARIMA predicts the next value from a linear combination of recent values (the AR - autoregressive part) and recent forecast errors (the MA - moving average part), after differencing the series to make it stationary (the I part).
Specifically, we are going to use forecast::auto.arima() to automatically select (hyperparameter optimisation) how many previous values to use for each part of this model (also known as the “order” of the model component).

For each window we treat the R-R series as a tachogram (interval value indexed by beat number), fit ARIMA on the first 80% of beats, forecast the rest, and record the mean absolute error.

arima_error <- function(window_df) {
  y   <- window_df$rr
  h   <- max(5, round(0.2 * length(y)))   # forecast horizon = last 20% of beats
  cut <- length(y) - h

  fit      <- forecast::auto.arima(y[seq_len(cut)])
  forecast <- forecast::forecast(fit, h = h)
  mean(abs(as.numeric(forecast$mean) - y[(cut + 1):length(y)]))   # forecast MAE
}

Here is the idea on one window of each rhythm - the red forecast tracks the regular rhythm closely but cannot keep up with AF:

demo <- rr_train |>
  group_by(rhythm) |>
  filter(window == last(window)) |>
  ungroup()

demo_forecast <- demo |>
  group_by(rhythm) |>
  group_modify(~ {
    y <- .x$rr
    h <- max(5, round(0.2 * length(y)))
    cut <- length(y) - h
    fit <- forecast::auto.arima(y[seq_len(cut)])
    fc  <- forecast::forecast(fit, h = h)
    tibble(beat     = seq_along(y),
           rr       = y,
           forecast = c(rep(NA_real_, cut), as.numeric(fc$mean)))
  }) |>
  ungroup()

# forecast error (mean absolute error over the forecast horizon) for each rhythm
demo_error <- demo_forecast |>
  filter(!is.na(forecast)) |>
  group_by(rhythm) |>
  summarise(mae = mean(abs(forecast - rr)), .groups = "drop") |>
  mutate(label = paste0("MAE = ", round(mae, 3), " s"))

ggplot(demo_forecast, aes(beat)) +
  geom_line(aes(y = rr), linewidth = 0.3) +
  geom_point(aes(y = rr), size = 0.5) +
  geom_line(aes(y = forecast), colour = "firebrick", linewidth = 0.6, na.rm = TRUE) +
  geom_text(data = demo_error, aes(x = Inf, y = Inf, label = label),
            hjust = 1.1, vjust = 1.5, size = 3.5, inherit.aes = FALSE) +
  facet_wrap(~rhythm, scales = "free_x") +
  labs(x = "Beat number", y = "R-R interval (s)",
       title = "ARIMA forecast (Red) vs Actual")

We build a per-window key table and a matching list of windows for the train and test sets separately. Each model below scores both sets - we look at the training scores to sanity-check that the score separates the classes, and keep the test scores untouched/unviewed for the final evaluation at the end. Now score every window. This fits one ARIMA model per window; it is cached so it only re-runs if the data changes.

train_grp  <- rr_train |> group_by(record, window, rhythm)
train_keys <- group_keys(train_grp)
train_list <- group_split(train_grp)

test_grp  <- rr_test |> group_by(record, window, rhythm)
test_keys <- group_keys(test_grp)
test_list <- group_split(test_grp)

train_keys$arima_err <- map_dbl(train_list, arima_error)
test_keys$arima_err  <- map_dbl(test_list, arima_error)

Now we can compare the forecast error in Normal vs AF across our training set:

train_keys |>
  group_by(rhythm) |>
  summarise(mean_forecast_error = mean(arima_err), .groups = "drop")

The mean error is higher for AF, but how well does the forecast error alone separate AF from Normal? We are treating arima_err directly as an AF score (larger = more AF-like). To make a prediction of AF vs Normal we to pick a threshold score i.e., if arima_err is greater than threshold we predict AF and below Normal.

To do this we look at the ROC curve, which traces sensitivity against 1 - specificity as we sweep the detection threshold across every value of arima_err. We then pick the threshold value that maximises the Area Under the Curve (or AUC).

train_keys |>
  roc_curve(truth = rhythm, arima_err, event_level = "first") |>
  autoplot() +
  labs(title = "ROC: ARIMA forecast error as an AF detector (train)")

train_keys |>
  roc_auc(truth = rhythm, arima_err, event_level = "first")

Q5 (1 marks): What other type of curve could we use to select the threshold? (Hint: think back to earlier labs)

7 Frequency Domain

The frequency-domain formulation we are going to use is does AF and Normal signals have distinguishable amounts of different frequencies in their signals?

We can infer the frequency distribution of a signal using the Fourier transform (base R’s fft()). A regular rhythm varies slowly (its R-R series wanders gently with breathing), so we expect its energy sits at low frequencies in a single concentrated peak. AF jumps from beat to beat, so its energy spreads out and up into higher frequencies. Rather than collapse this into one hand-picked score, we summarise each window’s spectrum with a handful of frequency-space features and then let a logistic regression learn how to weight them.

For each window we compute four interpretable frequency statistics:

  • hf_frac - fraction of power in the upper half of the band (AF pushes power up here),
  • spectral_centroid - the power-weighted mean frequency bin (the spectrum’s “centre of mass”; higher for AF),
  • spectral_entropy - how spread out / noisy the spectrum is, normalised to [0, 1] (a flat, broadband AF spectrum scores near 1; a single sharp peak scores near 0),
  • peak_frac - the share of power in the single biggest bin (a regular rhythm concentrates power in one peak, so this is high for Normal).
freq_features <- function(rr) {
  x     <- rr - mean(rr)                       # remove the DC offset
  n     <- length(x)
  power <- (Mod(fft(x))^2)[2:floor(n / 2)]     # one-sided power, dropping the DC term
  if (sum(power) == 0) {
    return(tibble(hf_frac = 0, spectral_centroid = 0,
                  spectral_entropy = 0, peak_frac = 0))
  }
  k    <- length(power)
  bins <- seq_len(k)
  p    <- power / sum(power)                    # spectrum as a probability distribution

  tibble(
    hf_frac           = sum(power[round(k / 2):k]) / sum(power),
    spectral_centroid = sum(bins * p),
    spectral_entropy  = -sum(p * log(p + 1e-12)) / log(k),
    peak_frac         = max(p)
  )
}

# attach the features to each window's key table. Assigning by column name
# overwrites any existing copies, so re-running this block never duplicates or
# mangles column names (unlike bind_cols).
train_feats <- map_dfr(train_list, ~ freq_features(.x$rr))
test_feats  <- map_dfr(test_list,  ~ freq_features(.x$rr))

train_keys[names(train_feats)] <- train_feats
test_keys[names(test_feats)]   <- test_feats

train_keys |>
  group_by(rhythm) |>
  summarise(across(c(hf_frac, spectral_centroid, spectral_entropy, peak_frac), mean),
            .groups = "drop")

Now we fit a simple logistic regression that predicts the probability a window is AF from those four features, using the training windows only. The fitted probability is the frequency model’s AF score (freq_prob).

freq_model <- glm(
  (rhythm == "AF") ~ hf_frac + spectral_centroid + spectral_entropy + peak_frac,
  family = binomial,
  data   = train_keys
)

summary(freq_model)
## 
## Call:
## glm(formula = (rhythm == "AF") ~ hf_frac + spectral_centroid + 
##     spectral_entropy + peak_frac, family = binomial, data = train_keys)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -9.9435     1.8149  -5.479 4.28e-08 ***
## hf_frac            -5.3487     0.9058  -5.905 3.53e-09 ***
## spectral_centroid   0.6396     0.1053   6.073 1.26e-09 ***
## spectral_entropy   11.1099     1.7850   6.224 4.84e-10 ***
## peak_frac           0.6069     1.4757   0.411    0.681    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1934.5  on 1562  degrees of freedom
## Residual deviance: 1473.2  on 1558  degrees of freedom
## AIC: 1483.2
## 
## Number of Fisher Scoring iterations: 5
# predicted probability of AF for every window
train_keys$freq_prob <- predict(freq_model, type = "response")
test_keys$freq_prob  <- predict(freq_model, newdata = test_keys, type = "response")

Q6 (2 marks) Based on the output from summary above which frequency-based features are the best predictors of AF in this dataset? Does this match our hypotheses (above) for what each of these features will capture?

Similarly to the time-based model above, we have to pick a threshold on the logistic regression output to decide when to make a prediction of AF. We can sweep possible thresholds from 0 to 1 to again trace out the ROC curve. Each point is one cut-off for calling a window AF; the dashed diagonal is chance, and the more the curve bows towards the top-left the better the model separates AF from Normal:

train_keys |>
  roc_curve(truth = rhythm, freq_prob, event_level = "first") |>
  autoplot() +
  labs(title = "ROC: logistic-regression frequency detector (train)")

train_keys |>
  roc_auc(truth = rhythm, freq_prob, event_level = "first")

Q7 (1 mark): What does it mean for the ROC curve to be below the dotted line?

8 State Space (Hidden Markov Model)

The state-space model formulation is: what sequence of hidden states best explains the observed intervals? where the rhythm moves between unobserved states that each emit specific intervals with different probabilities.

Rather than discretising each interval into a symbol, we model the interval values directly as continuous series of emitted values. We fit a two-state HMM on all the training windows using depmixS4. With two hidden states the model tends to discover one “regular” state (intervals tightly clustered around a mean) and one “irregular” state (a broad, high-variance distribution). Crucially we never show it the AF/Normal labels - the states in this case are learned purely from the distribution of R-R intervals.

# concatenate the training windows into one series, recording each window's length
# so depmixS4 treats the windows as separate sequences (the ntimes argument)
train_seq    <- bind_rows(train_list)
ntimes_train <- map_int(train_list, nrow)

hmm_spec <- depmix(rr ~ 1, data = train_seq, nstates = 2,
                   family = gaussian(), ntimes = ntimes_train)
hmm_fit  <- fit(hmm_spec, verbose = FALSE)
## converged at iteration 52 with logLik: 23365.62
summary(hmm_fit)
## Initial state probabilities model 
##   pr1   pr2 
## 0.167 0.833 
## 
## Transition matrix 
##         toS1  toS2
## fromS1 0.990 0.010
## fromS2 0.002 0.998
## 
## Response parameters 
## Resp 1 : gaussian 
##     Re1.(Intercept) Re1.sd
## St1           0.994  0.085
## St2           0.689  0.158

The two states differ in their emission distributions - one has a much larger standard deviation (the AF-like, irregular state). We use the training labels only to decide which hidden state that is: the state that accounts for most of the beats occurring in AF windows.

train_states <- posterior(hmm_fit, type = "viterbi")$state   # decoded state per beat

af_state <- tibble(rhythm = train_seq$rhythm, state = train_states) |>
  count(state, rhythm) |>
  group_by(state) |>
  summarise(af_frac = sum(n[rhythm == "AF"]) / sum(n), .groups = "drop") |>
  slice_max(af_frac, n = 1, with_ties = FALSE) |>
  pull(state)

Now to label a whole window we infer the underlying hidden states across the window and look at how much of it sits in the AF state. The fraction of a window’s beats in the AF state is a natural continuous AF score; a simple majority vote would threshold it at 0.5, but below we instead fit the best threshold from an ROC curve (below).

For new (test) data, depmixS4 predicts by building a model with the same structure, copying the fitted parameters into it, and decoding - it does not re-estimate anything:

# Viterbi-decode a list of windows under an already-fitted model
decode_states <- function(window_list, fitted) {
  seq_df <- bind_rows(window_list)
  spec   <- depmix(rr ~ 1, data = seq_df, nstates = 2,
                   family = gaussian(), ntimes = map_int(window_list, nrow))
  spec   <- setpars(spec, getpars(fitted))   # reuse the trained parameters
  posterior(spec, type = "viterbi")$state
}

# fraction of each window's beats that fall in the AF state
af_fraction <- function(states, window_list) {
  idx <- rep(seq_along(window_list), map_int(window_list, nrow))
  as.numeric(tapply(states == af_state, idx, mean))
}

train_keys$hmm_score <- af_fraction(train_states, train_list)
test_keys$hmm_score  <- af_fraction(decode_states(test_list, hmm_fit), test_list)

Now let’s fit the threshold we are going to use to determine how much of a window should be in the AF state to classify the window as AF.

roc_hmm <- train_keys |>
  roc_curve(truth = rhythm, hmm_score, event_level = "first")

hmm_threshold <- roc_hmm |>
  mutate(youden = sensitivity + specificity - 1) |>
  slice_max(youden, n = 1, with_ties = FALSE)

hmm_threshold
autoplot(roc_hmm) +
  geom_point(data = hmm_threshold,
             aes(x = 1 - specificity, y = sensitivity),
             colour = "firebrick", size = 3) +
  labs(title = "ROC: HMM AF-state fraction as a detector (train)",
       subtitle = paste0("Chosen threshold (max Youden's J) = ",
                         round(hmm_threshold$.threshold, 3)))

We then label every window by comparing its AF-state fraction to that fitted threshold:

label_hmm <- function(df) {
  df |> mutate(hmm_pred = factor(if_else(hmm_score >= hmm_threshold$.threshold,
                                         "AF", "Normal"),
                                 levels = c("AF", "Normal")))
}
train_keys <- label_hmm(train_keys)
test_keys  <- label_hmm(test_keys)

How well does that rule separate the rhythms on the training windows?

train_keys |>
  conf_mat(truth = rhythm, estimate = hmm_pred)
##           Truth
## Prediction   AF Normal
##     AF     1078    215
##     Normal    1    269

Q8 (2 marks): In this HMM, what is observed and what is hidden?

9 Comparing the Three Detectors

Each model gives every window a score where larger means more AF-like. To compare them fairly we rescale each score to [0, 1] (so they share an axis) and evaluate on the held-out test records. The ROC curve and its AUC summarise how well each score separates AF from normal across all possible cut-offs - the same yardstick both papers report.

test_scored <- test_keys |>
  mutate(across(c(arima_err, freq_prob, hmm_score), rescale)) |>   # higher = more AF
  pivot_longer(c(arima_err, freq_prob, hmm_score),
               names_to = "model", values_to = "score") |>
  mutate(
    model = recode(model,
                   arima_err = "Time (ARIMA)",
                   freq_prob = "Frequency (logistic regression)",
                   hmm_score = "State-space (HMM)"),
    truth = factor(rhythm, levels = c("AF", "Normal"))
  )

# ROC curves, all three models on one set of axes
test_scored |>
  group_by(model) |>
  roc_curve(truth, score, event_level = "first") |>
  autoplot() +
  labs(title = "ROC: three AF detectors on the test set")

test_scored |>
  group_by(model) |>
  roc_auc(truth, score, event_level = "first") |>
  dplyr::select(model, .estimate) |>   # dplyr::select - MASS (via depmixS4) also exports select()
  arrange(desc(.estimate))

Q9 (2 marks): How do these models compare to one another?

Q10 (3 marks): Suggest a way we could improve each of these models performance and suggest at least 3 additional analyses that it would be a good idea to perform before trying to use these in the real-world.

Reproducibility

sessionInfo()
## R version 4.6.0 (2026-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Arch Linux
## 
## Matrix products: default
## BLAS:   /usr/lib/libblas.so.3.12.0 
## LAPACK: /usr/lib/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Halifax
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] depmixS4_1.5-3  nlme_3.1-169    Rsolnp_2.0.1    MASS_7.3-65    
##  [5] nnet_7.3-20     scales_1.4.0    yardstick_1.4.0 HMM_1.0.2      
##  [9] fs_2.1.0        here_1.0.2      lubridate_1.9.5 forcats_1.0.1  
## [13] stringr_1.6.0   dplyr_1.2.1     purrr_1.2.2     readr_2.2.0    
## [17] tidyr_1.3.2     tibble_3.3.1    ggplot2_4.0.3   tidyverse_2.0.0
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1    timeDate_4052.112   farver_2.1.2       
##  [4] S7_0.2.2            fastmap_1.2.0       digest_0.6.39      
##  [7] rpart_4.1.27        timechange_0.4.0    lifecycle_1.0.5    
## [10] survival_3.8-6      magrittr_2.0.5      compiler_4.6.0     
## [13] rlang_1.2.0         sass_0.4.10         tools_4.6.0        
## [16] yaml_2.3.12         data.table_1.18.4   knitr_1.51         
## [19] labeling_0.4.3      DiceDesign_1.10     RColorBrewer_1.1-3 
## [22] parsnip_1.6.0       withr_3.0.2         numDeriv_2016.8-1.1
## [25] workflows_1.3.0     grid_4.6.0          stats4_4.6.0       
## [28] tune_2.1.0          colorspace_2.1-2    future_1.70.0      
## [31] globals_0.19.1      signal_1.8-1        cli_3.6.6          
## [34] rmarkdown_2.31      generics_0.1.4      rstudioapi_0.18.0  
## [37] future.apply_1.20.2 tzdb_0.5.0          cachem_1.1.0       
## [40] splines_4.6.0       dials_1.4.3         forecast_9.0.2     
## [43] parallel_4.6.0      urca_1.3-4          vctrs_0.7.3        
## [46] hardhat_1.4.3       Matrix_1.7-5        jsonlite_2.0.0     
## [49] hms_1.1.4           listenv_0.10.1      gower_1.0.2        
## [52] jquerylib_0.1.4     recipes_1.3.2       glue_1.8.1         
## [55] parallelly_1.47.0   codetools_0.2-20    rsample_1.3.2      
## [58] stringi_1.8.7       gtable_0.3.6        furrr_0.4.0        
## [61] pillar_1.11.1       htmltools_0.5.9     ipred_0.9-15       
## [64] lava_1.9.1          truncnorm_1.0-9     R6_2.6.1           
## [67] rprojroot_2.1.1     evaluate_1.0.5      lattice_0.22-9     
## [70] fracdiff_1.5-4      bslib_0.10.0        class_7.3-23       
## [73] Rcpp_1.1.1-1.1      prodlim_2026.03.11  xfun_0.57          
## [76] zoo_1.8-15          pkgconfig_2.0.3