#-----------------------------------------------------
# R-Script for:
# "Secondary (micro-)saccades: The influence of primary saccade end point
#  and target eccentricity on the process of postsaccadic fixation"
# by Sven Ohl, Stephan A. Brandt, & Reinhold Kliegl
# Vision Research (in press)
# 26.10.2011
#-----------------------------------------------------

#-----------------------------------------------------
# load packages
library(ggplot2)
library(lme4)
#-----------------------------------------------------

#-----------------------------------------------------
# load data
rm(list=ls())
load("OhlBrandtKliegl_visionresearch.rda")

	#-- Information on data variables --
	# subject:      subject number from 1 to 10
	# trialnum:     the trial number
	# orderevent:   the order of eye movements within one trial
	#		    ordered relative to the primary saccade.
	#		    e.g. "-1" is the last microsaccade before 
	#		    primary saccade onset
	# eccentricity: close vs. distant target condition
	# latency:	    when do eye movements occur with respect to
	#		    the end of the primary saccade
	# amplitude:    amplitude of eye movement
	# srt:	    saccadic reaction time of primary saccade
	# sacerror:     horizontal error of primary saccade in degree
	# direction:    0 = same direction as primary saccade
	#		    1 = opposite direction to primary saccade
	#		    2/ 3 = up and downward
	#		    direction has only been coded correctly for orderevent==1
#-----------------------------------------------------

#-----------------------------------------------------
# Select the first secondary (micro-)saccades
# orderevent == 0 is the primary saccade to the target
# orderevent == 1 is the first secondary (micro-)saccade

d <- data[which(data$orderevent==1),]

# create new variables:
d$over  <- ifelse(d$sacerror>0,d$sacerror,0)		#overshoot
d$under <- ifelse(d$sacerror<0,c(-1)*d$sacerror,0)	#undershoot

# square
d$over2 <- d$over^2
d$under2<- d$under^2

# eccentricity condition
# coding: close target == 0; distant target == 1
d$ecccond <- ifelse(d$eccentricity=="close target",0,1)
#----------------------------------------------------------

#----------------------------------------------------------
# create figure 1: distribution of secondary saccade latency and amplitude
fig1a <- ggplot(d,aes(x=latency)) + geom_histogram()
fig1a <- fig1a + theme_bw(base_size=16)+xlab("latency [in ms]")
fig1b <- ggplot(d,aes(x=amplitude))+geom_histogram()
fig1b <- fig1b + theme_bw(base_size=16) + xlab("amplitude [in deg]")

		#------------------------------------------------------
		# put together figure 1
		windows()
		grid.newpage()
		pushViewport(viewport(layout = grid.layout(2, 2)))
		vplayout <- function(x, y) 
		viewport(layout.pos.row = x, layout.pos.col = y)
		print(fig1a, vp = vplayout(1, 1:2))
		print(fig1b, vp = vplayout(2, 1:2))
		#------------------------------------------------------

# create figure 2: distribution of saccadic error 
# note saccadic error here is only for those events with secondary saccade
windows()
fig2 <- ggplot(d,aes(x=sacerror))+facet_grid(.~eccentricity)+geom_histogram()
fig2 <- fig2 + theme_bw(base_size=16)+ xlab("saccadic error [in deg]")
print(fig2)
#----------------------------------------------------------


###################################################################################################
#
#		START THE STATISTICAL ANALYSIS OF LATENCY; AMPLITUDE AND ORIENTATION
#
###################################################################################################

#----------------------------------------------------------
# Analysis of latency

# linear + quadratic model
summary(res1  <- lmer(latency~ecccond*(over+under+over2+under2)+(1|subject),data=d),REML=FALSE)

	#------------------------------------
	# get predictions of the model "res1"
	estimates <- coef(res1)$subject[1,]
	estimates[1] <- (coef(res1)$subject[,1]- ranef(res1)$subject)[1,1]
	
	e_hor <- seq(0,2,length.out=100)
	
	modelestimate<-function(sacerror,o,u,cdt)
	{
	latency <- estimates[[1]] + cdt*estimates[[2]] + sacerror*o*estimates[[3]] + 	
	sacerror*u*estimates[[4]] + sacerror^2*o*estimates[[5]] + sacerror^2*u*estimates[[6]] +
	sacerror*o*cdt*estimates[[7]] + sacerror*u*cdt*estimates[[8]] + sacerror^2*o*cdt*estimates[[9]] +
  	sacerror^2*u*cdt*estimates[[10]]
	
	if(u==1)
	{
	sacerror <- c(-1)*sacerror
	}
	eccentricity <- rep(cdt,times=length(sacerror))
	
	out <- cbind(latency,sacerror,eccentricity)
	}

	# predict:
	out1 <- modelestimate(e_hor,0,1,0)	# close target undershoot
	out2 <- modelestimate(e_hor,0,1,1)	# distant target undershoot
	out3 <- modelestimate(e_hor,1,0,0)	# close target overshoot
	out4 <- modelestimate(e_hor,1,0,1)	# distant target overshoot

	# put predictions together
	out  <- as.data.frame(rbind(out1,out2,out3,out4))
	out$eccentricity <- as.factor(out$eccentricity)
	levels(out$eccentricity) <- c("p: close target","p: distant target")
	#------------------------------------


# Plot the results of the lmer model along with the smooth of the raw data
figlatency <- ggplot(d,aes(x=sacerror,y=latency))+geom_smooth(aes(linetype=eccentricity,colour=eccentricity),size=1.2)
figlatency <- figlatency + geom_line(aes(x=sacerror,y=latency,linetype=eccentricity,colour=eccentricity),size=1.2,data=out)
figlatency <- figlatency + scale_colour_manual(values=c("red","blue","red","blue"))
figlatency <- figlatency + scale_linetype_manual(values=c(1,1,2,2))
figlatency <- figlatency + xlab("saccadic error [in deg]") + ylab("secondary saccade latency [ms]") + theme_bw(base_size=16)
figlatency <- figlatency + geom_vline(xintercept=0)
figlatency <- figlatency + xlim(quantile(d$sacerror,probs=c(0.025,0.975)))
figlatency <- figlatency + annotate("text", x = -1, y = 320,size=6, label = "undershoot")
figlatency <- figlatency + annotate("text", x = 0.5, y = 320,size=6, label = "overshoot")
figlatency <- figlatency + opts(panel.grid.major=theme_line(colour=NA))
figlatency <- figlatency + opts(panel.grid.minor=theme_line(colour=NA))
#----------------------------------------------------------

#----------------------------------------------------------
# Analysis of amplitude

# linear + quadratic model
summary(res2 <- lmer(amplitude~ecccond*(over+under+over2+under2) + (1|subject),data=d),REML=FALSE)

# plot the result of the lmer results along with smoothed raw data

	#------------------------------------
	# get predictions of the model "res2"
	estimates <- coef(res2)$subject[1,]
	estimates[1] <- (coef(res2)$subject[,1]	- ranef(res2)$subject)[1,1]
	
	e_hor <- seq(0,2,length.out=100)
	
	modelestimate<-function(sacerror,o,u,cdt)
	{
	amplitude <- estimates[[1]] + cdt*estimates[[2]] + sacerror*o*estimates[[3]] + 	
	sacerror*u*estimates[[4]] + sacerror^2*o*estimates[[5]] + sacerror^2*u*estimates[[6]] +
	sacerror*o*cdt*estimates[[7]] + sacerror*u*cdt*estimates[[8]] + sacerror^2*o*cdt*estimates[[9]] +
  	sacerror^2*u*cdt*estimates[[10]]
	
	if(u==1)
	{
	sacerror <- c(-1)*sacerror
	}
	eccentricity <- rep(cdt,times=length(sacerror))
	
	out <- cbind(amplitude,sacerror,eccentricity)
	}

	# predict:
	out1 <- modelestimate(e_hor,0,1,0)	# close target undershoot
	out2 <- modelestimate(e_hor,0,1,1)	# distant target undershoot
	out3 <- modelestimate(e_hor,1,0,0)	# close target overshoot
	out4 <- modelestimate(e_hor,1,0,1)	# distant target overshoot

	# put predictions together
	out  <- as.data.frame(rbind(out1,out2,out3,out4))
	out$eccentricity <- as.factor(out$eccentricity)
	levels(out$eccentricity) <- c("p: close target","p: distant target")
	#------------------------------------


figamp <- ggplot(d,aes(x=sacerror,y=amplitude,group=eccentricity))+geom_smooth(aes(linetype=eccentricity,colour=eccentricity),size=1.2)
figamp <- figamp + geom_line(aes(x=sacerror,y=amplitude,colour=eccentricity,linetype=eccentricity),size=1.2,data=out)
figamp <- figamp + scale_colour_manual(values=c("red","blue","red","blue"))
figamp <- figamp + scale_linetype_manual(values=c(1,1,2,2))
figamp <- figamp + xlab("saccadic error [in deg]") + ylab("secondary saccade amplitude [in deg]")
figamp <- figamp + theme_bw(base_size=16)
figamp <- figamp + geom_vline(xintercept=0)
figamp <- figamp + annotate("text", x = -1.5, y = 3,size=6, label = "undershoot")
figamp <- figamp + annotate("text", x = 1.5, y = 3,size=6, label = "overshoot")
figamp <- figamp + opts(panel.grid.major=theme_line(colour=NA))
figamp <- figamp + opts(panel.grid.minor=theme_line(colour=NA))
#----------------------------------------------------------


#----------------------------------------------------------
# Analysis of orientation

# coding for d$direction
# 0 = same direction as primary saccade
# 1 = opposite direction to primary saccade
# 2 = up
# 3 = down

table(d$direction)

# remove up and down movements
idx <- which(d$direction>=2)
d_logreg   <- d[-idx,]
table(d_logreg$direction)


# NESTED MODEL glmer
summary(res3 <- glmer(direction ~ ecccond*(over + under)+(1|subject),data=d_logreg,family=binomial))



	#------------------------------------
	# get predictions of the model "res6"
	estimates <- coef(res3)$subject[1,]
	estimates[1] <- (coef(res3)$subject[,1]	- ranef(res3)$subject)[1,1]
	
	e_hor <- seq(0,2,length.out=100)
	
	modelestimate<-function(sacerror,o,u,cdt)
	{
	direction <- estimates[[1]] + cdt*estimates[[2]] + sacerror*o*estimates[[3]] + 	
	sacerror*u*estimates[[4]] + sacerror*cdt*o*estimates[[5]] + sacerror*cdt*u*estimates[[6]]


	# logistic regression ->  1/(1+exp(-modelestimate))
	direction <- 1/(1+exp((-1)*direction))

	if(u==1)
	{
	sacerror <- c(-1)*sacerror
	}
	eccentricity <- rep(cdt,times=length(sacerror))
	
	out <- cbind(direction,sacerror,eccentricity)
	}

	# predict:
	out1 <- modelestimate(e_hor,0,1,0)	# close target undershoot
	out2 <- modelestimate(e_hor,0,1,1)	# distant target undershoot
	out3 <- modelestimate(e_hor,1,0,0)	# close target overshoot
	out4 <- modelestimate(e_hor,1,0,1)	# distant target overshoot

	# put predictions together
	out  <- as.data.frame(rbind(out1,out2,out3,out4))
	out$eccentricity <- as.factor(out$eccentricity)
	levels(out$eccentricity) <- c("close target","distant target")
	#------------------------------------


figdirection <- ggplot(d_logreg,aes(x=sacerror,y=direction))
figdirection <- figdirection + geom_line(aes(x=sacerror,y=direction,colour=eccentricity),size=1.2,data=out)
figdirection <- figdirection + geom_point(aes(x=sacerror,y=direction,colour=eccentricity),alpha=0.8,size=0.5,position=position_jitter(height=0.05))
figdirection <- figdirection + scale_colour_manual(values=c("red","blue"))
figdirection <- figdirection + xlab("saccadic error [in deg]") + ylab("p(secondary saccade orientation)") + theme_bw(base_size=16)
figdirection <- figdirection + geom_vline(xintercept=0)
figdirection <- figdirection + annotate("text", x = -1, y = 0.9,size=6, label = "undershoot")
figdirection <- figdirection + annotate("text", x = 0.5, y = 0.9,size=6, label = "overshoot")
figdirection <- figdirection + opts(legend.position=c(0.3,0.45))
figdirection <- figdirection + opts(panel.grid.major=theme_line(colour=NA))
figdirection <- figdirection + opts(panel.grid.minor=theme_line(colour=NA))
#----------------------------------------------------------