HeartRateLab.jl

interpolate(n::Array{Float64,1}; method::Symbol=:linear, fs::Int=10)

Interpolates the RR-intervals to fit a given sampling rate.

Arguments: n: the array that contains the RR-intervals method: the interpolation method, default=:linear (options: :constant, :linear, :quadratic, :cubic) fs: the sampling rate, default=10 Hz

Returns: the array with the interpolated values

interpolate_nans!(n::Array{Float64,1}; method::Symbol=:linear)

Interpolates nan values in a time series.

Arguments: n: the array that contains the RR-intervals method: the interpolation method, default=:linear (options: :constant, :linear, :quadratic, :cubic)

Returns: the array with the interpolated values

ms2bpm(n::Array{Float64,1})

Converts milliseconds to beats per minute (BPM).

Arguments: n: the array of RR-intervals in milliseconds

Returns: the array of RR-intervals in BPM

read_txt(infile::String)

Reads a text file with RR-intervals. The file should contain one RR-interval in milliseconds per line (the same way Kubios exports them).

Arguments: infile: the input file name

Returns: an array with the read data

read_wfdb(record::String,annotator::String)

Uses external ann2rr function from the WFDB Software Package to read the RR-intervals from a WFDB record. A possible alternative is to use the wfdb-python package.

Arguments: record: the record name annotator: the annotator of the record

Returns: an array with the read data

read_xdf(infile::String)

Reads an XDF file with RR-intervals. The file should contain one RR-interval stream, marked with "RR" in the name.

Arguments: infile: the input file name

Returns: an array with the read data

replace_bio_outliers(n::Array{T,1}) where T<:Real

Turns any rr-interval outside biologically plausible limits, i.e., less than 300 ms or greater than 2000 ms, into a NaN value.

Arguments: n: the array that contains the RR-intervals

Returns: the array without the intervals that are less than 300 ms or greater than 2000 ms

replace_ectopic_beats!(n::Array{T,1}; method::Symbol=:malik, threshold::Float64=0.2) where T<:Real

Replaces ectopic beats with NaN values.

Arguments: n: the array that contains the RR-intervals method: the method to replace the ectopic beats, default=:malik (options: :malik, :kamath, :acar, :karlsson, :custom) threshold: the threshold to replace the ectopic beats, default=0.2

Returns: the array without the ectopic beats

replace_statistical_outliers(n::Array{T,1}) where T<:Real

Turns any rr-interval outside the 2.5th and 97.5th percentiles into a NaN value.

Arguments: n: the array that contains the RR-intervals

Returns: the array without the intervals outside the 2.5th and 97.5th percentiles

replace_zeros(n::Array{T,1}) where T<:Real

Turns any zero value into a NaN value.

Arguments: n: the array that contains the RR-intervals

Returns: the array without the zero values

strip_extremes(n::Array{Float64,1})

Strips any NaN values from the extremes of the array of RR-intervals.

Arguments: n: the array that contains the RR-intervals

Returns: the array without the NaN values at the extremes

windowed(n::Array{Float64,1}; window_size::Int=60, stride::Int=1, time::Symbol=:beats, f::Function=identity)

This function returns views of the array of RR-intervals in a sliding window. The sliding window can be defined in beats or milliseconds.

Arguments: n: the array that contains the RR-intervals window_size: the size of the window, default=60 (beats) stride: the step size of the sliding window, default=1 (beat) time: the time unit of the window, default=:beats (options: :beats, :ms) f: reducing function to apply to the window, default=identity

Returns: the views of the array in a sliding window, with the reducing function applied to each window

Evaluation

Evaluation pipeline for comparing HRV data and model outputs.

This module provides functions for sliding-window analysis and feature comparison across different data modes (continuous, windowed, ensemble).

aic(loglik, k) -> Float64

Akaike Information Criterion: 2k - 2·loglik. Lower is better.

bic(loglik, k, n) -> Float64

Bayesian (Schwarz) Information Criterion: k·ln(n) - 2·loglik, where n is the number of observations. Lower is better; penalises parameters more than AIC.

eval_distance(real::DataFrame, ensemble::DataFrame; metric=:mahalanobis, features=nothing) -> NamedTuple

Compute feature-space distance between real and ensemble data.

Measures how far apart real observations and ensemble samples are in multi-dimensional feature space using different distance metrics.

Arguments

• real::DataFrame: Real feature observations (one or more rows)
• ensemble::DataFrame: Ensemble feature samples (typically many rows)
• metric=:mahalanobis: Distance metric: :mahalanobis, :euclidean
• features=nothing: Specific features to use. If nothing, uses common columns

Returns

• NamedTuple: (distance=Float64, metric=Symbol, feature_contributions=Dict)

Details

Distance Metrics:

• :euclidean — Standard L2 distance: √(Σ(x-y)²)

  • Simple, interpretable
  • Sensitive to scale differences

• :mahalanobis — Accounts for covariance structure: √((x-y)ᵀΣ⁻¹(x-y))

  • Accounts for feature correlations
  • Scale-invariant
  • More robust

Feature Contributions:

• Each feature's contribution to total distance
• Identifies which features drive the difference
• Useful for debugging model misspecification

Interpretation:

• Small distance (0-1) = good model fit
• Medium distance (1-5) = acceptable fit
• Large distance (>5) = poor model fit
• Depends on feature scaling and number of features

Examples

result = eval_distance(real_features, ensemble_features; metric=:mahalanobis)

println("Distance: $(result.distance)")
println("Metric: $(result.metric)")
sort!(result.feature_contributions, rev=true) |> display  # Top contributors

eval_distributional(real::DataFrame, ensemble::DataFrame; test=:ks, features=nothing) -> DataFrame

Compare feature distributions between real and model data using statistical tests.

For each feature, tests whether the ensemble distribution differs from real observations using Kolmogorov-Smirnov, Mann-Whitney U, or Anderson-Darling tests.

Arguments

• real::DataFrame: Feature observations (could be 1 row or many windows/subjects)
• ensemble::DataFrame: Ensemble feature samples (typically many rows from synthetic data)
• test=:ks: Statistical test to use: :ks (Kolmogorov-Smirnov), :mw (Mann-Whitney U), :ad (Anderson-Darling)
• features=nothing: Specific features to test. If nothing, uses all columns in both DataFrames

Returns

• DataFrame: One row per feature with columns: feature, statistic, pvalue, testname

Details

Statistical Tests:

• :ks — Kolmogorov-Smirnov test: compares empirical CDFs
• :mw — Mann-Whitney U test: non-parametric rank test
• :ad — Anderson-Darling test: goodness-of-fit test

Interpretation:

• Small p-value (< 0.05): Real and ensemble distributions are significantly different
• Large p-value (> 0.05): Cannot reject that distributions are equal
• statistic: Test-specific value (larger typically = more different)

Data Handling:

• NaN values are filtered out per feature
• Features with all NaN are skipped
• Single real observation compares point to ensemble distribution

Examples

# Compare real windowed data to synthetic ensemble
real_windows = windowed_feature_set(data; window_size=300, overlap=150)
ensemble = simulate_ensemble(model, params, 300; n_sim=100)
ensemble_features = extract_ensemble_features(ensemble)

# Test each feature
results = eval_distributional(real_windows, ensemble_features; test=:ks)

# View results
sort!(results, :p_value)  # Features most different from ensemble

eval_scalar(real::DataFrame, ensemble::DataFrame; features=nothing) -> DataFrame

Compare scalar statistics (means and errors) between real and ensemble data.

For each feature, computes mean of real observations and ensemble samples, then calculates absolute and relative errors.

Arguments

• real::DataFrame: Real feature observations (one or more rows)
• ensemble::DataFrame: Ensemble feature samples (typically many rows)
• features=nothing: Specific features to compare. If nothing, uses common columns

Returns

• DataFrame: One row per feature with columns: feature, realmean, simmean, abserror, relerror

Details

Metrics:

• real_mean: Mean of real feature observations
• sim_mean: Mean of ensemble feature samples
• abs_error: |realmean - simmean| (absolute difference)
• rel_error: abserror / (abs(realmean) + 1e-8) (relative error, epsilon-normalized)

Interpretation:

• Small errors (< 0.05) indicate excellent model fit
• Medium errors (0.05-0.15) indicate good model fit
• Large errors (> 0.15) indicate model needs improvement

Use Cases:

• Quick overall quality metric (complement to p-values)
• Identifying which features are biased
• Benchmarking against baseline models

Examples

real_windows = windowed_feature_set(data; window_size=300, overlap=150)
ensemble_features = extract_ensemble_features(ensemble)

errors = eval_scalar(real_windows, ensemble_features)

# Find most biased features
sort!(errors, :rel_error)
println(errors[1:5, :])  # Top 5 most biased

extract_ensemble_features(ensemble::Vector{Vector{Float64}}; features=nothing) -> DataFrame

Extract HRV features from all series in a synthetic ensemble.

This applies feature extraction to each independent series in an ensemble, producing a DataFrame with one row per series. This enables empirical feature distributions for Mode 3 (Sampled/Ensemble) analysis.

Arguments

• ensemble::Vector{Vector{Float64}}: Collection of IBI series (from simulate_ensemble())
• features=nothing: Optional list of specific features to extract. If nothing, uses valid_features() for signal length

Returns

• DataFrame: One row per series in ensemble, columns are feature names

Details

Signal Length:

• If all series have the same length, uses that for valid_features() filtering
• If series have different lengths, extracts all computable features (less strict)

Efficiency:

• Current: sequential (can be parallelized with Distributed.pmap())
• NaN handling: Features that fail to compute are marked as NaN

Integration with Evaluation:

# Typical workflow
lif = LIF()
result = fit(lif, data; method=:bayesian)
ensemble = simulate_ensemble(lif, result.params, length(data); n_sim=100)
ensemble_features = extract_ensemble_features(ensemble)

# Now use with comparison functions
eval_distributional(real_windows, ensemble_features; test=:ks)

Examples

# Generate ensemble and extract features
ensemble = simulate_ensemble(model, params, 300; n_sim=50)
features_df = extract_ensemble_features(ensemble)

# features_df has 50 rows (one per simulation) and ~40 columns (features)
# Can now compare against real data using eval_* functions

information_criteria(model, params, data; kwargs...) -> NamedTuple
information_criteria(result::ModelFitResult; kwargs...) -> NamedTuple

Compute the information-criterion summary for a single fitted model on data. Returns (model, n_params, n_obs, loglik, aic, bic). kwargs are forwarded to model_loglikelihood (e.g. n_sim_beats).

model_loglikelihood(model, params, data; n_sim_beats) -> Float64

Gaussian log-likelihood of the observed IBI series data under a model. The model contributes the mean μ_sim and variance σ²_sim of a simulated series; the predictive variance is σ² = max((1/n)·Σ(dataᵢ − μ_sim)², σ²_sim) and

logL = −½·(n·log 2π + n·log σ² + Σ(dataᵢ − μ_sim)² / σ²).

The max(·) combines two failure modes into one robust criterion:

• If the model under-predicts spread (σ²_sim too small — the near-deterministic mechanistic models), σ² falls back to the MLE residual variance Var(data) + (mean(data) − μ_sim)², which rewards matching the data's mean.
• If the model over-predicts spread (σ²_sim larger than the residual variance), σ² uses σ²_sim and the likelihood is penalised for being too wide.

Why this rather than a beat-indexed-residual or a raw marginal Normal(μ, σ) likelihood:

• Generative HRV models produce trajectories with their own phase, so dataᵢ − simᵢ residuals are meaningless even for a well-fitted model.
• A raw Normal(μ_sim, σ_sim) density is numerically ill-conditioned: the mechanistic models are near-deterministic, so σ_sim → 0 (dominated by solver noise), which makes logpdf diverge for reasons unrelated to fit quality.

This form is robust (σ² ≥ Var(data) > 0 for real data) and surfaces DMD's documented weakness — its reconstruction forces the mean toward a fixed ~800 ms, away from a typical recording's mean.

n_sim_beats controls the simulation length used to estimate μ_sim/σ²_sim; the default (500) gives stable moment estimates and is independent of the data length (some models — e.g. LIF's ODE solver — cannot cheaply produce tens of thousands of beats). If the model cannot simulate the given parameters at all, the likelihood is -Inf so that model ranks last rather than crashing a comparison.

model_n_params(model::AbstractHRVModel, params::NamedTuple) -> Int

Effective number of free parameters of a fitted model, the k in AIC/BIC.

For the mechanistic models the continuous fitted parameters live in params (LIF I ⇒ 1; VanDerPol μ,heart_rate ⇒ 2; Lorenz σ,ρ,β,threshold ⇒ 4), so the default counts the NamedTuple fields. DMD is special-cased: it stores no continuous parameters in params (it is data-driven) — its complexity is the number of retained dynamic modes (the truncation rank r), which is what must enter the criterion for an honest comparison.

rank_models(results::AbstractVector{ModelFitResult}; criterion=:bic, kwargs...) -> DataFrame

Rank a set of fitted models on their respective data by an information criterion.

Each ModelFitResult carries the fitted model, its parameters, and the data it was fitted to. The returned DataFrame is sorted best-first (lowest criterion) with columns:

• model — model type name
• n_params — effective free-parameter count k (model_n_params)
• n_obs — number of observed beats
• loglik — Gaussian-residual log-likelihood (model_loglikelihood)
• aic, bic— the two criteria
• delta — criterion - min(criterion) (0 for the best model)
• weight — Akaike/Schwarz weight exp(-Δ/2) / Σ exp(-Δ/2) (relative support)
• rank — 1 = best

criterion must be :aic or :bic (default :bic). kwargs are forwarded to model_loglikelihood.

simulate_ensemble(model::AbstractHRVModel, params::NamedTuple, n_beats::Int; n_sim::Int=100, rng=Random.default_rng()) -> Vector{Vector{Float64}}

Generate an ensemble of N independent synthetic IBI series from a model.

This enables Mode 3 (Sampled Windows) analysis: create a distribution of synthetic timeseries, then extract features from each to get an empirical feature distribution.

Arguments

• model::AbstractHRVModel: The HRV model to simulate from (LIF, VanDerPol, Lorenz, DMD, etc.)
• params::NamedTuple: Model parameters (typically from fit() result)
• n_beats::Int: Number of beats per simulated series
• n_sim::Int=100: Number of independent simulations to generate
• rng=Random.default_rng(): Random number generator for reproducibility

Returns

• Vector{Vector{Float64}}: Ensemble of nsim series, each is a Vector{Float64} of length ≈ nbeats

Details

Independence:

• Each series is generated independently from the same model and parameters
• Series will differ due to stochastic components in the model (noise, MCMC, etc.)
• With fixed RNG seed, results are reproducible

Stochastic vs Deterministic Models:

• Stochastic models (LIF, VDP with noise) will produce different series each run
• Deterministic models (Lorenz, DMD) may produce identical series without perturbation

Examples

# Fit model to data
lif = LIF()
result = fit(lif, data; method=:bayesian)

# Generate synthetic ensemble
ensemble = simulate_ensemble(lif, result.params, length(data); n_sim=100)
# Returns: Vector{Vector{Float64}} with 100 synthetic IBI series

# Extract features and compare
ensemble_features = extract_ensemble_features(ensemble)
eval_distributional(real_windows, ensemble_features; test=:ks)

windowed_feature_set(data::Vector{Float64}; window_size::Int=300, overlap::Int=150) -> DataFrame

Extract HRV features from overlapping time windows of IBI data.

This enables Mode 2 (Time Windows) analysis: convert one continuous timeseries into a distribution of feature vectors, one per window.

Arguments

• data::Vector{Float64}: IBI timeseries (milliseconds)
• window_size::Int=300: Number of beats per window (default 300 beats ≈ 4-5 minutes)
• overlap::Int=150: Number of overlapping beats between windows (default 150 = 50% overlap)

Returns

• DataFrame: One row per window, columns are valid features for window_size

Details

Window Computation:

• First window: data[1:window_size]
• Second window: data[(1 + window_size - overlap):(1 + window_size - overlap + window_size - 1)]
• Continue until insufficient beats remain for a complete window

Feature Selection:

• Only features valid for window_size beats are included (via valid_features())
• If a window has features that fail to compute, they are marked as NaN

Edge Cases:

• If data is shorter than window_size, returns empty DataFrame
• If overlap >= window_size, returns single window (or empty if too short)

Examples

# Load HRV data
data = read_txt("subject.txt")  # 1000 beats

# Extract features from 5-minute windows with 50% overlap
windows = windowed_feature_set(data; window_size=300, overlap=150)
# Returns DataFrame with ~7 rows (windows) × ~40 columns (features)

# Use in evaluation
model_features = extract_ensemble_features(ensemble)  # 100 rows
eval_distributional(windows, model_features; test=:ks)  # Compare

Performance

• Parallelizable: each window is independent
• Current implementation is sequential
• For large datasets (>10k beats), consider Distributed.pmap() wrapper