Extraction of effort-related features

Overview

In this script, we will extract effort-related features from the merged multimodal data we created for each trial in merging script.

Because we are dealing with time-varying data, we will extract number of statistics that characterize both instantaneous as well as cumulative nature of each (acoustic, motion, postural) timeseries in terms of effort.

Before we collect all relevant features, we normalize all time-varying features in the trial by minimum and maximum observed values for that feature for a participant. This is because we want to account for individual differences/strategies in using the whole range of motion and/or acoustic features. We therefore treat effort as a relative measure - but because we cannot access information about the minimum/maximum possible by a participant, we take the maximum and minimum that were observed within the whole experiment.

Note that this script collects more features than only the ones of current interest. The variables associated with the current study are selected in modeling script.

Code to prepare the environment 
import os
import glob
import pandas as pd
from scipy.signal import find_peaks
import numpy as np
import antropy as ent
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import seaborn as sns
import builtins
from pathlib import Path
from collections import defaultdict
from tdigest import TDigest
from math import e
from more_itertools import collapse
from scipy.stats import median_abs_deviation

curfolder = os.getcwd()

# Here we store our merged final timeseries data
mergedfolder = os.path.join(curfolder, "..", "03_TS_processing", "TS_merged")
filestotrack = glob.glob(os.path.join(mergedfolder, "merged_*.csv"))
filestotrack = [x for x in filestotrack if "merged_0" not in x] # Get rid of those from dyad 0 (pilot data)
filestotrack = [x for x in filestotrack if os.path.basename(x).split("_")[2] != "1"] # Get rid of those that have in the basename third element 1 (part 1, different experiment)
print(f"Number of files in total: {len(filestotrack)}")

# Here we store data related to exclusion
exclusionfolder = os.path.join(curfolder, "Cleaning")

# Here we store concept similarity data
conceptsimfolder = os.path.join(curfolder, "..", "04_ConceptSimilarity", "data")

# Here we store other potentially useful data
metadatafolder = os.path.join(curfolder, "..", "00_raw")

# Here we store the final data
datafolder = os.path.join(curfolder, "Datasets")
Number of files in total: 5238

Preparing additional metadata

To get additional metadata about the trials, we also need to load in and pre-process some more data, such as:

  • Similarity between answer and a target concept (computed with ConceptNet). This dataframe has information about the similarity between an answer and performed target concept, computed as cosine similarity of embeddings retrieved from ConceptNet.

  • Expressibility of a concept. This dataframe has information about expressibility of each concept in a stimuli list. The ratings were acquired from a separate group of people that participated in an online study (Ćwiek et al., 2026; Kadavá et al., 2024).

  • Response time. This dataframe has information about how long it took a guesser to provide an answer to a performance, measured from the end of the performance until pressing an Enter to confirm the answer.

This is the similarity dataframe

ID trial_order word answer exp English cosine_similarity sessionID
0 10 0 zwaaien zwaaien 2 to wave 1.000000 10_2
1 10 1 juichen dansen 2 to cheer 0.177888 10_2
2 10 2 juichen juichen 2 to cheer 1.000000 10_2
3 10 3 vis vis 2 fish 1.000000 10_2
4 10 4 hoog groot 2 high 0.448479 10_2
5 10 5 hoog rijken 2 high 0.233006 10_2
6 10 6 hoog hoog 2 high 1.000000 10_2
7 10 7 regen regen 2 rain 1.000000 10_2
8 10 8 verdrietig huilen 2 sad 0.389033 10_2
9 10 9 verdrietig verdrietig 2 sad 1.000000 10_2
10 10 10 langzaam slowmotion 2 slow 0.806700 10_2
11 10 11 langzaam joggen 2 slow 0.224116 10_2
12 10 12 langzaam langzaam 2 slow 1.000000 10_2
13 10 13 weten denken 2 to know 0.631044 10_2
14 10 14 weten idee 2 to know 0.291819 10_2


This is the expressibility dataframe

word modality fit
0 aanraken gesture 0.671
1 aanraken multimodal 0.695
2 aanraken vocal 0.207
3 aarde gesture 0.539
4 aarde multimodal 0.524
5 aarde vocal 0.195
6 ademen gesture 0.782
7 ademen multimodal 0.808
8 ademen vocal 0.785
9 alles gesture 0.351
10 alles multimodal 0.381
11 alles vocal 0.118
12 arm (lichaam) gesture 0.537
13 arm (lichaam) multimodal 0.575
14 arm (lichaam) vocal 0.127


This is the response time dataframe

TrialID response_time_sec prep_time_sec
3 10_2_3_p0 2.396 4.285
4 10_2_4_p0 1.965 3.144
5 10_2_5_p0 6.666 5.250
6 10_2_6_p0 2.039 4.844
7 10_2_7_p0 2.539 2.038
8 10_2_8_p0 2.860 2.142
9 10_2_9_p0 3.797 1.351
10 10_2_10_p0 4.073 2.930
11 10_2_11_p0 2.678 0.487
12 10_2_12_p0 5.783 6.776
13 10_2_13_p0 4.677 2.630
14 10_2_14_p0 1.540 3.616
15 10_2_15_p0 4.895 2.001
16 10_2_16_p0 2.307 3.057
21 10_2_21_p1 2.092 3.763

Excluding data

There are some trials that we need to exclude because of measurement errors or other issues. Some of them should be already excluded within processing pipeline, but we will double-check here.

Further, in previous scripts, we checked whether participants adhered to the conditions in the experiment, that is:

  • whether in gesture condition, they did not use any vocalizations
  • whether in vocal condition, they did not use any gestures/movements
  • whether in combined condition, they used both modalities

(see Condition adherence script)

  • and whether they did not use speech in any of the conditions

(see Speech check script)

Processing errors

First we exclude all errors related to (pre)processing

recording_errors = ['4_2', '21_2', '23_2', '25_1', '25_2', '35_1', '42_1', '42_2', '43_2', '44_2', '47_2', '71_1']
consent_errors = ['68_1', '68_2']
audio_errors = ['15_1', '16_2']
mt_errors = ['7_1', '7_2']

all_errors = recording_errors + consent_errors + audio_errors + mt_errors

# calculate all files in folderstotrack whose basename starts with amerged_+sessionID in all_errors
files_with_errors = [x for x in filestotrack if any(x.split(os.sep)[-1].startswith('merged_' + err) for err in all_errors)]

# print length and percentage of files with errors
print(f"Number of files with errors: {len(files_with_errors)} out of {len(filestotrack)} ({len(files_with_errors)/len(filestotrack)*100:.2f}%)")

# get rid of them
filestotrack = [x for x in filestotrack if x not in files_with_errors]
Number of files with errors: 78 out of 5238 (1.49%)

Now we exclude trials that had high error in calculating inverse kinematics (joint angles)

# Load in OpenSim error
IK_error = pd.read_csv(os.path.join(curfolder, "..", "02_MotionTracking_processing", "errors", "opensim_above_threshold.csv"))

# How many trials are affected from each part
IK_error['part'] = IK_error['pcnfolder'].str.split('_').str[1]
print(f"Number of trials affected from part 2: {IK_error[IK_error['part'] == '2'].shape[0]}")

ik_to_exclude = IK_error['pcnfolder'].tolist()
print(f"Which is {len(ik_to_exclude)/len(filestotrack)*100:.2f}% of all trials.")

files_with_ik_errors = [x for x in filestotrack if any(x.split(os.sep)[-1].startswith('merged_' + err) for err in ik_to_exclude)]
filestotrack = [x for x in filestotrack if x not in files_with_ik_errors]
Number of trials affected from part 2: 628
Which is 16.69% of all trials.

Exclusion of trials with condition violation

Now we exclude trials where modality rules were not met (e.g., voice in gesture condition) and where speech was used to explain concepts

def extract_trial_id_from_path(path):
    """Extract trial ID from file path by processing filename components.

    Args:
        path (str): File path containing trial information in filename

    Returns:
        str: Extracted trial ID in format 'X_X_X_X_X' or 'X_X_X_X_X_X' depending on presence of 'rein'
    """
    name = Path(path).stem          # merged_10_1_11_p1
    name = name.replace("merged_", "")
    parts = name.split("_")
    if "rein" in name:
        return "_".join(parts[:6])
    else:
        return "_".join(parts[:5])

# load in
condition = pd.read_csv(os.path.join(exclusionfolder, "condition_violation_annotation_filled.csv"))
speech = pd.read_csv(os.path.join(exclusionfolder, "control_for_speech_check_cleaned.csv"))

#### CONDITION

# in condition, keep only those that has CUT == 0
condition_clean = condition[condition['CUT'] == 0].copy()
print(f"Number of trials after condition clean: {condition_clean.shape[0]}")

condition_clean['TrialID'] = condition_clean['FILE'].apply(lambda x: '_'.join(x.split('_')[:5]) if 'rein' not in x else '_'.join(x.split('_')[:6]))
condition_clean['TrialID'] = condition_clean['TrialID'].str.replace('trial_', '', regex=False)
exclude_condition = condition_clean['TrialID'].tolist()

#### SPEECH

# in speech, get rid of all rows where there is 'no speech' in column check
speech_clean = speech[speech['check'] != 'no speech'].copy()
speech_clean = speech_clean[~speech_clean['check'].str.contains('cut|condition broke', na=False)].copy()
print(f"Number of trials after speech clean: {speech_clean.shape[0]}")

speech_clean['TrialID'] = speech_clean['filename'].apply(lambda x: '_'.join(x.split('_')[:4] + [x.split('_')[7]]) if 'rein' not in x else '_'.join(x.split('_')[:5] + [x.split('_')[8]]))
speech_clean['TrialID'] = speech_clean['TrialID'].str.replace('trial_', '', regex=False)
exclude_speech = speech_clean['TrialID'].tolist()

# find trials to exclude
trials_to_exclude = set(exclude_condition).union(set(exclude_speech))
print(f"Total number of trials to exclude: {len(trials_to_exclude)}")

# how many percent it is
print(f"Which is {len(trials_to_exclude)/len(filestotrack)*100:.2f}% of all trials.")

file_trial_map = {
    f: extract_trial_id_from_path(f)
    for f in filestotrack
}

trials_to_exclude = set(exclude_condition).union(set(exclude_speech))

filestotrack_clean = [
    f for f, trial_id in file_trial_map.items()
    if trial_id not in trials_to_exclude
]

excluded = set(file_trial_map.values()) & trials_to_exclude
missing  = trials_to_exclude - set(file_trial_map.values())

print(f"Excluded {len(excluded)} trials from filestotrack.")
print(f"Trials not found in filestotrack: {len(missing)}")
print("Total amount of files after final cleaning: ", len(filestotrack_clean))
Number of trials after condition clean: 165
Number of trials after speech clean: 52
Total number of trials to exclude: 206
Which is 4.54% of all trials.
Excluded 134 trials from filestotrack.
Trials not found in filestotrack: 72
Total amount of files after final cleaning:  4408

Extracting features: rationale

In the following part, we will be extracting summaries from the existing timeseries that we pre-processed in the previous scripts. The goal is to capture the effort-related features that might be relevant in a repair, i.e., when a participant has to correct themselves to regain understanding with their partner who previously did not understand the performed concept.

Because we are dealing with time-varying data, we will extract number of statistics that characterize both instantaneous as well as cumulative nature of each (acoustic, motion, postural) timeseries in terms of effort. These are:

  • global mean and standard deviation
  • number, mean and standard deviation of the peaks (using function find_peaks from scipy.signal). For acceleration and moment change, we also calculate the measurements relating to negative peaks.
  • range of the values
  • integral

Additionally, we will compute sample entropy to measure the complexity of a timeseries.

Furthermore, we will utilize measurements that capture characteristics beyond the statistics sketched above. These include:

  • intermittency (as dimensionless jerk); used in Pouw et al. (2021);
  • bounding box of movement volume (i.e., gesture space); used in Żywiczyński et al. (2024);
  • vowel space area (VSA); used in Berisha et al. (2014);
  • motor complexity; computed as the slope of PCA; similar to Yan et al. (2020);
  • number of submovements; used in Pouw et al. (2021); Trujillo et al. (2018);
  • number of moving articulators.

Overal mean+std per variable

Basic statistics capture properties of the timeseries in several dimensions: - globally (mean, standard deviation of the timeseries, range of the values, rate of the feature) - locally (number, mean and standard deviation of the peaks) - cummulatively (integral of the timeseries)

We also collect the time stamps of the peaks, as we will use them to compute the synchronization between the modalities.

First, to be able to apply the find_peaks function properly, we will need for each column a mean and SD of the timeseries across all trials.

if not os.path.exists(os.path.join(datafolder, "general_means_stds_nonzscored.csv")):

    # Running stats for z-scores
    running_sum = defaultdict(float)
    running_sum_sq = defaultdict(float)
    running_count = defaultdict(int)

    for file in filestotrack_clean:   
        print(f"Processing file: {file}")
        df = pd.read_csv(file)

        # Get all numerical cols (except time)
        numcols = df.select_dtypes(include=[np.number]).columns
        numcols = [x for x in numcols if x != 'time']  

        for col in numcols:
            #print(f"  Processing column: {col}")
            unnormed_col = df[col]

            # Update online stats
            valid = unnormed_col.notna()
            running_sum[col] += unnormed_col[valid].sum()
            running_sum_sq[col] += (unnormed_col[valid] ** 2).sum()
            running_count[col] += valid.sum()

    # Create the general mean and std
    general_means = {col + '_general_mean': running_sum[col] / running_count[col] for col in running_sum}
    general_stds = {
        col + '_general_std': np.sqrt((running_sum_sq[col] / running_count[col]) - (running_sum[col] / running_count[col]) ** 2)
        for col in running_sum
    }

    # Combine into one dataframe
    df_means_nonz = pd.DataFrame([{**general_means, **general_stds}])

    # save it
    df_means_nonz.to_csv(os.path.join(datafolder, 'general_means_stds_nonzscored.csv'), index=False)

else:
    # load it
    df_means_nonz = pd.read_csv(os.path.join(datafolder, 'general_means_stds_nonzscored.csv'))

df_means_nonz.head(20)
left_back_general_mean right_forward_general_mean right_back_general_mean left_forward_general_mean COPXc_general_mean COPYc_general_mean COPc_general_mean audio_general_mean envelope_general_mean envelope_norm_general_mean ... f1_clean_general_std f2_clean_general_std f3_clean_general_std f1_clean_vel_general_std f2_clean_vel_general_std f3_clean_vel_general_std lowerbody_power_general_std leg_power_general_std head_power_general_std arm_power_general_std
0 1.18748 0.973936 1.547235 1.390193 0.000002 0.000004 0.584969 -0.000002 0.002621 0.039264 ... 414.496998 525.487178 474.195414 37654.390067 50577.277891 53458.115579 64.481032 1.883556 3.323033 21.473082

1 rows × 1052 columns

Basic/descriptive statistics

Now we can collect basic statistics.

Note that when we look for peaks, we take the mean value for that signal over all trials. Additionally, in timeseries that have both positive and negative values (e.g., acceleration), we will look for both positive and negative peaks.

def get_statistics(cols, df, dictionary, means_df, plot=False):
    """Calculate statistics for specified columns in a dataframe.

    Args:
        cols: List of column names to process
        df: Input dataframe containing the signals
        dictionary: Dictionary to store results
        means_df: Dataframe containing mean and standard deviation values
        plot: Boolean flag to enable plotting

    Returns:
        Dictionary containing calculated statistics for each column
    """

    fs = 500  # Sampling frequency

    for col in cols:

        # exclude edge artefacts
        edge_beg = int(0.15 * fs)  # 150 ms
        edge_end = int(0.2 * fs)   # 200 ms
        full_signal = df[col].values
        trimmed_signal = df[col].values[edge_beg:-edge_end]

        # access sd from means_df
        col_sd = means_df[col + '_general_std'].iloc[0]

        scale = col_sd # global scale, taking into consideration variability across all trials
        
        if 'acc' in col or 'Acc' in col or 'jerk' in col:
            distance = int(0.08 * fs)   # 80 ms
            prominence = 0.4 * scale
        elif 'moment' in col:
            distance = int(0.20 * fs)   # 200 ms
            prominence = 0.4 * scale
        elif 'speed' in col or 'power' in col:
            distance = int(0.20 * fs)   # 200 ms
            prominence = 0.4 * scale
        elif 'COPc' in col:
            distance = int(0.25 * fs)   # 200 ms
            prominence = 0.4 * scale
        elif 'f0' in col or 'f1' in col or 'f2' in col or 'f3' in col:
            distance = int(0.20 * fs)   # 200 ms
            prominence = 0.1 * scale
        elif 'envelope' in col:
            distance = int(0.15 * fs)   # 150 ms
            prominence = 0.4 * scale
        else:
            distance = int(0.10 * fs)
            prominence = 0.5 * scale

        # for 2nd derivative of position data, joint kinematics, or moments, we calculate peaks on absolute effort
        if 'acc' in col or 'Acc' in col or 'moment_sum_change' in col:
            signal_for_peaks = np.abs(full_signal - np.median(full_signal))
        else:
            signal_for_peaks = full_signal    

        min_height = col_sd * 0.3
        peaks_full, _ = find_peaks(signal_for_peaks, distance=distance, prominence=prominence, height=min_height)
        valid_mask = (peaks_full >= edge_beg) & (peaks_full < len(full_signal) - edge_end)
        peaks_valid = peaks_full[valid_mask]
        peaks_values = signal_for_peaks[peaks_valid]
        peaks_times = df['time'].values[peaks_valid].tolist()
        peaks_sign = np.sign(df[col].values[peaks_valid]).tolist() # negative or positive peaks
        if len(peaks_values) > 0:
            peak_mean = peaks_values.mean()
            peak_std = peaks_values.std()
            peak_n = len(peaks_values)
        else:
            # Get global maximum within the same valid region used for peak detection
            valid_signal = signal_for_peaks[edge_beg : len(full_signal) - edge_end]
            peak_mean = valid_signal.max() if len(valid_signal) > 0 else np.nan
            peak_std = np.nan
            peak_n = 0

        # Plot if needed
        if plot:
            plt.plot(df['time'], signal_for_peaks, label="Signal")
            plt.scatter(peaks_times, peaks_values, color='r', label="Pos Peaks", marker="x")
            plt.axvline(df['time'].iloc[edge_beg], color='k', linestyle='--', label='Edge Artefact Start')
            plt.axvline(df['time'].iloc[-edge_end-1], color='k', linestyle='--', label='Edge Artefact End')
            plt.title(col)
            plt.legend()
            plt.legend(fontsize='small')
            plt.show()

        # Calculate remaining features on trimmed signal + excluding NAN values
        signal = trimmed_signal[~np.isnan(trimmed_signal)]
        Gmean = np.mean(signal)   # General mean
        Gstd = np.std(signal)     # General std
        integral = np.trapz(signal, dx=2)  # Integral
        duration = len(signal) * 2
        integral_norm = integral / duration if duration > 0 else np.nan  # Normalized integral
        max_val = np.nanmax(trimmed_signal)
        min_val = np.nanmin(trimmed_signal)
        range_val = max_val - min_val       # Range                  
        
        # Save all in dictionary
        dictionary[col] = [Gmean, Gstd, peak_mean, peak_std, peak_n, peaks_times, peaks_sign, integral, integral_norm, range_val, max_val, min_val]

    return dictionary

# Function to adapt row
def adapt_row(row_to_process):
    """Process a row of data by extracting and organizing statistical features from signal data.

    This function takes a row containing signal data and statistical measures, extracts the relevant
    features, and organizes them into a standardized format. It handles both peak-related and general
    statistical measures, and filters out columns that don't meet certain criteria.

    Args:
        row_to_process: A pandas Series or DataFrame row containing signal data and pre-calculated measures

    Returns:
        A processed row containing only the extracted statistical features, with NaN columns removed
    """
    for col in row_to_process.columns:
        # Calculate for all expcet some already calculated measures
        if 'inter' not in col and 'sampEn' not in col and 'bbmv' not in col and 'duration' not in col:
            row_to_process[col + '_Gmean'] = row_to_process[col].apply(lambda x: x[0])
            row_to_process[col + '_Gstd'] = row_to_process[col].apply(lambda x: x[1])
            row_to_process[col + '_peak_mean'] = row_to_process[col].apply(lambda x: x[2])
            row_to_process[col + '_peak_std'] = row_to_process[col].apply(lambda x: x[3])
            row_to_process[col + '_peak_n'] = row_to_process[col].apply(lambda x: x[4])
            row_to_process[col + '_peak_times'] = row_to_process[col].apply(lambda x: x[5])
            row_to_process[col + '_peak_signs'] = row_to_process[col].apply(lambda x: x[6])
            row_to_process[col + '_integral'] = row_to_process[col].apply(lambda x: x[7])
            row_to_process[col + '_integral_norm'] = row_to_process[col].apply(lambda x: x[8])
            row_to_process[col + '_range'] = row_to_process[col].apply(lambda x: x[9])
            row_to_process[col + '_max'] = row_to_process[col].apply(lambda x: x[10])
            row_to_process[col + '_min'] = row_to_process[col].apply(lambda x: x[11])

    # Now keep only this newly created cols
    row_final = row_to_process[[col for col in row_to_process.columns if any(x in col for x in ['Gmean', 'Gstd', 'peak_mean', 'peak_std', 'peak_n', 'sampen', 'inter', 'integral', 'integral_norm', 'peak_times', 'peak_signs', 'bbmv', 'range', 'duration', 'negpeak_mean', 'negpeak_std', 'negpeak_n', 'negpeak_times', 'max', 'min']) or col in ['TrialID', 'condition']]]

    # Get rid of cols with NaNs
    row_final = row_final.dropna(axis=1, how='all')
    
    return row_final

###### USAGE

sample = filestotrack_clean[200]
print(f"Processing sample file: {sample}")
df_sample = pd.read_csv(sample)
print(f"FileInfo: {df_sample['FileInfo'].iloc[0]}")

sample_dict = {}
sample_col = ['arm_moment_sum_change']
#sample_col = ['envelope_norm']
#sample_col = ['COPc'] 
#sample_col = ['f0']
#sample_col = ['RWrist_speed'] #'arm_power', 'RWrist_speed'
#sample_col = ['RWrist_acc']

if df_sample.shape[0] > 0:
    sample_dict = get_statistics(sample_col, df_sample, sample_dict, df_means_nonz, plot=True)

# Convert dictionary to dataframe
sample_features = pd.DataFrame({key: [value] for key, value in sample_dict.items()})
sample_features = adapt_row(sample_features)

# Here we see the results
sample_features.head()
Processing sample file: f:\FLESH_ContinuousBodilyEffort\05_TS_featureExtraction\..\03_TS_processing\TS_merged\merged_12_2_69_p1.csv
FileInfo: p1_hoorn_geluiden_c1

arm_moment_sum_change_Gmean arm_moment_sum_change_Gstd arm_moment_sum_change_peak_mean arm_moment_sum_change_peak_n arm_moment_sum_change_peak_times arm_moment_sum_change_peak_signs arm_moment_sum_change_integral arm_moment_sum_change_integral_norm arm_moment_sum_change_range arm_moment_sum_change_max arm_moment_sum_change_min
0 2.110741 1.135988 2.087184 0 [] [] 10976.181168 2.109993 4.037998 4.173805 0.135807

Bounding box of movement volume

def get_bbmv(df, group, kp_dict):
    """
    Calculate the bounding box movement volume (BBMV) for a given group of keypoints.

    Parameters:
    -----------
    df : pandas.DataFrame
        DataFrame containing coordinate data
    group : str
        Name of the group to calculate BBMV for
    kp_dict : dict
        Dictionary mapping group names to keypoint prefixes

    Returns:
    --------
    float
        Natural logarithm of the sum of BBMVs for all keypoints in the group
    """
    coordinates = [col for col in df.columns if any(x in col for x in ['_x', '_y', '_z'])]

    # Prepare columns that belong to a group (e.g., arm)
    kp = kp_dict[group]
    colstoBBMV = [col for col in coordinates if any(x in col for x in kp)]

    # Keep only unique names without coordinates
    kincols = list(set([col.split('_')[0] for col in colstoBBMV]))

    bbmvs = {}
    for col in kincols:
        # Span of x, y, z
        x_span = df[col + '_x'].max() - df[col + '_x'].min()
        y_span = df[col + '_y'].max() - df[col + '_y'].min()
        z_span = df[col + '_z'].max() - df[col + '_z'].min()

        # Calculate BBMV
        bbmv = x_span * y_span * z_span
        bbmvs[col] = bbmv
        
    # Get the sum for the whole group
    bbmv_sum = sum(bbmvs.values())

    # Natural logarithm
    bbmv_sum = np.log(bbmv_sum) 

    return bbmv_sum

#### USAGE

# These are groups for BBMV
kp_arms = ['RWrist_z', 'RElbow_z', 'RShoulder_z', 'LWrist_z', 'LElbow_z', 'LShoulder_z']
kp_lower = ['RAnkle_z', 'RKnee_z', 'LAnkle_z', 'LKnee_z']
kp_legs = ['RAnkle_z', 'RKnee_z', 'LAnkle_z', 'LKnee_z']
kp_head = ['Head_z']
kp_keys = {'arm': kp_arms, 'lowerbody': kp_lower, 'leg': kp_legs, 'head': kp_head}
subdf = df_sample

# Calculate BBMV for sample
bbmv = get_bbmv(subdf, 'arm', kp_keys)
print('BBMV of the current trial: ', bbmv)
BBMV of the current trial:  3.1627629746370705

Symmetry of arm movement / multi-movement asymmetry of submovements

Symmetry of arm movement is a measure of how much the participant uses both arms. Symmetrical movement is expected to be less effortful than asymmetrical movement, as asymmetrical movement might require more cognitive resources to coordinate the two arms. We will compute symmetry as the correlation between left and right arm trajectories, similar to Xiong et al. (2002). If it’s close to 1, it means that the participant uses both arms in a symmetrical way.

def compute_ArmCoupling(df, keypoints=['Wrist'], dimensions=['x', 'y', 'z'], absolute=True, z=False):
    """
    Calculate the average correlation between left and right arm movements across specified keypoints and dimensions.

    Parameters:
    - df: DataFrame containing the movement data
    - keypoints: List of body keypoints to analyze (default: ['Wrist'])
    - dimensions: List of spatial dimensions to consider (default: ['x', 'y', 'z'])
    - absolute: If True, returns absolute correlation values (default: True)
    - z: If True, uses z-scored columns (default: False)

    Returns:
    - float: Average correlation between left and right movements
    """
    correlations = []
    
    for kp in keypoints:
        for dim in dimensions:
            if z:
                left_col = f"L{kp}_{dim}_z"
                right_col = f"R{kp}_{dim}_z"
            else:
                left_col = f"L{kp}_{dim}"
                right_col = f"R{kp}_{dim}"
                
            if left_col in df.columns and right_col in df.columns:
                corr = np.corrcoef(df[left_col], df[right_col])[0, 1]  # Pearson correlation
                if absolute:
                    corr = np.abs(corr)
                correlations.append(corr)
    
    return np.nanmean(correlations)  # Average correlation across keypoints and dimensions

def compute_multimovementAssymetry(df, z=False):
    """
    Calculate the asymmetry in movement between left and right arms across multiple joints.

    Parameters:
    - df: DataFrame containing the movement data
    - z: If True, uses z-scored columns (default: False)

    Returns:
    - dict: Dictionary containing:
        - mean_integral_difference: Mean difference in integrated velocities across joints
        - joint_integral_differences: List of differences for each joint
        - total_left_integral_velocity: Total integrated velocity for left joints
        - total_right_integral_velocity: Total integrated velocity for right joints
    """
    joints = ['elbow_flex', 'wrist_flex', 'wrist_dev']
    
    integral_diff_left = 0
    integral_diff_right = 0
    joint_integral_differences = []

    for joint in joints:
        if z:
            joint_r_speed = f"{joint}_r_speed_z"
            joint_l_speed = f"{joint}_l_speed_z"
        else:
            joint_r_speed = f"{joint}_r_speed"
            joint_l_speed = f"{joint}_l_speed"
            
        if joint_r_speed in df.columns and joint_l_speed in df.columns:
            # Compute the integral of velocity (sum of speed over time) for left and right joints
            integral_velocity_left = np.sum(np.abs(df[joint_l_speed]))  # Sum of absolute speed
            integral_velocity_right = np.sum(np.abs(df[joint_r_speed]))  # Sum of absolute speed

            # Compute the difference in integral of angular velocity between left and right 
            # the more +, the more distance is left is travelling), 
            # but note that because od z-scoring, 0 does not reflect perfect symmetry in distance travelled by both hands
            joint_integral_diff = integral_velocity_left - integral_velocity_right
            joint_integral_differences.append(joint_integral_diff)

            # Accumulate for total integral difference calculation
            integral_diff_left += integral_velocity_left
            integral_diff_right += integral_velocity_right

    # Compute the overall difference in integrated velocities
    total_integral_difference = np.mean(joint_integral_differences)  # Mean of differences across joints

    return {
        'mean_integral_difference': total_integral_difference,
        'joint_integral_differences': joint_integral_differences,
        'total_left_integral_velocity': integral_diff_left,
        'total_right_integral_velocity': integral_diff_right
    }

### USAGE

coupling_score = compute_ArmCoupling(df_sample, keypoints=['Wrist', 'Elbow'], dimensions=['x', 'y', 'z'], z=False)
arm_asymmetry = compute_multimovementAssymetry(df_sample, z=False)

# Plot the trajectories
plt.figure(figsize=(10, 8))
ax = plt.axes(projection='3d')

# Plot the left arm
ax.plot3D(df_sample['LWrist_x'], df_sample['LWrist_y'], df_sample['LWrist_z'], 'r')
ax.plot3D(df_sample['LElbow_x'], df_sample['LElbow_y'], df_sample['LElbow_z'], 'r', linestyle='--')

# Plot the right arm
ax.plot3D(df_sample['RWrist_x'], df_sample['RWrist_y'], df_sample['RWrist_z'], 'b')
ax.plot3D(df_sample['RElbow_x'], df_sample['RElbow_y'], df_sample['RElbow_z'], 'b', linestyle='--')

plt.show()

# Print the symmetry score
print("Coupling Score:", coupling_score)
print("Arm Asymmetry:", arm_asymmetry['mean_integral_difference'])

Coupling Score: 0.5357169102055633
Arm Asymmetry: -45.318783240188054

Intermittency

Intermittency characterizes the ‘unsmoothess’ of the movement. We follow Pouw et al. (2021)’s adaptation of Hogan & Sternad (2009) and compute dimensionless jerk measure. This measure is computed as the integral of the second derivative of the speed (i.e., jerk) squared multiplied by the duration cubed over the maximum squared velocity.

def get_intermittency(jerk_values, speed_values):
    """
    Calculate intermittency based on jerk and speed values.

    Args:
        jerk_values (np.ndarray): Array of jerk values
        speed_values (np.ndarray): Array of speed values

    Returns:
        float: Calculated intermittency value or np.nan if all speeds are zero
    """
    smoothed_jerk = jerk_values
    speed = speed_values
    
    if not np.all(speed == 0):
        integrated_squared_jerk = np.sum(smoothed_jerk ** 2)
        max_squared_speed = np.max(speed ** 2)
        D3 = len(speed) ** 3
        jerk_dimensionless = integrated_squared_jerk * (D3 / max_squared_speed)
        smoothness = jerk_dimensionless
    else:
        smoothness = np.nan

    return smoothness

### USAGE

# Calculate on sample
subdf = df_sample
group = 'arm'

# Get the columns that contain group and sum/power (these are cols with aggregated measures)
cols = [col for col in subdf.columns if group in col and 'sum' in col or 'power' in col]

# Calculate intermittency for kinematic measures (from Pose2sim)
speed = subdf[group + "_speedKin_sum"].values
jerk = subdf[group + "_jerkKin_sum"].values
intermittency_kin = get_intermittency(jerk, speed)
intermittency_kin = np.log(intermittency_kin) # natural log
print('Intermittency, kinematics: ', intermittency_kin)

# Calculate intermittency for inverse kinematic measures (from OpenSim)
speed = subdf[group + "_angSpeed_sum"].values
jerk = subdf[group + "_angJerk_sum"].values
intermittency_ik = get_intermittency(jerk, speed)
intermittency_ik = np.log(intermittency_ik) # natural log
print('Intermittency, joint angles: ', intermittency_ik)
Intermittency, kinematics:  21.50655846100739
Intermittency, joint angles:  17.559605779397916

Vowel space area

Vowel space area (VSA) is defined as the are of the quadrilateral formed by the four corner vowels /i/, /a/, /u/ and /ɑ/ in the F1-F2 space. It is a measure of the articulatory ‘working’ space and researchers often use it to characterize speech motor control (Berisha et al., 2014).

Despite the fact that we are not dealing with speech data, we can still compute VSA for the first two formants. Following Daniel R. McCloy’s tutorial for R, we will compute VSA for the sample timeseries as the area of the convex hull encompassing all tokens.

def getVSA(f1, f2, plot=False):
    """
    Calculate the Vowel Space Area (VSA) from formant frequencies.

    Parameters:
    f1 (array-like): First formant frequencies in Hz
    f2 (array-like): Second formant frequencies in Hz
    plot (bool): If True, generates a plot of the convex hull

    Returns:
    float: Natural logarithm of the 2D convex hull volume
    """
   
    # 2d convex hull
    points = np.array([f1, f2]).T
    hull_2d = ConvexHull(points)
    volume_2d = hull_2d.volume
    volume_2d = np.log(volume_2d) # natural log

    if plot:
        plt.figure()
        plt.plot(points[:, 0], points[:, 1], 'o', markersize=3, label='Formant Points')
        
        for simplex in hull_2d.simplices:
            plt.plot(points[simplex, 0], points[simplex, 1], 'r-', linewidth=1)
        
        plt.xlabel('F1 (Hz)')
        plt.ylabel('F2 (Hz)')
        plt.gca().invert_yaxis()
        plt.legend()
        plt.title('2D Convex Hull of Vowel Space Area')
        plt.show()
    
    return volume_2d

### USAGE

f1_clean = df_sample['f1_clean'].dropna()
f2_clean = df_sample['f2_clean'].dropna()
vsa = getVSA(f1_clean, f2_clean, plot=True)

Motor complexity

As a part of exploratory analysis, we also want to assess the motor complexity of the movements as a proxy of coordinative - therefore cognitive - effort relating to the signal.

Using OpenSim (Seth et al., 2018), we have been able to extract joint angle measurements for our data (see here). We can use them to assess a motor complexity or alternatively, motor performance of a participant’s movements. We will use Principal Component Analysis (PCA) to assess the motor complexity. We will look at the number of principal components that explain a certain amount of variance in the data and the slope of the explained variance. This will give us a measure of how many (uncorrelated) dimensions/modes cover the most (80 or 95%) of the features of a signal, therefore how complex the movement is in terms of its coordination patterns (Daffertshofer et al., 2004; Yan et al., 2020).

def get_PCA(df, plot=False):
    """
    Perform Principal Component Analysis (PCA) on input DataFrame.

    Parameters:
    -----------
    df : pandas.DataFrame
        Input data for PCA analysis
    plot : bool, optional
        Whether to generate and display the explained variance plot (default: False)

    Returns:
    --------
    tuple
        Contains:
        - n_components_for_80_variance: Number of components explaining 80% variance
        - slope_80: Slope of cumulative variance up to 80% threshold
        - n_components_for_95_variance: Number of components explaining 95% variance
        - slope_95: Slope of cumulative variance up to 95% threshold
    """
    # Step 1: Standardize the Data
    scaler = StandardScaler()
    standardized_data = scaler.fit_transform(df)
    
    # Step 2: Apply PCA
    pca = PCA()
    pca.fit(standardized_data)
    
    # Step 3: Explained Variance
    explained_variance = pca.explained_variance_ratio_
    cumulative_explained_variance = explained_variance.cumsum()

    # Step 4: Find the number of components that explain 95% variance
    n_components_for_80_variance = np.argmax(cumulative_explained_variance >= 0.8) + 1
    
    # Step 5: Compute the slope for the first n_components_for_95_variance
    if n_components_for_80_variance > 1:
        slope_80 = (cumulative_explained_variance[n_components_for_80_variance-1] - cumulative_explained_variance[0]) / (n_components_for_80_variance - 1)
    else:
        slope_80 = cumulative_explained_variance[0]  # Only one component case

    # 95% variance
    n_components_for_95_variance = np.argmax(cumulative_explained_variance >= 0.95) + 1
    
    # Step 5: Compute the slope for the first n_components_for_95_variance
    if n_components_for_95_variance > 1:
        slope_95 = (cumulative_explained_variance[n_components_for_95_variance-1] - cumulative_explained_variance[0]) / (n_components_for_95_variance - 1)
    else:
        slope_95 = cumulative_explained_variance[0]  # Only one component case

    sns.set_style("whitegrid")  

    main_color = "black"   # Black for main cumulative variance line
    red_color = "red"      # Basic red for 80% variance threshold
    blue_color = "blue"    # Basic blue for 95% variance threshold

    if plot:
        plt.grid(False)
        plt.gca().spines['top'].set_visible(False)
        
        # Adjust x-axis indexing: Components start from 1, not 0
        x_values = range(1, len(cumulative_explained_variance) + 1)
        
        # Plot cumulative explained variance (black line with markers)
        plt.plot(x_values, cumulative_explained_variance, marker='o', linestyle='-', markersize=4, linewidth=2, color=main_color, alpha=0.9, label='Cumulative Explained Variance')
        
        # Add vertical dashed lines for 80% and 95% variance thresholds
        plt.axvline(x=n_components_for_80_variance, color=red_color, linestyle='--', 
                    linewidth=2, alpha=0.9)
        
        plt.axvline(x=n_components_for_95_variance, color=blue_color, linestyle='--', 
                    linewidth=2, alpha=0.9)

        plt.text(n_components_for_80_variance + 0.3, 0.8, "80% Variance", 
                color=red_color, fontsize=12, fontweight='bold')
        
        plt.text(n_components_for_95_variance + 0.3, 0.95, "95% Variance", 
                color=blue_color, fontsize=12, fontweight='bold')
        plt.xlabel('Number of Principal Components', fontsize=16, color='black')
        plt.ylabel('Cumulative Explained Variance', fontsize=16, color='black')
        plt.xticks(x_values)  # Explicitly setting x-ticks
        
        plt.grid(True, alpha=0.3)  
        plt.gcf().set_size_inches(6, 4)
        plt.tight_layout()
        plt.xlim(0, len(cumulative_explained_variance) +1)  # Adjust x-axis limits
        plt.show()

    return n_components_for_80_variance, slope_80, n_components_for_95_variance, slope_95

# Calculate on sample
PCAcolls_all = ['pelvis_tilt', 'pelvis_list', 'pelvis_rotation', 'pelvis_tx',
       'pelvis_ty', 'pelvis_tz', 'hip_flexion_r', 'hip_adduction_r',
       'hip_rotation_r', 'knee_angle_r', 'knee_angle_r_beta', 'ankle_angle_r',
       'subtalar_angle_r', 'hip_flexion_l', 'hip_adduction_l',
       'hip_rotation_l', 'knee_angle_l', 'knee_angle_l_beta', 'ankle_angle_l',
       'subtalar_angle_l', 'L5_S1_Flex_Ext', 'L5_S1_Lat_Bending',
       'L5_S1_axial_rotation', 'L4_L5_Flex_Ext', 'L4_L5_Lat_Bending',
       'L4_L5_axial_rotation', 'L3_L4_Flex_Ext', 'L3_L4_Lat_Bending',
       'L3_L4_axial_rotation', 'L2_L3_Flex_Ext', 'L2_L3_Lat_Bending',
       'L2_L3_axial_rotation', 'L1_L2_Flex_Ext', 'L1_L2_Lat_Bending',
       'L1_L2_axial_rotation', 'L1_T12_Flex_Ext', 'L1_T12_Lat_Bending',
       'L1_T12_axial_rotation', 'neck_flexion', 'neck_bending',
       'neck_rotation', 'arm_flex_r', 'arm_add_r', 'arm_rot_r', 'elbow_flex_r',
       'pro_sup_r', 'wrist_flex_r', 'wrist_dev_r', 'arm_flex_l', 'arm_add_l',
       'arm_rot_l', 'elbow_flex_l', 'pro_sup_l', 'wrist_flex_l',
       'wrist_dev_l']

PCAcolls_arm = ['arm_flex_r', 'arm_add_r', 'arm_rot_r', 'elbow_flex_r', 'pro_sup_r',
       'wrist_flex_r', 'wrist_dev_r', 'arm_flex_l', 'arm_add_l', 'arm_rot_l',
       'elbow_flex_l', 'pro_sup_l', 'wrist_flex_l', 'wrist_dev_l']

# Calculate PCA for all
subdf = df_sample
PCA_all = get_PCA(subdf[PCAcolls_all], plot=True)
print('PCA for whole-body movement complexity', PCA_all)

# Calculate PCA for arms
subdf = df_sample
PCA_arms = get_PCA(subdf[PCAcolls_arm], plot=True)
print('PCA for arm movement complexity', PCA_arms)

PCA for whole-body movement complexity (5, 0.12029050407842116, 10, 0.06644314091273001)

PCA for arm movement complexity (4, 0.14837689758332867, 7, 0.09488107077615593)

Extracting features for all trials

import warnings
warnings.filterwarnings("ignore")

# Initialize the dataframe
features_df = pd.DataFrame()
features_errors = []

##########################
##### preparations #######
##########################

# These are our four dimensions for which we have aggregated movement data
groups = ['arm', 'lowerbody', 'leg', 'head']

# These are groups for BBMV
kp_arms = ['RWrist', 'RElbow', 'RShoulder', 'LWrist', 'LElbow', 'LShoulder']
kp_lower = ['RAnkle', 'RKnee', 'LAnkle', 'LKnee']
kp_legs = ['RAnkle', 'RKnee', 'LAnkle', 'LKnee']
kp_head = ['Head']
kp_keys = {'arm': kp_arms, 'lowerbody': kp_lower, 'leg': kp_legs, 'head': kp_head}

# These are cols associated with balance
balancecols = ['COPc']

# These are cols associated with acoustics
voccols = ['envelope_norm', 'envelope_change', 'f0', 'f1_clean', 'f1_clean_vel', 'f2_clean', 'f2_clean_vel', 'f3_clean', 'f3_clean_vel', 'CoG']

# These are cols for motor complexity (PCA)
PCAcolls_all = ['pelvis_tilt', 'pelvis_list', 'pelvis_rotation', 'pelvis_tx',
       'pelvis_ty', 'pelvis_tz', 'hip_flexion_r', 'hip_adduction_r',
       'hip_rotation_r', 'knee_angle_r', 'knee_angle_r_beta', 'ankle_angle_r',
       'subtalar_angle_r', 'hip_flexion_l', 'hip_adduction_l',
       'hip_rotation_l', 'knee_angle_l', 'knee_angle_l_beta', 'ankle_angle_l',
       'subtalar_angle_l', 'L5_S1_Flex_Ext', 'L5_S1_Lat_Bending',
       'L5_S1_axial_rotation', 'L4_L5_Flex_Ext', 'L4_L5_Lat_Bending',
       'L4_L5_axial_rotation', 'L3_L4_Flex_Ext', 'L3_L4_Lat_Bending',
       'L3_L4_axial_rotation', 'L2_L3_Flex_Ext', 'L2_L3_Lat_Bending',
       'L2_L3_axial_rotation', 'L1_L2_Flex_Ext', 'L1_L2_Lat_Bending',
       'L1_L2_axial_rotation', 'L1_T12_Flex_Ext', 'L1_T12_Lat_Bending',
       'L1_T12_axial_rotation', 'neck_flexion', 'neck_bending',
       'neck_rotation', 'arm_flex_r', 'arm_add_r', 'arm_rot_r', 'elbow_flex_r',
       'pro_sup_r', 'wrist_flex_r', 'wrist_dev_r', 'arm_flex_l', 'arm_add_l',
       'arm_rot_l', 'elbow_flex_l', 'pro_sup_l', 'wrist_flex_l',
       'wrist_dev_l']

PCAcolls_arm = ['arm_flex_r', 'arm_add_r', 'arm_rot_r', 'elbow_flex_r', 'pro_sup_r',
       'wrist_flex_r', 'wrist_dev_r', 'arm_flex_l', 'arm_add_l', 'arm_rot_l',
       'elbow_flex_l', 'pro_sup_l', 'wrist_flex_l', 'wrist_dev_l']

fs = 500
df_means = df_means_nonz

##########################
####### main loop ########
##########################

for file in filestotrack_clean:

    # if third element of basename is 1, we skip this file (it's from part 1 of the experiment
    if os.path.basename(file).split('_')[2] == '1':
        continue
    
    df = pd.read_csv(file)
    print('working on file: ', file)
    
    # Get metadata
    df["concept"] = df["FileInfo"].apply(lambda x: x.split("_")[1])
    df["modality"] = df["FileInfo"].apply(lambda x: x.split("_")[2])
    # if gebaren, change to gesture
    df["modality"] = df["modality"].replace({'gebaren': 'gesture', 'gebaren': 'gesture'})
    # if combinatie, change to multimodal
    df["modality"] = df["modality"].replace({'combinatie': 'multimodal', 'combinatie': 'multimodal'})
    # if geluiden, change to vocal
    df["modality"] = df["modality"].replace({'geluiden': 'vocal', 'geluiden': 'vocal'})

    df["correction"] = df["FileInfo"].apply(lambda x: x.split("_")[3])
    df.drop(columns=["FileInfo"], inplace=True)

    df_normed = df.copy()
    
    ##########################
    ##### movement feat ######
    ##########################

    if df_normed.shape[0] > 0:

        # Prepare dictionary to store the values
        body_sumfeat = {}

        # Loop over our four dimensions-groups and collect features
        for group in groups:

            # If there is no movement, we set lal values to 0 and continue
            if df_normed.shape[0] == 0:
                body_sumfeat_row = pd.DataFrame()
                continue

            else:
                # Get the columns that contain group and sum/power (these are cols with aggregated measures)
                cols = [col for col in df_normed.columns if group in col and 'sum' in col or 'power' in col]

                # Get statistics
                body_sumfeat = get_statistics(cols, df_normed, body_sumfeat, df_means)
    
                # Calculate intermittency for kinematic measures (from Pose2sim)
                # create edges for removing edge artefacts

                edge_beg = int(0.15 * fs)  # 100 ms
                edge_end = int(0.2 * fs)  # 200 ms
                df_normed_trimmed = df_normed.iloc[edge_beg:-edge_end].reset_index(drop=True)  

                speed = df_normed_trimmed[group + "_speedKin_sum"].values
                jerk = df_normed_trimmed[group + "_jerkKin_sum"].values
                intermittency = get_intermittency(jerk, speed)
                intermittency = np.log(intermittency) # natural log
                body_sumfeat[group + "_inter_Kin"] = intermittency

                # Calculate intermittency for inverse kinematic measures (from OpenSim)
                speed = df_normed_trimmed[group + "_angSpeed_sum"].values
                jerk = df_normed[group + "_angJerk_sum"].values
                intermittency = get_intermittency(jerk, speed)
                intermittency = np.log(intermittency) # natural log
                body_sumfeat[group + "_inter_IK"] = intermittency
                
                # Get bounding box of movement volume (BBMV)
                bbmv_sum = get_bbmv(df_normed_trimmed, group, kp_keys)          
                body_sumfeat[group + "_bbmv"] = bbmv_sum
            
    else:
        body_sumfeat = {}

    # Convert dictionary to dataframe
    body_sumfeat_row = pd.DataFrame({key: [value] for key, value in body_sumfeat.items()})

    # Adapt
    body_sumfeat_row = adapt_row(body_sumfeat_row)

    # Get BBMV for the whole body (so concatenate all the groups)
    bbmv_total = body_sumfeat_row[[col for col in body_sumfeat_row.columns if 'bbmv' in col and 'rate' not in col]].sum(axis=1)
    body_sumfeat_row['bbmv_total'] = bbmv_total

    ####################
    ###### balance #####
    ####################

    # Collect features (note that here we collect features for the segment of the general movement, i.e., movement of any body part)
    balance_sumfeat = {}
    balance_sumfeat = get_statistics(balancecols, df_normed, balance_sumfeat, df_means)

    # Convert dictionary to dataframe
    balance_sumfeat_row = pd.DataFrame({key: [value] for key, value in balance_sumfeat.items()})

    # Adapt
    balance_sumfeat_row = adapt_row(balance_sumfeat_row)

    ####################
    #### acoustics #####
    ####################

    if df_normed.shape[0] > 0:

        # Collect features
        vocfeat_dict = {}
        vocfeat_dict = get_statistics(voccols, df_normed, vocfeat_dict, df_means)

    else:
        vocfeat_dict = {}

    # Convert dictionary to dataframe
    vocfeat_row = pd.DataFrame({key: [value] for key, value in vocfeat_dict.items()})

    # Adapt
    vocfeat_row = adapt_row(vocfeat_row)

    # If f1-f3 are not only nan, calculate the volume of the convex hull (i.e., vowel space area)
    try:
        if df_normed_trimmed['f1_clean'].isnull().all() and df_normed_trimmed['f2_clean'].isnull().all():
            print("No f1, f2 in the dataframe")
            vocfeat_row['VSA_f1f2'] = 0
        else:
            f1_clean = df_normed_trimmed['f1_clean'].dropna()
            f2_clean = df_normed_trimmed['f2_clean'].dropna()
            VSA2d = getVSA(f1_clean, f2_clean) # get VSA
            vocfeat_row['VSA_f1f2'] = VSA2d
    except Exception as e:
        print("Error in calculating VSA: ", e)
        vocfeat_row['VSA_f1f2'] = 0
        features_errors.append({'file': file, 'VSA error': str(e)})

    # Concatenate the rows for movement, balance, and acoustics
    newrow_concat = pd.concat([body_sumfeat_row, balance_sumfeat_row, vocfeat_row], axis=1)

    ###################################
    ###### arms submovements ##########
    ###################################

    # Prepare an empty row to store the features
    newrow_arms = pd.DataFrame()

    if df_normed_trimmed.shape[0] == 0:
        print("No movement in arms")
        arms_submov = 0 # number of submovements   
        wrists_tempVar = 0 # temporal variability
    else:
        # Count submovements in wrist
        submovcols = [col for col in df_normed_trimmed.columns if "Wrist" in col and 'speed' in col]
        arms_submov = []
        wrists_tempVar = []
        for col in submovcols:
            # Calculate the number of peaks that delimit unique submovements
            scale = median_abs_deviation(df_normed_trimmed[col], scale='normal')
            prominence = scale * 0.5  # set prominence to half the MAD
            height = np.percentile(df_normed_trimmed[col], 60)  # set height to 60th percentile
            
            fullpeaks, _ = find_peaks(df_normed[col], height=height, prominence=prominence, distance=60)
            valid_mask = (fullpeaks >= edge_beg) & (fullpeaks < (df_normed.shape[0] - edge_end))  # only keep peaks in the trimmed area
            peaks_valid = fullpeaks[valid_mask]
            peaks_time = df_normed['time'].values[peaks_valid].tolist()
            arms_submov.append(len(peaks_valid))
            
            # Calculate the temporal variability of the submovements
            intervals = np.diff(peaks_time)
            tempVar = np.std(intervals)
            wrists_tempVar.append(tempVar)

    # Append to the row
    newrow_arms['arms_submovements'] = [arms_submov]
    newrow_arms['Wrist_tempVar'] = [wrists_tempVar]

    ###################################
    ######### arm symmetry ############
    ###################################

    # Calculate symmetry for arms
    symmetry_score = compute_ArmCoupling(df_normed_trimmed, keypoints=['Wrist', 'Elbow'], dimensions=['x', 'y', 'z'])
    arm_asymmetry = compute_multimovementAssymetry(df_normed_trimmed)
    newrow_arms['arm_coupling'] = symmetry_score
    newrow_arms['arm_asymmetry'] = arm_asymmetry['mean_integral_difference']


    ###################################
    ###### Motor complexity (PCA) #####
    ###################################

    if df_normed_trimmed.shape[0] == 0:
        body_nComp_80 = 0
        body_slope_80 = 0
        body_nComp_95 = 0
        body_slope_95 = 0
    else:
        bodydf = df_normed_trimmed.loc[:, PCAcolls_all]
        body_nComp_80, body_slope_80, body_nComp_95, body_slope_95 = get_PCA(bodydf)

    # Append to the row
    newrow_concat['body_nComp_80'] = body_nComp_80
    newrow_concat['body_slope_80'] = body_slope_80
    newrow_concat['body_nComp_95'] = body_nComp_95
    newrow_concat['body_slope_95'] = body_slope_95

    # Arm complexity

    if df_normed_trimmed.shape[0] == 0:
        arm_nComp_80 = 0
        arm_slope_80 = 0
        arm_nComp_95 = 0
        arm_slope_95 = 0
    else:
        armdf = df_normed_trimmed.loc[:, PCAcolls_arm]
        arm_nComp_80, arm_slope_80, arm_nComp_95, arm_slope_95 = get_PCA(armdf)

    # Append to the row
    newrow_concat['arm_nComp_80'] = arm_nComp_80
    newrow_concat['arm_slope_80'] = arm_slope_80
    newrow_concat['arm_nComp_95'] = arm_nComp_95
    newrow_concat['arm_slope_95'] = arm_slope_95

    # Concatenate newrow and newrow_arms
    newrow_final = pd.concat([newrow_concat, newrow_arms], axis=1)
    #newrow_final = newrow_concat

    # Add duration of the trial
    newrow_final['duration_trial'] = df_normed_trimmed["time"].iloc[-1] - df_normed_trimmed["time"].iloc[0]
    
    # Add metadata
    # if there is rein in TrialID, put it away
    df_normed["TrialID"] = df_normed["TrialID"].apply(lambda x: x.replace("_rein", ""))
    newrow_final["TrialID"] = df_normed["TrialID"].values[0]
    newrow_final["concept"] = df_normed["concept"].values[0]
    newrow_final["modality"] = df_normed["modality"].values[0]
    newrow_final["correction"] = df_normed["correction"].values[0]

    # Add the row the the overal df
    features_df = pd.concat([features_df, newrow_final], ignore_index=True)

    ###################################
    ########### Answer info ###########
    ###################################
    
    index = int(df_normed["TrialID"].values[0].split("_")[2])
    concept = df_normed["concept"].values[0]
    session = df_normed["TrialID"].values[0].split("_")[0] + "_" + df_normed["TrialID"].values[0].split("_")[1]

    # make subsimilaritydata which has only the rows for the current session
    subsimilaritydata = similaritydata[similaritydata["sessionID"] == session]

    # Check what was the answer and cosine similarity for the index
    try:
        answer_fol = subsimilaritydata[subsimilaritydata["trial_order"] == index]["answer"].values[0]
        answer_fol_sim = subsimilaritydata[subsimilaritydata["trial_order"] == index]["cosine_similarity"].values[0]

        # Check what was the answer and cosine similarity for index before - but only if the previous index has the same concept as the current one
        if index - 1 in subsimilaritydata["trial_order"].values and subsimilaritydata[subsimilaritydata["trial_order"] == index - 1]["word"].values[0] == concept:
            answer_prev = subsimilaritydata[subsimilaritydata["trial_order"] == index - 1]["answer"].values[0]
            answer_prev_sim = subsimilaritydata[subsimilaritydata["trial_order"] == index - 1]["cosine_similarity"].values[0]
        else:
            print("previous index doesn't have the same concept")
            answer_prev = None
            answer_prev_sim = None

        # Append to the features_df
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "answer_fol"] = answer_fol
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "answer_fol_dist"] = answer_fol_sim
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "answer_prev"] = answer_prev
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "answer_prev_dist"] = answer_prev_sim
    except IndexError:
        print("IndexError: No answer for the index", index)
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "answer_fol"] = None
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "answer_fol_dist"] = None
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "answer_prev"] = None
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "answer_prev_dist"] = None
        features_errors.append('no answer for the index ' + df_normed["TrialID"].values[0])

    # Check what is the expressibility for the concept & modality
    try:
        express = express_df[(express_df["word"] == concept) & (express_df["modality"] == df_normed["modality"].values[0])]["fit"].values[0]
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "expressibility"] = express
    except IndexError:
        #print("IndexError: No expressibility for the concept", concept, "and modality", df_normed["modality"].values[0])
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "expressibility"] = None
        features_errors.append('no expressibility for the concept ' + concept + ' and modality ' + df_normed["modality"].values[0])

    # Add response time from answerdata
    try:
        response_time = answerdata[answerdata["TrialID"] == df_normed["TrialID"].values[0]]["response_time_sec"].values[0]
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "response_time_sec"] = response_time
    except IndexError:
        print("IndexError: No response time for TrialID", df_normed["TrialID"].values[0])
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "response_time_sec"] = None
        features_errors.append('no response time for TrialID ' + df_normed["TrialID"].values[0])
    try:
        prep_time = answerdata[answerdata["TrialID"] == df_normed["TrialID"].values[0]]["prep_time_sec"].values[0]
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "prep_time_sec"] = prep_time
    except IndexError:
        print("IndexError: No prep time for TrialID", df_normed["TrialID"].values[0])
        features_df.loc[features_df["TrialID"] == df_normed["TrialID"].values[0], "prep_time_sec"] = None
        features_errors.append('no prep time for TrialID ' + df_normed["TrialID"].values[0])


# Create a concept id to group together corrections of one concept
for row in features_df.iterrows():
    row = row[1]
    features_df.loc[row.name, "concept_id"] = row["concept"] + "_" + row["modality"] + "_" + row["TrialID"].split("_")[-1]

# Get rid of columns that have only NA values
features_df = features_df.dropna(axis=1, how='all')

# Round all num cols to 3
numcols = features_df.select_dtypes(include=[np.number]).columns
features_df[numcols] = features_df[numcols].round(5)
features_errors
[{'file': 'f:\\FLESH_ContinuousBodilyEffort\\07_TS_featureExtraction\\..\\03_TS_processing\\TS_merged\\merged_1_2_47_p1.csv',
  'VSA error': "QH6421 qhull internal error (qh_maxsimplex): qh.MAXwidth required for qh_maxsimplex.  Used to estimate determinate\n\nWhile executing:  | qhull i Qt\nOptions selected for Qhull 2019.1.r 2019/06/21:\n  run-id 538753744  incidence  Qtriangulate  _pre-merge  _zero-centrum\n  _max-width  0  Error-roundoff 1.5e-12  _one-merge 7.5e-12\n  _near-inside 3.8e-11  Visible-distance 3e-12  U-max-coplanar 3e-12\n  Width-outside 6e-12  _wide-facet 1.8e-11  _maxoutside 9e-12\n\nA Qhull internal error has occurred.  Please send the input and output to\nqhull_bug@qhull.org. If you can duplicate the error with logging ('T4z'), please\ninclude the log file.\n"}]
# save
features_df.to_csv(os.path.join(datafolder, "effort_features_preliminary_nonzscored.csv"), index=False)

Aftermath

Now we want to check whether all sequences have all trials that they should have. We except that the answer is no, because we excluded some trials (see above)

# load in effort_features_preliminary.csv
features_df = pd.read_csv(os.path.join(datafolder, "effort_features_preliminary_nonzscored.csv"))

# correct concept_id such that it contains first eleemnt of TrialID (session) as a first thing
features_df["concept_id"] = features_df["TrialID"].apply(lambda x: x.split("_")[0]) + "_" + features_df["concept_id"]

# Identify how are which trials corrected
c0 = features_df[features_df["correction"] == "c0"]["concept_id"].unique()
c1 = features_df[features_df["correction"] == "c1"]["concept_id"].unique()
c2 = features_df[features_df["correction"] == "c2"]["concept_id"].unique()

only_c0 = [x for x in c0 if x not in c1 and x not in c2] # Guessed on first try
print("Guessed on first try: ", len(only_c0))
c0_c1 = [x for x in c0 if x in c1 and x not in c2] # Guessed on second try
print("Guessed on second try: ", len(c0_c1))
all_c = [x for x in c0 if x in c1 and x in c2] # Guessed on third try or never
print("Guessed on third try (or never): ", len(all_c))

# Check whether we have some data that are outside of these three categories (could be that we lost some trials in the processing pipeline)
all_concepts = features_df["concept_id"].unique()
missing = [x for x in all_concepts if x not in only_c0 and x not in c0_c1 and x not in all_c]
if len(missing) > 0:
    print("Missing: ", missing) # Currently we miss no trials
else:   
    print("No missing trials")

orphaned = [x for x in features_df["concept_id"].unique() if x not in c0]
if orphaned:
    print("Concepts with no c0:", orphaned)

# Add column correction_info to features_df
features_df["correction_info"] = features_df["correction"]

# If the concept_id is in only_c0, then correction_info is c0_only
features_df.loc[features_df["concept_id"].isin(only_c0), "correction_info"] = "c0_only"

# Save it
features_df.to_csv(os.path.join(datafolder, "features_df_nonzscored_final.csv"), index=False)
Guessed on first try:  943
Guessed on second try:  405
Guessed on third try (or never):  790
Missing:  ['10_lopen_vocal_p1', '10_langzaam_gesture_p0', '10_man_multimodal_p1', '10_vechten_vocal_p1', '11_kauwen_gesture_p0', '12_onweer_gesture_p0', '13_vrouw_multimodal_p1', '14_man_multimodal_p0', '14_ei_multimodal_p0', '15_piepen_gesture_p1', '15_vlieg_multimodal_p1', '15_stil_multimodal_p1', '15_huilen_multimodal_p1', '15_vallen_vocal_p0', '15_luidruchtig_multimodal_p0', '15_berg_vocal_p1', '18_vuur_gesture_p0', '18_kotsen_gesture_p0', '18_ei_gesture_p0', '19_ziek_multimodal_p0', '19_onweer_multimodal_p1', '19_blij_vocal_p0', '19_staan_gesture_p0', '1_verbranden_multimodal_p1', '1_vlieg_gesture_p0', '1_vogel_gesture_p1', '20_vuur_vocal_p0', '20_vrouw_multimodal_p0', '20_vechten_multimodal_p1', '24_kotsen_multimodal_p0', '24_hoorn_multimodal_p1', '24_bitter_vocal_p1', '24_verdrietig_vocal_p1', '24_springen_gesture_p1', '24_geven_multimodal_p0', '26_grommen_gesture_p1', '26_bliksem_vocal_p1', '26_vrouw_multimodal_p0', '26_misschien_multimodal_p1', '26_bang_multimodal_p1', '26_kind_vocal_p0', '26_man_gesture_p0', '26_piepen_gesture_p0', '26_weten_gesture_p0', '27_piepen_multimodal_p0', '27_kruipen_multimodal_p1', '27_luidruchtig_vocal_p0', '28_rennen_vocal_p0', '28_staan_vocal_p0', '28_klein_gesture_p0', '28_regen_gesture_p1', '29_heet_gesture_p1', '29_kind_multimodal_p0', '2_wind_gesture_p1', '2_hoog_multimodal_p1', '2_man_vocal_p0', '2_vlieg_vocal_p1', '2_niet_gesture_p1', '30_staart_multimodal_p1', '30_geven_vocal_p1', '30_hoorn_gesture_p0', '30_kat_gesture_p1', '30_man_multimodal_p0', '31_walgen_multimodal_p1', '31_weten_multimodal_p1', '31_ver_multimodal_p0', '32_horen_multimodal_p1', '33_kotsen_multimodal_p1', '33_luidruchtig_gesture_p0', '33_man_gesture_p1', '33_vis_vocal_p1', '33_vlieg_vocal_p1', '33_dik_multimodal_p0', '34_kauwen_gesture_p0', '34_ver_vocal_p0', '34_vlieg_multimodal_p0', '35_kruipen_vocal_p0', '35_jagen_vocal_p0', '36_klein_vocal_p0', '36_boos_multimodal_p0', '36_onweer_multimodal_p0', '36_ziek_gesture_p0', '37_hoog_gesture_p0', '37_scherp_gesture_p1', '37_horen_vocal_p1', '37_vlieg_multimodal_p0', '38_zout_gesture_p0', '38_vallen_multimodal_p0', '38_kind_multimodal_p1', '38_ver_multimodal_p1', '38_luidruchtig_vocal_p0', '38_dik_vocal_p0', '38_lopen_vocal_p0', '39_zoemen_gesture_p1', '39_hond_gesture_p1', '39_niet_gesture_p0', '3_dik_vocal_p0', '3_horen_vocal_p0', '3_vuur_gesture_p1', '40_vrouw_gesture_p1', '40_bliksem_gesture_p1', '40_niet_multimodal_p0', '40_vlieg_vocal_p0', '40_geven_vocal_p0', '41_zoemen_gesture_p1', '45_misschien_multimodal_p1', '45_bliksem_gesture_p0', '45_ik_multimodal_p0', '46_blij_gesture_p1', '46_hoog_multimodal_p1', '48_water_vocal_p1', '48_hoorn_multimodal_p0', '49_jagen_multimodal_p0', '49_kind_multimodal_p1', '49_springen_multimodal_p1', '50_geur_multimodal_p1', '50_onweer_multimodal_p1', '50_vangen_multimodal_p1', '50_man_gesture_p1', '51_walgen_multimodal_p1', '52_vangen_vocal_p0', '52_kauwen_vocal_p0', '52_klein_vocal_p1', '52_zout_multimodal_p0', '52_kind_multimodal_p1', '52_zoet_multimodal_p1', '52_staan_multimodal_p1', '53_vrouw_gesture_p1', '53_geur_vocal_p0', '53_ik_multimodal_p0', '54_niet_gesture_p1', '54_zoemen_gesture_p0', '54_staan_vocal_p1', '54_bliksem_vocal_p1', '55_zoemen_gesture_p1', '55_kind_gesture_p1', '55_jagen_gesture_p1', '56_jagen_gesture_p0', '56_zoemen_gesture_p1', '56_piepen_multimodal_p0', '58_geur_gesture_p1', '58_regen_vocal_p0', '58_huilen_gesture_p0', '59_slang_gesture_p0', '5_vallen_vocal_p1', '5_horen_multimodal_p0', '5_zoemen_multimodal_p1', '5_misschien_multimodal_p1', '5_zoet_multimodal_p1', '5_water_multimodal_p1', '5_vuur_vocal_p0', '5_regen_vocal_p0', '5_bang_vocal_p0', '60_berg_gesture_p0', '60_dood_multimodal_p0', '60_walgen_multimodal_p0', '60_verdrietig_gesture_p0', '61_bitter_vocal_p1', '61_luidruchtig_vocal_p1', '61_hoog_multimodal_p1', '61_geur_gesture_p0', '61_niet_multimodal_p0', '61_zout_vocal_p0', '62_koud_multimodal_p0', '62_zoemen_gesture_p0', '62_kauwen_multimodal_p1', '62_ei_vocal_p0', '62_misschien_vocal_p0', '62_klein_gesture_p0', '63_geven_multimodal_p0', '63_heet_multimodal_p1', '63_zout_multimodal_p1', '63_weten_multimodal_p1', '63_luidruchtig_vocal_p1', '64_bijten_multimodal_p1', '64_staan_multimodal_p1', '64_sterk_vocal_p0', '64_vuur_multimodal_p0', '65_ver_vocal_p0', '66_berg_gesture_p1', '69_staan_multimodal_p0', '69_jagen_multimodal_p1', '6_zoet_gesture_p1', '6_zuigen_gesture_p1', '6_staart_vocal_p0', '6_zoemen_gesture_p0', '6_walgen_gesture_p0', '6_onweer_gesture_p0', '70_ei_gesture_p1', '71_berg_gesture_p0', '71_likken_vocal_p1', '72_kind_gesture_p1', '8_geur_vocal_p1', '8_geven_vocal_p1', '8_verbranden_multimodal_p1', '8_zoemen_gesture_p0', '8_luidruchtig_gesture_p0', '8_zoet_multimodal_p0', '9_niet_gesture_p1', '9_weten_vocal_p1', '9_vogel_multimodal_p1']
Concepts with no c0: ['10_lopen_vocal_p1', '10_vechten_vocal_p1', '12_onweer_gesture_p0', '13_vrouw_multimodal_p1', '14_man_multimodal_p0', '15_piepen_gesture_p1', '15_stil_multimodal_p1', '15_huilen_multimodal_p1', '15_vallen_vocal_p0', '15_berg_vocal_p1', '18_vuur_gesture_p0', '18_ei_gesture_p0', '19_ziek_multimodal_p0', '19_onweer_multimodal_p1', '19_staan_gesture_p0', '1_vlieg_gesture_p0', '20_vuur_vocal_p0', '20_vrouw_multimodal_p0', '20_vechten_multimodal_p1', '24_kotsen_multimodal_p0', '24_hoorn_multimodal_p1', '24_verdrietig_vocal_p1', '24_springen_gesture_p1', '24_geven_multimodal_p0', '26_grommen_gesture_p1', '26_bliksem_vocal_p1', '26_misschien_multimodal_p1', '26_bang_multimodal_p1', '26_kind_vocal_p0', '26_piepen_gesture_p0', '26_weten_gesture_p0', '28_staan_vocal_p0', '28_klein_gesture_p0', '28_regen_gesture_p1', '29_heet_gesture_p1', '29_kind_multimodal_p0', '2_wind_gesture_p1', '2_vlieg_vocal_p1', '2_niet_gesture_p1', '30_staart_multimodal_p1', '30_hoorn_gesture_p0', '30_kat_gesture_p1', '30_man_multimodal_p0', '31_ver_multimodal_p0', '33_kotsen_multimodal_p1', '33_man_gesture_p1', '33_vis_vocal_p1', '33_dik_multimodal_p0', '34_vlieg_multimodal_p0', '35_kruipen_vocal_p0', '35_jagen_vocal_p0', '36_klein_vocal_p0', '36_boos_multimodal_p0', '36_ziek_gesture_p0', '37_hoog_gesture_p0', '37_scherp_gesture_p1', '37_horen_vocal_p1', '37_vlieg_multimodal_p0', '38_vallen_multimodal_p0', '38_ver_multimodal_p1', '38_dik_vocal_p0', '38_lopen_vocal_p0', '39_zoemen_gesture_p1', '39_hond_gesture_p1', '39_niet_gesture_p0', '3_dik_vocal_p0', '3_vuur_gesture_p1', '40_vrouw_gesture_p1', '40_bliksem_gesture_p1', '40_niet_multimodal_p0', '40_geven_vocal_p0', '45_bliksem_gesture_p0', '45_ik_multimodal_p0', '46_blij_gesture_p1', '46_hoog_multimodal_p1', '48_hoorn_multimodal_p0', '49_jagen_multimodal_p0', '49_kind_multimodal_p1', '50_geur_multimodal_p1', '50_onweer_multimodal_p1', '50_vangen_multimodal_p1', '50_man_gesture_p1', '51_walgen_multimodal_p1', '52_vangen_vocal_p0', '52_zout_multimodal_p0', '52_kind_multimodal_p1', '52_staan_multimodal_p1', '53_vrouw_gesture_p1', '54_niet_gesture_p1', '54_staan_vocal_p1', '54_bliksem_vocal_p1', '55_zoemen_gesture_p1', '56_jagen_gesture_p0', '56_zoemen_gesture_p1', '58_geur_gesture_p1', '58_regen_vocal_p0', '58_huilen_gesture_p0', '5_vallen_vocal_p1', '5_horen_multimodal_p0', '5_zoemen_multimodal_p1', '5_misschien_multimodal_p1', '5_zoet_multimodal_p1', '5_water_multimodal_p1', '5_vuur_vocal_p0', '5_bang_vocal_p0', '60_dood_multimodal_p0', '60_verdrietig_gesture_p0', '61_bitter_vocal_p1', '61_luidruchtig_vocal_p1', '61_hoog_multimodal_p1', '61_geur_gesture_p0', '61_niet_multimodal_p0', '61_zout_vocal_p0', '62_koud_multimodal_p0', '62_zoemen_gesture_p0', '62_kauwen_multimodal_p1', '62_misschien_vocal_p0', '62_klein_gesture_p0', '63_geven_multimodal_p0', '63_heet_multimodal_p1', '63_zout_multimodal_p1', '64_staan_multimodal_p1', '64_sterk_vocal_p0', '64_vuur_multimodal_p0', '65_ver_vocal_p0', '69_staan_multimodal_p0', '69_jagen_multimodal_p1', '6_zoet_gesture_p1', '6_zoemen_gesture_p0', '6_walgen_gesture_p0', '70_ei_gesture_p1', '71_likken_vocal_p1', '72_kind_gesture_p1', '8_geven_vocal_p1', '8_verbranden_multimodal_p1', '8_zoemen_gesture_p0', '8_luidruchtig_gesture_p0', '9_niet_gesture_p1', '9_vogel_multimodal_p1']

Now we also want to add metadata for modelling

# create pcn_id
features_df['pcn_id'] = features_df['TrialID'].apply(lambda x: x.split('_')[0] + '_' + x.split('_')[-1])

# this is all demodata
demo = pd.read_csv(os.path.join(metadatafolder, 'Demographics_all', 'all_demodata.csv'))

# leave only cols with pcn_ID, BFI_ and Familiarity
demo = demo[['pcn_ID', 'BFI_extra', 'BFI_agree', 'BFI_consc', 'BFI_negemo', 'BFI_open', 'Familiarity']]
demo['pcn_ID'] = demo['pcn_ID'].apply(lambda x: x.split('_')[0] + '_p' + x.split('_')[1])

# correct 19_p7 to 19_p1
demo['pcn_ID'] = demo['pcn_ID'].replace('19_p7', '19_p1')
# 32_p1 to 32_p0
demo['pcn_ID'] = demo['pcn_ID'].replace('32_p1', '32_p0')
# and 32_p2 to 32_p1
demo['pcn_ID'] = demo['pcn_ID'].replace('32_p2', '32_p1')

# merge on pcn_ID
features_demo = pd.merge(features_df, demo, left_on='pcn_id', right_on='pcn_ID', how='left')

# save
features_demo.to_csv(os.path.join(datafolder, 'features_df_nonzscored_final_with_demo.csv'), index=False)

We will be also dealing with some loss of the data because some features was not possible to extract. Let’s see for which features it’s the most of the data (>50%)

lowerbody_power_peak_std    0.763550
leg_power_peak_std          0.627145
dtype: float64

Creating z-scored dataframe for relative effort

For exploratory analysis that utilizes PCA and XGBoost (not part of the current study), we also want to create dataframe that contains z-scored features (normalized per participant)

features_df_zscored = features_df.copy()

# create pcnID which is first and last element of TrialID
features_df_zscored["pcnID"] = features_df_zscored["TrialID"].apply(lambda x: x.split("_")[0] + "_" + x.split("_")[-1])
numcols = features_df_zscored.select_dtypes(include=[np.number]).columns

# exclude these 
excludecols = ['answer_fol_dist', 'answer_prev_dist', 'expressibility', 'response_time_sec', 'prep_time_sec']
numcols = [col for col in numcols if col not in excludecols]

# z-score within pcnID
features_df_zscored[numcols] = features_df_zscored.groupby("pcnID")[numcols].transform(lambda x: (x - x.mean()) / x.std())

# save it
features_df_zscored.to_csv(os.path.join(datafolder, "features_df_zscored_feat_final.csv"), index=False)

# save it also separately per modality
features_df_zscored_gesture = features_df_zscored[features_df_zscored["modality"] == "gebaren"]
features_df_zscored_combination = features_df_zscored[features_df_zscored["modality"] == "combinatie"]
features_df_zscored_vocal = features_df_zscored[features_df_zscored["modality"] == "geluiden"]
features_df_zscored_gesture.to_csv(os.path.join(datafolder, "features_df_zscored_feat_gesture.csv"), index=False)
features_df_zscored_combination.to_csv(os.path.join(datafolder, "features_df_zscored_feat_combination.csv"), index=False)
features_df_zscored_vocal.to_csv(os.path.join(datafolder, "features_df_zscored_feat_vocal.csv"), index=False)

References

Berisha, V., Sandoval, S., Utianski, R., Liss, J., & Spanias, A. (2014). Characterizing the distribution of the quadrilateral vowel space area. The Journal of the Acoustical Society of America, 135(1), 421–427. https://doi.org/10.1121/1.4829528
Ćwiek, A., Fuchs, S., Pouw, W., & Kadavá, Š. (2026). Self-reported expressibility predicts communicative success: Open dataset, validation, and simulation. Behavior Research Methods. https://doi.org/10.3758/s13428-026-03071-4
Daffertshofer, A., Lamoth, C. J. C., Meijer, O. G., & Beek, P. J. (2004). PCA in studying coordination and variability: A tutorial. Clinical Biomechanics (Bristol, Avon), 19(4), 415–428. https://doi.org/10.1016/j.clinbiomech.2004.01.005
Hogan, N., & Sternad, D. (2009). Sensitivity of Smoothness Measures to Movement Duration, Amplitude, and Arrests. Journal of Motor Behavior, 41(6), 529–534. https://doi.org/10.3200/35-09-004-RC
Kadavá, Š., Ćwiek, A., Fuchs, S., & Pouw, W. (2024). What do we mean when we say gestures are more expressive than vocalizations? An experimental and simulation study. Proceedings of the Annual Meeting of the Cognitive Science Society, 46(0). https://escholarship.org/uc/item/2mp1v3v5
Pouw, W., Dingemanse, M., Motamedi, Y., & Özyürek, A. (2021). A Systematic Investigation of Gesture Kinematics in Evolving Manual Languages in the Lab. Cognitive Science, 45(7), e13014. https://doi.org/10.1111/cogs.13014
Seth, A., Hicks, J. L., Uchida, T. K., Habib, A., Dembia, C. L., Dunne, J. J., Ong, C. F., DeMers, M. S., Rajagopal, A., Millard, M., Hamner, S. R., Arnold, E. M., Yong, J. R., Lakshmikanth, S. K., Sherman, M. A., Ku, J. P., & Delp, S. L. (2018). OpenSim: Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement. PLOS Computational Biology, 14(7), e1006223. https://doi.org/10.1371/journal.pcbi.1006223
Trujillo, J., Simanova, I., Bekkering, H., & Özyürek, A. (2018). Communicative intent modulates production and comprehension of actions and gestures: A Kinect study. Cognition, 180, 38–51. https://doi.org/10.1016/j.cognition.2018.04.003
Xiong, Y., Quek, F., & Mcneill, D. (2002). Hand gesture symmetric behavior detection and analysis in natural conversation. 179–184. https://doi.org/10.1109/ICMI.2002.1166989
Yan, Y., Goodman, J. M., Moore, D. D., Solla, S. A., & Bensmaia, S. J. (2020). Unexpected complexity of everyday manual behaviors. Nature Communications, 11(1), 3564. https://doi.org/10.1038/s41467-020-17404-0
Żywiczyński, P., Placiński, M., Sibierska, M., Boruta-Żywiczyńska, M., Wacewicz, S., Meina, M., & Gärdenfors, P. (2024). Praxis, demonstration and pantomime: A motion capture investigation of differences in action performances. Language and Cognition, 1–28. https://doi.org/10.1017/langcog.2024.8