# 21 January, 2008
# Control analyses
# MM

library(ggplot2)
library(reshape)
library(lme4)

rm(list=ls())
load("MM.rda")
d$TYP <- relevel(d$TYP, ref="single MS")

# Look at density
# ... main effect: MTC x TYP
m <- ggplot(d, aes(x=mt, group=MTC, colour=MTC)) 
m <- m + geom_density(aes(x=mt)) + facet_grid( CV+EXP ~ TYP)  + scale_x_continuous(limits=c(0,1000), breaks=seq(from=0,to=1000,100))

# Select data: not Exp 5; not zero and not first MS
ix <- which(d$EXP == "VVb" | d$TYP == "no MS" | d$TYP == "first MS" )
d <- d[-ix, ]
d[] <- lapply(d,function(x) x[drop=TRUE])

# Add new variables
# ...  convert condition labels also to factors
d <- cbind(d, colsplit(d$EXP, names=c("Cue", "Target")))
d$Cue <- relevel(d$Cue, 2)
d$Target <- relevel(d$Target, 2)

# Assign relevant subset of data
table(d$TYP)

# VARIOUS CONTROL ANALYSES

# Interesting characteristics of microsaccades; most of them approach significance if between-id is ignored
cor(cbind(d$rt, d$mt,d$mt_r, d$sacAmp, d$sacVPeak, d$mn))

# Look at microsaccade timing, mt & mt_r (r = 0.93)
d$EXP2 <- C(d$EXP, matrix(c(-.5, -.5,   +.5,   +.5,   # (1)  type of cue
                            -.5, +.5,   -.5,   +.5,   # (2)  type of target
                            +.5, -.5,   -.5,   +.5    # (3)  cue x target 
                                                   ), 4, 3), 3)  
contrasts(d$EXP2)

# ... relative to cue onset	
d$mt.c <- scale(d$mt, center=TRUE, scale=FALSE)
d$mt.sq <- d$mt.c*d$mt.c

print(lmer(rt ~  mt.c + mt.sq + CV + mtc.e + EXP2 + CV:mt.c + CV:mt.sq + (1|ID), d), cor=FALSE)  # significant

qplot(mt, rt, data=d,  group=CV, linetype=CV, geom=c("smooth"),  #, method=lm, formula=y ~ poly(x,1) , 
      xlab="Cue-locked Onset of Microsaccade [ms]", ylab="Reaction Time [ms]") +
      scale_x_continuous(breaks=seq(from=0,to=1500,100)) + scale_y_continuous(breaks=seq(from=500,to=800,50)) +
      coord_cartesian(xlim=c(0, 1500), ylim=c(500,800) ) 

qplot(mt, rt, data=d, group=CV, linetype=CV, facets = EXP ~ ., geom=c("smooth"), #method=lm, formula=y ~ poly(x,2) , 
      xlab="Cue-locked Onset of Microsaccade [ms]", ylab="Reaction Time [ms]") +
      scale_x_continuous(breaks=seq(from=0,to=1500,100)) + scale_y_continuous(breaks=seq(from=500,to=800,50)) +
      coord_cartesian(xlim=c(0, 1500), ylim=c(500,800) ) 

# ... relative to target onset
d$mt_r.c <- scale(d$mt_r, center=TRUE, scale=FALSE)
d$mt_r.sq <- d$mt_r.c*d$mt_r.c
print(lmer(rt ~  mt_r.c + mt_r.sq + CV + mtc.e + EXP2 + CV:mt_r.c + CV:mt_r.sq + (1|ID), d), cor=FALSE)  # significant

qplot(mt_r, rt, data=d, colour=CV, geom=c("smooth"), #method=lm, formula=y ~ poly(x,1) , 
      xlab="Target-locked Onset of Microsaccade [ms]", ylab="Reaction Time") +
      scale_x_continuous(breaks=seq(-1500, 0, 100)) + scale_y_continuous(breaks=seq(from=500,to=800,50)) +
      coord_cartesian(xlim=c(-1500, 0), ylim=c(500,800) ) 
     
qplot(mt, rt, data=d, colour=MTC, geom=c("smooth"), #method=lm, formula=y ~ poly(x,1) , 
      xlab="Cue-locked Onset of Microsaccade [ms]", ylab="Reaction Time [ms]") +
      scale_x_continuous(breaks=seq(from=0,to=1500,100)) + scale_y_continuous(breaks=seq(from=500,to=800,50)) +
      coord_cartesian(xlim=c(0, 1500), ylim=c(500,800) ) 
 
qplot(mt_r, rt, data=d, colour=MTC, geom=c("smooth"), #method=lm, formula=y ~ poly(x,1) , 
     xlab="Target-locked Onset of Microsaccade [ms]", ylab="Reaction Time") +
     scale_x_continuous(breaks=seq(-1500, 0, 100)) + scale_y_continuous(breaks=seq(from=500,to=800,50)) +
     coord_cartesian(xlim=c(-1500, 0), ylim=c(500,800) ) 
     
qplot(mt, rt, data=d, colour=EXP, geom=c("smooth"), #method=lm, formula=y ~ poly(x,1) , 
     xlab="Cue-locked Onset of Microsaccade [ms]", ylab="Reaction Time [ms]") +
     scale_x_continuous(breaks=seq(from=0,to=1500,100)) + scale_y_continuous(breaks=seq(from=500,to=800,50)) +
     coord_cartesian(xlim=c(0, 1500), ylim=c(500,800) ) 
     
qplot(mt_r, rt, data=d, colour=EXP, geom=c("smooth"), #method=lm, formula=y ~ poly(x,1) , 
      xlab="Target-locked Onset of Microsaccade [ms]", ylab="Reaction Time") +
      scale_x_continuous(breaks=seq(-1500, 0, 100)) + scale_y_continuous(breaks=seq(from=500,to=800,50)) +
      coord_cartesian(xlim=c(-1500, 0), ylim=c(500,800) ) 

qplot(mt, rt, data=d, facets=MTC ~ CV, colour=EXP, geom=c("smooth"), method=lm, formula=y ~ poly(x,2) , 
      xlab="Cue-locked Onset of Microsaccade [ms]", ylab="Reaction Time") +
      scale_x_continuous(breaks=seq(from=0,to=1500,100)) + scale_y_continuous(breaks=seq(from=500,to=800,50)) +
      coord_cartesian(xlim=c(0, 1500), ylim=c(500,800) ) 
       
lmer(rt ~  mt_r+cv.e+mtc.e+mt_r:cv.e+(1|ID), d)  # significant

# Look at microsaccade magnitude, sacAmp & sacVPeak (r = 0.93)
# ... amplitude 
qplot(sacAmp, rt, data=d, facets=MTC ~ CV, colour=EXP, geom=c("smooth"), method=lm, formula=y ~ poly(x,2) , 
      xlab="Amplitude of Microsaccade", ylab="Reaction Time") +
      scale_y_continuous(breaks=seq(from=400,to=1000,100)) +
      coord_cartesian(xlim=quantile(d$sacAmp, c(.05, .95) ), ylim=c(400,1000) ) 

lmer(rt ~  sacAmp+cv.e+mtc.e+sacAmp:cv.e+(1|ID), d)  # significant

# ... peak velocity 
qplot(sacVPeak, rt, data=d, facets=MTC ~ CV, colour=EXP, geom=c("smooth"), method=lm, formula=y ~ poly(x,2) , 
      xlab="Velocity Peak of Microsaccade", ylab="Reaction Time") +
      scale_x_continuous(breaks=seq(10, 50, 10)) +
      scale_y_continuous(breaks=seq(400, 1000, 100)) +
      coord_cartesian(xlim=quantile(d$sacVPeak, c(.05, .95) ), ylim=c(400, 1000) )

lmer(rt ~  sacVPeak+cv.e+mtc.e+sacVPeak:cv.e+(1|ID), d)  # significant

# Look at number of microsaccades
qplot(mn, rt, data=d, facets=MTC ~ CV, colour=EXP, geom=c("smooth"), method=lm, formula=y ~ poly(x,2) , 
     xlab="Number of Microsaccades in CTI", ylab="Reaction Time") +
     scale_y_continuous(breaks=seq(400, 1000, 100)) +
     coord_cartesian(xlim=c(2,4), ylim=c(400, 1000) )

lmer(rt ~  mn+cv.e+mtc.e+mn:cv.e+(1|ID), d)  # significant