Analysis of Pupillary Data

1. Analysis of Pupillary Data Society for Psychophysiological Research Workshop, Boston, September 14, 2011

2. Approaches to Pupillary Data Analysis • Scale data to mm (can be performed at several points up to statistical analysis, discussed earlier) • Identify and correct for blinks and other artifacts • Smooth data through filtering algorithms • Segregate trials according to conditions (for repeated measures designs) • Identify stimulus onset for stimulus-locked analyses • Identify time of RT for response-locked analyses • Apply signal averaging for repeated measures designs as appropriate • Determine critical variables to be used for analyses (absolute diameters, amplitude of change, peak responses, time of critical responses, frequency analyses, etc. • Consider more sophisticated approaches including: Principal Components Analysis, Fourier/Wavelet Analyses, waveform comparison models

3. Notes on Pupillary Signal Averaging • Signal averaging approaches may be applied to pupillary data in repeated measures experiments just as they are in event-related potential experiments, with all of the same caveats • The advantage is that small changes on the order of 0.01 mm or better (sometimes as detailed as several one-thousandths of a mm) can be reliable obtained. • The pupil is affected by many types of ongoing changes in visual stimulation, CNS activity, etc., which make signal averaging advantageous, but has a higher signal-to-noise ration than the EEG, so that fewer trials are needed to extract a good averaged response, perhaps with minimums of: • 3-5 trials for light reactions • 5-10 trials for motor-initiated dilations • 5-10 or more trials for other cognitive and emotional events • Note that Orienting or highly novel events are limited in the number of trials that can be obtained because of the nature of the task

4. Pupillary Response at Visual ThresholdG. Hakerem & S. Sutton, Nature, 1966

5. Pupillary Light Reactions • A single pupillary light reaction is not likely to be a good representation; usually at least 3-5 responses should be averaged • The light reaction in darkness is a function of both stimulus characteristics and time of dark adaptation • When a series of light stimuli are presented, the first response will start at a larger diameter and show a stronger contraction than subsequent stimuli. • Unless standard stimuli and interstimulus intervals are employed, data are likely to be interpretable only through complex non-linear analyses, not simple averaging across all stimuli. • At high rates of visual stimulation, pupillary responses become integrated (i.e., individual light reactions not identifiable).

6. (for description of multiple measures, see Steinhauer et al., Psychophysiology, 1992)

7. Baseline Start of Contraction Maximum Velocity Minimum Diameter Redilation

8. Correcting for Artifacts • Identify changes in diameter between subsequent points that are too large to be physiologically meaningful, or where pupil diameter is out of range (e.g., scaled value is 0 or 11 mm) • If the blink/closure/artifact occurs at a critical time (during a brief visual stimulus, at the minimum of constriction, at the peak of a dilation interval) then it is usually best to discard the trial entirely • To apply a linear correction, identify both the last good data point before the artifact (p1) and the first good data point after the artifact (pL). Take the difference between them (D = pL - p1), divide by the number of points in the artifact plus one (D / (nbad +1)), and add this amount to p1 incrementally for each of the substituted points. • In our laboratories, we have found linear interpretation to be an appropriate and easily implemented correction to most types of data, although several researchers have employed more complex curve fitting to individual data.

9. Filtering Data (1) • Why filter if the pupil is a (relatively) slowly changing measure? Except for some light reactions, factors such as magnification, measurement, and signal processing introduce noise aspects that are unrelated to true pupillary changes. • Filtering eliminates small artifacts that can be seen in the electronically derived continuous pupil waveform. • Should filtering be performed before or after signal averaging? Actually, it does not matter computationally, but makes subsequent artifact correction and component identification easier if performed first. • Do not filter if your are concerned about very precise time transients (e.g., start of light reaction to nearest 0.016 msec) as filtering smooths any sharp peaks and may shift some peak activity slightly.

10. Filtering Data (2) • There are a number of different approaches to filtering. We will describe a simple procedure described by Ruchkin & Glaser (1978). • Each point is recalculated beginning with the original data set, using a band of points centered on each original point. • The larger the number of points, the greater the effective filter. • Assume an initial frequency of 60 Hz (16.7 msec between samples) • The frequency of the filter f0 is characterized as 1/[(2L+1)T], where L is the number of points below and above the center point, and T is the sampling interval (16.7 ms). • Ruchkin and Glaser noted that the transfer function has a high secondary peak, so that the “first pass” is repeated again on the same data, thus compromising a “two pass” filter. • The half power frequency (or -3 dB reduction) yields a gain of 0.707, so that the effective bandpass is 0.44 * f0 . Ruchkin, D.S. & Glaser, E.M.,1978

11. Filtering Data (3) To calculate filter characteristic: • a) multiply sampling rate by 0.44; for 60 Hz, this equals 60 Hz * 0.44 = 26.4 Hz • b) divide by the number of points (2L+1) • c) Thus, for data recorded at 60 Hz, a 3 point filter -> Low pass of 8.8 Hz a 5 point filter -> Low pass of 5.3 Hz a 7 point filter -> Low pass of 3.8 Hz a 9 point filter -> Low pass of 2.9 Hz

12. Filtering Data (4) To operationalize filter (here for L=1, 3 point filter) 1) set up 3 arrays a, b, c 2) read data into array a 3) for each point x in array a, substitute the mean of xa-1, xa, and xa+1 for xb in array b; that is, xb = (xa-1+ xa+ xa+1)/3 (repeat leading and trailing values for first and final points) 4) repeat, calculating array c values from array b This filter results in no phase shifts. Remember that high and low amplitude values will be attenuated

13. Tonic Data • In some limited cases, only overall average diameter is needed (is pupil larger for different sustained tasks, is it smaller at end of a session?). Define a specific interval and take the average of points, possibly report variation.

14. Phasic Data (data of D. Friedman et al., EEG Cl. Neuro., 1973) Prestimulus or Baseline: Often use average of at least 1000 msec, minimum 200 msec, prior to stimulus. Peak Dilation: Often seen 1200-1400 msec after critical stimulus. Use either peak or averageof several points surrounding peak. During recording in light, earlier peaks or increases in diameter related to inhibition of the parasympathetic system may be seen. Often, integrated activity over much longer durations may be desired

15. Digit Span Task Often, activity over much longer durations may need to be analyzed

16. Probability Effects, Auditory Counting Task

17. Probability Effects, Auditory Choice Reaction Task

18. Principal Components Analysis • Data for entire waveforms are entered into PCA, typically each average, each condition, each subject (can also be performed for individual subjects) • Factors are based on similar variance effects across conditions, may or may not correspond to peaks or sustained activity periods • Factors are extracted in decreasing order of variance explained • Factor loadings reflect influence of variables across time points • A factor score – representing time by loading effects – is obtained for each factor for each waveform • Factor scores may be used in repeated measures analyses

21. Steinhauer & Hakerem, 1992

23. Pupil: Time Frequency Analysis Pupillary Sleepiness Test Franzen et al (in prep) seconds

24. Wavelet analysis of Oscillatory Activity S1a Hypothesis: Decreased oscillatory activity in patient groups would be indicative of disruption of central parasympathetic control. Pupil diameter was assessed in darkness for 11 minutes as subjects fixated three small red LEDs. The pupil was digitized 60 times/sec with a resolution of 0.05 mm using an ISCAN ETL-400 system. Initial Data Reduction: The pupil was first corrected for blinks (S1a), then smoothed and detrended (correcting for baseline, S1b). PUI can be computed as the absolute change in diameter across time, essentially a string measure of activity, but this includes all frequencies of change. However, Ludtke et al. (1998) suggested using low frequency slow changes only. A wavelet analysis implemented in MatLab assessed frequencies at each sampling point, expressing PUI as the sum of power in the 0-0.6 Hz frequency range (S1c). S1b S1c

25. Power Derived from Wavelet Analysis

26. MATLAB IMPLEMENTATIONS Thanks to Greg Siegle, Ph.D. Depts. Of Psychiatry and Psychology University of Pittsburgh

27. The environment • Matlab is a general purpose language for mathematical and graphical operations • All operations can be done interactively or by calling stored functions • The basic computational unit is a row x col matrix x=[1 2 3;4 5 6] x = 1 2 3 4 5 6 Scalers (e.g., 3) and vectors (e.g., [3 4 5]) are supported

28. Plotting • pts=0:pi/20:2.*pi; x=cos(pts); y=sin(pts); • plot(sin(pts)); plot(pts,y); plot(y); • plot (x,y); • axis on/off, clf • xlabel(‘time’); ylabel(‘dilation’); title(‘my graph’); • plot([x; y]’); • plot([x; y]); • legend (‘cos’,’sin’); • Specialized plots, e.g., errorbar(pts,y,x); • hist(pts) • figure; figure(2); • subplot(2,2,1); plot(x,y);

29. Functions • Each function, accessible as a command, is stored in its own file e.g., storing a function called parabola.m function y=parabola(x) % squares its input % usage: y=parabola(x) y=x.^2; Allows you to say “plot(parabola (1:10))” and help (parabola)

30. Operations in functions For ct=1:5 b=b+1; end while ct<5 b=b.^2; ct=ct+1; end if a>2 b=4; else b=5; end Scientific Functions Trig: sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh Rounding: floor, ceil, round, fix Modular: rem, mod Expon.: exp, log, log2, log10, sqrt Primes: factor, primes Matrix: det, inv, pinv, eig, svd, fft and many more Polynomials: roots, polyfit, polyval

31. The Pupil Toolkit The Pupil Toolkit Greg Siegle, Ph.D. University of Pittsburgh, School of Medicine

32. The Pupil Toolkit Goals • Be able to read in and process a single subject’s pupillary data in 1 step. • Have average graphs to show subjects directly after running them • Have relevant diagnostic information produced • Have a stock set of statistics automatically calculated on pupil waveforms when they are read in. • Make it easy to analyze new experiments

33. The Pupil Toolkit Standard flow for 1-subject’s pupil data • Read in the data • Clean the data (i.e., eliminate blinks) • Segment the data into trials (from markers or eprime) • Drop bad trials • Make condition-related averages • Graph the condition-related averages

34. Single subject diagnostic quality control graphs available 1 minute after testing The Pupil Toolkit via p=readiscan(‘1002kvid.isc’); Or run a stored procedure unique to your experiment, e.g. p=procsilkkvidexample(fname); Which: 1) reads the data file via p=readiscan(fname) 2) preprocesses it via p=stublinks(p) 3) segments it into trials via p=segmentpupiltrials(p) 4) calculates statistics via p.stats=pupiltrialstats(p)

35. The Pupil Toolkit >> p=procsilkkvidexample(1002) Raw data >> p=procsilkkvidexample(1002) found 34280 records in data dropped 1 of 51 trials p = FileName: '1002KVID.isc' header: [1x1 struct] EventTrain: [34280x1 double] RescaleData: [1x34280 double] BlinkTimes: [34280x1 logical] NoBlinks: [34280x1 double] NoBlinksUnsmoothed: [1x34280 double] NoBlinksDetrend: [34280x1 double] EventTicks: [156x1 double] EventCodes: [156x1 double] RescaleFactor: 60 EventTimes: [156x1 double] AllSeconds: [34280x1 double] TrialStarts: [51x1 double] TrialEnds: [51x1 double] TrialLengths: [51x1 double] NumTrials: 51 StimLatencies: [51x1 double] StimTypes: [51x1 double] TrialTypesNoDrops: [51x1 double] PupilTrials: [51x661 double] EventTrials: [51x661 double] BlinkTrials: [51x661 double] DetrendPupilTrials: [51x661 double] NormedPupTrials: [51x661 double] NormedDetrendPupTrials: [51x661 double] TrialSeconds: [1x661 double] Suspect: [1x51 logical] stats: [1x1 struct] drops: [1x51 logical] TrialTypes: [51x1 double] Conditions: [1 2 4] ConditionMeans: [3x661 double] ConditionSds: [3x661 double] condstats: [1x1 struct] After trial segmentation Valence Identification Task

36. p=procvfg(5107,'numorder',2); The Pupil Toolkit Digit Sorting Task (sort 3, 4 or 5 numbers)

37. The Pupil Toolkit Reading data from different pupilometers • ISCAN • readiscan • readiscan05text • readiscanbehavobstext • ASL – to use must first install ASL matlab drivers • readasl2006 • readasl2007 • readasl2008 • readasltext • readasltextlunalab

38. The Pupil Toolkit How the blink elimination routine works data = stublinks(data,graphics,lrtask,manualblinks,lowthresh) • Smooths the data using a 3 pt flat filter applied twice • Identifies blinks via many criteria (e..g, >.5mm change in 1 sample) • Kills single good values between blinks, • Interpolates values between blinks • Fixes blinks at beginning and end of data collection

39. The Pupil Toolkit Graphs of a single subject’s data – all trials showpupiltrials(p) plotpupiltrialmatrix(p); rqplottrialmatrix(p) pupilautocorrgraph(p) threedpupilgraph(p);

40. The Pupil Toolkit Single subject’s data – aggregate waveforms available 1 minute after testing via rqplotaggpupilcondmeans(p) Start with plotpupilcondmeans(p)

41. The Pupil Toolkit .4 .3 .2 .1 GRPPR 0.0 Control - nonpers Control - persrel -.1 Depressed - nonpers -.2 Depressed - persrel 0 288 432 864 144 576 720 1152 1296 1440 1584 1728 2016 2160 2304 2592 2736 1008 1872 2448 Aggregate data – available after some work

42. The Pupil Toolkit Statistics available on every trialvia p.stats=pupiltrialstats(p) • Whole trial: % blinks, blink_in baseline, mean amplitude, slope • Trial Segments: mean baseline amplitude • Relative to peak: peak, peak latency, , Slope post peak • Within a window: mean amplitude, peak, peak latency, max, min, slope • Relative to stimulus: Peak amplitude post stimulus, peak latency post stimulus, • Relative to rt: peak amplitude post rt, slope post rt

43. Example – the face mask study • Question – do different masks cause differential light reflexes following faces? • Data in facemask\data\fmpall.mat

45. pentagons swishes blank Themasks shuffled fractal faceblur

46. 1000 – Aggie – grand average face mask

47. 1000 – Aggie – condition means face mask p=procfacemask(1000)

48. 1001 – Greg – condition means face mask p=procfacemask(1001)

49. Mean of 13 subjects mask face Metrics: flat & low dip peak-trough and similar median to the face. Winner: swishes fmaggmeans(pall,1)

50. The Pupil Toolkit Comparing waveforms • Comparing waveforms is not trivial • We have implemented functions for computing tests at every sample along the waveform. • Unless these comparisons are a-priori I recommend using these only in the context of a group x time or condition x time interaction done on a dimension-reduced dataset. • Controlling type I error • Guthrie & Buchwald’s (1991) technique • Guthrie D, Buchwald JS (1991): Significance testing of difference potentials. Psychophysiology 28:240-244. • Blair & Karniski’s (1993) technique • Blair RC, Karniski W (1993): An alternative method for significance testing of waveform difference potentials. Psychophysiology 30:518-524