Using timecorr

Using timecorr#

Using timecorr: A Comprehensive Guide#

The timecorr toolbox is designed to compute dynamic correlations and explore higher-order correlational structure in timeseries data. This tutorial will walk you through the core concepts and functionality step by step.

Key Concepts#

timecorr operates in two main steps:

  1. Calculate dynamic correlations: Compute moment-by-moment correlations using kernel-based weighting

  2. Dimensionality reduction: Project correlations back to original data space to prevent “dimension explosion”

By repeating these steps, you can approximate higher-order correlations (correlations between correlations) in a computationally tractable way.

Applications#

This approach is particularly useful for:

  • Neuroscience: Dynamic brain connectivity analysis

  • Finance: Time-varying market correlations

  • Climate Science: Environmental variable relationships

  • Social Sciences: Dynamic network analysis

Load in required libraries#

[ ]:

# Load Required Libraries

import timecorr as tc
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import warnings
warnings.simplefilter("ignore")
%matplotlib inline

# Set style for better plots
plt.style.use('default')
sns.set_palette("husl")

Simulate some data#

Data Simulation#

The simulate_data function generates synthetic timeseries data for testing and demonstration. Let’s explore different simulation options and understand the data structure.

Simulation Parameters#

  • ``S``: Number of subjects/datasets

  • ``T``: Number of timepoints

  • ``K``: Number of features (variables)

  • ``datagen``: Type of data generation (‘ramping’, ‘block’, ‘random’, ‘constant’)

  • ``set_random_seed``: For reproducible results

[ ]:

# Simulate single subject's timeseries (reduced size for tutorial)
sim_1 = tc.simulate_data(S=1, T=100, K=20, set_random_seed=100)

print(f'Single subject data:')
print(f'  Shape: {np.shape(sim_1)}')
print(f'  Type: {type(sim_1)}')
print(f'  Data format: (timepoints, features) = ({sim_1.shape[0]}, {sim_1.shape[1]})')

# Visualize the simulated data
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(sim_1[:, :5])  # Plot first 5 features
plt.title('Sample Features Over Time')
plt.xlabel('Timepoints')
plt.ylabel('Signal Value')
plt.legend([f'Feature {i+1}' for i in range(5)])

plt.subplot(1, 2, 2)
plt.imshow(sim_1.T, aspect='auto', cmap='RdBu_r')
plt.title('All Features Heatmap')
plt.xlabel('Timepoints')
plt.ylabel('Features')
plt.colorbar()
plt.tight_layout()
plt.show()

[ ]:

# Simulate multiple subjects' timeseries
sim_3 = tc.simulate_data(S=3, T=100, K=20, set_random_seed=100)

print(f'Multi-subject data:')
print(f'  Type: {type(sim_3)}')
print(f'  Number of subjects: {len(sim_3)}')
print(f'  Each subject shape: {np.shape(sim_3[0])}')
print(f'  First subject type: {type(sim_3[0])}')

# Visualize data from all subjects
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for i, subject_data in enumerate(sim_3):
    axes[i].imshow(subject_data.T, aspect='auto', cmap='RdBu_r')
    axes[i].set_title(f'Subject {i+1}')
    axes[i].set_xlabel('Timepoints')
    if i == 0:
        axes[i].set_ylabel('Features')

plt.tight_layout()
plt.show()

print(f'\\nData Structure Summary:')
print(f'- Single subject: Returns numpy array of shape (T, K)')
print(f'- Multiple subjects: Returns list of {len(sim_3)} arrays, each of shape (T, K)')

[ ]:

# specify kernel:
width = 5
gaussian = {'name': 'Gaussian', 'weights': tc.gaussian_weights, 'params': {'var': width}}

# calcuate the dynamic correlations use a gaussian kernel and width of 5 for 1 simulate subject
vec_corrs = tc.timecorr(sim_1, weights_function=gaussian['weights'], weights_params=gaussian['params'])

[ ]:

# Compare different kernel functions
print("Kernel Comparison:")
T = 100  # Number of timepoints

# Define different kernels - note: Mexican Hat uses 'sigma' not 'var'
kernels = {
    'Gaussian': {'weights': tc.gaussian_weights, 'params': {'var': 5}},
    'Laplace': {'weights': tc.laplace_weights, 'params': {'scale': 5}},
    'Mexican Hat': {'weights': tc.mexican_hat_weights, 'params': {'sigma': 5}}
}

# Visualize kernel shapes
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
axes = axes.ravel()

for i, (name, kernel) in enumerate(kernels.items()):
    # Kernel functions expect params as a dictionary argument, not keyword arguments
    weights = kernel['weights'](T, kernel['params'])
    axes[i].plot(weights)
    axes[i].set_title(f'{name} Kernel')
    axes[i].set_xlabel('Timepoints')
    axes[i].set_ylabel('Weight')
    axes[i].grid(True)

# Compute dynamic correlations with Gaussian kernel
width = 5
gaussian = {'name': 'Gaussian', 'weights': tc.gaussian_weights, 'params': {'var': width}}

print(f"\\nComputing dynamic correlations with {gaussian['name']} kernel (var={width})...")
vec_corrs = tc.timecorr(sim_1, weights_function=gaussian['weights'], weights_params=gaussian['params'])

print(f"Input data shape: {sim_1.shape}")
print(f"Vectorized correlations shape: {vec_corrs.shape}")
print(f"Explanation: {sim_1.shape[1]} features -> {vec_corrs.shape[1]} unique correlations")
print(f"Formula: K*(K+1)/2 = {sim_1.shape[1]}*({sim_1.shape[1]}+1)/2 = {sim_1.shape[1]*(sim_1.shape[1]+1)//2}")

# Show kernel weights in last subplot
axes[3].plot(gaussian['weights'](T, gaussian['params']))
axes[3].set_title('Gaussian Kernel (Used)')
axes[3].set_xlabel('Timepoints')
axes[3].set_ylabel('Weight')
axes[3].grid(True)

plt.tight_layout()
plt.show()

Understanding Vectorized vs. Matrix Format#

The timecorr function returns correlations in vectorized format - the upper triangle of correlation matrices is flattened into a vector. This is memory-efficient and avoids redundancy since correlation matrices are symmetric.

Format Conversion#

  • ``vec2mat()``: Convert vectorized correlations to full matrices

  • ``mat2vec()``: Convert matrices back to vectorized format

[ ]:

# Understanding Vectorized vs. Matrix Format

print("Format Conversion Examples:")

# Show vectorized format (what timecorr returns)
print(f"1. Vectorized correlations shape: {vec_corrs.shape}")
print(f"   Interpretation: {vec_corrs.shape[0]} timepoints × {vec_corrs.shape[1]} unique correlations")

# Convert to matrix format for visualization
matrix_corrs = tc.vec2mat(vec_corrs)
print(f"\\n2. Matrix format shape: {matrix_corrs.shape}")
print(f"   Interpretation: {matrix_corrs.shape[0]} × {matrix_corrs.shape[1]} correlation matrix for each of {matrix_corrs.shape[2]} timepoints")

# Verify the conversion
print(f"\\n3. Consistency check:")
print(f"   Original features: {sim_1.shape[1]}")
print(f"   Matrix dimensions: {matrix_corrs.shape[0]} × {matrix_corrs.shape[1]}")
print(f"   Matrix is square: {matrix_corrs.shape[0] == matrix_corrs.shape[1]}")

# Visualize correlation matrices at different timepoints
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
timepoints = [25, 50, 75]

for i, t in enumerate(timepoints):
    im = axes[i].imshow(matrix_corrs[:, :, t], cmap='RdBu_r', vmin=-1, vmax=1)
    axes[i].set_title(f'Correlation Matrix at t={t}')
    axes[i].set_xlabel('Features')
    if i == 0:
        axes[i].set_ylabel('Features')

plt.colorbar(im, ax=axes, fraction=0.046, pad=0.04)
plt.tight_layout()
plt.show()

# Show temporal evolution of specific correlations
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
# Plot evolution of a few feature pairs
feature_pairs = [(0, 1), (0, 5), (0, 10), (5, 10)]
for i, (f1, f2) in enumerate(feature_pairs):
    plt.plot(matrix_corrs[f1, f2, :], label=f'Features {f1+1}-{f2+1}')
plt.title('Temporal Evolution of Feature Correlations')
plt.xlabel('Timepoints')
plt.ylabel('Correlation')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
# Show distribution of correlations at a single timepoint
correlations_t50 = matrix_corrs[:, :, 50]
# Extract upper triangle (unique correlations)
upper_triangle = correlations_t50[np.triu_indices_from(correlations_t50, k=1)]
plt.hist(upper_triangle, bins=20, alpha=0.7, edgecolor='black')
plt.title('Distribution of Correlations at t=50')
plt.xlabel('Correlation Value')
plt.ylabel('Frequency')
plt.grid(True)

plt.tight_layout()
plt.show()

# Convert back to vectorized format to verify round-trip
vec_corrs_roundtrip = tc.mat2vec(matrix_corrs)
print(f"\\n4. Round-trip conversion:")
print(f"   Original vectorized shape: {vec_corrs.shape}")
print(f"   Round-trip vectorized shape: {vec_corrs_roundtrip.shape}")
print(f"   Are they equal? {np.allclose(vec_corrs, vec_corrs_roundtrip)}")

print(f"\\n5. Memory efficiency:")
print(f"   Vectorized format: {vec_corrs.nbytes:,} bytes")
print(f"   Matrix format: {matrix_corrs.nbytes:,} bytes")
print(f"   Matrix is {matrix_corrs.nbytes / vec_corrs.nbytes:.1f}x larger")

[ ]:

# calcuate the dynamic isfc correlations use a Laplace kernel
# and width of 10 for 3 simulated subjects, and take the element-wise average correlations across matrices.
width = 10
laplace = {'name': 'Laplace', 'weights': tc.laplace_weights, 'params': {'scale': width}}

dyna_corrs = tc.timecorr(sim_3, combine=tc.corrmean_combine,
                         weights_function=laplace['weights'], weights_params=laplace['params'])

Higher order correlations#

Ok, now that we’ve gone over how to calculate dynamic correlations, let’s walk through reducing the correlations back to the original size of the data using the rfun parameter. Again, you have several options. If you want more information, please checkout the API documentation.

The default for rfun is None, which we used for calculating the dynamic correlations, but in this example we’ll use PCA.

[ ]:

# Compare different multi-subject approaches
width = 10
laplace = {'name': 'Laplace', 'weights': tc.laplace_weights, 'params': {'scale': width}}

print("Multi-Subject Analysis Comparison:")
print(f"Using {laplace['name']} kernel with scale={width}")

# 1. Individual subject correlations (no combination)
individual_corrs = []
for i, subject in enumerate(sim_3):
    corr = tc.timecorr(subject, weights_function=laplace['weights'], weights_params=laplace['params'])
    individual_corrs.append(corr)
    print(f"  Subject {i+1} correlations shape: {corr.shape}")

# 2. ISFC with mean combination
print(f"\\n2. ISFC Analysis:")
dyna_corrs = tc.timecorr(sim_3,
                         combine=tc.corrmean_combine,
                         weights_function=laplace['weights'],
                         weights_params=laplace['params'])
print(f"  ISFC combined shape: {dyna_corrs.shape}")
print(f"  Interpretation: Shared correlation patterns across all {len(sim_3)} subjects")

# 3. Compare individual vs. group patterns
print(f"\\n3. Pattern Comparison:")
# Convert to matrices for visualization
individual_mats = [tc.vec2mat(corr) for corr in individual_corrs]
group_mat = tc.vec2mat(dyna_corrs)

# Visualize at a specific timepoint
t = 50
fig, axes = plt.subplots(1, 4, figsize=(16, 4))

# Individual subjects
for i in range(3):
    im = axes[i].imshow(individual_mats[i][:, :, t], cmap='RdBu_r', vmin=-1, vmax=1)
    axes[i].set_title(f'Subject {i+1} (t={t})')
    axes[i].set_xlabel('Features')
    if i == 0:
        axes[i].set_ylabel('Features')

# Group average (ISFC)
im = axes[3].imshow(group_mat[:, :, t], cmap='RdBu_r', vmin=-1, vmax=1)
axes[3].set_title(f'ISFC Group (t={t})')
axes[3].set_xlabel('Features')

plt.colorbar(im, ax=axes, fraction=0.046, pad=0.04)
plt.tight_layout()
plt.show()

# Show temporal dynamics comparison
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
for i in range(3):
    plt.plot(individual_mats[i][0, 5, :], alpha=0.7, linewidth=1, label=f'Subject {i+1}')
plt.plot(group_mat[0, 5, :], 'k-', linewidth=3, label='ISFC Group')
plt.title('Feature Pair (1,6) Correlation Over Time')
plt.xlabel('Timepoints')
plt.ylabel('Correlation')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
# Correlation between individual and group patterns
individual_avg = np.mean(individual_mats, axis=0)
correlation_with_group = np.corrcoef(individual_avg.flatten(), group_mat.flatten())[0, 1]
plt.scatter(individual_avg.flatten(), group_mat.flatten(), alpha=0.5)
plt.plot([-1, 1], [-1, 1], 'r--', linewidth=2)
plt.title(f'Individual vs. ISFC Correlations\\nr = {correlation_with_group:.3f}')
plt.xlabel('Average Individual Correlations')
plt.ylabel('ISFC Group Correlations')
plt.grid(True)

plt.tight_layout()
plt.show()

Higher-Order Correlations: Preventing Dimension Explosion#

A key innovation of timecorr is the ability to compute higher-order correlations - correlations between correlation patterns - in a computationally tractable way.

The Problem#

Without dimensionality reduction, computing correlations on correlation data leads to exponential growth:

  • Original data: K features → K(K+1)/2 correlations

  • 2nd order: K(K+1)/2 features → (K(K+1)/2)((K(K+1)/2)+1)/2 correlations

  • This quickly becomes computationally intractable!

The Solution: Dimensionality Reduction (rfun)#

The rfun parameter applies dimensionality reduction to project correlations back to the original feature space:

  • ``’PCA’``: Principal Component Analysis (most common)

  • ``’ICA’``: Independent Component Analysis

  • ``’FactorAnalysis’``: Factor analysis

  • Graph measures: 'pagerank_centrality', 'eigenvector_centrality'

  • Custom functions: You can define your own reduction method

Summary and Next Steps#

What We’ve Covered#

In this tutorial, you’ve learned:

  1. Data Simulation: How to generate synthetic timeseries data for testing

  2. Dynamic Correlations: Computing time-varying correlation patterns using kernel functions

  3. Kernel Functions: Different weighting approaches (Gaussian, Laplace, Mexican Hat)

  4. Data Formats: Converting between vectorized and matrix correlation formats

  5. Multi-Subject Analysis: ISFC approaches for finding shared patterns

  6. Higher-Order Correlations: Using dimensionality reduction to explore correlations of correlations

  7. Multiple Levels: Building hierarchical correlation structures

Key Parameters Summary#

  • ``weights_function``: Kernel type (tc.gaussian_weights, tc.laplace_weights, etc.)

  • ``weights_params``: Kernel parameters ({‘var’: width} or {‘scale’: width})

  • ``cfun``: Correlation function (None, tc.isfc, tc.wisfc, tc.autofc)

  • ``rfun``: Dimensionality reduction (‘PCA’, ‘ICA’, ‘FactorAnalysis’, graph measures)

  • ``combine``: How to combine multi-subject results (tc.corrmean_combine, tc.mean_combine)

Real-World Applications#

The techniques demonstrated here apply to:

  • Neuroscience: fMRI dynamic functional connectivity, EEG/MEG source connectivity

  • Finance: Time-varying market correlations, risk assessment

  • Climate: Environmental variable relationships over seasons/years

  • Social Networks: Dynamic community structure, behavioral synchrony

Advanced Topics#

For more advanced applications, explore:

  • Statistical Testing: Permutation tests for significance

  • Cross-Validation: Using decoder functions for validation

  • Custom Kernels: Defining your own weighting functions

  • Graph Theory: Network measures and centrality metrics

  • Optimization: Memory and performance considerations for large datasets

Next Steps#

  1. Try the Applications Tutorial for domain-specific examples

  2. Explore the API Documentation for complete function references

  3. Experiment with your own data using the patterns shown here

  4. Check out the research paper for theoretical background

Happy correlating! 🎯

[ ]:

# Demonstrate dimensionality reduction and higher-order correlations
print("Higher-Order Correlation Analysis:")

# Step 1: Compute correlations with PCA reduction
print("\\n1. Dynamic correlations with PCA reduction:")
dyna_corrs_reduced = tc.timecorr(sim_3,
                                 rfun='PCA',
                                 weights_function=laplace['weights'],
                                 weights_params=laplace['params'])

print(f"  Original data shape: {np.shape(sim_3)} (list of {len(sim_3)} arrays)")
print(f"  Reduced correlations shape: {np.shape(dyna_corrs_reduced)}")
print(f"  Key insight: Same shape as original, but now represents correlation patterns!")

# Step 2: Compare different reduction methods
print("\\n2. Different reduction methods:")
reduction_methods = ['PCA', 'ICA']  # Removed FactorAnalysis to avoid potential issues
reduced_results = {}

for method in reduction_methods:
    try:
        result = tc.timecorr(sim_3,
                           rfun=method,
                           weights_function=laplace['weights'],
                           weights_params=laplace['params'])
        reduced_results[method] = result
        print(f"  {method}: {result.shape}")
    except Exception as e:
        print(f"  {method}: Failed ({str(e)})")

# Visualize the reduced correlation patterns
if reduced_results:
    fig, axes = plt.subplots(len(reduced_results), 3, figsize=(15, 4*len(reduced_results)))
    if len(reduced_results) == 1:
        axes = axes.reshape(1, -1)

    for i, (method, result) in enumerate(reduced_results.items()):
        for j in range(3):  # For each subject
            axes[i, j].imshow(result[j].T, aspect='auto', cmap='RdBu_r')
            axes[i, j].set_title(f'{method} - Subject {j+1}')
            axes[i, j].set_xlabel('Timepoints')
            if j == 0:
                axes[i, j].set_ylabel('Reduced Features')

    plt.tight_layout()
    plt.show()

# Step 3: Demonstrate the dimension problem without reduction
print("\\n3. Dimension explosion without reduction:")
K = sim_1.shape[1]  # Number of features
correlations_level1 = K * (K + 1) // 2
correlations_level2 = correlations_level1 * (correlations_level1 + 1) // 2

print(f"  Original features (K): {K}")
print(f"  1st order correlations: {correlations_level1}")
print(f"  2nd order correlations (without reduction): {correlations_level2:,}")
print(f"  Reduction ratio: {correlations_level2 / K:.1f}x larger!")
print(f"  \\n  With rfun='PCA': Always maintains {K} features per subject")
Contents