# Combining data from five EXPeriments
# May 14, 2007
# revision: June 13, 2007
# revision: January 21, 2008
# revision: February 11, 2008: correct early onset MS problem

# setwd("/Users/kliegl/documents/projects/EyeMoves/Posner2006.2/")
setwd("/Users/jochen/Projects/work/eyelinkTest/Posner2006_2/analysis/R")
rm(list=ls())

libpath <- '/Users/jochen/Projects/R library/potsdam'
# d.all <- read.table("data/Posner2006_2.sac", header=TRUE)
# d.all <- read.table("../sac/Posner2006_2.sac", header=TRUE)
d.all <- read.table("/Users/jochen/Projects/papers/notYetKnown/analyses/martin/Posner2006_2.sac", header=TRUE)
d.all$CTI <- d.all$tTarEin - d.all$tCueEin

# ... keep only relevant vars
d.all <- d.all[,-c(3, 5, 6, 9:12, 14:16, 20, 36, 39)]

######### EXPerimental codes ##################
d.all$ID <- as.factor(d.all$VPCode %/% 1000)   # integer division
d.all$sess <- d.all$VPCode %% 1000             # modulo

# ... rename tarTyp to EXP
d.all$exp <- ifelse(d.all$tarTyp==1, -1, 1)
d.all$EXP <- as.factor(d.all$tarTyp)
levels(d.all$EXP) <- c("manual RT", "saccadic RT")

d.all$CV <- as.factor(d.all$cueTyp)
levels(d.all$CV) <- c("invalid", "neutral",  "valid")

# cv  = cueTyp
d.all$cv <- d.all$cueTyp
d.all$CV <- as.factor(d.all$cv)
levels(d.all$CV) <- c("invalid", "neutral", "valid")


d.all$SACTYPE <- as.factor(d.all$sacType)
levels(d.all$SACTYPE) <- c("MS", "larger sacc", "manuell Antwortzeile")

# cue direction: left = -1, neutral = 0, right = 1
d.all$cd <- d.all$cv
d.all$cd[d.all$tarPos == -1 & d.all$cv == -1] <-  1  #left target & invalid cue  -> right cue
d.all$cd[d.all$tarPos == -1 & d.all$cv ==  1] <- -1  #left target & valid cue    -> left cue
d.all$cd[d.all$tarPos ==  1 & d.all$cv == -1] <- -1  #right target & invalid cue -> left cue
d.all$cd[d.all$tarPos ==  1 & d.all$cv ==  1] <-  1  #right target & valid cue   -> right cue
d.all$CD <- as.factor(d.all$cd)
levels(d.all$CD) <- c("left",  "neutral", "right")

# target position: left = -1, right = 1
d.all$tpos <- d.all$tarPos
d.all$TPOS <- as.factor(d.all$tpos)
levels(d.all$TPOS) <- c("left",  "right")

d.all$mt <- NA
d.all$mt <- d.all$sacOnset
d.all$mt_r <- NA
# dies war in submission falsch--tTarEin bezieht sich auf trial-Onset,  mt dagegen auf tCueEin ...
d.all$mt_r <- d.all$mt - (d.all$tTarEin - d.all$tCueEin)

d.all$md <- ifelse(d.all$sacAngle1 <= -pi/2 | d.all$sacAngle1 > pi/2, -1, 1)
d.all$MD <- as.factor(d.all$md)
levels(d.all$MD) <- c("left",  "right")

# ... ms congruency with cue:   congruent = 1, incongruent = -1
d.all$mcc <- NA
d.all[which(d.all$md==d.all$cd),  "mcc"] <- 1
d.all[which(d.all$md==-d.all$cd),  "mcc"] <- -1
d.all$MCC <- as.factor(d.all$mcc)
levels(d.all$MCC) <- c("incongruent", "congruent")  # check levels(d$MCC)!

# ... ms congruency with target: congruent = +.5, congruent = -.5
d.all$mtc <- NA
d.all[which(d.all$md==d.all$tpos), "mtc"] <- 1
d.all[which(d.all$md==-d.all$tpos), "mtc"] <- -1
d.all$MTC <- as.factor(d.all$mtc)
levels(d.all$MTC) <- c("incongruent", "congruent")   # check levels(d$MTC)!

################################################
# dependent variables and transformations
d.all$rt <- d.all$reaRT
d.all$lrt <- log(d.all$rt)

d.all$sacBeforeResponseIx <- d.all$sacOnset < (d.all$tTarEin - d.all$tCueEin)
sacMax100MsAfterResponseIx <- d.all$sacOnset <= (d.all$tTarEin - d.all$tCueEin + 100) 

d.all$sacDx <- d.all$sacxOffset - d.all$sacxOnset
d.all$SACDIR[d.all$sacDx > 0] <- 1
d.all$SACDIR[d.all$sacDx < 0] <- -1

######## filter criteria #######################
# prior: blinks after cue, sac > 1°, no response
# now:
# TRIALS WITH MS
# altes Kriterium (Submission),  schließt MS nach Reaktion ein
# d <- d.all[d.all$sacType==1
#          & d.all$sacOnset >= 1     
#          & d.all$sacOnset <= d.all$tTarEin + 100
#          & d.all$reaRT <= 1000
#          & d.all$reaCor == 1, ]
# modifiziertes altes Kriterium (Submission),  schließt MS nach nur bis 100 ms nach Reaktion ein
# d <- d.all[d.all$sacType==1
#          & d.all$sacOnset >= 1     
#          & d.all$sacOnset + d.all$tCueEin <= d.all$tTarEin + 100
#          & d.all$reaRT <= 1000
#          & d.all$reaCor == 1, ]
# neues Kriterium,  schließt MS nach Reaktion aus
d <- d.all[d.all$sacType==1
         & d.all$sacOnset >= 1     
		 & d.all$sacBeforeResponseIx
         & d.all$reaRT <= 1000
         & d.all$reaCor == 1, ]

# TRIALS WITHOUT MS
d0 <- d.all[d.all$numBinSac==1 & d.all$sacType>=2 & d.all$reaRT <= 1000 & d.all$reaCor == 1, ]

#############################################################
# Define exclusive sets: single, first and last of several ms
# source("functions/parseq_fast.R")
# ix <- parseq_fast(d$trial)
source(paste(libpath,  "/parseq.R",  sep=""))
ix <- parseq(d$trial)

d$typ <- NA  # variable containing typ
d0$typ <- 0

d$mn <- NA   # number of ms in trial
d0$mn <- 0

# ... single-ms trials
ix.s <- which(ix[,3]==1)
sngl <- ix[ix.s, 1]
d$typ[sngl] <- -1
d$mn[sngl] <- 1

# ... first ms of trials with several ms
ix.f <- which(ix[,3]!=1)
frst <- ix[ix.f, 1]
frst_of_n <- ix[ix.f, 3]
d$typ[frst] <- 1
d$mn[frst] <- frst_of_n

# ... last ms of trials with several ms
ix.l <- which(ix[,3]!=1) # unique(ix.l == ix.f) == TRUE
last <- ix[ix.l, 2]
last_of_n <- ix[ix.l, 3]
d$typ[last] <- 2
d$mn[last] <- last_of_n

# ... remove MS between first and last MS
ix <- which(is.na(d$mn))
d <- d[-ix, ]

d$TYP <- factor(d$typ)
levels(d$TYP) <- c("single MS", "first MS", "last MS")

d0$TYP <- factor(d0$typ)
levels(d0$TYP) <- "no MS"

ix <- parseq(d$trial)
cor(d[ix[ix.f,1], "SACDIR"], d[ix[ix.f,2], "SACDIR"], use="pairwise.complete.obs")
ixvp1 <- intersect(which(d$ID=="1"), ix[ix.f,1])
ixvp2 <- intersect(which(d$ID=="1"), ix[ix.f,2])
cortab <- as.data.frame(array(dim = c(length(unique(d$ID)), 6)))
names(cortab) <- c("ID", "cor", "mtFirst", "mtLast", "deltaMt", "N")
cortab$ID <- unique(d$ID)
for (vp in unique(d$ID)) {
	ixvp1 <- intersect(which(d$ID==vp), ix[ix.f,1])
	ixvp2 <- intersect(which(d$ID==vp), ix[ix.f,2])
	# correlation between saccade direction (left vs. right, symbolic) of first and second saccade
	cortab[cortab$ID==vp,"cor"] <- cor(d[ixvp1, "SACDIR"], d[ixvp2 , "SACDIR"], use="pairwise.complete.obs")
	cortab[cortab$ID==vp,"mtFirst"] <- mean(d[ixvp1, "mt"])
	cortab[cortab$ID==vp,"mtLast"] <- mean(d[ixvp2, "mt"])
	cortab[cortab$ID==vp,"deltaMt"] <- mean(d[ixvp2, "mt"] - d[ixvp1, "mt"])
	cortab[cortab$ID==vp,"N"] <- length(ixvp1)
}
cortab

# some stats
apply(cortab[,2:6], 2, mean)
apply(cortab[,2:6], 2, range)
apply(cortab[,2:6], 2, sd)/sqrt(18)
mean(tapply(d$mt[d$TYP=="single MS"], d$ID[d$TYP=="single MS"], mean))
sd(tapply(d$mt[d$TYP=="single MS"], d$ID[d$TYP=="single MS"], mean))

##################################################
# generate compatible data frames
d0 <- d0[, c("ID", "trial",  "EXP", "TYP", "CV", "MTC", "MCC", "rt", "lrt", "mt", "mt_r", "sacAmp", "sacVPeak", "sacxOnset", "sacxOffset", "sacDx", "SACDIR", "mn", "exp", "cd", "cv", "mtc", "tpos", "sess", "CTI")]
d <-   d[, c("ID", "trial", "EXP", "TYP", "CV", "MTC", "MCC", "rt", "lrt", "mt", "mt_r", "sacAmp", "sacVPeak", "sacxOnset", "sacxOffset", "sacDx", "SACDIR", "mn", "exp", "cd", "cv", "mtc", "tpos", "sess", "CTI")]
d <- rbind(d, d0)
d[] <- lapply(d,function(x) x[drop=TRUE])

# Save MMVV
# save(d, file="data/Posner2006_2.rda")
save(list=c("d", "d.all", "cortab"), file="./Posner2006_2.rda")

# EXTRAS
library(reshape)
library(ggplot2)
library(lme4)

# cue validity effect
rm(list=ls())
# load("data/Posner2006_2.rda")
load("./Posner2006_2.rda")

d.rs <- melt(d, id=c("ID", "EXP", "TYP", "CV", "MTC", "sess"), measure=c("rt", "lrt", "mn"), subset=TYP !="first MS")
(table.id <- cast(d.rs, ID + sess + EXP + CV ~ variable, function(x) c(M=mean(x), N=length(x))))
summary(aov(rt_M ~ EXP*sess*CV, data=table.id))
summary(aov(mn_M ~ EXP*sess*CV, data=table.id))


# > summary(aov(rt_M ~ EXP*sess*CV, data=table.id))
#              Df Sum Sq Mean Sq F value    Pr(>F)    
# EXP           1   3.32    3.32  6.8774 0.0089272 ** 
# sess          1   6.91    6.91 14.3345 0.0001667 ***
# CV            2   0.03    0.01  0.0308 0.9696245    
# EXP:sess      1   0.06    0.06  0.1202 0.7288845    
# EXP:CV        2   0.05    0.02  0.0475 0.9536254    
# sess:CV       2   0.09    0.04  0.0887 0.9151182    
# EXP:sess:CV   2   0.10    0.05  0.0997 0.9050975    
# Residuals   672 324.17    0.48                      
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

# EXP:sess: deutlichere Verbesserung über Sessions in manueller EXP
# EXP:CV: stärkere CV-Effekte bei manueller als bei sakkadischer EXP (bei sakkadischer kaum ein Effekt invalide vs. neutral)
# sess:CV (n.s.) keine Veränderung der CV Effekte über die Sitzungen


# summary(aov(mn_M ~ EXP*sess*CV, data=table.id))
#              Df  Sum Sq Mean Sq   F value    Pr(>F)    
# EXP           1   3.32    3.32  6.8774 0.0089272 ** 
# sess          1   6.91    6.91 14.3345 0.0001667 ***#
# über die Sitzungen abnehmende Anzahl MS
# mehr MS in Trials mit sakkadischen Antworten !
# alles andere n.s.


# PART 1: EFFECT OF N OF MS IN TRIAL ON RT
# include "No MS" trials, i.e. leave out MTC, use only first of several MS
ix <- which(d$TYP !="last MS")
d.sub <- d[ix, ]
d.sub[] <- lapply(d.sub,function(x) x[drop=TRUE])

d.sub$TYP <- relevel(d.sub$TYP, ref="no MS")
d.sub$TYP2 <- C(d.sub$TYP, matrix(c(-2, +1, +1,  0, -1, 1), 3, 2), 2)
# c1: no ms vs. ms; c2: single vs. first
contrasts(d.sub$TYP2)
d.sub$CV <-  relevel(d.sub$CV, ref="neutral")          # ref is "neutral Cue"
d.sub$MTC <- relevel(d.sub$MTC, ref="congruent")       # ref is "target-congruent MS"
d.sub$sess <- d.sub$sess-3.5                           # ref is middle of Exp


print(m0 <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+TYP2)^2 + (1|ID), data=d.sub), cor=FALSE)
print(m0.man <- lmer(rt ~ poly(sess,2) + CV*TYP2 + (1|ID), data=d.sub, subset=EXP=="manual RT"), cor=FALSE)
print(m0.sac <- lmer(rt ~ poly(sess,2) + CV*TYP2 + (1|ID), data=d.sub, subset=EXP=="saccadic RT"), cor=FALSE)
# Main message: MS for neutral cues slightly increase RT, for invalid cues they slightly reduce cost, and for valid cues they increase benefit!
# also: for saccadic responding,  invalid cues produce no cost (relative to neutral)

# wiederholung der analysen,  beschränkt auf MS, die maximal 400 ms nach Cue auftauchen
subix <- d.sub$TYP=="no MS" | (d.sub$mt>0 & d.sub$mt<400)
print(m0 <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+TYP2)^2 + (1|ID), data=d.sub,  subset=subix), cor=FALSE)
print(m0.man <- lmer(rt ~ poly(sess,2) + CV*TYP2 + (1|ID), data=d.sub, subset=subix & EXP=="manual RT"), cor=FALSE)
print(m0.sac <- lmer(rt ~ poly(sess,2) + CV*TYP2 + (1|ID), data=d.sub, subset=subix & EXP=="saccadic RT"), cor=FALSE)

d.sub.rs <- melt(d.sub, id=c("ID", "EXP", "TYP", "CV", "MTC", "sess"), measure=c("rt", "lrt", "mn"))
(table <- cast(d.sub.rs, EXP + CV + TYP ~ variable, function(x) c(M=mean(x), N=length(x))))
table$CV <- relevel(table$CV, ref="invalid")
levels(table$TYP) <- c("0 MS", "1 MS", "2+ MS")
qplot(x=CV, y=rt_M, data=table, colour=TYP, group=TYP, facets= . ~ EXP, geom=c("point", "line"))
ggsave(file="figures/rt.CV.TYP.EXP.P2.pdf")
# längere RT bei MS--(Abwesenheit von) MS als Vigilanzindikator?

qplot(x=CV, y=lrt_M, data=table, colour=TYP, group=TYP, facets= . ~ EXP, geom=c("point", "line"))

# PART 2: EFFECT OF MTC IN TRIAL ON RT
ix <- which(d$TYP == "no MS")
d.sub12 <- d[-ix, ]
d.sub12[] <- lapply(d.sub12,function(x) x[drop=TRUE])

# include MTC as factor, leave out "No MS" trials; combine for nesting within TYP?, separate LMEs for MTC of first and last MS
print(m1 <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+TYP)^2 + (1|ID), data=d.sub12, subset=TYP !="last MS" ), cor=FALSE)
print(m2 <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+TYP)^2 + (1|ID), data=d.sub12, subset=TYP !="first MS"), cor=FALSE)
# muss identisch sein, da first und last Sakkaden aus denselben Trials sind

print(m1a <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+MTC+TYP)^2 + (1|ID), data=d.sub12, subset=TYP !="last MS" ), cor=FALSE)
print(m2a <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+MTC+TYP)^2 + (1|ID), data=d.sub12, subset=TYP !="first MS"), cor=FALSE)

print(m1b <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+MTC)^3 + (1|ID), data=d.sub12, subset=TYP !="last MS"), cor=FALSE)
print(m2b <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+MTC)^3 + (1|ID), data=d.sub12, subset=TYP !="first MS"), cor=FALSE)

print(m1c <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+MTC+TYP)^3 + (1|ID), data=d.sub12, subset=TYP !="last MS"), cor=FALSE)
print(m2c <- lmer(rt ~ poly(sess,2)*EXP + (EXP+CV+MTC+TYP)^3 + (1|ID), data=d.sub12, subset=TYP !="first MS"), cor=FALSE)


contrasts(d.sub12$TYP) <- cbind(c(-1, 0, 1), c( -1, 2, -1))
contrasts(d.sub12$TYP)
# contrasts(d.sub12$TYP) <- contr.SAS(3)
# contrasts(d.sub12$TYP)
print(m3c.man <- lmer(rt ~ poly(sess,2) + (CV+MTC+TYP)^3 + (1|ID), data=d.sub12, subset= EXP == "manual RT"), cor=FALSE)
print(m3c.sac <- lmer(rt ~ poly(sess,2) + (CV+MTC+TYP)^3 + (1|ID), data=d.sub12, subset= EXP == "saccadic RT"), cor=FALSE)

d.sub12.rs <- melt(d.sub12, id=c("ID", "EXP", "TYP", "CV", "MTC", "sess"), measure=c("rt", "lrt", "mn"))
(table12 <- cast(d.sub12.rs, EXP + CV + TYP + MTC ~ variable, function(x) c(M=mean(x), N=length(x))))
table12$CV <- relevel(table12$CV, ref="invalid")
levels(table12$TYP) <- c("single MS", "first MS", "last MS")
qplot(x=CV, y=rt_M, colour=MTC, group=MTC,  facets= EXP ~ TYP, data=table12, geom = c("point", "line"))
ggsave(file="figures/rt.CV.MTC.TYP.EXP.P2.pdf")

ix <- which(d$TYP == "first MS")
d.sub13 <- d[ix, ]
d.sub13[] <- lapply(d.sub13,function(x) x[drop=TRUE])

d.sub13.rs <- melt(d.sub13, id=c("ID", "EXP", "CV", "MTC", "sess"), measure=c("rt", "lrt", "mn"))
(table13 <- cast(d.sub13.rs, EXP + CV + MTC ~ variable, function(x) c(M=mean(x), N=length(x))))

qplot(x=CV, y=rt_M, colour=MTC, group=MTC, facets= EXP ~ ., data=table13, geom = c("point", "line"))
qplot(x=CV, y=lrt_M, colour=MTC, group=MTC, facets= EXP ~ ., data=table13, geom = c("point", "line"))
qplot(x=CV, y=rt_N, colour=MTC, group=MTC, facets= EXP ~ ., data=table13, geom = c("point", "line"))

#

# Extra for plots
#	   function(x) c(M=mean(x), SE=1.96*sd(x)/sqrt(length(x)), N=length(x) ))
#p <- p +geom_errorbar(aes(max=M+SE, min=M-SE), width=0.1)
#p


# korrolaranalyse: MTC-EFfekte ohne "Wissen" um CV
tapply(d.sub12$rt, list(d.sub12$TYP, d.sub12$MTC, d.sub12$EXP), mean)
ix <- d.sub12$mt < 500
tapply(d.sub12$rt[ix], list(d.sub12$TYP[ix], d.sub12$EXP[ix], d.sub12$MTC[ix]), mean)
print(m1aFirst <- lmer(rt ~ poly(sess,2)*EXP + (EXP*MTC) + (1|ID), data=d.sub12, subset=TYP =="first MS" & mt>200 & mt<400), cor=FALSE)
print(m1aSingle <- lmer(rt ~ poly(sess,2)*EXP + (EXP*MTC) + (1|ID), data=d.sub12, subset=TYP =="single MS" ), cor=FALSE)
print(m2aLast <- lmer(rt ~ poly(sess,2)*EXP + (EXP*MTC) + (1|ID), data=d.sub12, subset=TYP =="last MS" ), cor=FALSE)

# Fazit: für erste MS haben wir Kongruenzeffekt,  für letzte umgekehrten,  für single erst mal keinen,  aber s.u.: letzteres liegt an gegenläufigen Effekten früher und später MS
msplit <- ix
msplit <- NA
singleIx <- d.sub12$TYP=="single MS"
saccIx <- d.sub12$EXP=="saccadic RT"
manIx  <- d.sub12$EXP=="manual RT"
mdsacc <- median(d.sub12[singleIx & saccIx,"rt"])
mdman <- median(d.sub12[singleIx & manIx,"rt"])
msplit[singleIx & saccIx] <- d.sub12[singleIx & saccIx,  "rt"] < mdsacc
msplit[singleIx & manIx] <- d.sub12[singleIx & manIx,  "rt"] < mdman

print(m1aSingleMdSplit <- lmer(rt ~ poly(sess,2)*EXP + (EXP*MTC*as.factor(msplit)) + (1|ID), data=d.sub12, subset=TYP =="single MS" ), cor=FALSE)


# aber: bei Mediansplit zwischen frühen und späten MS finden wir auch single MS einen Effekt: frühe mit MTC,  und späte mit negativer MTC() führen zu schnelleren Reaktionen