# Analyses for Angele, Slattery, Yang, Kliegl, & Rayner (2008): Parafoveal processing in reading: Manipulating n + 1 and n + 2 previews simultaneously. Visual Cognition, 16(6), 697–707.

rm(list=ls())
require(reshape)

workingdir <- choose.dir()

#Alternative: Set working directory directly:

# workingdir <- "/working directory"
setwd(workingdir)

# directory must contain "Item by subject data.csv"

isd <- read.csv2("Item by subject data.csv", header=TRUE, skip=1)

isd <- isd[1:4552,]
a <- melt(isd, id=c("subjectID", "ItemID", "frequency", "preview", "Cond.No.in.the.scipt"))

a$value[a$value==0] <- NA

levels(a$variable) <- c("ffd.n", "ffd.n1", "ffd.n2", "del.X", "sfd.n", "sfd.n1", "sfd.n2", "del.X1", "gzd.n", "gzd.n1", "gzd.n2", "del.X2", "lpo.n", "lpo.n1", "lpo.n2")
a <- cbind(a,colsplit(a$variable, names=c("measure","twx"),split="\\."))

a<-subset(a,measure!="del")

a$id <- factor(a$subjectID)
a$sn <- factor(a$ItemID)

# make data frame for analysis

d.all <- with(a,data.frame(id = id, sn = sn, twx = as.numeric(twx), prevw = as.numeric(preview), freqc = as.numeric(frequency), measure=measure, value=value))
d.all$twx <- as.factor(d.all$twx)
levels(d.all$twx) <- c("N", "N+1", "N+2")
d.all$prevw <- as.factor(d.all$prevw)
levels(d.all$prevw) <- c("both different", "both identical", "n+1", "n+2")
d.all$freqc <- as.factor(d.all$freqc)
levels(d.all$freqc) <- c("high", "low")

# contrast coding
# x1: both different vs. N+1 different (2 vs 3): N+2 preview effect when N+1 was masked
# x2: both identical vs. N+2 different (1 vs 4): N+2 preview effect when N+1 was not masked
# x3: 2 and 3 vs 1 and 4: N+1 preview effect

d.all$x1<-ifelse(d.all$prevw=="both different" | d.all$prevw=="n+1",ifelse(d.all$prevw=="both different",.5,-.5),0)
d.all$x2<-ifelse(d.all$prevw=="both identical" | d.all$prevw=="n+2",ifelse(d.all$prevw=="both identical",-.5,.5),0)
d.all$x3<-ifelse(d.all$prevw=="both identical" | d.all$prevw=="n+2",-.25,.25)

d.m <- d.all[!is.na(d.all$value),] 

d.m$value <- log(d.m$value) # log transform all values
library(lme4)

ffd.n <- lmer(data=d.m, value ~ x1 + x2 + x3 + freqc + (1|id) + (1|sn), subset= measure=="ffd" & twx == "N")
ffd.n1 <- lmer(data=d.m, value ~ x1 + x2 + x3 + freqc + (1|id) + (1|sn), subset= measure=="ffd" & twx == "N+1")
ffd.n2 <- lmer(data=d.m, value ~ x1 + x2 + x3 + freqc + (1|id) + (1|sn), subset= measure=="ffd" & twx == "N+2")

gzd.n <- lmer(data=d.m, value ~ x1 + x2 + x3 + freqc + (1|id) + (1|sn), subset= measure=="gzd" & twx == "N")
gzd.n1 <- lmer(data=d.m, value ~ x1 + x2 + x3 + freqc + (1|id) + (1|sn), subset= measure=="gzd" & twx == "N+1")
gzd.n2 <- lmer(data=d.m, value ~ x1 + x2 + x3 + freqc + (1|id) + (1|sn), subset= measure=="gzd" & twx == "N+2")

#extract coefficients and variance estimates for power analysis

get.s <- function(x) { c(s.sn = sqrt(as.numeric(VarCorr(x)$sn)), s.id = sqrt(as.numeric(VarCorr(x)$id)), s.y =  as.numeric(attr(VarCorr(x),"sc"))) }
allef <- function(x) { c(fixef(x), get.s(x), N = length(summary(x)@frame$value)
) }

all.ffd.f <- rbind(allef(ffd.n),allef(ffd.n1),allef(ffd.n2))
rownames(all.ffd.f) <- c("FFD n", "FFD n+1", "FFD n+2")

all.gzd.f <- rbind(allef(gzd.n),allef(gzd.n1),allef(gzd.n2))
rownames(all.gzd.f) <- c("GZD n", "GZD n+1", "GZD n+2")

all.f <- rbind(all.ffd.f,all.gzd.f)

save(all.f, file="lmer estimates for power analysis log.rData")