<jupyter_start><jupyter_text>Analysis of pupil data======================There are many valid approaches to the analysis of pupillometry data, but the optimal approach will depend on the type of experiment being run and the research question in mind. [Kelbsch et al. (2019)](https://www.frontiersin.org/articles/10.3389/fneur.2019.00129/full) provide an informative view on standards in pupillometry, and there are [some](https://github.com/samhforbes/PupillometryR) [great](https://github.com/beOn/cili) [examples](https://github.com/ihrke/pypillometry) of community-developed packages for streamlining pupillometry data analysis. Whilst it is always worth exploring the wider options that are available, *PyPlr* includes a set of pandas-reliant (you won't regret [learning pandas!](https://bitbucket.org/hrojas/learn-pandas/src/master/)) scripting tools for implementing a standard data processing pipeline that is optimised to work with some of the idiosynchrasies of Pupil Labs data. So far we have found these tools to be adequate for analysing our own data, but we welcome suggestions for improvements. See [here](https://github.com/pupil-labs/pupil-tutorials/blob/master/01_load_exported_data_and_visualize_pupillometry.ipynb) for a general introduction to analysing Pupil Labs pupillometry data. What follows is a typical *PyPlr* analysis pipeline with tools from `pyplr.utils` and `pyplr.preproc`. Export with Pupil Player------------------------The first step in our pipeline is to export the data using [Pupil Player](https://docs.pupil-labs.com/core/software/pupil-player/), making sure the required plugins (e.g., Annotation Player plugin) are enabled. Currently this must be done individually for each recording. Below is printed the file structure of a typical recording after export (more info on recording format [here](https://docs.pupil-labs.com/developer/core/recording-format/recording-format)).<jupyter_code>from pyplr import utils

# Pupil Labs recording directory
rec_dir = '/Users/jtm/OneDrive - Nexus365/protocols/pipr_protocol/JTM'
utils.print_file_structure(rec_dir)<jupyter_output>JTM/
    .DS_Store
    annotation.pldata
    annotation_player.pldata
    annotation_player_timestamps.npy
    annotation_timestamps.npy
    blinks.pldata
    blinks_timestamps.npy
    eye1.intrinsics
    eye1.mp4
    eye1_lookup.npy
    eye1_timestamps.npy
    fixations.pldata
    fixations_timestamps.npy
    gaze.pldata
    gaze_timestamps.npy
    info.player.json
    notify.pldata
    notify_timestamps.npy
    pupil.pldata
    pupil_timestamps.npy
    surface_definitions_v01
    user_info.csv
    world.intrinsics
    world.mp4
    world_lookup.npy
    world_timestamps.npy
    analysis/
        plr_extraction.png
        plr_metrics.csv
        plr_overall_metrics.png
        processed.csv
        pupil_processing.png
    exports/
        000/
            annotations.csv
            export_info.csv
            gaze_positions.csv
            pupil_gaze_positions_info.txt
            pupil_positions.csv
            world.mp4
            world_timestamps.csv
            world_timestamps[...]<jupyter_text>Load exported data------------------Now we can set up some constants and use `pyplr.utils` to get things moving. `new_subject(...)` returns a dictionary with the root directory, recording id, data directory and a newly made output directory for the analysis results. Then, passing the data directory to `load_pupil(...)`, we load the specified columns from the *pupil_positions.csv* exported data.<jupyter_code># Columns to load
use_cols = ['confidence',
            'method',
            'pupil_timestamp',
            'eye_id',
            'diameter_3d',
            'diameter']

# Get a handle on a subject
s = utils.new_subject(
    rec_dir, export='000', out_dir_nm='pyplr_analysis')

# Load pupil data
samples = utils.load_pupil(
    s['data_dir'], eye_id='best', cols=use_cols)
samples<jupyter_output>************************************************************
*************************** JTM ****************************
************************************************************
Loaded 48552 samples<jupyter_text>Preprocessing-------------Preprocessing pupillometry data is relatively straight-forward and typically involves dealing with signal loss due to eye blinks and smoothing out any high frequency noise. There are [sophisticated algorithms](https://osf.io/jyz43/) for detecting blinks in a pupil time course, but Pupil Player has its own [Blink Detection](https://docs.pupil-labs.com/core/software/pupil-capture/blink-detection) plugin that detects blinks based on rapid changes in the [confidence metric](https://docs.pupil-labs.com/core/terminology/confidence). Whilst you can export these blink events and use the timestamps to index and mask the pupil timecourse, we find that it is effective simply to mask the pupil data with thresholds on the first derivative and confidence metric, and then to follow up with linear interpolation and smoothing.<jupyter_code>from pyplr import graphing
from pyplr import preproc

# Sampling frequency
SAMPLE_RATE = 120

# Pupil columns to analyse
pupil_cols = ['diameter_3d', 'diameter']

# Make figure for processing
f, axs = graphing.pupil_preprocessing(nrows=5, subject='Example')

# Plot the raw data
samples[pupil_cols].plot(title='Raw', ax=axs[0], legend=True)
axs[0].legend(loc='center right', labels=['mm', 'pixels'])

# Mask first derivative
samples = preproc.mask_pupil_first_derivative(
    samples, threshold=3.0, mask_cols=pupil_cols)
samples[pupil_cols].plot(
    title='Masked 1st deriv (3*SD)', ax=axs[1], legend=False)

# Mask confidence
samples = preproc.mask_pupil_confidence(
    samples, threshold=0.8, mask_cols=pupil_cols)
samples[pupil_cols].plot(
    title='Masked confidence (<0.8)', ax=axs[2], legend=False)

# Interpolate
samples = preproc.interpolate_pupil(
    samples, interp_cols=pupil_cols)
samples[pupil_cols].plot(
    title='Linear interpolation', ax=axs[3], legend=False)

# Smooth  
samples = preproc.butterworth_series(
    samples, fields=pupil_cols, filt_order=3, 
    cutoff_freq=4/(SAMPLE_RATE/2))
samples[pupil_cols].plot(
    title='3rd order Butterworth filter with 4 Hz cut-off', 
    ax=axs[4], legend=False);<jupyter_output><empty_output><jupyter_text>Extraction----------Having cleaned the timecourse we need to extract the events of interest. This requires the annotation events sent during the recording, which have the timestamps needed for extraction. These can be loaded with `utlis.load_annotations(...)`.<jupyter_code>events = utils.load_annotations(s['data_dir'])
events<jupyter_output>Loaded 6 events<jupyter_text>Now we can pass our samples and events to `utils.extract(...)` along with some parameters defining the number of samples to extract and how much to offset the data relative to the event. Here we extract-65 second epochs and offset by a 5 second baseline period.<jupyter_code># Number of samples to extract and which sample
# should mark the onset of the event
DURATION = 7800
ONSET_IDX = 600

# Extract the event ranges
ranges = utils.extract(
    samples, 
    events, 
    offset=-ONSET_IDX, 
    duration=DURATION, 
    borrow_attributes=['color'])
ranges<jupyter_output>Extracted ranges for 6 events<jupyter_text>``utils.extract(...)`` returns a pandas DataFrame with a two-level (event, onset) [hierarchichal index](https://pandas.pydata.org/pandas-docs/stable/user_guide/advanced.html). The original index to the data is preserved in a column along with any attributes borrowed from the annotations. Now is a good time to create additional columns expressing the pupil data as percent signal change from baseline.<jupyter_code># Calculate baselines
baselines = ranges.loc[:, range(0, ONSET_IDX), :].mean(level=0)

# New columns for percent signal change
ranges = preproc.percent_signal_change(
    ranges, baselines, pupil_cols)
ranges<jupyter_output><empty_output><jupyter_text>Plotting and parametrisation----------------------------It is common practice when analysing the PLR to describe it in terms of paremeters relating to time, velocity and acceleration. So once you have an array of data representing a pupil's response to light, be it from a single trial or an average of multiple trials, simply pass it to `pyplr.plr.PLR` along with some basic parameters.<jupyter_code>from pyplr.plr import PLR

average_plr = ranges.mean(level=1)['diameter_3d'].to_numpy()
plr = PLR(average_plr, 
          sample_rate=SAMPLE_RATE, 
          onset_idx=ONSET_IDX, 
          stim_duration=1)<jupyter_output><empty_output><jupyter_text>Now you can do a quick plot with options to display the velocity and acceleration profiles and the derived parameters:<jupyter_code>fig = plr.plot(vel=True, acc=True, print_params=True)<jupyter_output><empty_output><jupyter_text>Or simply access the parameters as a DataFrame.<jupyter_code>params = plr.parameters()
params<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11184357541899441, 'fraction_numerical': 0.009832402234636871, 'mean_word_length': 4.441337386018237, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 33, 'https://': 12, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '426566', 'license': 'MIT', 'n_tokens_mistral': 2775, 'n_tokens_neox': 2637, 'n_words': 828, 'path': 'docs/04d_analysis.ipynb', 'repo_name': 'willpatera/cvd_pupillometry', 'size': '8950'}