-
Notifications
You must be signed in to change notification settings - Fork 0
/
exp1_saccades.Rmd
231 lines (180 loc) · 12.1 KB
/
exp1_saccades.Rmd
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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
---
title: "Exp. 1 - saccade analysis"
author: "Matteo Lisi"
date: "8/1/2018"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.path="./Figs/",message=FALSE,warning=FALSE)
```
This file documents the analysis of saccadic eye movements of Experiment 1. The analysis of saccadic eye movements for the other experiments proceeded through the same steps reported here. (Note that this document may contain some additional analysis that, due to space limitations of the journal's format, were not included in the published article.)
The preprocessing steps are reported exactly in the article, here I report the statistical analyses starting from a dataset containing saccades informations (onsets, offesets, latency, etc.). It requires some packages (all available on CRAN): lme4, optimx, numDeriv, ggplot2, plus some functions available [here](https://github.com/mattelisi/mlisi).
The Matlab code for running the experiment that generated the data, as well as the analysis of raw gaze recordings is available [here](https://github.com/mattelisi/gaussianblobnoise-saccade).
```{r preliminaries, echo=FALSE, results="hide"}
rm(list=ls())
# load custom functions
# library(devtools)
# install_github("mattelisi/mlisi") # if not already installed
library(mlisi)
library(ggplot2) # nicer theme
nice_theme <- theme_bw()+theme(text=element_text(family="Helvetica",size=9),panel.border=element_blank(),strip.background = element_rect(fill="white",color="white",size=0),strip.text=element_text(size=rel(0.8)),panel.grid.major.x=element_blank(),panel.grid.major.y=element_blank(),panel.grid.minor=element_blank(),axis.line.x=element_line(size=.4),axis.line.y=element_line(size=.4),axis.text.x=element_text(size=7,color="black"),axis.text.y=element_text(size=7,color="black"),axis.line=element_line(size=.4), axis.ticks=element_line(color="black"))
# load data
d <- read.table("./data/exp1_saccade.txt",header=T,sep="\t")
# exclude secondary saccades and micro-saccades
d <- d[which(d$sacType==1),] # this include only the first large saccade leaving fixation area in each trial
d <- d[which(d$blink==0),] # remove trials with blinks
# prepare data
#d$tarX <- abs(d$ecc * d$side)
# any(is.na(d$tarX)) # sanity check
d$sacXresp <- abs(d$sacxOffset - d$sacxOnset)
ppd <- 50.1079 # pixel per degree conversion factor
d$sigma <- d$sigma / ppd
d$gain <- d$sacXresp/abs(d$ecc * d$side - d$sacxOnset)
d$sigmaf <- factor(round(d$sigma,digits=1))
#remove "DC component" as in van Opstal & van Gisbergen, 1989, Vision Research, "Scatter in the metrics of saccades and properties of the collicular map"
for(i in 1:nrow(d)){
d$sacxOffset[i] <- d$sacxOffset[i] - mean(d$sacxOnset[which(d$id==d$id[i] & d$session==d$session[i])],na.rm=T)
d$sacxOnset[i] <- d$sacxOnset[i] - mean(d$sacxOnset[which(d$id==d$id[i] & d$session==d$session[i])],na.rm=T)
d$sacyOffset[i] <- d$sacyOffset[i] - mean(d$sacyOnset[which(d$id==d$id[i] & d$session==d$session[i])],na.rm=T)
d$sacyOnset[i] <- d$sacyOnset[i] - mean(d$sacyOnset[which(d$id==d$id[i] & d$session==d$session[i])],na.rm=T)
}
# signed saccade error
d$tarX_s <- (d$ecc * d$side)
d$sacXresp_s <- (d$sacxOffset - d$sacxOnset)
d$sacErr <- sqrt((d$sacxOffset-d$tarX_s)^2 + d$sacyOffset^2)
## rt cut off
# sacOnse is measured with respect to target onset (this is what I want)
# sacRT wrt is with respect to the previous (micro)saccade
d$sacRT <- d$sacOnset
round(mean(d$sacRT<100,na.rm=T) * 100, digits=4) # this compute the percentage of excluded trials
d <- d[which(d$sacRT>=100),]
# # mistakes (i.e. saccades in the wrong direction)
# mean(sign(d$sacXresp_s) != sign(d$tarX_s),na.rm=T) * 100 # total prop. of directional errors
# d$dir_errors <- sign(d$sacXresp_s) != sign(d$tarX_s)
# tapply(d$dir_errors,list(d$sigmaf,d$id),mean,na.rm=T) * 100 # prop. of directional errors for each subject
# absolute error
round(mean(d$sacErr>(d$tarX-2),na.rm=T) * 100, digits=2)
d <- d[which(d$sacErr<(d$tarX-2)),]
# remove unrealistic values of peak velocities (about 0.02%) likely due to blinks
mean(d$sacVPeak>5000,na.rm=T) * 100
d <- d[which(d$sacVPeak<5000),]
# adaptive amplitude filter
filter_group <- paste(d$id, d$sigmaf, d$tarX, sep="_")
out_index <- outfilter(d$sacXresp, filter_group, nsd=3)
round(length(out_index)/nrow(d) * 100,digits=2)
d <- d[-out_index,]
# normalize saccade responses; again as in van Opstal & van Gisbergen, 1989
d$sacXresp <- d$tarX * d$gain
# label sessions according to range of the targets
sess <- d$session
for(i in 1:nrow(d)){
c_i <- with(d[d$vp==d$vp[i] & d$id==d$id[i] & d$session==d$session[i],], (max(tarX)+min(tarX))/2)
sess[i] <- ifelse(c_i>10,"large",ifelse(c_i<10,"small",NA))
}
d$session_number <- d$session
d$session <- sess
# drop unnecessary columns
drops <- c("vp","sameDir","dir_errors","tedfFixOff","gap","sacRT_fd","vpcode","largebef","sacType","sacNumber","sacNumberAfterOnset","sacNumberAfterPrimary","tarx","tary","cxm","cym","tedfFix","tFix","tOn","tOff","tSac","sacAngle1","sacAmp","sacAngle2")
d <- d[ , !(names(d) %in% drops)]
```
ALl the cut-off and preliminary operations are recorded the .Rmd script (commands not shown in the final output). The dataset ready for the analysis is the following. It includes only primary saccades, secondary saccades will be analysed later.
```{r showdataset}
str(d)
```
## Explorative plots
(Some of these plots may not be included in the final version of the article.)
This plot show the cumulative distribution of saccade amplitudes (horizontal).
The vertical dashed line indicate the target distance.
```{r cumulative, echo=FALSE, fig.width=4, fig.height=5}
ggplot(d, aes(x=sacXresp,group=sigmaf,color=sigmaf))+geom_vline(aes(xintercept=tarX),color="black",lty=2,size=0.4)+ stat_ecdf(geom = "step",size=0.6)+facet_grid(tarX~.)+nice_theme+scale_color_manual(values=c("black","dark grey","blue"),name=expression(paste(sigma[blob], " [deg]")))+labs(y="cumulative probability",x="saccade amplitude [deg]")+scale_x_continuous(breaks=seq(5,13,1))+coord_cartesian(xlim=c(5.5,12.5))
```
Average saccade amplitudes
```{r ampli, echo=FALSE, fig.width=2.8, fig.height=2.1}
dag <- aggregate(sacXresp ~ sigmaf+id +tarX , d, mean)
dag2 <- aggregate(sacXresp ~ sigmaf +tarX , dag , mean)
dag2$sacXresp_se <- aggregate(sacXresp ~ sigmaf+tarX, dag, bootMeanSE)$sacXresp
#pdf(file="ampli_exp1.pdf",height=2.1,width=2.8)
ggplot(dag2, aes(x=tarX,y=sacXresp,ymin=sacXresp-sacXresp_se,ymax=sacXresp+sacXresp_se,color=sigmaf,group=sigmaf))+geom_abline(intercept=0,slope=1,lty=2,size=0.2)+geom_errorbar(width=0,alpha=1)+geom_line()+geom_point(size=1.5)+nice_theme+labs(x="target distance [deg]", y="saccade amplitude [deg]")+scale_color_manual(values=c("black","dark grey","blue"),name=expression(paste(sigma[blob], " [deg]")))+scale_x_continuous(limits=c(7.5,12.5))+scale_y_continuous(limits=c(6.8,11.5))
#dev.off()
```
Variability of saccade amplitude gain
```{r spread, echo=FALSE, fig.width=1.6, fig.height=1.6}
dag <- aggregate(gain ~ sigmaf+sigma+id , d, sd)
dag2 <- aggregate(gain ~ sigmaf+sigma , dag , mean)
dag2$gain_se <- aggregate(gain ~ sigmaf+sigma, dag, bootMeanSE)$gain
#pdf(file="gain_SD.pdf",height=1.6,width=1.4)
ggplot(dag2, aes(x=sigma,y=gain,ymin=gain-gain_se,ymax=gain+gain_se,color=sigmaf,group=1))+geom_line(color="black")+geom_errorbar(width=0,alpha=1,color="black")+geom_point(size=2.3)+nice_theme+labs(x=expression(paste(sigma," [deg]")), y="variability [gain SD]")+scale_color_manual(values=c("black","dark grey","blue"),guide=F)+scale_x_continuous(limits=c(0.2,1.61),breaks=c(0.3,0.9,1.5))+scale_y_continuous(limits=c(0.1, 0.18))
#dev.off()
```
Saccade gain as a function of saccadic latency. Horizontal coordinates are computed by splitting the data according to the quartile of the latency distributions of each observer (for each observer I compute the mean latency in each quartile, then average the results across observers).
```{r}
# This assign each observation to a quartile of individual latency distribution
d$RTbin <- NA
d$RTbin_mean <- NA
for(i in unique(d$id)){
for(s_i in unique(d$sigmaf)){
d$RTbin[d$id==i & d$sigmaf==s_i]<-cut(d$sacRT[d$id==i & d$sigmaf==s_i],breaks=quantile(d$sacRT[d$id==i & d$sigmaf==s_i]),labels=c("1st quartile","2nd quartile","3rd quartile","4th quartile"),names=T)
for(q_i in unique(d$RTbin[d$id==i & d$sigmaf==s_i])){
d$RTbin_mean[d$id==i & d$RTbin==q_i & d$sigmaf==s_i] <- mean(d$sacRT[d$id==i & d$RTbin==q_i & d$sigmaf==s_i],na.rm=T)
}
}
}
```
```{r latency-gain, echo=FALSE, fig.width=2.8, fig.height=2.1}
dag <- aggregate(cbind(gain,RTbin_mean) ~ sigmaf+id+ RTbin , d, mean)
dag2 <- aggregate(cbind(gain, RTbin_mean) ~ sigmaf + RTbin , dag , mean)
dag2$se <- aggregate(gain ~ sigmaf + RTbin, dag, bootMeanSE)$gain
dag2$rt_se <- aggregate(RTbin_mean ~ sigmaf + RTbin, dag, bootMeanSE)$RTbin_mean
#pdf(file="ampli-RT_exp1.pdf",height=2.1,width=2.8)
ggplot(dag2, aes(x=RTbin_mean,y=gain,ymin=gain-se,ymax=gain+se,color=sigmaf,group=sigmaf))+geom_hline(aes(yintercept=1),lty=2,size=0.4)+geom_errorbar(width=0,alpha=1)+geom_line()+geom_point(size=1.5)+nice_theme+labs(x="saccade latency [ms]", y="saccade gain")+scale_color_manual(values=c("black","dark grey","blue"),name=expression(paste(sigma[blob], " [deg]")))+scale_y_continuous(breaks=seq(0,1,0.05),limits=c(0.75,1))+scale_x_continuous(breaks=seq(100,600,50))
#dev.off()
```
## Statistical analyses
### Positional uncertainty increases saccade variability and hypometria
Same as reported in article.
Repeated-measures ANOVA on average saccade variability.
```{r}
summary(aov(gain~sigmaf+Error(id/sigmaf), aggregate(gain~sigmaf+sigma+id, d, sd)))
```
Linear mixed-effect model fit on saccade amplitudes.
```{r}
library(lme4)
m1 <- lmer(sacXresp ~ tarX * sigmaf + (tarX + sigmaf |id), d)
summary(m1)
```
Correlation between changes in variability and undershoot bias across observers.
```{r, corr1}
# the change in undershoot bias is measured by the change in the LMM linear coefficient for each participants
gain_change <- fixef(m1)[6]+ranef(m1)$id[,4]
# sanity check: OK
#all(aggregate(gain ~ sigmaf+id, d[d$sigmaf=="1.5",], sd)$id==row.names(ranef(m1)$id))
# the difference in variability is measured by the difference in variance between sigma=1.5 and sigma=0.3
var_increase <- aggregate(gain ~ sigmaf+id, d[d$sigmaf=="1.5",], var)$gain - aggregate(gain ~ sigmaf+id, d[d$sigmaf=="0.3",], var)$gain
# Pearson correlation
cor.test(var_increase, gain_change)
```
Analysis of saccade latency: first split the data according to the quartiles of individual saccade latencies, then run ANOVA test.
```{r}
# ANOVA with gain as dependent variable
dag <- aggregate(gain ~ sigmaf+id + RTbin , d, mean)
dag$RTbin <- factor(dag$RTbin)
summary(aov(gain~sigmaf*RTbin + Error(id/(sigmaf*RTbin)), dag ))
# ANOVA with gain SD as dependent variable
dag <- aggregate(gain ~ sigmaf+id + RTbin , d, sd)
dag$RTbin <- factor(dag$RTbin)
summary(aov(gain~sigmaf*RTbin + Error(id/(sigmaf*RTbin)), dag ))
```
### Saccadic range-effect
Before modelling the saccadic range effect as a function of uncertainty (see script `range_effect.R`) we check whether saccade amplitudes toward the intermediate target (at 10deg eccentricity, which was present in both sessions) are influenced by the range. This test is motivated specifically by recent studies that reported no range effect for dot-like, clearly visible targets. We do that with a two-tailed t-test.
```{r range-test-1}
dag <- aggregate(sacXresp ~ sigmaf+sigma + session +tarX + id, d, mean)
l_sx <- with(dag,sacXresp[sigmaf=="0.3" & session=="large" & tarX==10])
s_sx <- with(dag,sacXresp[sigmaf=="0.3" & session=="small" & tarX==10])
t.test(l_sx-s_sx, var.equal=T)
```
The test is not significant, indicating that a range-effect is not present when uncertainty is small (i.e., when $sigma=0.3$).
We can check whether the range-effect, as measured by the differences across sessions in the amplitudes of saccades made to the 10 deg target, varies with $sigma$. We use a repeated-measures ANOVA.
```{r}
summary(aov(sacXresp ~ sigma*session + Error(id/(sigmaf*session)), dag[dag$tarX==10,]))
```
The ANOVA indicates an interaction between $sigma$ (three levels, 0.3, 0.9 and 1.5 deg) and the session (large vs. small eccentricity range). The same analysis was run on the data from Experiment 3 (not shown here).