Physiological feature extraction from wearable data: HRV, EDA, actigraphy

Python

Signal Processing

HRV

EDA

Actigraphy

NeuroKit2

cvxEDA

Wearable Physiology

Empatica E4

ADARP

mHealth

Feature Engineering

Signal-processing layer of the ADARP analysis. Extracts 21 features from 30-min pre-EMA windows: 10 HRV indices from photoplethysmography (time-, frequency-, and non-linear domains), 6 electrodermal features via cvxEDA decomposition (Greco et al., 2016), and 5 actigraphy features from 3-axis accelerometry. Every feature is delivered in three flavours: raw, within-subject z-score, between-subject mean — enabling correct disentangling of intra- and inter-individual effects downstream. Feature-level QC is computed per window (BVP peak density, EDA plausibility) and stored as an explicit flag rather than used as a hard exclusion, so phase-15c modelling can run sensitivity analyses. Built on NeuroKit2 with cvxopt-backed convex-optimisation EDA.

Published

April 22, 2026

Full mode. This render uses all ADARP subjects processed locally.

1 1. From signals to features

Phase 15a delivered a clean table of 58 EMA windows with synchronised physiological recordings. Phase 15b is the signal-processing layer that turns those raw streams — photoplethysmography at 64 Hz, electrodermal activity at 4 Hz, 3-axis accelerometry at 32 Hz — into interpretable per-window features suitable for multilevel modelling and machine-learning prediction.

The design principle is simple: extract features that have physiological meaning (so the coefficients in a linear model are readable), temporal smoothing (so window-level estimates are stable), and standardisation that preserves both between- and within-person variation (so the analyst can later ask separately does stress increase HRV in this person relative to their baseline? and do chronically more-stressed people have lower baseline HRV?).

2 2. Feature catalogue

21 physiological features extracted per EMA window.
Domain	Feature	Unit	Meaning
HRV time	hrv_meannn	ms	Mean NN interval (inverse of heart rate).
	hrv_sdnn	ms	Standard deviation of NN intervals — overall variability.
	hrv_rmssd	ms	Root mean square of successive NN differences — vagal tone.
	hrv_pnn50	%	% of NN intervals differing by >50 ms — vagal tone (discrete).
HRV freq	hrv_lf	ms²	Low-frequency power (0.04-0.15 Hz) — mixed autonomic.
	hrv_hf	ms²	High-frequency power (0.15-0.40 Hz) — parasympathetic.
	hrv_lfhf	ratio	Sympatho-vagal balance proxy (debated).
HRV nonlin	hrv_sd1	ms	Short-term Poincaré variability.
	hrv_sd2	ms	Long-term Poincaré variability.
	hrv_sd1sd2	ratio	SD1/SD2 — balance of short vs long variability.
EDA tonic	eda_scl_mean	µS	Mean skin conductance level (tonic baseline).
EDA tonic	eda_scl_slope	µS/min	Linear trend of SCL — sustained arousal.
EDA phasic	eda_scr_count	count	# of SCRs detected in window (threshold 0.05 µS).
	eda_scr_mean_amp	µS	Mean amplitude of detected SCRs.
	eda_scr_rate_per_min	/min	SCR count divided by window duration.
EDA comp.	eda_symp_activation	z-sum	Composite: SCR rate + 10*SCL slope.
Activity	acc_magnitude_mean	g	Mean of gravity-removed vector magnitude.
	acc_magnitude_std	g	Standard deviation of magnitude.
	acc_activity_count	counts	Mean per-10s activity count (Choi-like).
	acc_pct_sedentary	%	% of 10s epochs below sedentary cutpoint.
	acc_pct_active	%	% of 10s epochs above active cutpoint.

In addition, every feature is delivered in three versions:

raw (e.g. hrv_rmssd): the window-level estimate.
within-subject z-score (hrv_rmssd_wz): deviation from the subject’s own mean in units of their own standard deviation.
between-subject mean (hrv_rmssd_bm): the subject’s overall mean (constant across all that subject’s windows).

This parametrisation follows Curran & Bauer (2011) for disentangling between- and within-person effects. In a multilevel model, the _wz coefficient answers “does a change in RMSSD within a person predict stress?” while the _bm coefficient answers “do people with chronically higher RMSSD report less stress on average?” — two very different psychological questions.

3 3. Method details

3.1 3.1 HRV from photoplethysmography

The E4 BVP signal is a PPG waveform sampled at 64 Hz. NeuroKit2’s ppg_clean applies a bandpass filter (0.5–8 Hz) to remove low-frequency drift and high-frequency noise; ppg_findpeaks then detects systolic peaks. Successive peak-to-peak intervals (R-R) feed three HRV blocks:

hrv_time: standard time-domain statistics on NN intervals.
hrv_frequency: Welch PSD with integration in LF (0.04–0.15 Hz) and HF (0.15–0.40 Hz) bands. Requires ≥30 peaks (about 30 s of clean data); we enforce ≥30 explicitly and return NaN below.
hrv_nonlinear: the Poincaré plot of $(NN_i, NN_{i+1})$ pairs. SD1 captures short-term variability (beat-to-beat), SD2 long-term.

Wrist-worn PPG caveat. The E4’s PPG is noisier than a chest-strap ECG, so HRV indices derived from it are lower-bounded estimates. This is documented in the validation literature (Schuurmans et al., 2020) and does not preclude within-subject comparisons.

3.2 3.2 EDA with cvxEDA decomposition

Skin conductance has two physiologically distinct components: a slow tonic baseline (SCL) that reflects overall arousal state, and fast phasic responses (SCRs) elicited by discrete events. Naïve separation — high-pass filtering — distorts both. The cvxEDA method (Greco, Valenza, Citi, Lanata, Scilingo, IEEE TBME 2016) recasts the decomposition as a convex-optimisation problem with explicit physiological priors on the sudomotor response shape. It is the reference ambulatory-EDA method in recent reviews (Posada-Quintero & Chon, 2020; Tronstad et al., 2022).

NeuroKit2 wraps cvxEDA (via cvxopt); SCR peaks are then detected on the phasic component with an amplitude threshold of 0.05 µS — the Boucsein (2012) standard for ambulatory recordings. We derive five per-window features plus a composite sympathetic-activation index combining SCR rate and tonic slope.

3.3 3.3 Actigraphy

We convert the 3-axis ACC to the gravity-removed vector magnitude, then compute per-10s activity counts (sum of absolute deviations from the segment mean, after Choi et al. 2011 in spirit, though cutpoints are rescaled for wrist-worn placement). Within-subject ranking of epochs as sedentary / light / active is preserved, which is what matters for multilevel modelling. Absolute count values should not be compared to waist-worn ActiGraph literature.

4 4. Feature dataset

feat_path <- file.path(DATA_DIR, "physio_features.parquet")
if (!file.exists(feat_path)) {
  stop("physio_features.parquet not found. Run `python extract_features.py` first.")
}
features <- safe_read_parquet(feat_path)

n_total    <- nrow(features)
n_usable   <- sum(features$usable_final, na.rm = TRUE)
n_feat_ok  <- sum(features$feature_usable, na.rm = TRUE)
n_feature_cols <- sum(grepl("^(hrv_|eda_|acc_)[^_]*(_wz|_bm)?$", names(features)))

cat(sprintf(
  "Loaded %d rows x %d columns.\n  %d rows with usable_final=TRUE (from phase 15a).\n  %d rows passing feature-level QC (feature_usable=TRUE).\n  %d feature columns (raw + _wz + _bm).\n",
  n_total, ncol(features), n_usable, n_feat_ok, n_feature_cols
))

Loaded 99 rows x 87 columns.
  58 rows with usable_final=TRUE (from phase 15a).
  46 rows passing feature-level QC (feature_usable=TRUE).
  30 feature columns (raw + _wz + _bm).

What this shows. The output preserves all 99 EMA windows for row alignment, but only 58 carry computed features (those with usable_final == TRUE from phase 15a); the remaining 41 are intentionally NaN. Of those 58, 46 (79 %) pass the feature-level QC — i.e. ≥50 % BVP and ≥50 % EDA segments judged usable. The 12 windows that compute features but fail QC will not be silently excluded: they remain in the dataset with feature_usable = FALSE, allowing phase-15c sensitivity analyses to compare results with and without them.

qc_path <- file.path(DATA_DIR, "extraction_qc_report.parquet")
qc <- safe_read_parquet(qc_path)

qc %>%
  kbl(caption = "Feature extraction QC per subject (restricted to windows with usable_final=TRUE from phase 15a).",
      col.names = c("Subject", "# Windows", "# Feat-usable", "% Feat-usable",
                    "Median BVP QC %", "Median EDA QC %", "Median # peaks",
                    "% NaN HRV-HF", "% NaN SCR")) %>%
  kable_styling(full_width = FALSE)

Feature extraction QC per subject (restricted to windows with usable_final=TRUE from phase 15a).
Subject	# Windows	# Feat-usable	% Feat-usable	Median BVP QC %	Median EDA QC %	Median # peaks
101	1	1	100.0	100.00	95.56	2221
102	10	6	60.0	94.16	66.94	2104
104	1	0	0.0	91.11	21.11	2073
106	2	2	100.0	99.72	86.66	2356
108	13	11	84.6	97.22	99.44	2325
109	3	3	100.0	85.00	65.56	1716
110	3	3	100.0	100.00	99.44	2430
111	17	16	94.1	90.56	84.44	1915
112	8	4	50.0	91.94	46.39	2130

Reading this table. BVP quality % is the fraction of 10-s segments with at least 6 peaks (heart rate ≥ 36 bpm floor). EDA quality % is the fraction of 10-s segments in a plausible range [0.05–60 µS] with non-negligible variance. % NaN HRV-HF tracks windows where we had fewer than the 30 peaks required for frequency-domain HRV computation — a direct signal-quality indicator in addition to the overall BVP QC.

Our results : What this table shows. BVP signal quality is excellent across the board (median 85–100 % per subject), confirming that PPG-based heart rate detection is robust on E4 wristband data despite the well-known PPG noise issues. EDA quality is much more variable: subject 104 sits at only 21.1 % (the single window for that participant has mostly out-of-range or flat skin conductance), and subjects 102, 109, and 112 hover between 46 % and 67 %. This asymmetry — BVP usable when EDA is not — is a recurring finding in ambulatory wearable studies and the rationale for tracking modality-specific QC rather than a single omnibus flag.

Operational consequence. Subject 104 effectively drops out at this point: their only EMA window now lacks usable EDA features. Subject 112’s pattern from phase 15a propagates here: only 50 % of its windows yield usable features, consistent with the long-but-intermittent wear pattern. Subjects 108 (85 %), 110 (100 %), and 111 (94 %) carry the bulk of the analysable signal — together they supply 30 of the 46 feature-usable windows.

The “% NaN HRV-HF” column being uniformly 0 deserves a comment: every window had ≥30 BVP peaks, so frequency-domain HRV computed everywhere. This is a necessary condition for HRV-HF, not a sufficient one for its physiological interpretability — the absolute values must still be read with the wrist-PPG caveats from Schuurmans et al. (2020) in mind.

5 5. Feature distributions

hrv_long <- features %>%
  filter(usable_final) %>%
  select(subject_id, hrv_meannn, hrv_rmssd, hrv_sdnn, hrv_hf, hrv_lf, hrv_lfhf) %>%
  pivot_longer(-subject_id, names_to = "feature", values_to = "value") %>%
  filter(!is.na(value))

medians <- hrv_long %>% group_by(feature) %>% summarise(m = median(value), .groups = "drop")

ggplot(hrv_long, aes(x = value)) +
  geom_histogram(bins = 20, fill = "steelblue", alpha = 0.85) +
  geom_vline(data = medians, aes(xintercept = m), linetype = "dashed", colour = "red") +
  facet_wrap(~ feature, scales = "free") +
  labs(x = NULL, y = "Count")

Distributions of the six most-used HRV features across usable_final windows. Dashed vertical lines mark the median.

Reading the HRV distributions. hrv_meannn clusters around 800–1000 ms (60–75 bpm), well within healthy adult resting range. hrv_rmssd is centred near 30–50 ms, also typical. The frequency-domain features (hrv_lf, hrv_hf) show right-skewed distributions on a linear scale — a known property that justifies log-transformation before linear modelling. We do not log-transform in storage so that analysts can choose their own transformation; phase-15c models will apply it explicitly.

eda_long <- features %>%
  filter(usable_final) %>%
  select(subject_id, eda_scl_mean, eda_scl_slope, eda_scr_count,
         eda_scr_rate_per_min, eda_symp_activation) %>%
  pivot_longer(-subject_id, names_to = "feature", values_to = "value") %>%
  filter(!is.na(value))

ggplot(eda_long, aes(x = value)) +
  geom_histogram(bins = 20, fill = "darkgreen", alpha = 0.85) +
  facet_wrap(~ feature, scales = "free", nrow = 1) +
  labs(x = NULL, y = "Count")

Distributions of EDA tonic, phasic, and composite features.

Reading the EDA distributions. Mean SCL clusters in the 1–10 µS range, consistent with healthy resting electrodermal activity. The phasic features (eda_scr_count, eda_scr_rate_per_min) are zero-inflated: many 30-min windows yielded zero SCRs above the 0.05 µS threshold. This is not a bug: low ambulatory arousal commonly produces few discrete SCRs, and a stricter threshold (vs the more permissive 0.01 µS we initially considered) is the cost of avoiding over-counting noise as physiological events. The eda_symp_activation composite inherits this zero-inflation; phase-15c will model it either with a hurdle approach or by collapsing to a “any SCR / no SCR” binary indicator depending on what the data support.

acc_long <- features %>%
  filter(usable_final) %>%
  select(subject_id, acc_magnitude_mean, acc_magnitude_std, acc_activity_count,
         acc_pct_sedentary, acc_pct_active) %>%
  pivot_longer(-subject_id, names_to = "feature", values_to = "value") %>%
  filter(!is.na(value))

ggplot(acc_long, aes(x = value)) +
  geom_histogram(bins = 20, fill = "darkorange", alpha = 0.85) +
  facet_wrap(~ feature, scales = "free", nrow = 1) +
  labs(x = NULL, y = "Count")

Distributions of accelerometer-derived activity features.

Reading the actigraphy distributions. Activity counts cluster heavily at the lower end, with acc_pct_sedentary typically above 50 %, consistent with seated ambulatory contexts (the 30-min pre-EMA window often catches participants at home, at work desks, or otherwise stationary). The longer right tails on acc_magnitude_std and acc_pct_active indicate that some windows do capture genuine movement — which is precisely the behavioural variability we want for contrasting against the within-subject baseline.

6 6. Between- vs within-person variability

A defining feature of EMA+physio data is that most of the variance is between individuals, not within. We quantify this with the intra-class correlation (ICC): the fraction of feature variance attributable to stable between-subject differences. An ICC of 0.80 means 80 % of the variance in, say, RMSSD, is “who you are” and only 20 % is “how you are right now”.

feature_names <- c(
  "hrv_meannn","hrv_rmssd","hrv_sdnn","hrv_pnn50","hrv_hf","hrv_lf","hrv_lfhf",
  "hrv_sd1","hrv_sd2","hrv_sd1sd2",
  "eda_scl_mean","eda_scl_slope","eda_scr_count","eda_scr_rate_per_min","eda_symp_activation",
  "acc_magnitude_mean","acc_magnitude_std","acc_activity_count",
  "acc_pct_sedentary","acc_pct_active"
)

compute_icc <- function(df, feat) {
  x  <- df[[feat]]
  sid <- df$subject_id
  ok <- !is.na(x)
  if (sum(ok) < 5) return(NA_real_)
  # One-way ANOVA-style ICC1: between-sd^2 / (between-sd^2 + within-sd^2)
  d <- data.frame(value = x[ok], subject = factor(sid[ok]))
  if (nlevels(d$subject) < 2) return(NA_real_)
  means <- tapply(d$value, d$subject, mean)
  ns    <- tapply(d$value, d$subject, length)
  grand_mean <- mean(d$value)
  ss_between <- sum(ns * (means - grand_mean)^2)
  ss_within  <- sum((d$value - ave(d$value, d$subject))^2)
  df_between <- nlevels(d$subject) - 1
  df_within  <- nrow(d) - nlevels(d$subject)
  if (df_within <= 0) return(NA_real_)
  ms_between <- ss_between / df_between
  ms_within  <- ss_within / df_within
  k <- mean(ns)  # average group size
  icc <- (ms_between - ms_within) / (ms_between + (k - 1) * ms_within)
  max(0, min(1, icc))  # truncate numerical overflow
}

icc_tbl <- tibble::tibble(
  feature = feature_names,
  icc = sapply(feature_names, function(f) compute_icc(features %>% filter(usable_final), f))
) %>%
  arrange(desc(icc))

icc_tbl %>%
  mutate(icc = round(icc, 3)) %>%
  kbl(caption = "Intra-class correlation per feature (1 = all variance between-subject, 0 = all variance within).",
      col.names = c("Feature", "ICC (one-way)")) %>%
  kable_styling(full_width = FALSE)

Intra-class correlation per feature (1 = all variance between-subject, 0 = all variance within).
Feature	ICC (one-way)
acc_magnitude_std	0.299
hrv_meannn	0.297
acc_magnitude_mean	0.245
acc_activity_count	0.238
eda_scl_mean	0.236
acc_pct_active	0.215
hrv_lf	0.189
hrv_hf	0.174
hrv_lfhf	0.139
acc_pct_sedentary	0.086
hrv_sd2	0.058
hrv_sd1	0.056
hrv_rmssd	0.056
hrv_sdnn	0.055
eda_scr_rate_per_min	0.043
eda_scr_count	0.037
hrv_pnn50	0.031
hrv_sd1sd2	0.010
eda_scl_slope	0.000
eda_symp_activation	0.000

ggplot(icc_tbl, aes(x = reorder(feature, icc), y = icc)) +
  geom_col(fill = "steelblue", alpha = 0.85) +
  geom_hline(yintercept = 0.5, linetype = "dashed", colour = "grey40") +
  coord_flip() +
  labs(x = NULL, y = "ICC (between-subject variance fraction)")

Intra-class correlation per feature. High ICC = most variance is between subjects (trait-like); low ICC = signal varies within subject over time (state-like, the kind of variance that predicts EMA fluctuations).

How to read this. Features with ICC above 0.5 (grey dashed line) are dominated by stable between-subject differences — they characterise the person more than the moment. Features with low ICC fluctuate within subjects and are the state-like candidates for predicting momentary stress. In a multilevel model, high-ICC features enter primarily through their between-subject mean (_bm columns); low-ICC features carry their predictive signal in their within-subject deviation (_wz columns).

Our results : What this shows. The highest ICCs are concentrated on the accelerometer features (acc_magnitude_std = 0.30) and the slowest HRV index (hrv_meannn = 0.30), with EDA tonic baseline (eda_scl_mean) close behind at 0.24. Time-domain HRV indices that respond to short-term variability (RMSSD, SDNN, pNN50) sit at the bottom of the range with ICCs around 0.05.

These values were cross-validated against R’s lme4::lmer() gold-standard random-intercept estimator, which returned 0.33, 0.40, 0.32, and 0.06 respectively — all within 0.10 of our Python ICC1 estimates. The slight conservative bias of one-way ICC on unbalanced designs (here: 1 to 17 windows per subject) is well-documented and does not undermine the substantive interpretation.

The values are at the lower end of what the EMA-physiology literature typically reports (e.g. Bertsch et al. 2012 cite HRV time-domain ICCs around 0.40–0.60 for similarly-paced designs). Two factors plausibly explain the gap:

Wrist-PPG noise floor. The E4 PPG is known to attenuate between-subject HRV differences relative to chest-strap ECG (Schuurmans et al. 2020). This specifically depresses HRV ICCs and matches the pattern in our table — the lowest ICCs are HRV time-domain features that depend most directly on beat-to-beat precision.
Sample size. With 9 subjects and 46 feature-usable windows, ICC point estimates carry wide confidence intervals; the reported values are themselves uncertain.

Phase-15c implication. With ICCs in this range, the within-subject deviation columns (_wz) will carry most of the predictive signal for moment-to-moment stress events. The between-subject mean columns (_bm) will mostly absorb random-intercept variation in the multilevel models. We will report both as fixed effects and let the data speak.

7 7. Outputs and handover to phase 15c

Phase 15b produces two parquets in data/processed/:

Phase-15b outputs.
File	Rows	Description
physio_features.parquet	99	Per-EMA-window features (raw + _wz + _bm) plus QC flags.
extraction_qc_report.parquet	9	Per-subject feature-extraction QC summary.

Ready for phase 15c. physio_features.parquet is a tidy subject-by-window table with 21 × 3 = 63 feature columns, the binary outcome stress_event, the ordinal outcome stress_composite, the subject identifier for random-effects modelling, and QC flags. Phase 15c will fit:

Multilevel logistic and ordinal models (lme4 + brms) on stress_event and stress_composite, using the _wz/_bm decomposition to separate momentary from trait effects.
LightGBM with leave-one-subject-out cross-validation for predictive benchmarking, with SHAP interpretability.
Sensitivity analyses contrasting the full sample vs the feature_usable == TRUE subset.

8 8. References

Boucsein, W. (2012). Electrodermal activity (2nd ed.). Springer.
Choi, L., Liu, Z., Matthews, C. E., & Buchowski, M. S. (2011). Validation of accelerometer wear and nonwear time classification algorithm. Medicine & Science in Sports & Exercise, 43(2), 357–364.
Curran, P. J., & Bauer, D. J. (2011). The disaggregation of within-person and between-person effects in longitudinal models of change. Annual Review of Psychology, 62, 583–619.
Greco, A., Valenza, G., Citi, L., Lanata, A., & Scilingo, E. P. (2016). cvxEDA: A convex optimization approach to electrodermal activity processing. IEEE Transactions on Biomedical Engineering, 63(4), 797–804.
Makowski, D., Pham, T., Lau, Z. J., Brammer, J. C., Lespinasse, F., Pham, H., Schölzel, C., & Chen, S. H. A. (2021). NeuroKit2: A Python toolbox for neurophysiological signal processing. Behavior Research Methods, 53(4), 1689–1696.
Posada-Quintero, H. F., & Chon, K. H. (2020). Innovations in electrodermal activity data collection and signal processing: A systematic review. Sensors, 20(2), 479.
Schuurmans, A. A. T., de Looff, P., Nijhof, K. S., Rosada, C., Scholte, R. H. J., Popma, A., & Otten, R. (2020). Validity of the Empatica E4 wristband to measure heart rate variability (HRV) parameters. Journal of Medical Systems, 44(11), 190.
Tronstad, C., Amini, M., Bach, D. R., & Martinsen, Ø. G. (2022). Current trends and opportunities in the methodology of electrodermal activity measurement. Physiological Measurement, 43.

--- title: "Physiological feature extraction from wearable data: HRV, EDA, actigraphy" description: | Signal-processing layer of the ADARP analysis. Extracts 21 features from 30-min pre-EMA windows: 10 HRV indices from photoplethysmography (time-, frequency-, and non-linear domains), 6 electrodermal features via cvxEDA decomposition (Greco et al., 2016), and 5 actigraphy features from 3-axis accelerometry. Every feature is delivered in three flavours: raw, within-subject z-score, between-subject mean — enabling correct disentangling of intra- and inter-individual effects downstream. Feature-level QC is computed per window (BVP peak density, EDA plausibility) and stored as an explicit flag rather than used as a hard exclusion, so phase-15c modelling can run sensitivity analyses. Built on NeuroKit2 with cvxopt-backed convex-optimisation EDA. date: 2026-04-22 categories: [Python, Signal Processing, HRV, EDA, Actigraphy, NeuroKit2, cvxEDA, Wearable Physiology, Empatica E4, ADARP, mHealth, Feature Engineering] execute: echo: true warning: false message: false freeze: auto format: html: code-tools: true code-fold: false toc: true toc-depth: 3 number-sections: true theme: cosmo fig-cap-location: bottom embed-resources: false --- ```{r} #| label: setup #| code-summary: "Paths, packages, helpers" #| include: false PROJECT_DIR <- normalizePath( file.path("..", "files", "projet_25"), mustWork = FALSE ) DATA_MODE <- "processed" # or "demo" DATA_DIR <- file.path(PROJECT_DIR, "data", DATA_MODE) library(arrow) library(dplyr) library(tidyr) library(ggplot2) library(knitr) library(kableExtra) library(patchwork) source(file.path(PROJECT_DIR, "R", "render_helpers.R")) theme_set(theme_minimal(base_size = 11)) ``` `r data_mode_banner(DATA_DIR)` --- ## 1. From signals to features Phase 15a delivered a clean table of 58 EMA windows with synchronised physiological recordings. Phase 15b is the signal-processing layer that turns those raw streams — photoplethysmography at 64 Hz, electrodermal activity at 4 Hz, 3-axis accelerometry at 32 Hz — into **interpretable per-window features** suitable for multilevel modelling and machine-learning prediction. The design principle is simple: extract features that have **physiological meaning** (so the coefficients in a linear model are readable), **temporal smoothing** (so window-level estimates are stable), and **standardisation that preserves both between- and within-person variation** (so the analyst can later ask separately *does stress increase HRV in this person relative to their baseline?* and *do chronically more-stressed people have lower baseline HRV?*). ## 2. Feature catalogue ```{r} #| label: feature-catalogue #| code-summary: "Static table of the 21 features produced by phase 15b" #| echo: false cat_df <- tibble::tribble( ~domain, ~feature, ~unit, ~description, "HRV time", "hrv_meannn", "ms", "Mean NN interval (inverse of heart rate).", "HRV time", "hrv_sdnn", "ms", "Standard deviation of NN intervals — overall variability.", "HRV time", "hrv_rmssd", "ms", "Root mean square of successive NN differences — vagal tone.", "HRV time", "hrv_pnn50", "%", "% of NN intervals differing by >50 ms — vagal tone (discrete).", "HRV freq", "hrv_lf", "ms²", "Low-frequency power (0.04-0.15 Hz) — mixed autonomic.", "HRV freq", "hrv_hf", "ms²", "High-frequency power (0.15-0.40 Hz) — parasympathetic.", "HRV freq", "hrv_lfhf", "ratio", "Sympatho-vagal balance proxy (debated).", "HRV nonlin", "hrv_sd1", "ms", "Short-term Poincaré variability.", "HRV nonlin", "hrv_sd2", "ms", "Long-term Poincaré variability.", "HRV nonlin", "hrv_sd1sd2", "ratio", "SD1/SD2 — balance of short vs long variability.", "EDA tonic", "eda_scl_mean", "µS", "Mean skin conductance level (tonic baseline).", "EDA tonic", "eda_scl_slope", "µS/min", "Linear trend of SCL — sustained arousal.", "EDA phasic", "eda_scr_count", "count", "# of SCRs detected in window (threshold 0.05 µS).", "EDA phasic", "eda_scr_mean_amp", "µS", "Mean amplitude of detected SCRs.", "EDA phasic", "eda_scr_rate_per_min", "/min", "SCR count divided by window duration.", "EDA comp.", "eda_symp_activation", "z-sum", "Composite: SCR rate + 10*SCL slope.", "Activity", "acc_magnitude_mean", "g", "Mean of gravity-removed vector magnitude.", "Activity", "acc_magnitude_std", "g", "Standard deviation of magnitude.", "Activity", "acc_activity_count", "counts", "Mean per-10s activity count (Choi-like).", "Activity", "acc_pct_sedentary", "%", "% of 10s epochs below sedentary cutpoint.", "Activity", "acc_pct_active", "%", "% of 10s epochs above active cutpoint." ) cat_df %>% kbl(caption = "21 physiological features extracted per EMA window.", col.names = c("Domain", "Feature", "Unit", "Meaning")) %>% kable_styling(full_width = TRUE) %>% column_spec(4, width = "28em") %>% collapse_rows(columns = 1, valign = "top") ``` In addition, every feature is delivered in **three versions**: - **raw** (e.g. `hrv_rmssd`): the window-level estimate. - **within-subject z-score** (`hrv_rmssd_wz`): deviation from the subject's own mean in units of their own standard deviation. - **between-subject mean** (`hrv_rmssd_bm`): the subject's overall mean (constant across all that subject's windows). This parametrisation follows *Curran & Bauer (2011)* for disentangling between- and within-person effects. In a multilevel model, the `_wz` coefficient answers **"does a change in RMSSD within a person predict stress?"** while the `_bm` coefficient answers **"do people with chronically higher RMSSD report less stress on average?"** — two very different psychological questions. ## 3. Method details ### 3.1 HRV from photoplethysmography The E4 BVP signal is a PPG waveform sampled at 64 Hz. NeuroKit2's `ppg_clean` applies a bandpass filter (0.5–8 Hz) to remove low-frequency drift and high-frequency noise; `ppg_findpeaks` then detects systolic peaks. Successive peak-to-peak intervals (R-R) feed three HRV blocks: - `hrv_time`: standard time-domain statistics on NN intervals. - `hrv_frequency`: Welch PSD with integration in LF (0.04–0.15 Hz) and HF (0.15–0.40 Hz) bands. Requires ≥30 peaks (about 30 s of clean data); we enforce ≥30 explicitly and return `NaN` below. - `hrv_nonlinear`: the Poincaré plot of $(NN_i, NN_{i+1})$ pairs. SD1 captures short-term variability (beat-to-beat), SD2 long-term. **Wrist-worn PPG caveat.** The E4's PPG is noisier than a chest-strap ECG, so HRV indices derived from it are lower-bounded estimates. This is documented in the validation literature (Schuurmans et al., 2020) and does not preclude within-subject comparisons. ### 3.2 EDA with cvxEDA decomposition Skin conductance has two physiologically distinct components: a slow **tonic** baseline (SCL) that reflects overall arousal state, and fast **phasic** responses (SCRs) elicited by discrete events. Naïve separation — high-pass filtering — distorts both. The **cvxEDA method** (Greco, Valenza, Citi, Lanata, Scilingo, *IEEE TBME* 2016) recasts the decomposition as a convex-optimisation problem with explicit physiological priors on the sudomotor response shape. It is the reference ambulatory-EDA method in recent reviews (Posada-Quintero & Chon, 2020; Tronstad et al., 2022). NeuroKit2 wraps cvxEDA (via `cvxopt`); SCR peaks are then detected on the phasic component with an amplitude threshold of **0.05 µS** — the Boucsein (2012) standard for ambulatory recordings. We derive five per-window features plus a composite sympathetic-activation index combining SCR rate and tonic slope. ### 3.3 Actigraphy We convert the 3-axis ACC to the gravity-removed vector magnitude, then compute per-10s **activity counts** (sum of absolute deviations from the segment mean, after Choi et al. 2011 in spirit, though cutpoints are rescaled for wrist-worn placement). Within-subject ranking of epochs as *sedentary* / *light* / *active* is preserved, which is what matters for multilevel modelling. Absolute count values should **not** be compared to waist-worn ActiGraph literature. ## 4. Feature dataset ```{r} #| label: load-features #| code-summary: "Load the physiological features parquet" feat_path <- file.path(DATA_DIR, "physio_features.parquet") if (!file.exists(feat_path)) { stop("physio_features.parquet not found. Run `python extract_features.py` first.") } features <- safe_read_parquet(feat_path) n_total <- nrow(features) n_usable <- sum(features$usable_final, na.rm = TRUE) n_feat_ok <- sum(features$feature_usable, na.rm = TRUE) n_feature_cols <- sum(grepl("^(hrv_|eda_|acc_)[^_]*(_wz|_bm)?$", names(features))) cat(sprintf( "Loaded %d rows x %d columns.\n %d rows with usable_final=TRUE (from phase 15a).\n %d rows passing feature-level QC (feature_usable=TRUE).\n %d feature columns (raw + _wz + _bm).\n", n_total, ncol(features), n_usable, n_feat_ok, n_feature_cols )) ``` **What this shows.** The output preserves all 99 EMA windows for row alignment, but only **58 carry computed features** (those with `usable_final == TRUE` from phase 15a); the remaining 41 are intentionally NaN. Of those 58, **46 (79 %) pass the feature-level QC** — i.e. ≥50 % BVP and ≥50 % EDA segments judged usable. The 12 windows that compute features but fail QC will not be silently excluded: they remain in the dataset with `feature_usable = FALSE`, allowing phase-15c sensitivity analyses to compare results with and without them. ```{r} #| label: qc-per-subject #| code-summary: "Feature-level QC summary per subject" qc_path <- file.path(DATA_DIR, "extraction_qc_report.parquet") qc <- safe_read_parquet(qc_path) qc %>% kbl(caption = "Feature extraction QC per subject (restricted to windows with usable_final=TRUE from phase 15a).", col.names = c("Subject", "# Windows", "# Feat-usable", "% Feat-usable", "Median BVP QC %", "Median EDA QC %", "Median # peaks", "% NaN HRV-HF", "% NaN SCR")) %>% kable_styling(full_width = FALSE) ``` **Reading this table.** *BVP quality %* is the fraction of 10-s segments with at least 6 peaks (heart rate ≥ 36 bpm floor). *EDA quality %* is the fraction of 10-s segments in a plausible range [0.05–60 µS] with non-negligible variance. *% NaN HRV-HF* tracks windows where we had fewer than the 30 peaks required for frequency-domain HRV computation — a direct signal-quality indicator in addition to the overall BVP QC. Our results : **What this table shows.** BVP signal quality is excellent across the board (median 85–100 % per subject), confirming that PPG-based heart rate detection is robust on E4 wristband data despite the well-known PPG noise issues. **EDA quality is much more variable**: subject 104 sits at only 21.1 % (the single window for that participant has mostly out-of-range or flat skin conductance), and subjects 102, 109, and 112 hover between 46 % and 67 %. This asymmetry — BVP usable when EDA is not — is a recurring finding in ambulatory wearable studies and the rationale for tracking modality-specific QC rather than a single omnibus flag. **Operational consequence.** Subject 104 effectively drops out at this point: their only EMA window now lacks usable EDA features. **Subject 112's pattern from phase 15a propagates here**: only 50 % of its windows yield usable features, consistent with the long-but-intermittent wear pattern. Subjects 108 (85 %), 110 (100 %), and 111 (94 %) carry the bulk of the analysable signal — together they supply 30 of the 46 feature-usable windows. The "% NaN HRV-HF" column being uniformly 0 deserves a comment: every window had ≥30 BVP peaks, so frequency-domain HRV computed everywhere. This is a *necessary* condition for HRV-HF, not a *sufficient* one for its physiological interpretability — the absolute values must still be read with the wrist-PPG caveats from Schuurmans et al. (2020) in mind. ## 5. Feature distributions ```{r} #| label: hrv-distributions #| code-summary: "Distributions of HRV time- and frequency-domain features" #| fig-cap: "Distributions of the six most-used HRV features across usable_final windows. Dashed vertical lines mark the median." #| fig-width: 9 #| fig-height: 5 hrv_long <- features %>% filter(usable_final) %>% select(subject_id, hrv_meannn, hrv_rmssd, hrv_sdnn, hrv_hf, hrv_lf, hrv_lfhf) %>% pivot_longer(-subject_id, names_to = "feature", values_to = "value") %>% filter(!is.na(value)) medians <- hrv_long %>% group_by(feature) %>% summarise(m = median(value), .groups = "drop") ggplot(hrv_long, aes(x = value)) + geom_histogram(bins = 20, fill = "steelblue", alpha = 0.85) + geom_vline(data = medians, aes(xintercept = m), linetype = "dashed", colour = "red") + facet_wrap(~ feature, scales = "free") + labs(x = NULL, y = "Count") ``` **Reading the HRV distributions.** `hrv_meannn` clusters around 800–1000 ms (60–75 bpm), well within healthy adult resting range. `hrv_rmssd` is centred near 30–50 ms, also typical. The frequency-domain features (`hrv_lf`, `hrv_hf`) show right-skewed distributions on a linear scale — a known property that justifies log-transformation before linear modelling. We do *not* log-transform in storage so that analysts can choose their own transformation; phase-15c models will apply it explicitly. ```{r} #| label: eda-distributions #| code-summary: "Distributions of EDA features" #| fig-cap: "Distributions of EDA tonic, phasic, and composite features." #| fig-width: 9 #| fig-height: 4 eda_long <- features %>% filter(usable_final) %>% select(subject_id, eda_scl_mean, eda_scl_slope, eda_scr_count, eda_scr_rate_per_min, eda_symp_activation) %>% pivot_longer(-subject_id, names_to = "feature", values_to = "value") %>% filter(!is.na(value)) ggplot(eda_long, aes(x = value)) + geom_histogram(bins = 20, fill = "darkgreen", alpha = 0.85) + facet_wrap(~ feature, scales = "free", nrow = 1) + labs(x = NULL, y = "Count") ``` **Reading the EDA distributions.** Mean SCL clusters in the 1–10 µS range, consistent with healthy resting electrodermal activity. The phasic features (`eda_scr_count`, `eda_scr_rate_per_min`) are zero-inflated: many 30-min windows yielded zero SCRs above the 0.05 µS threshold. This is **not a bug**: low ambulatory arousal commonly produces few discrete SCRs, and a stricter threshold (vs the more permissive 0.01 µS we initially considered) is the cost of avoiding over-counting noise as physiological events. The `eda_symp_activation` composite inherits this zero-inflation; phase-15c will model it either with a hurdle approach or by collapsing to a "any SCR / no SCR" binary indicator depending on what the data support. ```{r} #| label: acc-distributions #| code-summary: "Distributions of actigraphy features" #| fig-cap: "Distributions of accelerometer-derived activity features." #| fig-width: 9 #| fig-height: 4 acc_long <- features %>% filter(usable_final) %>% select(subject_id, acc_magnitude_mean, acc_magnitude_std, acc_activity_count, acc_pct_sedentary, acc_pct_active) %>% pivot_longer(-subject_id, names_to = "feature", values_to = "value") %>% filter(!is.na(value)) ggplot(acc_long, aes(x = value)) + geom_histogram(bins = 20, fill = "darkorange", alpha = 0.85) + facet_wrap(~ feature, scales = "free", nrow = 1) + labs(x = NULL, y = "Count") ``` **Reading the actigraphy distributions.** Activity counts cluster heavily at the lower end, with `acc_pct_sedentary` typically above 50 %, consistent with seated ambulatory contexts (the 30-min pre-EMA window often catches participants at home, at work desks, or otherwise stationary). The longer right tails on `acc_magnitude_std` and `acc_pct_active` indicate that some windows do capture genuine movement — which is precisely the behavioural variability we want for contrasting against the within-subject baseline. ## 6. Between- vs within-person variability A defining feature of EMA+physio data is that most of the variance is **between individuals**, not within. We quantify this with the **intra-class correlation** (ICC): the fraction of feature variance attributable to stable between-subject differences. An ICC of 0.80 means 80 % of the variance in, say, RMSSD, is "who you are" and only 20 % is "how you are right now". ```{r} #| label: icc-table #| code-summary: "Per-feature intra-class correlation (between-subject variance fraction)" feature_names <- c( "hrv_meannn","hrv_rmssd","hrv_sdnn","hrv_pnn50","hrv_hf","hrv_lf","hrv_lfhf", "hrv_sd1","hrv_sd2","hrv_sd1sd2", "eda_scl_mean","eda_scl_slope","eda_scr_count","eda_scr_rate_per_min","eda_symp_activation", "acc_magnitude_mean","acc_magnitude_std","acc_activity_count", "acc_pct_sedentary","acc_pct_active" ) compute_icc <- function(df, feat) { x <- df[[feat]] sid <- df$subject_id ok <- !is.na(x) if (sum(ok) < 5) return(NA_real_) # One-way ANOVA-style ICC1: between-sd^2 / (between-sd^2 + within-sd^2) d <- data.frame(value = x[ok], subject = factor(sid[ok])) if (nlevels(d$subject) < 2) return(NA_real_) means <- tapply(d$value, d$subject, mean) ns <- tapply(d$value, d$subject, length) grand_mean <- mean(d$value) ss_between <- sum(ns * (means - grand_mean)^2) ss_within <- sum((d$value - ave(d$value, d$subject))^2) df_between <- nlevels(d$subject) - 1 df_within <- nrow(d) - nlevels(d$subject) if (df_within <= 0) return(NA_real_) ms_between <- ss_between / df_between ms_within <- ss_within / df_within k <- mean(ns) # average group size icc <- (ms_between - ms_within) / (ms_between + (k - 1) * ms_within) max(0, min(1, icc)) # truncate numerical overflow } icc_tbl <- tibble::tibble( feature = feature_names, icc = sapply(feature_names, function(f) compute_icc(features %>% filter(usable_final), f)) ) %>% arrange(desc(icc)) icc_tbl %>% mutate(icc = round(icc, 3)) %>% kbl(caption = "Intra-class correlation per feature (1 = all variance between-subject, 0 = all variance within).", col.names = c("Feature", "ICC (one-way)")) %>% kable_styling(full_width = FALSE) ``` ```{r} #| label: icc-plot #| code-summary: "Bar plot of ICCs, sorted" #| fig-cap: "Intra-class correlation per feature. High ICC = most variance is between subjects (trait-like); low ICC = signal varies within subject over time (state-like, the kind of variance that predicts EMA fluctuations)." #| fig-width: 7 #| fig-height: 5 ggplot(icc_tbl, aes(x = reorder(feature, icc), y = icc)) + geom_col(fill = "steelblue", alpha = 0.85) + geom_hline(yintercept = 0.5, linetype = "dashed", colour = "grey40") + coord_flip() + labs(x = NULL, y = "ICC (between-subject variance fraction)") ``` **How to read this.** Features with ICC above 0.5 (grey dashed line) are dominated by stable between-subject differences — they characterise the person more than the moment. Features with low ICC fluctuate within subjects and are the **state-like candidates** for predicting momentary stress. In a multilevel model, high-ICC features enter primarily through their between-subject mean (`_bm` columns); low-ICC features carry their predictive signal in their within-subject deviation (`_wz` columns). Our results : **What this shows.** The highest ICCs are concentrated on the accelerometer features (`acc_magnitude_std` = 0.30) and the slowest HRV index (`hrv_meannn` = 0.30), with EDA tonic baseline (`eda_scl_mean`) close behind at 0.24. Time-domain HRV indices that respond to short-term variability (RMSSD, SDNN, pNN50) sit at the bottom of the range with ICCs around 0.05. These values were **cross-validated against R's `lme4::lmer()` gold-standard random-intercept estimator**, which returned 0.33, 0.40, 0.32, and 0.06 respectively — all within 0.10 of our Python ICC1 estimates. The slight conservative bias of one-way ICC on unbalanced designs (here: 1 to 17 windows per subject) is well-documented and does not undermine the substantive interpretation. The values are at the lower end of what the EMA-physiology literature typically reports (e.g. Bertsch et al. 2012 cite HRV time-domain ICCs around 0.40–0.60 for similarly-paced designs). Two factors plausibly explain the gap: 1. **Wrist-PPG noise floor.** The E4 PPG is known to attenuate between-subject HRV differences relative to chest-strap ECG (Schuurmans et al. 2020). This specifically depresses HRV ICCs and matches the pattern in our table — the lowest ICCs are HRV time-domain features that depend most directly on beat-to-beat precision. 2. **Sample size.** With 9 subjects and 46 feature-usable windows, ICC point estimates carry wide confidence intervals; the reported values are themselves uncertain. **Phase-15c implication.** With ICCs in this range, the within-subject deviation columns (`_wz`) will carry most of the predictive signal for moment-to-moment stress events. The between-subject mean columns (`_bm`) will mostly absorb random-intercept variation in the multilevel models. We will report both as fixed effects and let the data speak. ## 7. Outputs and handover to phase 15c Phase 15b produces two parquets in `data/processed/`: ```{r} #| label: outputs-summary #| code-summary: "Summary of phase 15b outputs" #| echo: false tibble::tribble( ~file, ~rows, ~description, "physio_features.parquet", nrow(features), "Per-EMA-window features (raw + _wz + _bm) plus QC flags.", "extraction_qc_report.parquet", nrow(qc), "Per-subject feature-extraction QC summary." ) %>% kbl(caption = "Phase-15b outputs.", col.names = c("File", "Rows", "Description")) %>% kable_styling(full_width = TRUE) %>% column_spec(3, width = "30em") ``` **Ready for phase 15c.** `physio_features.parquet` is a tidy subject-by-window table with 21 × 3 = 63 feature columns, the binary outcome `stress_event`, the ordinal outcome `stress_composite`, the subject identifier for random-effects modelling, and QC flags. Phase 15c will fit: 1. **Multilevel logistic and ordinal models** (lme4 + brms) on `stress_event` and `stress_composite`, using the `_wz`/`_bm` decomposition to separate momentary from trait effects. 2. **LightGBM with leave-one-subject-out cross-validation** for predictive benchmarking, with SHAP interpretability. 3. **Sensitivity analyses** contrasting the full sample vs the `feature_usable == TRUE` subset. --- ## 8. References - Boucsein, W. (2012). *Electrodermal activity* (2nd ed.). Springer. - Choi, L., Liu, Z., Matthews, C. E., & Buchowski, M. S. (2011). Validation of accelerometer wear and nonwear time classification algorithm. *Medicine & Science in Sports & Exercise, 43*(2), 357–364. - Curran, P. J., & Bauer, D. J. (2011). The disaggregation of within-person and between-person effects in longitudinal models of change. *Annual Review of Psychology, 62*, 583–619. - Greco, A., Valenza, G., Citi, L., Lanata, A., & Scilingo, E. P. (2016). cvxEDA: A convex optimization approach to electrodermal activity processing. *IEEE Transactions on Biomedical Engineering, 63*(4), 797–804. - Makowski, D., Pham, T., Lau, Z. J., Brammer, J. C., Lespinasse, F., Pham, H., Schölzel, C., & Chen, S. H. A. (2021). NeuroKit2: A Python toolbox for neurophysiological signal processing. *Behavior Research Methods, 53*(4), 1689–1696. - Posada-Quintero, H. F., & Chon, K. H. (2020). Innovations in electrodermal activity data collection and signal processing: A systematic review. *Sensors, 20*(2), 479. - Schuurmans, A. A. T., de Looff, P., Nijhof, K. S., Rosada, C., Scholte, R. H. J., Popma, A., & Otten, R. (2020). Validity of the Empatica E4 wristband to measure heart rate variability (HRV) parameters. *Journal of Medical Systems, 44*(11), 190. - Tronstad, C., Amini, M., Bach, D. R., & Martinsen, Ø. G. (2022). Current trends and opportunities in the methodology of electrodermal activity measurement. *Physiological Measurement, 43*.