Add collective behavior metrics

[khan-u]

Summary

Analyzing collective behavior in multi-animal tracking data requires quantifying group-level patterns beyond individual kinematics.

This issue proposes implementing 4 core collective behavior metrics + 2 supporting metrics grounded in foundational frameworks:

Core Metrics quantify collective state and hierarchical structure:

• Polarization, Milling: Distinguish swarm, torus, and dynamic/highly parallel groups [2]
• Group Spread: Quantify spatial dispersion via radius of gyration
• Leadership: Reveal influence hierarchies via velocity cross-correlation [3]

Supporting Metrics combine with Core Metrics to reveal pairwise coordination mechanisms:

• Egocentric Angle: Conspecific position in self-centered reference frame [5]
• Radial/Tangent Velocity: Decompose relative motion into radial and tangential components [5]

Design properties:

1. Framework-agnostic: Works with SLEAP, DeepLabCut, LightningPose, Anipose outputs
2. Species-agnostic: Generalizes across mice, fish, birds, insects, primates
3. N-individuals: Supports collective behavior analysis of 2 to N individuals

Proposed Metrics

Core Metrics

Group-level metrics quantifying collective state (polarization, milling, spread) and structure (leadership).

Metric	Question Answered	Formula
Polarization	Are their body orientations or movement directions aligned?	Magnitude of summed unit direction vectors (normalized)
Milling	Are they rotating around a center?	Magnitude of summed angular momentum (normalized)
Group Spread	How dispersed is the group?	Radius of gyration: RMS distance from centroid
Leadership	Who influences whom?	Pairwise velocity cross-correlation at time lag τ

Polarization quantifies alignment of direction vectors:

$$\Phi = \frac{1}{N} \left\| \sum_{i=1}^{N} \hat{u}_i \right\|$$

where $\hat{u}_i$ is the unit direction vector for individual $i$, and $N$ is the number of valid individuals at each time point.

Output range:

• 1.0 = perfect alignment (all direction vectors point the same way)
• 0.0 = no alignment (random directions or exact cancellation)

       High Polarization (≈1)          Low Polarization (≈0)
        Aligned directions             Random/canceling directions

            ● →                           ● →
            ● →                         ← ●
            ● →                           ● ↓
            ● →                           ● ↑


     Magnitude of sum ≈ N            Magnitude of sum ≈ 0

Polarization can be computed in two modes:

(1) Orientation polarization (body-axis mode): Alignment of body orientations (which way they're facing) from body_axis_keypoints=(origin, target) like ("tail_base", "neck"). Computes orientation as target - origin within a single frame, so it works from the first frame (index 0). Optional validate_ap=True runs prior-free anterior-posterior axis validation to verify or suggest keypoint pairs (see #945).

(2) Heading polarization (displacement mode): Alignment of movement directions (which way they're traveling) when body_axis_keypoints=None. Computes direction as position[t] - position[t - displacement_frames], so frames at indices 0 through displacement_frames - 1 are NaN (no prior frame to compare).

Parameters

Parameter	Type	Default	Description
data	xarray.DataArray	—	Position data with time, space, and individuals dims. Requires x, y space coordinates. Must include keypoints dim when using body_axis_keypoints; for heading polarization, pre-select a keypoint via .sel(keypoints="...") else the first (index 0) will be used.
body_axis_keypoints	tuple[Hashable, Hashable]	None	(origin, target) keypoint pair defining body axis for orientation polarization. When None, computes heading polarization from displacement.
displacement_frames	int	1	Frame interval for displacement-based heading polarization. Ignored when body_axis_keypoints is set.
return_angle	bool	False	Return mean angle as second output. Returns mean body orientation angle (orientation polarization) or mean movement direction angle (heading polarization).
in_degrees	bool	False	If True, mean angle is returned in degrees; otherwise radians. Only relevant when return_angle=True.
validate_ap	bool	False	If True, run anterior-posterior axis validation when body_axis_keypoints is provided. Results stored in polarization.attrs["ap_validation_result"].
ap_validation_config	dict	None	Configuration overrides for AP validation. See movement.kinematics.body_axis.ValidateAPConfig for options.

Returns

Condition	Return Value
return_angle=False	DataArray named "polarization", dims ("time",), values clipped to [0, 1]
return_angle=True	Tuple (polarization, mean_angle) where mean_angle is a DataArray named "mean_angle", dims ("time",)

Edge Case Handling

Condition	Behavior
Missing data (NaN)	Excluded per individual, per frame
Zero-length direction vector	Excluded (stationary individuals in heading mode, or coincident keypoints in orientation mode)
All direction vectors invalid	Returns NaN for that frame
Exact vector cancellation	Polarization = 0, angle = NaN
3D position data	Z-coordinate ignored; computed in x/y plane
Heading polarization, first N frames	NaN (no prior frame for displacement)

---

Milling quantifies rotational motion around the group center:

        High Milling (≈1)              Low Milling (≈0)
     Coordinated rotation           Random/linear motion

           ●→                            ● →
                                       ↗
        ↑●  ⊙  ●↓                    ●       ● ↓
                                         ↘
           ←●                          ← ●

     ⊙ = group centroid           Arrows = direction vectors
     Rotating clockwise            No consistent rotation
     (viewed from above)

Each individual contributes angular momentum relative to the centroid. When all rotate in the same direction (clockwise or counter-clockwise), angular momenta sum constructively → high milling. Random or linear motion → angular momenta cancel → low milling.

Like polarization, milling can be computed in two modes:

(1) Displacement mode (default, body_axis_keypoints=None): Uses velocity vectors to measure actual rotational motion. This is the standard physics definition of angular momentum (r × v).

(2) Body-axis mode: Uses body orientation vectors to measure rotational posture (tangential alignment). Useful when position data is noisy or to detect milling readiness before movement begins.

Displacement mode is typically preferred since milling fundamentally describes rotational motion, not just rotational orientation.

---

Group Spread quantifies spatial dispersion via radius of gyration:

        High Spread                     Low Spread
      (large Rg value)               (small Rg value)

    ●               ●                     ● ●
            ⊙                              ⊙
    ●               ●                     ● ●

     ⊙ = centroid                    ⊙ = centroid
     Large RMS distance              Small RMS distance

Collective State Space

Polarization and milling define the collective state; spread quantifies spatial extent:

Polarization	Milling	Spread	Collective State
Low	Low	High	Dispersed swarm
High	Low	Low	Tight school
Low	High	Low	Compact torus/mill
High	Low	High	Loose migrating flock

Note: Couzin et al. 2002 defines four states based on polarization and milling alone: swarm, torus, dynamic parallel group, and highly parallel group. The table above extends this framework by adding spread as a third dimension to distinguish spatial configurations within these states.

Leadership

Leadership is computed via velocity cross-correlation at time lag τ:

    Time →

    Individual i (Leader):     →   →   ↗   ↗   ↑   ↑   ↖   ↖   ←   ←
                               t=0 t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9

    Individual j (Follower):   ↓   →   →   ↗   ↗   ↑   ↑   ↖   ↖   ←
                               t=0 t=1 t=2 t=3 t=4 t=5 t=6 t=7 t=8 t=9
                                   └───────────────────────────────┘
                                    j copies i's velocity with τ=1 lag

    
     j's velocity at time t resembles i's velocity from 1 frame earlier → i leads j by τ=1 frame                              

The algorithm computes correlation between velocity_i(t) and velocity_j(t + τ) across a range of lags. The lag τ with maximum correlation indicates the temporal delay in influence.

Interpreting τ:

Cross-correlation peak	Influence direction	Meaning
Peak at +τ	i leads j	i's velocity predicts j's future velocity
Peak at −τ	j leads i	j's velocity predicts i's future velocity
Peak at τ≈0	Simultaneous	No leader-follower relationship

Leadership reveals hierarchical structure that may interact with collective state.

---

Leadership + Collective State Interaction

In pigeon flocks, leaders consistently occupy front positions, as spatial position strongly correlates with hierarchical rank. [3] Whether leadership varies across different collective states (polarized vs. unpolarized) is an open empirical question.

Hypothesized interactions (to be tested with combined metrics):

Leadership	+ State Context	Interpretation
Clear leader (+τ peak)	High polarization	Coordinated migration following a leader
No clear leader (τ≈0)	High polarization	Synchronized movement without hierarchy
Leader initiates	State transition	Leader drives shifts (e.g., swarm → school)

---

Supporting + Core Metrics Interaction

Focal-other pairwise metrics (computed for all N×(N-1) dyads) designed to combine with Core metrics, revealing how pairs coordinate to achieve collective states.

Metric	Questions Answered	Formula
Egocentric Angle	Where are others relative to my direction? • During polarization: V-formation or single-file? • During milling: tracking neighbor or monitoring center? • Social attention allocation across roles	Angle between focal's direction vector and vector to conspecific (see below)
Radial/Tangent Velocity	Are they approaching or passing by? • Leader approaching followers or vice versa? • Chase/flee vs. parallel travel? • Approach-avoidance dynamics during state transitions	Decompose relative velocity into radial and tangential components (see below)

Egocentric Angle measures the angular position of a conspecific relative to the focal individual's direction:

              Conspecific B
                   ●
                  /
                 /
                /  θ = egocentric angle
               /
    Focal A  ●───────────→  (A's direction vector)

    θ = 0°   → B is directly ahead of A
    θ = 90°  → B is to A's side
    θ = 180° → B is directly behind A

This reveals spatial arrangements: V-formation (neighbors at ~45°), single-file (leader ahead at 0°, follower behind at 180°), or side-by-side travel (~90°).

---

Radial/Tangent Velocity decomposes the relative velocity between two individuals into two orthogonal components:

• Radial component: Velocity along the line connecting the pair, defined as -d(distance)/dt.

  • Positive → closing distance (approaching)
  • Negative → increasing distance (separating)

• Tangential component: Velocity perpendicular to the connecting line.

  • Captures lateral/sideways motion ("passing by")
  • High tangential + low radial → orbiting or parallel travel

                B
                ●
                │
                │  ← radial axis (A─B line)
                │     radial = -d(distance)/dt
                │     + = closing, − = separating
                │
                ●─────────  ← tangential axis
                A              (perpendicular to A─B)

Together they answer: "Are A and B approaching, separating, or passing by each other?"

---

Relationship to Existing ROI Module

The movement.roi module already implements egocentric angle and approach vector computations for individual-to-region relationships:

Existing (ROI module)	Proposed (Collective module)
compute_egocentric_angle_to_nearest_point()	compute_egocentric_angle()
compute_approach_vector()	compute_radial_tangent_velocity()
Reference: static region (polygon/line)	Reference: other moving individuals
Output: (time, individuals, keypoints)	Output: (time, individuals, individuals_other)
Use case: "Where is the nest relative to me?"	Use case: "Where is my neighbor relative to me?"

The underlying math is similar—both rely on compute_signed_angle_2d from movement.utils.vector. The collective functions can reuse this utility while extending the pattern to pairwise individual relationships across all N×(N-1) dyads.

See existing examples for API patterns to follow:

• boundary_angles.py: compute_velocity(positions) as forward vector, in_degrees=True for angle output
• compute_kinematics.py: compute_norm() for vector magnitudes
• compute_head_direction.py: bilateral keypoint patterns (left_keypoint, right_keypoint) for perpendicular forward vector

Already reused in compute_polarization:

• compute_norm from movement.utils.vector for vector magnitude calculations
• compute_signed_angle_2d from movement.utils.vector for mean angle computation
• convert_to_unit from movement.utils.vector for unit vector normalization

---

Identity Tracking Requirements

Metric	Identity Tracking Needed?
Leadership	Required across frames (temporal cross-correlation)
Supporting metrics	Required across frames (pairwise relationships)
Polarization/Milling	Displacement mode: required (position across frames). Body-axis mode: not required (single-frame orientation).
Group Spread	Not required (instantaneous spatial measure)

Methods: SLEAP top-down-id for identity tracking, with EKS for pose refinement [5] or post-hoc Hungarian matching

---

API Design

# Proposed function signatures
from movement.kinematics import (
    # Core: State & Structure
    compute_polarization,
    compute_milling,
    compute_group_spread,
    compute_leadership,
    # Supporting: Pairwise
    compute_egocentric_angle,
    compute_radial_tangent_velocity,
)

# Core metrics - collective state

# Orientation polarization: body orientation alignment from keypoint pair
polarization = compute_polarization(
    ds.position,
    body_axis_keypoints=("tail_base", "neck"),
)

# Orientation polarization with mean body orientation angle (radians)
polarization, mean_angle = compute_polarization(
    ds.position,
    body_axis_keypoints=("tail_base", "neck"),
    return_angle=True,
)

# Orientation polarization with AP axis validation
polarization = compute_polarization(
    ds.position,
    body_axis_keypoints=("tail_base", "neck"),
    validate_ap=True,  # run anterior-posterior validation
)
# Validation results stored in polarization.attrs["ap_validation_result"]

# AP validation with custom config overrides
polarization = compute_polarization(
    ds.position,
    body_axis_keypoints=("tail_base", "neck"),
    validate_ap=True,
    ap_validation_config={
        "lateral_var_weight": 0.5,  # reduce penalty for side-to-side motion (default: 1.0)
        "confidence_floor": 0.2,    # stricter confidence warning (default: 0.1)
    },
)

# Heading polarization: movement direction alignment from displacement
polarization = compute_polarization(
    ds.position.sel(keypoints="thorax"),  # pre-select keypoint, or first is used
    displacement_frames=2,  # compare position over 2-frame window
)

# Heading polarization: first keypoint (index 0) used if none selected
polarization = compute_polarization(ds.position)

# Heading polarization with mean movement direction angle (degrees)
polarization, mean_angle = compute_polarization(
    ds.position.sel(keypoints="thorax"),
    return_angle=True,
    in_degrees=True,
)

milling = compute_milling(
    ds.position,
    body_axis_keypoints=("tail_base", "neck"),  # uses displacement mode if None
)

spread = compute_group_spread(ds.position, method="radius_of_gyration")

# Core metrics - hierarchical structure
leadership = compute_leadership(
    ds.position,  # velocity computed internally via differentiation
    max_lag=30,  # frames to search for optimal lag τ
)

# Supporting metrics - pairwise coordination
ego_angle = compute_egocentric_angle(
    ds.position,
    body_axis_keypoints=("tail_base", "neck"),
)

radial_tangent = compute_radial_tangent_velocity(
    ds.position,  # velocity computed internally via differentiation
)

# Return dimensions:
# compute_polarization()               -> (time,)
# compute_polarization(return_angle=True) -> tuple: (time,), (time,)
# compute_milling()                    -> (time,)
# compute_group_spread()               -> (time,)
# compute_leadership()                 -> (individuals, individuals_other, 2) where [:,:,0]=correlation, [:,:,1]=optimal_lag
# compute_egocentric_angle()           -> (time, individuals, individuals_other)
# compute_radial_tangent_velocity()    -> (time, individuals, individuals_other, 2) where [:,:,:,0]=radial, [:,:,:,1]=tangent

# Computational complexity:
# Core state metrics: O(N) per frame, O(N*T) total
# Leadership: O(N²*T*L) where L=max_lag
# Supporting: O(N²*T) for all pairwise relationships

Future Work: 3D Polarization Support

PR #948 extends the vector.py utilities to support 3D coordinate transforms. A subsequent PR can extend compute_polarization() to support 3D data:

Planned changes:

• Replace return_angle/in_degrees with return_direction parameter:
  • False (default): return only polarization
  • "unit_vector": also return mean heading/orientation as unit vector
  • "angle_radians": also return angle(s) in radians
  • "angle_degrees": also return angle(s) in degrees
• Use cart2pol/cart2sph for angle conversion (replaces compute_signed_angle_2d)
• 2D data returns (polarization, azimuth)
• 3D data returns (polarization, azimuth, elevation)

References

1. Movement zebras example
   • Demonstrates manual polarization calculation for 44 zebras
   • This proposal packages that approach into reusable functions
2. Couzin et al. 2002 - "Collective memory and spatial sorting in animal groups" Journal of Theoretical Biology 218(1): 1-11
   • Foundational model defining collective states (swarm, torus, dynamic/highly parallel groups) via polarization and angular momentum metrics
3. Nagy et al. 2010 - "Hierarchical group dynamics in pigeon flocks" Nature 464, 890-893
   • Leader-follower dynamics and group alignment metrics
4. Ballerini et al. 2008 - "Interaction ruling animal collective behavior depends on topological rather than metric distance" PNAS 105(4): 1232-1237
   • Starling flocks; topological interaction (6-7 nearest neighbors); uses α-shape algorithm for border detection
5. Cheng et al. 2025 - "Asymmetric Social Representations in the Prefrontal Cortex for Cooperative Behavior" bioRxiv
   • Source for egocentric angle, approach/tangent velocity used for behavioral classification (Sharp, Track, Sync, Join); validates SLEAP + EKS workflow for identity tracking.

Related Issues

• Explore existing packages for relevant metrics, features  #109 - Explore existing packages for relevant metrics (mentions swaRm for collective behavior: https://github.com/swarm-lab/swaRm)
• Write a guide on spatial coordinate systems #374 - Spatial coordinate systems (attachment: Spatial.coordinate.systems.pdf)
• Support re-IDing data #792 - Support re-IDing data (related to multi-animal identity tracking)
• Zulip Meeting Minutes - 2026-03-06, meeting notes on collective behavior metrics and mosaic behavior classification, potential identity tracking difficulties, and Niko's mention of polarization metric.
• feat: add compute_polarization() to collective behavior metrics #875 - compute_polarization() implementation
• ROI statistics - Entry / Exits #879 - ROI statistics - Entry/Exit, co-occupancy in ROIs.
• feat(body_axis): add skeleton-agnostic AP inference #945 - Prior-free body-axis AP inference (validate_ap parameter)
• feat(vector): add 3D support for coordinate transforms and vector ops #948 - 3D support for coordinate transforms & vector operations (extends vector.py utilities used by compute_polarization)

Checklist for Implementation

Core Metrics

• compute_polarization() feat: add compute_polarization() to collective behavior metrics #875
• validate_ap + ap_validation_config parameters feat(body_axis): add skeleton-agnostic AP inference #945
• compute_milling()
• compute_group_spread() Adding function: compute_group_spread #909
• compute_leadership()

Supporting Metrics

• compute_egocentric_angle()
• compute_radial_tangent_velocity()

Infrastructure

• 3D support for vector.py utilities feat(vector): add 3D support for coordinate transforms and vector ops #948 (extends compute_norm, convert_to_unit, cart2pol, pol2cart; adds cart2sph, sph2cart)
• 3D support for compute_polarization() (future work, depends on feat(vector): add 3D support for coordinate transforms and vector ops #948)

Testing & Documentation

• Unit tests for compute_polarization() - 49 tests across 5 classes (see feat: add compute_polarization() to collective behavior metrics #875):

Test Class	# Tests	Coverage
TestComputePolarizationValidation	13	DataArray type, required dims, space coord labels, empty keypoints, keypoint validation, displacement_frames type/range
TestComputePolarizationBehavior	18	Aligned→1, opposite→0, perpendicular→0, partial alignment, single individual, NaN exclusion, stationary exclusion, zero-length body-axis exclusion, permutation invariance, translation/scaling/rotation invariance (displacement mode), body-axis mode invariance
TestHeadingSourceSelection	5	Body-axis vs displacement mode, first-keypoint fallback, explicit .sel() selection, z-coord ignored, first-frame validity with angle
TestDisplacementFrames	4	First-N-frames NaN, reference-frame NaN propagation, multi-frame smoothing, oversized displacement window
TestReturnAngle	9	Output naming/dims, cardinal directions, diagonal motion, cancellation→NaN, angle rotation under global rotation, wraparound near ±π, in_degrees conversion

• Unit tests for validate_ap - TestValidateAPConfig (see feat(body_axis): add skeleton-agnostic AP inference #945)
• Unit tests for 3D vector operations - 7 tests + 2 fixtures (see feat(vector): add 3D support for coordinate transforms and vector ops #948)
• Example scripts for multi-animal datasets:
  • compute_polarization_viz.py - Polarization dynamics visualization (see feat: add compute_polarization() to collective behavior metrics #875)
  • compute_polarization_AP_inference.py - AP validation demo (see feat(body_axis): add skeleton-agnostic AP inference #945)
  • best_individual_selection_metrics.py - Reference individual selection metrics approach (see feat(body_axis): add skeleton-agnostic AP inference #945)
  • grid_search_ap_thresholds.py - AP config grid search optimization (see feat(body_axis): add skeleton-agnostic AP inference #945)