# RTs library(lme4) library(reshape) library(ggplot2) rm(list=ls()) parseq_fast<-function(a) { # PARSEQ parses a column vector according to its incremental sequence of unchanging # values. Use for identifying sequential blocks of data. # b1: begin b2: end b3: length of data block. # J. Laubrock, 04-10-2004 # Transformation von matlab in R by Sven Ohl 12-03-2006 # speed-up by jochen laubrock, 05-30-06 la = length(a) dim(a) <- c(la,1) diff <- a[2:la,]-a[1:la-1,] tmp <- which(diff!=0) b <- cbind(c(1,tmp+1), c(tmp, la)) b <- cbind(b, b[,2]-b[,1]+1) } load("MM.rda") # Select data: not Exp 5; not "first MS" ix <- which(d$EXP == "VVb" | d$TYP == "first MS") d <- d[-ix, ] d[] <- lapply(d,function(x) x[drop=TRUE]) d$MN <- relevel(d$TYP, 3) levels(d$MN) <- c("0", "1", "2+") # 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) ################################# d.mn <- melt(d, id=c("ID", "EXP", "Cue", "Target", "MN", "CV"), measure=c("mn")) cast(d.mn, EXP ~ ., mean) # Frequency of microsaccade numbers table(d$mn)/nrow(d) ################################## # LME ################################## # Generate means and figure based on all trials, allow for missing design cells # Table 1 (used in manuscript) d.rt <- melt(d, id=c("ID", "EXP", "Cue", "Target", "MN", "CV", "mn"), measure=c("rt")) cast(d.rt, mn ~., function(x) c( N=length(x), p=round(length(x)/10043,digits=3))) tab0 <- cast(d.rt, EXP + MN ~ CV, function(x) c(M=round(mean(x)), SE= round(sd(x)/sqrt(length(x))), N=length(x))) # Figure 2 # ... set up data tab <- cast(d.rt, EXP + MN + CV ~ ., function(x) c(M=round(mean(x)), SE=round(sd(x)/sqrt(length(x))), N=length(x), lower=round(mean(x)-(sd(x)/sqrt(length(x)))), upper=round(mean(x)+(sd(x)/sqrt(length(x)))) )) # names(tab)[[4]] <- "MeanRT" tab # ... plot p <- qplot(MN, MeanRT, data=tab, group=CV, colour=CV, shape=CV, facets = . ~ EXP, xlab="Number of MS", ylab = "Reaction Time [ms]") p <- p + geom_point(size=3) + geom_line(size=1) + geom_errorbar(aes(max=upper, min=lower, width=0.1), size=1) p ################################# # LMEs # ... (0) Overall LME on number of microsaccades and RTs d.mn <- melt(d, id=c("ID", "EXP", "Cue", "Target", "CV", "cv.e"), measure=c("mn")) d.mn.ID <- melt(cast(d.mn, ID+Cue+Target+CV ~ ., mean )) cast(d.mn, EXP ~ ., mean, margins=TRUE ) cast(d.mn, Cue+Target ~ ., mean, margins=TRUE) # ... ... LME: basically VV: 1.5, VA and AV: 2.0, AA: 2.5, that is two main effects print(m0 <- lmer(mn ~ EXP*CV + (1|ID), d), cor=FALSE) source("contrasts.R") print(m0.a <- lmer(mn ~ EXP0*CV + (1|ID), d), cor=FALSE) # should be identical with m0 print(m1 <- lmer(mn ~ EXP1*CV + (1|ID), d), cor=FALSE) print(m2 <- lmer(mn ~ EXP2*cv.e + (1|ID), d), cor=FALSE) # Use the effect coding for manuscript # ... ... ID-level ANOVA m0.id <- coef(m0)[[1]][1] names(d.mn.ID)[[5]] <- "NofMS" summary(aov(NofMS ~ Cue*Target*CV+Error(ID/CV), d.mn.ID)) # ... (1) CV crossed with repeated contrasts of MN d$design1 <- factor(paste(d$MN, d$CV, sep="")) d$design1 <- C(d$design1, matrix(c(-.5, +.5, -.5, +.5, -.5, +.5, # (1) CV -.5, -.5, +.5, +.5, 0, 0, # (2) MN1: 0 vs 1 0, 0, -.5, -.5, +.5, +.5, # (3) MN2: 1 vs 2 +.5, -.5, -.5, +.5, 0, 0, # (4) CV x MN1 0, 0, +.5, -.5, -.5, +.5 # (5) CV x MN2 ), 6, 5), 6) contrasts(d$design1) print(m1.VV <- lmer(rt ~ design1 + (1|ID), d, subset=EXP=="VV"), cor=FALSE) print(m1.VA <- lmer(rt ~ design1 + (1|ID), d, subset=EXP=="VA"), cor=FALSE) print(m1.AV <- lmer(rt ~ design1 + (1|ID), d, subset=EXP=="AV"), cor=FALSE) print(m1.AA <- lmer(rt ~ design1 + (1|ID), d, subset=EXP=="AA"), cor=FALSE) # ... (2) CV nested within levels of MN, used for post-hoc follow-up levels(d$MN)[3] <- "2" d$design2 <- factor(paste(d$MN, d$CV, sep="")) d$design2 <- C(d$design2, matrix(c(-.5, +.5, 0, 0, 0, 0, # (1) CV| MN = 0 0, 0, -.5, +.5, 0, 0, # (2) CV| MN = 1 0, 0, 0, 0, -.5, +.5, # (3) CV| MN = 2 -.5, -.5, +.5, +.5, 0, 0, # (4) MN1: 0 vs 1 0, 0, -.5, -.5, +.5, +.5 # (5) MN2: 1 vs 2 ), 6, 5), 6) contrasts(d$design2) print(m2.VV <- lmer(rt ~ design2 + (1|ID), d, subset=EXP=="VV"), cor=FALSE) print(m2.VA <- lmer(rt ~ design2 + (1|ID), d, subset=EXP=="VA"), cor=FALSE) print(m2.AV <- lmer(rt ~ design2 + (1|ID), d, subset=EXP=="AV"), cor=FALSE) print(m2.AA <- lmer(rt ~ design2 + (1|ID), d, subset=EXP=="AA"), cor=FALSE) # ... (3) LME for Cue x MN interaction, probably not quite correct yet. cast(d.rt, Cue ~ MN, mean ) d$Cue <- relevel(d$Cue, ref=2) d$design3 <- factor(paste(d$Cue, d$MN, sep="")) d$design3 <- C(d$design3, matrix(c(-.5, +.5, 0, 0, 0, 0, # (1) MN1| Cue=A 0, -.5, +.5, 0, 0, 0, # (2) MN2| Cue=A 0, 0, 0, -.5, +.5, 0, # (3) MN1| Cue=V 0, 0, 0, 0, -.5, +.5, # (4) MN2| Cue=V -.5, -.5, -.5, +.5, +.5, +.5 # (5) Cue ), 6, 5), 6) d$design3a <- C(d$design3, matrix(c(-.5, +.5, 0, -.5, +.5, 0, # (1) MN1 0, -.5, +.5, 0, -.5, +.5, # (2) MN2 0, 0, 0, 1, 1, 1, # (3) Cue 0, 0, 0, -.5, +.5, 0, # (3) MN1 x Cue 0, 0, 0, 0, -.5, +.5 # (4) MN2 x Cue ), 6, 5), 6) contrasts(d$design3a) print(m4 <- lmer(rt ~ design3a + (1|ID), d), cor=FALSE) # ... (4) LME for trials without microsaccades, follow-up cast(d.rt, Target ~ CV, mean, subset=MN=="0") d$design4 <- factor(paste(d$Target, d$CV, sep="")) d$design4 <- C(d$design4, matrix(c(-.5, +.5, 0, 0, # (1) CV| Target=V 0, 0, -.5, +.5, # (2) CV| Target=A -.5, -.5, +.5, +.5 # (3) Target ), 4, 3), 3) contrasts(d$design4) print(m4.NoMS <- lmer(rt ~ Cue*design4 + (1|ID), d, subset=MN=="0"), cor=FALSE) # ... (5) LME for overall analyses 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) d$MN0 <- C(d$MN, matrix(c( 1, -.5, -.5, # (1) 0 vs. 1 or 2+ 0, -.5, +.5 # (2) type of target ), 3, 2), 2) print(m5 <- lmer(rt ~ EXP2*cv.e*MN0 + (1|ID), d), cor=FALSE) # Use the effect coding for manuscript