BCI Classification

Introduction

Brain-Computer Interfaces (BCIs) enable direct communication between the brain and external devices by translating neural activity into control signals. This vignette demonstrates how to analyze and classify BCI data using PhysioEEG, focusing on motor imagery and steady-state visual evoked potential (SSVEP) paradigms.

Motor Imagery Paradigm

Motor imagery BCI relies on event-related desynchronization (ERD) and event-related synchronization (ERS) in the mu (8-12 Hz) and beta (13-30 Hz) frequency bands. When users imagine left or right hand movements, characteristic power decreases occur over contralateral sensorimotor cortex.

SSVEP Paradigm

SSVEP BCIs use flickering visual stimuli at specific frequencies to elicit steady-state responses in visual cortex. Different flicker frequencies enable multi-class control by detecting which frequency dominates the EEG spectrum.

Classification Pipeline

A typical BCI classification workflow involves: 1. Preprocessing (filtering, artifact removal) 2. Feature extraction (CSP, bandpower, Riemannian geometry) 3. Classification (LDA, SVM with cross-validation) 4. Performance evaluation (accuracy, kappa, information transfer rate)

Setup

library(PhysioEEG)
library(SummarizedExperiment)
library(ggplot2)
library(dplyr)

Create simulated BCI data with motor imagery trials (left vs. right hand):

# Generate 80 trials: 40 left hand, 40 right hand
# make_eeg_bci() creates n_trials PER class, so n_trials = 40 yields 80 total
# 16 channels, 250 Hz sampling rate, 4-second trials
eeg <- make_eeg_bci(
  n_trials = 40,
  n_channels = 16,
  sr = 250,
  trial_sec = 4
)

# Trial condition labels ("left"/"right") are stored in metadata
condition <- metadata(eeg)$labels

# Examine data structure
eeg
dim(assay(eeg, "raw"))  # 1000 timepoints x 16 channels x 80 trials

# Check trial labels
table(condition)  # Should show balanced left/right classes

Motor Imagery Analysis

Motor imagery produces ERD during movement imagination and ERS after task completion. The eegMotorImagery() function computes power changes in motor-related frequency bands.

# Compute ERD/ERS in mu and beta bands
eeg <- eegMotorImagery(
  eeg,
  bands = list(
    mu = c(8, 12),
    beta = c(13, 30)
  ),
  assay_name = "raw"
)

# ERD/ERS results stored in metadata as an n_trials x (n_channels * n_bands)
# matrix with column names "<channel>_<band>" (e.g. "C3_mu")
erd_results <- metadata(eeg)$erd_ers

# Examine mu-band power changes (one column per channel)
mu_cols <- grep("_mu$", colnames(erd_results))
mu_power <- erd_results[, mu_cols, drop = FALSE]
head(mu_power)

# Laterality index LI = (R - L) / (|R| + |L|) from C4 (right) vs C3 (left).
# Positive values indicate right-hemisphere dominance (left hand imagery),
# negative values indicate left-hemisphere dominance (right hand imagery).
lat <- eegLateralization(
  eeg,
  left_ch = "C3",
  right_ch = "C4",
  band = c(8, 12),
  method = "power"
)
laterality <- lat$per_trial   # per-trial LI (n_trials x 1)
lat$summary                   # biomarker summary (mean LI with CI)

Visualize ERD/ERS patterns:

# Plot per-channel mu-band ERD/ERS for left vs. right conditions
library(ggplot2)

# Channel labels for the mu-band columns
ch_names <- sub("_mu$", "", colnames(mu_power))

# Average mu-band power change per channel for each condition
mu_left <- colMeans(mu_power[condition == "left", , drop = FALSE])
mu_right <- colMeans(mu_power[condition == "right", , drop = FALSE])

df_plot <- data.frame(
  channel = factor(rep(ch_names, 2), levels = ch_names),
  power = c(mu_left, mu_right),
  condition = rep(c("Left", "Right"), each = length(ch_names))
)

ggplot(df_plot, aes(x = channel, y = power, fill = condition)) +
  geom_col(position = "dodge") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    title = "Motor Imagery ERD/ERS (Mu Band)",
    x = "Channel",
    y = "Power Change (%)",
    fill = "Condition"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Common Spatial Patterns (CSP)

CSP is a powerful technique for extracting spatial filters that maximize variance differences between two classes. It finds projections where one class has high variance and the other has low variance.

# Apply bandpass filter to isolate the mu+beta band.
# eegFilter() operates on continuous (2D) data, so band-pass each epoch's
# channels (8-30 Hz, zero-phase Butterworth) into a new "filtered" assay.
sr <- samplingRate(eeg)
bf <- signal::butter(4, c(8, 30) / (sr / 2), type = "pass")
raw3d <- assay(eeg, "raw")
assay(eeg, "filtered") <- array(
  apply(raw3d, c(2, 3), function(sig) signal::filtfilt(bf, sig)),
  dim = dim(raw3d)
)

# Extract CSP features
# n_filters = 3 selects the 3 top + 3 bottom eigenvectors = 6 components
eeg <- eegCSP(
  eeg,
  labels = condition,
  n_filters = 3,
  assay_name = "filtered"
)

# CSP spatial filters stored in metadata as (n_components x n_channels);
# transpose to view as channels x components
csp_filters <- t(metadata(eeg)$csp$filters)
dim(csp_filters)  # 16 channels x 6 components

# Projected data (log-variance features for classification)
csp_features <- metadata(eeg)$csp_features
dim(csp_features)  # 80 trials x 6 features

Understanding CSP components:

# First components maximize variance for class 1 (left hand)
# Last components maximize variance for class 2 (right hand)

# Visualize first CSP filter (most discriminative for left)
filter_1 <- csp_filters[, 1]
names(filter_1) <- colData(eeg)$label

# Create topographic visualization data
topo_data <- data.frame(
  channel = names(filter_1),
  weight = filter_1
)

print(topo_data)

# Examine CSP feature distributions
boxplot(
  csp_features[, 1] ~ condition,
  main = "CSP Component 1 (Left-specific)",
  xlab = "Condition",
  ylab = "Log Variance"
)

boxplot(
  csp_features[, 6] ~ condition,
  main = "CSP Component 6 (Right-specific)",
  xlab = "Condition",
  ylab = "Log Variance"
)

Feature Extraction

The eegBCIfeatures() function provides multiple feature extraction methods that can be combined for robust classification.

# eegBCIfeatures() extracts one method per call and returns a matrix.
# Collect the three feature families into a named list.
features_all <- list(
  bandpower  = eegBCIfeatures(eeg, method = "bandpower",  labels = condition,
                              assay_name = "filtered"),
  csp        = eegBCIfeatures(eeg, method = "csp",        labels = condition,
                              assay_name = "filtered"),
  riemannian = eegBCIfeatures(eeg, method = "riemannian", labels = condition,
                              assay_name = "filtered")
)

# Each list element is a feature matrix
names(features_all)

# Bandpower features: power in specified frequency bands per channel
bp_features <- features_all$bandpower
dim(bp_features)  # 80 trials x (16 channels * n_bands)

# CSP features: log-variance of CSP-projected signals
csp_features <- features_all$csp
dim(csp_features)  # 80 trials x 6 components

# Riemannian features: tangent space projection of covariance matrices
riem_features <- features_all$riemannian
dim(riem_features)  # 80 trials x n_tangent_features

Bandpower-specific feature extraction:

# Extract power in the default mu/beta frequency bands
features_bp <- eegBCIfeatures(
  eeg,
  method = "bandpower",
  labels = condition,
  assay_name = "raw"
)

# Bandpower features are log band power per channel and band
# (default bands: mu 8-13 Hz, beta 13-30 Hz)

# Normalize features (z-score per feature dimension)
features_norm <- scale(features_bp)

# Examine feature importance via correlation with labels
label_numeric <- ifelse(condition == "left", 0, 1)
correlations <- cor(features_norm, label_numeric)
head(sort(abs(correlations), decreasing = TRUE))

Riemannian geometry features:

# Riemannian approach treats covariance matrices as points
# on a manifold and projects to tangent space

features_riem <- eegBCIfeatures(
  eeg,
  method = "riemannian",
  labels = condition,
  assay_name = "filtered"
)

# These features often outperform traditional methods
# for BCI classification due to better handling of
# non-Euclidean structure of covariance matrices

riem_data <- features_riem
dim(riem_data)

SSVEP Detection

For frequency-tagged visual stimuli, SSVEP analysis identifies the target frequency in the EEG spectrum.

# Detect SSVEP across candidate flicker frequencies using CCA.
# Each candidate uses 3 harmonics in its reference signals.
ssvep_results <- eegSSVEP(
  eeg,
  frequencies = c(12, 15, 20),
  n_harmonics = 3,
  assay_name = "raw"
)

# Results: one row per candidate frequency with CCA correlation and SNR
ssvep_results

# SNR at each candidate frequency (named by frequency)
snr_values <- ssvep_results$snr
names(snr_values) <- ssvep_results$frequency

# High SNR indicates a strong SSVEP response
mean(snr_values)
sd(snr_values)

# Binary detection: SSVEP present vs. absent (threshold at SNR = 3)
ssvep_detected <- snr_values > 3
table(ssvep_detected)

Multi-frequency SSVEP (for multi-class BCI):

# In a multi-class SSVEP BCI, test several candidate target frequencies at once
target_freqs <- c(7, 9, 11, 13)  # 4-class BCI

ssvep_all <- eegSSVEP(eeg, frequencies = target_freqs, n_harmonics = 2,
                      assay_name = "raw")

# The candidate frequency with the highest SNR is the detected class
predicted_freq <- ssvep_all$frequency[which.max(ssvep_all$snr)]
predicted_freq

Classification

The eegBCIclassify() function performs classification with cross-validation to estimate generalization performance.

# eegBCIclassify() returns a per-trial data.frame of predictions
# (columns: trial, predicted_class, confidence, true_class). Derive
# accuracy, Cohen's kappa, and the confusion matrix from those predictions.
bci_metrics <- function(res) {
  lv <- sort(unique(c(res$true_class, res$predicted_class)))
  cm <- table(
    true = factor(res$true_class, levels = lv),
    predicted = factor(res$predicted_class, levels = lv)
  )
  acc <- sum(diag(cm)) / sum(cm)
  pe <- sum(rowSums(cm) * colSums(cm)) / sum(cm)^2
  kappa <- if (pe < 1) (acc - pe) / (1 - pe) else NA_real_
  list(accuracy = acc, kappa = kappa, confusion_matrix = cm)
}

# Classify using CSP features and LDA
# 5-fold cross-validation for robust performance estimation
results_lda <- eegBCIclassify(
  eeg,
  features = csp_features,
  labels = condition,
  method = "lda",
  cv_folds = 5
)

# Extract performance metrics
m_lda <- bci_metrics(results_lda)
accuracy_lda <- m_lda$accuracy
kappa_lda <- m_lda$kappa
confusion_lda <- m_lda$confusion_matrix

cat(sprintf("LDA Accuracy: %.2f%%\n", accuracy_lda * 100))
cat(sprintf("Cohen's Kappa: %.3f\n", kappa_lda))
print(confusion_lda)

Compare multiple classifiers:

# LDA: Linear Discriminant Analysis (assumes Gaussian distributions)
results_lda <- eegBCIclassify(
  eeg,
  features = csp_features,
  labels = condition,
  method = "lda",
  cv_folds = 10
)

# Shrinkage LDA: regularized covariance estimate, robust for high-dimensional
# feature spaces with limited trials
results_slda <- eegBCIclassify(
  eeg,
  features = csp_features,
  labels = condition,
  method = "shrinkage_lda",
  cv_folds = 10
)

# Compare performance
m_lda <- bci_metrics(results_lda)
m_slda <- bci_metrics(results_slda)
comparison <- data.frame(
  Method = c("LDA", "Shrinkage-LDA"),
  Accuracy = c(m_lda$accuracy, m_slda$accuracy),
  Kappa = c(m_lda$kappa, m_slda$kappa)
)

print(comparison)

Feature comparison for classification:

# Classify using different feature sets
# 1. CSP features
acc_csp <- bci_metrics(eegBCIclassify(
  eeg,
  features = features_all$csp,
  labels = condition,
  method = "lda",
  cv_folds = 5
))$accuracy

# 2. Bandpower features
acc_bp <- bci_metrics(eegBCIclassify(
  eeg,
  features = features_all$bandpower,
  labels = condition,
  method = "lda",
  cv_folds = 5
))$accuracy

# 3. Riemannian features
acc_riem <- bci_metrics(eegBCIclassify(
  eeg,
  features = features_all$riemannian,
  labels = condition,
  method = "lda",
  cv_folds = 5
))$accuracy

# Compare feature sets
feature_comparison <- data.frame(
  Features = c("CSP", "Bandpower", "Riemannian"),
  Accuracy = c(acc_csp, acc_bp, acc_riem)
)

print(feature_comparison)

Performance Evaluation

Classification Accuracy

Accuracy is the proportion of correctly classified trials. For balanced two-class problems, chance level is 50%.

# Overall accuracy
results <- eegBCIclassify(
  eeg,
  features = csp_features,
  labels = condition,
  method = "lda",
  cv_folds = 5
)

m <- bci_metrics(results)
overall_acc <- m$accuracy

# Per-class accuracy (sensitivity/recall)
conf_mat <- m$confusion_matrix
class_acc_left <- conf_mat["left", "left"] / sum(conf_mat["left", ])
class_acc_right <- conf_mat["right", "right"] / sum(conf_mat["right", ])

cat(sprintf("Overall Accuracy: %.2f%%\n", overall_acc * 100))
cat(sprintf("Left Hand Accuracy: %.2f%%\n", class_acc_left * 100))
cat(sprintf("Right Hand Accuracy: %.2f%%\n", class_acc_right * 100))

Cohen’s Kappa

Kappa accounts for chance agreement and is more informative than raw accuracy for imbalanced datasets.

kappa <- m$kappa

# Interpretation:
# < 0.00: Poor agreement (worse than chance)
# 0.00-0.20: Slight agreement
# 0.21-0.40: Fair agreement
# 0.41-0.60: Moderate agreement
# 0.61-0.80: Substantial agreement
# 0.81-1.00: Almost perfect agreement

cat(sprintf("Cohen's Kappa: %.3f\n", kappa))

if (kappa < 0.4) {
  cat("Classification performance: Fair\n")
} else if (kappa < 0.6) {
  cat("Classification performance: Moderate\n")
} else if (kappa < 0.8) {
  cat("Classification performance: Substantial\n")
} else {
  cat("Classification performance: Excellent\n")
}

Information Transfer Rate (ITR)

ITR measures the communication speed of a BCI in bits per minute. It depends on accuracy, number of classes, and selection time.

# ITR formula for binary classification:
# ITR = (log2(N) + P*log2(P) + (1-P)*log2((1-P)/(N-1))) * (60/T)
# where N = number of classes, P = accuracy, T = trial duration (sec)

calculate_itr <- function(accuracy, n_classes, trial_duration) {
  if (accuracy <= 1/n_classes) return(0)  # Below chance

  P <- accuracy
  N <- n_classes

  # Bits per trial
  bits_per_trial <- log2(N) + P * log2(P) +
                    (1 - P) * log2((1 - P) / (N - 1))

  # Bits per minute
  itr <- bits_per_trial * (60 / trial_duration)
  return(itr)
}

# For our 2-class BCI with 4-second trials
itr <- calculate_itr(
  accuracy = overall_acc,
  n_classes = 2,
  trial_duration = 4
)

cat(sprintf("Information Transfer Rate: %.2f bits/min\n", itr))

# Compare different scenarios
scenarios <- expand.grid(
  accuracy = seq(0.6, 0.95, by = 0.05),
  trial_duration = c(3, 4, 5)
)

scenarios$itr <- mapply(
  calculate_itr,
  accuracy = scenarios$accuracy,
  trial_duration = scenarios$trial_duration,
  MoreArgs = list(n_classes = 2)
)

# Find optimal configuration
best <- scenarios[which.max(scenarios$itr), ]
cat(sprintf("\nOptimal configuration:\n"))
cat(sprintf("  Accuracy: %.2f%%\n", best$accuracy * 100))
cat(sprintf("  Trial Duration: %.1f s\n", best$trial_duration))
cat(sprintf("  ITR: %.2f bits/min\n", best$itr))

Cross-Validation Strategies

Proper cross-validation is critical for reliable performance estimation.

# 5-fold CV: Standard for small datasets (n < 100)
cv5 <- eegBCIclassify(
  eeg,
  features = csp_features,
  labels = condition,
  method = "lda",
  cv_folds = 5
)

# 10-fold CV: More reliable estimate, higher computational cost
cv10 <- eegBCIclassify(
  eeg,
  features = csp_features,
  labels = condition,
  method = "lda",
  cv_folds = 10
)

# Leave-one-out CV: Maximum variance reduction (n = n_trials)
# Only for small datasets due to computational cost
cv_loo <- eegBCIclassify(
  eeg,
  features = csp_features,
  labels = condition,
  method = "lda",
  cv_folds = nrow(csp_features)  # 80 folds = leave-one-out
)

cv_comparison <- data.frame(
  Method = c("5-Fold", "10-Fold", "Leave-One-Out"),
  Accuracy = c(bci_metrics(cv5)$accuracy, bci_metrics(cv10)$accuracy,
               bci_metrics(cv_loo)$accuracy),
  Kappa = c(bci_metrics(cv5)$kappa, bci_metrics(cv10)$kappa,
            bci_metrics(cv_loo)$kappa)
)

print(cv_comparison)

Summary

This vignette demonstrated a complete BCI classification pipeline:

Data preparation: Simulated motor imagery EEG with left/right hand classes
Motor imagery analysis: ERD/ERS computation in mu and beta bands
Spatial filtering: CSP for optimal discrimination between classes
Feature extraction: Bandpower, CSP, and Riemannian geometry methods
SSVEP detection: Frequency-domain analysis for visual BCIs
Classification: LDA and SVM with cross-validation
Performance metrics: Accuracy, kappa, and information transfer rate

For real BCI applications, additional considerations include: - Online vs. offline analysis pipelines - Adaptive classifiers that update with user feedback - Multi-session training and transfer learning - Real-time artifact rejection - User-specific calibration and optimization

The PhysioEEG package provides flexible tools for exploring these advanced topics in BCI research.

PhysioEEG Development Team

2026-08-25