Transcript PICAN MATLAB training - University of Pittsburgh Department of

Analysis of Pupillary Data
Society for Psychophysiological Research Workshop, Boston, September 14, 2011
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
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
Pupillary Response at Visual Threshold
G. Hakerem & S. Sutton, Nature, 1966
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).
Average Pupillary Light Reaction
Start
Pupil Diamter (mm)
7
6.5
Amplitude of
Light Reaction
End
6
Light
3400
3000
2600
2200
1800
1400
1000
600
200
-200
5.5
Time (ms)
(for description of multiple measures, see Steinhauer et al., Psychophysiology, 1992)
Am plitude
Baseline
Start of Contraction
Pupil Diameter (mm)
6.1
5.9
5.7
5.5
Maximum
Velocity
Minimum
Diameter
Redilation
5.3
5.1
Time
104
424
744
1064
1384
1704
2024
2344
2664
2344
2664
2984
3304
3624
Tim e (m sec)
Diameter Diff.
(mm)
Velocity
0.06
0.04
0.02
0.00
-0.02
Time
104
424
744
1064
1384
1704
2024
2984
3304
3624
Tim e (m sec)
Velocity Diff
(mm)
Acceleration
0.010
0.005
0.000
-0.005
Time
104
424
744
1064
1384
1704
2024
Time (msec)
2344
2664
2984
3304
3624
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.
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.
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
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
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
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.
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 average of 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
Digit Span Task
Sevens
Control
0.3
Schiz
0.2
0.1
0
0
1 2 3 4 5
6 7 8
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Tens
Control
0.3
Schiz
0.2
0.1
0
0 1 2
3 4 5
6 7 8
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Often, activity over much longer durations may need to be analyzed
Probability Effects, Auditory Counting Task
.33
.67
1.00
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
-0.05
-0.1
-0.15
0
0.5
1
1.5
2
2.5
Probability Effects, Auditory Choice Reaction Task
0.7
.33
.67
1.00
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
0
0.5
1
1.5
2
2.5
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
Recording in Light
Recording in Darkness
0.14
0.14
0.12
0.1
Sympathetic
0.08
Parasympathetic
0.06
Summation
0.04
0.02
Change in Diameter (mm)
Change in Diameter (mm)
0.12
0.1
Sympathetic
0.08
Parasympathetic
0.06
Summation
0.04
0.02
0
0
300
600
900
1200
Time (msec)
Steinhauer & Hakerem, 1992
1500
300
600
900
Time (msec)
1200
1500
Pupil: Time Frequency Analysis
Pupillary Sleepiness Test
Franzen et al (in prep)
seconds
Wavelet analysis of
Oscillatory Activity
S1a
Pupil diameter (mm)
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.
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).
7
6
2
S1b
Change in pupil diameter (mm)
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.
8
5
0
4
6
Time (minutes)
8
10
12
8
10
12
8
10
12
Detrended Pupil Diameter
2
1
0
-1
-2
0
2
4
S1c
6
Time (minutes)
Running PUI
250
PUI (power, 0-0.6 Hz)
Initial Data Reduction: The pupil was
first corrected for blinks (S1a), then
smoothed and detrended (correcting
for baseline, S1b).
Absolute Pupil Diameter
9
200
150
100
50
0
0
2
4
6
Time (minutes)
Detrended Pupil Diameter
Change in pupil diameter (mm)
Change in pupil diameter (mm)
Detrended Pupil Diameter
2
1
0
-1
-2
0
2
4
6
Time (minutes)
8
10
12
2
1
0
-1
-2
0
2
4
6
Time (minutes)
8
10
8
10
12
Running PUI
Running PUI
250
PUI (power, 0-0.6 Hz)
200
150
100
50
0
0
2
4
6
Time (minutes)
8
10
200
150
100
50
12
0
0
2
4
6
Time (minutes)
Power Derived from Wavelet Analysis
lnPUI by Subject Group
1
lnPUI (0-0.6 Hz)
PUI (power, 0-0.6 Hz)
250
0.8
0.6
0.4
Control
0.2
Alcohol
Schizophrenic
0
Depressed
-0.2
-0.4
1
2
3
Total
First Min
Last 5
12
MATLAB
IMPLEMENTATIONS
Thanks to Greg Siegle, Ph.D.
Depts. Of Psychiatry and Psychology
University of Pittsburgh
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
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);
1
0.5
0
-0.5
1
-1
0
10
20
30
40
50
0.5
0
-0.5
-1
-1
-0.5
0
0.5
1
1
0.5
0
1
-0.5
0.5
-1
0
10
20
30
40
50
0
-0.5
-1
1
1.2
1.4
1.6
1.8
2
2
1
0
10
– hist(pts)
-1
8
6
4
2
0
-1
-0.5
0
0.5
1
• figure; figure(2);
• subplot(2,2,1); plot(x,y);
-2
-2
0
2
4
6
8
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)
Operations in functions
if a>2
b=4;
else
b=5;
end
For ct=1:5
b=b+1;
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
while ct<5
b=b.^2;
ct=ct+1;
end
The
The Pupil Toolkit
Pupil
Toolkit
Greg Siegle, Ph.D.
University of Pittsburgh,
School of Medicine
Goals
The
Pupil
Toolkit
• 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
Standard flow for 1subject’s pupil data
The
Pupil
Toolkit
• 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
Single subject diagnostic quality control graphs
available 1 minute after testing
The
Pupil
Toolkit
via p=readiscan(‘1002kvid.isc’);
1002KVID.isc diagnostic plot
8
pupil
rblnk
serial
detrend
blinks
7
6
iscan units
5
4
3
2
1
0
0.5
1
1.5
2
tick
2.5
3
4
x 10
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)
>> p=procsilkkvidexample(1002)
Raw data
After trial segmentation
1002KVID.isc diagnostic plot
>> p=procsilkkvidexample(1002)
found 34280 records in data
dropped 1 of 51 trials
8
0
0.5
1
1.5
2
2.5
4
2
1
0
0
5 10
t
5 10
t
18
2
1
0
0
5 10
t
mm
mm
0
2
1
0
mm
mm
mm
5 10
t
17
5 10
t
13
0
5 10
t
14
0
5 10
t
19
2
1
0
2
1
0
0
0
10
20
30
40
50
0
10
20
30
40
50
0
10
20
30
40
50
0.3
0.5
1
0.5
0
5 10
t
10
2
1
0
0
5 10
t
15
0
5 10
t
20
0
5 10
t
2
1
0
2
1
0
positive
negative
neutral
0.4
1
0
5 10
t
p=
0
0.5
dilation
manual
suspect
0
2
1
0
0
5 10
t
9
2
1
0
0.5
0
bad rt
5 10
t
12
2
1
0
0
2
1
0
x 10
1
blink at baseline> 50% blinks
5 10
t
0
5 10
t
8
5
mm
mm
0
3
tick
2
1
0
mm
mm
1
0
2
1
0
mm
5 10
t
16
2
1
0
4
mm
0
2
1
0
2
1
0
5 10
t
7
mm
mm
2
0
mm
5 10
t
11
3
mm
0
3
2
1
0
5 10
t
6
2
1
0
4
mm
iscan units
mm
mm
0
5
2
mm
6
2
1
0
mm
1
pupil
rblnk
serial
detrend
blinks
7
0.2
0.1
1
0
0.5
0
0
10
20
30
40
50
60
0
10
20
30
40
50
60
1
0.5
0
-0.1
-0.2
0
2
4
6
seconds
Valence Identification Task
8
10
12
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]
The
Pupil
Toolkit
p=procvfg(5107,'numorder',2);
(sort 3, 4 or 5 numbers)
510702.502 diagnostic plot
5
manual
pupil
rblnk
event
4.5
suspect
4
3.5
2.5
1.5
1
0.5
0
-1
10
0
10
blink at baseline
> 50% blinks
2
1
5
10
15
20
25
30
35
40
5
10
15
20
25
30
35
40
5
10
15
20
25
30
35
40
5
10
15
20
25
30
35
40
5
10
15
20
25
30
35
40
0.5
bad rt
3
0
-1
10
0
-1
10
0.5
0
0.5
1
1.5
2
tick
2.5
3
3.5
4
0
0
4
x 10
harmonic mean rts for numorder-5107-2.txt
1500
510702.502
1400
0.3
ms
3
4
5
0.35
1300
1200
1100
0.25
1000
0.2
0
1
2
# stimuli
3
4
0.15
1400
0.1
1300
ms
pupil diameter (mm)
mm diameter
The
Pupil
Toolkit
Digit Sorting Task
0.05
0
1200
1100
-0.05
1000
0.5
-0.1
0
2
4
6
8
10
seconds
12
14
16
18
1
1.5
2
Absent/Present
2.5
Reading data from
different pupilometers
The
Pupil
Toolkit
• ISCAN
– readiscan
– readiscan05text
– readiscanbehavobstext
• ASL – to use must first install ASL matlab drivers
–
–
–
–
–
readasl2006
readasl2007
readasl2008
readasltext
readasltextlunalab
How the blink elimination
routine works
The
Pupil
Toolkit
data =
stublinks(data,graphics,lrtask,manualblink
s,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
Graphs of a single subject’s data – all trials
showpupiltrials(p)
2
0
-1
1011 20
t
0
mm
1
0
-1
1012 20
t
1016 20
t
0
0
1017 20
t
0
-1
0
10 20
t
0
10 20
t
109 20
t
0
0
1
0
0
1014 20
t
0
10 20
t
0
1015 20
t
0
1020 20
t
0
10 20
t
0
1
0
-1
0
1
1019 20
t
0
-1
0
1010 20
t
-1
1
1018 20
t
0
1
-1
-1
0
mm
mm
1013 20
t
-1
1
mm
0
-1
0
1
-1
0
mm
0
0
1
-1
1
108 20
t
-1
mm
0
1
mm
0
1
0
-1
mm
0
-1
107 20
t
rqplottrialmatrix(p)
1
0
-1
mm
0
1
mm
mm
106 20
t
-1
mm
0
1
0
plotpupiltrialmatrix(p);
5
mm
-1
1
mm
1
0
mm
-1
4
mm
0
mm
1
mm
mm
1
3
1
mm
1
The
Pupil
Toolkit
0
-1
0
10 20
t
threedpupilgraph(p);
pupilautocorrgraph(p)
Single subject’s data – aggregate waveforms available
1 minute after testing via rqplotaggpupilcondmeans(p)
Start with plotpupilcondmeans(p)
Stimulus Aligned Data - All Conditions
The
Pupil
Toolkit
Stimulus Aligned Data - Collapsed on Personal Relevance
+
N
+p
-p
Np
0.8
0.6
+
N
0.6
0.4
mm
mm
0.4
0.8
0.2
0.2
0
0
-0.2
-0.2
0
5
10
t (s)
15
20
0
Mean +/- St Dev
5
10
t (s)
15
20
Response Aligned Data - Collapsed on Personal Relevance
0.8
+
N
0.6
mm
1
0.5
0
5
mm
0.4
6
0.2
4
0
5
0
3
10
15
t (s)
2
20
1
valence
-0.2
-5
0
5
10
t (s)
15
20
The
Pupil
Toolkit
Aggregate data – available after some work
.4
.3
.2
.1
GRPPR
0.0
Control - nonpers
Control - persrel
-.1
Depressed - nonpers
-.2
Depressed - persrel
Statistics available on every trial
via p.stats=pupiltrialstats(p)
The
Pupil
Toolkit
• 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
Example – the face mask study
• Question – do different masks cause
differential light reflexes following faces?
• Data in facemask\data\fmpall.mat
blank
pentagons
swishes
faceblur
shuffled
fractal
The
masks
1000 – Aggie – grand average
face
mask
0.1
0.05
0
-0.05
-0.1
0
100
200
300
400
500
600
700
0.6
time (s)
1000 – Aggie –
condition means
rts
0.65
0.55
0.5
0.45
0.4
0
1
2
3
4
5
6
7
face
mask
0.3
1.
2.
3.
4.
5.
6.
0.2
dilation
0.1
0
-0.1
-0.2
-0.3
-0.4
0
p=procfacemask(1000)
2
4
6
seconds
8
10
12
pentagons
swishes
faceblur
shuffled
blank
fractal
1001 – Greg –
condition means
rts
0.8
time (s)
0.6
0.4
0.2
0
0
1
2
3
4
5
6
face
mask
0.25
pentagons
swishes
faceblur
shuffled
blank
fractal
0.2
0.15
dilation
0.1
0.05
0
-0.05
-0.1
-0.15
-0.2
0
p=procfacemask(1001)
2
4
6
seconds
8
10
12
7
Mean of 13 subjects
mask
face
0.25
1.
2.
3.
4.
5.
6.
0.2
0.15
0.1
pentagons
swishes
faceblur
shuffled
blank
fractal
Metrics:
flat & low dip
peak-trough
and similar median
to the face.
Winner:
swishes
dilation
0.05
0
-0.05
-0.1
-0.15
-0.2
-0.25
0
2
4
fmaggmeans(pall,1)
6
seconds
8
10
12
Comparing waveforms
The
Pupil
Toolkit
• 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 dimensionreduced 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:518524
T-tests at every sample, marking intervals
significantly long enough to care about via
Guthrie & Buchwald (1991)
What’s the emotion?
UGLY
The
Pupil
Toolkit
Put the digits
in numerical order
7 4 3 1 5
Never Depressed (41)
Unmedicated
Depressed (47)
0.3
0.2
mm
0.1
0
5
10
15
20
25
seconds
Regions of significant differences. Yellow is p<threshold you
pass as an argument (p<0.1 by default) and red is less than p<.05. The
black line is intervals longer than the patch length you pass in.
Pupil toolkit functions to
implement Guthrie &
Buchwald’s (1991) technique
The
Pupil
Toolkit
• gutautocorr
– gives the autocorrelation (acorr) of waveforms in matrix X (subjects x
samples) after removing k principal components
– [acorr,kmin]=gutautocorr(X,k)
• gsgutsims
– Runs simulations for Guthrie and Buchwald’s (1991) technique to yield
minimum # of consecutive tests necessary for a difference to be
considered significant at p<.05
– [minlen]=gsgutsims(N,T,sig,auto,numsims)
• N=# subs, T= sampling interval, sig = target waveform-wise significance
(usually 0.05), auto = autocorrelation in the data, numsims = # of simulations
(default = 1000)
• gsgutsimsbetween(Ng1,Ng2,T,sig,ro,numsims)
Pupil toolkit functions to
implement Blair & Karnitski’s
technique
The
Pupil
Toolkit
• getblairkarniskitmax
– [tmaxthresh,p05tmax] =
getblairkarniskitmax(data,group,sigthresh)
– generates all permutations of data to conditions and
for each permutation does a t-test at each time-point.
– We then select the tmax for which 95% (or other
threshold) of the permutations are rejected, such that
95% of the permutations have NO significant t-tests
– We then apply that threshold to the successive t-tests
in our waveform of interest.
Functions to compare waveforms:
Within Subjects
The
Pupil
Toolkit
• diffwavgraph
– Contrasts 2 conditions via t-tests at each sample
– [s,h]=diffwavgraph(wavcond1,wavcond2,samprate,res
amprate,alpha,outliers,patchlen,bw,pscale,wavcond3,
pscalemag,xax,linewidth)
– NOTE: This is the only function in the set which is
well documented… If you get this you’ll get the
rest…
• condwavgraph
– contrasts all conditions within subjects via anova at
each sample
– [s,h]=condwavgraph(condwavs,samprate,resamprate,
alpha,outliers,patchlen,bw,pscale,xax)
diffwavgraph help file
The
Pupil
Toolkit
•
•
•
•
graphs cond1 v. cond2
expects 2 matrices, each with subjects in rows and wavs in columns
and significance of difference for each time point
usage: s=diffwavgraph(wavcond1,wavcond2,samprate,resamprate,alpha,outliers,patchlen,bw,pscale,wavcond3,pscalemag,xax,linewidth)
•
•
•
•
•
•
wavecond1: this is the matrix of N rows for condition 1
wavecond2: this is the matrix of N rows for condition 2
samprate: the sampling rate, in Hz
resamprate: The rate at which to resample the data - usually the same as samprate
–
You can leave the samprate and resamprate the same and the routine will run fastest, and in the most principled way. The reason to consider resampling is to
decrease the autocorrelation in the data. The more you resample, the "rougher" the data will be, and thus the less points you'll need in a row to get
significance. So, in case you play with resampling in getting the autocorrelation, I let you throw that in as a parameter...
alpha: threshold for significance, usually set to .05 or .1
–
And like Guthrie and Buchwald, I like the .1 threshold. That said, if regions I like are not coming out, I often like to see what the actual significance of patches
"would be" so that I can know whether it's a power issue. Towards that end, I'll often play, in the privacy of my darkened office with the door locked, with
thresholds of .3 or .5...
outliers:
This is an easy way to recompute patches with specific people taken out, just to see if things change. By default it should be a vector with N
rows and one column of all zeros. If you put ones on any row, those people are not counted in computing mean waveforms or significance tests.
patchlen: the length of consecutive data points required for sig
–
note: patchlen does NOT account for resamprate. So, even if you downsample to 1hz, a patchlen of 17 refers to 17 points in a row in the original sampled
space. Thus you must change it yourself in the calling routine...
bw: whether or not graphs should be in black and white
–
For display on the screen, set bw=0 and graphs will appear in color. For publications in which color is costly, set bw=1 and it will make your graphs in black
and white, with dotted lines as necessary.
pscale:
This should be zero for most applications. That will set the significance bars to a uniform height. If pscalemag is not zero, the significance bars are of the height
of the p-value, scaled by pscale.
wavcond3:
This is an optional 3rd condition, which is plotted but not included in tests. Set it to zero if there is no third condition.
pscalemag:
–
This is how large the bars for significance should appear under the x axis (if pscale =0). So if the y axis goes from -10 to 10, you might make pscale = 1. But if
the y axis goes from -.1 to .1 you might make pscale = 0.02. making pscalemag negative puts the significance bars below the x axis
xax:
–
This is the units you want on the x-axis. By default it puts the x-axis in seconds. But if you want it in ticks, pass a vector which counts from 1:wavelen.
linewidth:
This is how wide the lines are on the plots. 0.5 by default. If you want to thicken then up, e.g., for a poster, pass in values > 0.5
•
•
•
•
•
Function by Greg Siegle, Ph.D.
cite as Siegle, G. J. (2003) The Pupil Toolkit, Available directly from the author, as used in, e.g.,
Siegle GJ, Steinhauer SR, Carter CS, Ramel W, Thase ME (2003): Do the seconds turn into hours? Relationships
between sustained pupil dilation in response to emotional information and self-reported rumination.
Cognitive Therapy and Research 27:365-382.
•
•
•
•
•
•
•
•
•
•
•
•
•
Are the conditions different?
0.1
0.05
dilation
0
-0.05
1.
2.
3.
4.
5.
6.
-0.1
-0.15
-0.2
-0.25
0
2
4
6
seconds
8
pentagons
swishes
faceblur
shuffled
blank
fractal
10
0.57 to 6.08: F(5,8)=12.57, p=0.00
fmaggmeans(pall,2)
12
Toolboxes
The
Pupil
Toolkit
• Toolboxes provide lots of functionality in a
specific domain
• Toolboxes you might be interested in
– Database
– Statistics
– Signal processing
– Wavelet
Database toolbox
• Database toolbox (from
Start)
– Works w/ all odbc datasources
and builds sql queries
– Click on generate-m-file under
query
– Suggests making query with
query variables as input
variable
The
Pupil
Toolkit
Stats toolbox
•
•
•
Distribution fitting tool
Can have nominal data with nominal class via
nominal command – reduces space greatly and
lets you select, e.g., by category, and make
plots by category.
W/ curve fitting (cftool) there’s a goodness of fit
metric and can get quantitative measures of fit
for hyp testing
– And can use custom equation by selecting new
fit->type of fit->custom
– Fit-options = linear least squares or other as I
choose
– Can plot confidence interval around curve
function
– Can output fit params via save-to-workspace
•
Statsregressionnonlinearmixed effects
The
Pupil
Toolkit
Pupil Dilation as a continuous
measure of cognitive load
proportion of
maximum
dilation
Pupil dilation
outside fMRI
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
-0.05
-0.1
Pupil dilation
during fMRI
0.3
0.25
3
4
5
0.2
0.15
0.1
0.05
0
0
2
4
6
8
10 12 14 16 18
0
2
4
6
8
10 12 14 16 18
seconds
seconds
# digits
3
4
5
0.2
0.1
0
4 8 12 16
seconds
0.4
0.3
0.2
0.1
0
% change
0
0.2
0.1
0
0
0
4
8 12 16
seconds
Siegle et al, 2003, Neuroimage
4
8
12 16
seconds
eye-tracking colored by pupil
Courtesy of Greg Siegle, PhD.
To Play or Download these
presentations:
http://www.wpic.pitt.edu/research/biometrics
Go to the Lab Publications link
Click on the Biometrics Archives
Look under SPR 2011 Workshop
Select Presentations or References
http://www.wpic.pitt.edu/research/biometrics
(select Lab Publications)
Download Presentations or select Reference List
For additional information:
http://www.wpic.pitt.edu/research/biometrics
[email protected]
412-954-5366