##### Potsdam Mind Research Repository #####

## This script is part of the publication package for the article:
# Kliegl, R., Hohenstein, S., Yan, M., & McDonald, S. (in press).
# How preview space/time translates into preview cost/benefit for fixation durations during reading.
# Quarterly Journal of Experimental Psychology. 

## The data is from:
# McDonald, S. Parafoveal preview benefit in reading is only obtained from
# the saccade goal. Vision Research (2006) vol. 46 (26) pp. 4416-24
# http://dx.doi.org/10.1016/j.visres.2006.08.027


# History:
#  ... plus reanalyses, 14-11-2010, Reinhold Kliegl 
#  ... corrected and modified, 09-02-2011, Sven Hohenstein
#  ... corrected and modified, 11-02-2011, Reinhold Kliegl & Sven Hohenstein
#  ... modified and updated, 23-04-2012, Sven Hohenstein
#  ... updated, 19-06-2013, Sven Hohenstein


rm(list =ls())

library(ggplot2) # the script works for ggplot 0.9.3.1
library(lme4)
# library(grid)
library(gridExtra)
source("remef.v0.6.7.R")


data <- read.table("data.Mid_post_bound_first-pass_16SUs.V1.dat", header=TRUE)
# 12 pixels per letter - 96 in total for 8-letter (including spc before) target words.

data$insa <- data$insa/12
data$lp <- data$lp/12
data$wlpix <- data$wlpix/12
names(data)[7] <- "wl"


ix <- which(data$TB=="T")

dat <- data.frame(cbind(data[ix,  c(1, 3:8, 10:14)], data[ix-1, c(5:8, 10:14) ]))
#names(dat)[11] <- "sfd"                  # rename ffd to sfd 
dat$sfd = ifelse(dat$nfix==1, dat$ffd, NA)
dat$sfd.1 = ifelse(dat$nfix.1==1, dat$ffd.1, NA)

levels(dat$prev) <- c("correct", "incorrect")
levels(dat$bound) <- c("mid-word", "post-word")

#dat <- dat[ , -c(5, 6, 12 )]

# Center covariates
dat$cfrq    <- scale(dat$lnfreq,   center=TRUE, scale=FALSE)
dat$cinsa   <- scale(dat$insa,     center=TRUE, scale=FALSE)
dat$cfrq.1  <- scale(dat$lnfreq.1, center=TRUE, scale=FALSE)
dat$cinsa.1 <- scale(dat$insa.1,   center=TRUE, scale=FALSE)
dat$cgaze.1   <- scale(log(dat$gaze.1),   center=TRUE, scale=FALSE)
dat$cffd.1   <- scale(log(dat$ffd.1),   center=TRUE, scale=FALSE)
dat$csfd.1   <- scale(log(dat$sfd.1),   center=TRUE, scale=FALSE)


# Delete 21 trials w/o eye-movement info about pretarget word
dat <- subset(dat, !is.na(cgaze.1))  

# Set up contrasts as nested within boundary condition
dat$cond <- factor(paste(dat$bound, dat$prev, sep="_"), labels=c("M-C", "M-I", "P-C", "P-I"))
cmat <- matrix(c( -1/2, -1/2, +1/2, +1/2,            # Main effect Bound
	              -1,    +1,    0,    0,             # Nested effect 1: Prev | Bound == mid-word
	               0,     0,   -1,   +1),  4,  3)    # Nested effect 2: Prev | Bound == post-word  

rownames(cmat) <- c("M-C", "M-I", "P-C", "P-I")
colnames(cmat) <- c(".Bound", ".Prev|Bound==M", ".Prev|Bound==P")

(contrasts(dat$cond) <- cmat)

######## SFD with preview nested in boundary #######
print(m0.lSFD.ML  <- lmer(log(sfd) ~ qlp + cfrq + cinsa*cond + (1 | subj) + (1 | word), dat, REML=FALSE), cor=FALSE)       
print(m1.lSFD.ML  <- lmer(log(sfd) ~ qlp + cfrq + (cinsa + cgaze.1)*cond + (1 | subj) + (1 | word), dat, REML=FALSE), cor=FALSE)     
anova(m0.lSFD.ML, m1.lSFD.ML)

# ... final baseline model
print(m0.lSFD  <- lmer(log(sfd) ~ qlp + cfrq + cinsa*cond + (1 | subj) + (1 | word), dat), cor=FALSE)       

# ... final expanded model
print(m1.lSFD  <- lmer(log(sfd) ~ qlp + cfrq + (cinsa + cgaze.1)*cond + (1 | subj) + (1 | word), dat), cor=FALSE)     

# other combinations
print(m1.lFFD  <- lmer(log(ffd) ~ qlp + cfrq + (cinsa + cgaze.1)*cond + (1 | subj) + (1 | word), dat), cor=FALSE)     
print(m1.lGD   <- lmer(log(gaze)~ qlp + cfrq + (cinsa + cgaze.1)*cond + (1 | subj) + (1 | word), dat), cor=FALSE)     

print(m2.lSFD  <- lmer(log(sfd) ~ qlp + cfrq + (cinsa + cffd.1)*cond + (1 | subj) + (1 | word), dat), cor=FALSE)     
print(m2.lFFD  <- lmer(log(ffd) ~ qlp + cfrq + (cinsa + cffd.1)*cond + (1 | subj) + (1 | word), dat), cor=FALSE)     
print(m2.lGD   <- lmer(log(gaze)~ qlp + cfrq + (cinsa + cffd.1)*cond + (1 | subj) + (1 | word), dat), cor=FALSE)     

# ... ... check the two slopes in post-word condition (NOT OPTIMAL; TOO MUCH LOSS OF POWER)
ixPC <- which(dat$bound=="post-word" & dat$prev=="correct")
print(m1.lSFD.PC  <- lmer(log(sfd) ~ qlp + cfrq + cinsa + cgaze.1 + (1 | subj) + (1 | word), dat[ixPC, ]), cor=FALSE)       

ixPI <- which(dat$bound=="post-word" & dat$prev=="incorrect")
# extended models (for cgaze.1)
print(m1.lSFD.PI  <- lmer(log(sfd) ~ qlp + cfrq + cinsa + cgaze.1 + (1 | subj) + (1 | word), dat[ixPI, ]), cor=FALSE)     
print(m1.lFFD.PI  <- lmer(log(ffd) ~ qlp + cfrq + cinsa + cgaze.1 + (1 | subj) + (1 | word), dat[ixPI, ]), cor=FALSE)    
print(m1.lGD.PI  <- lmer(log(gaze) ~ qlp + cfrq + cinsa + cgaze.1 + (1 | subj) + (1 | word), dat[ixPI, ]), cor=FALSE)    
# base models (for insa)
print(m1.lSFD.PI  <- lmer(log(sfd) ~ qlp + cfrq + cinsa + (1 | subj) + (1 | word), dat[ixPI, ]), cor=FALSE)     
print(m1.lFFD.PI  <- lmer(log(ffd) ~ qlp + cfrq + cinsa + (1 | subj) + (1 | word), dat[ixPI, ]), cor=FALSE)    
print(m1.lGD.PI  <- lmer(log(gaze) ~ qlp + cfrq + cinsa + (1 | subj) + (1 | word), dat[ixPI, ]), cor=FALSE)    

  

#################### FIGURES 5 (SFD; effect plots) ##########################
# log SFD 
dat$sfd.insa[!is.na(dat$sfd)] <- exp(remef(m1.lSFD, fix = c(2:3, 5, 12:14), ran = list(1,1)))

# ... 3-factor insa interaction (for log SFD) --> FIGURE 5a
xbreaks <- 2:11
ybreaks <- seq(150, 350, 50)
p <- ggplot(data=dat, aes(x = insa, y = sfd.insa, linetype = prev)) + facet_grid( . ~ bound)
p <- p + stat_smooth(method = "lm", formula = y ~ poly(x, 1), size = 1.2, colour = "black") + scale_linetype("Preview")
p <- p + scale_x_reverse("Incoming saccade amplitude [letters]", breaks = xbreaks, labels = as.character(xbreaks)) 
p <- p + scale_y_continuous("Target single-fixation duration [ms]", breaks = ybreaks, labels = as.character(ybreaks), trans = "log")
p <- p + coord_cartesian(ylim = c(min(ybreaks), max(ybreaks)), xlim = c(1.5,12)) + theme(legend.position = c (.62, .28))
#p <- p + geom_vline(xintercept = mean(dat$insa), legend = FALSE, linetype = 1, size = 0.5, colour = "red")
(sfd.insa <- p)

# ... 3-factor gaze.1 interaction (for log SFD) --> FIGURE 5b
# Note: there are always exactly two fixations for mid-word boundary!
dat$sfd.gaze.1[!is.na(dat$sfd)] <- exp(remef(m1.lSFD, fix = c(2:3, 4, 9:11), ran = list(1,1)))

dat$gaze.1b <- ifelse(dat$bound=="mid-word", dat$gaze.1, dat$gaze.1 * 2)
xbreaks <- 2 * seq(150, 350, 50)
# dat$gaze.1b <- ifelse(dat$bound=="mid-word", dat$gaze.1/2, dat$gaze.1)
# xbreaks <- seq(150, 350, 50)
ybreaks <- seq(150, 350, 50)
p <- ggplot(data=dat, aes(x = gaze.1b, y = sfd.gaze.1, linetype = prev)) 
p <- p + stat_smooth(method = "lm", formula = y ~ poly(x, 1), size = 1.2, colour = "black") + scale_linetype("Preview")
p <- p + scale_x_continuous("Gaze duration on pretarget [ms]", breaks = xbreaks, labels = as.character(xbreaks), trans = "log") 
p <- p + scale_y_continuous("Target single-fixation duration [ms]", breaks = ybreaks, labels = as.character(ybreaks), trans = "log")    
p <- p + facet_grid( . ~ bound)
p <- p + coord_cartesian(ylim = c(min(ybreaks), max(ybreaks)), xlim = 2 * c(150, 400)) + theme(legend.position = c (.62, .28))
# p <- p + coord_cartesian(ylim = c(min(ybreaks), max(ybreaks)), xlim = c(150, 400)) + theme(legend.position = c (.62, .28))
#p <- p + geom_vline(xintercept = mean(dat$gaze.1b*(3/4)), legend = FALSE, linetype = 1, size = 0.5, colour = "red")
(sfd.gaze1 <- p)


# modify x-axis tick labels of right-hand facet
grob <- ggplotGrob(sfd.gaze1)
grob[["grobs"]][[8]][["children"]][["axis"]][["grobs"]][[2]][["label"]] <- as.character(xbreaks / 2)

grid.arrange(sfd.insa, grob, ncol = 1)