# GRAPHICS for:
# Wotschack, C., & Kliegl, R. (2013). 
# Reading strategy modulates parafoveal-on-foveal effects in sentence reading.
# The Quarterly Journal of Experimental Psychology, 66(3), 548-562.

# Christiane Wotschack, November 2011

# adapted to current version of ggplot2.
# Reinhold Kliegl, March 2013

rm(list=ls())

# LIBRARIES
library(grid)
library(ggplot2)
library(gtools)

#--------------------------------------------
# FIGURE 1: FREQUENCY EFFECTS (YOUNG ADULTS)
#--------------------------------------------
# 09/2009

load("yng_bd3.rda")

# VARIABLES
CND <- "cnd" # experimental condition (0 = easy-question, 1 = hard-question)
y    <- "dur" # single fixation duration
VAR   <- c("f","f2") # frequency word n, frequency word n+1
NB    <- c( 9, 9) 	#Number of Bins per VAR
			 
					
# for each condition seperately
bo <- b[which(b$cnd==0),]
bp <- b[which(b$cnd==1),]

# create categorized varibales
VARq <- paste(VAR,".q" ,sep="") 
VARp <- paste(VAR,".qm",sep="")

for (i in 1:length(VAR)){   
	bo[,VARq[i]] <- quantcut(bo[,VAR[i]], q=seq(0,1, length.out=NB[i]+1) ) 
	bp[,VARq[i]] <- quantcut(bp[,VAR[i]], q=seq(0,1, length.out=NB[i]+1) )
	bo <- merge(bo, data.frame(tmp =tapply(bo[,VAR[i]], bo[,VARq[i]], mean)), by.x=VARq[i], by.y=0, sort=FALSE)
	bp <- merge(bp, data.frame(tmp =tapply(bp[,VAR[i]], bp[,VARq[i]], mean)), by.x=VARq[i], by.y=0, sort=FALSE)
	colnames(bo)[colnames(bo)=="tmp" ] <- VARp[i]
	colnames(bp)[colnames(bp)=="tmp" ] <- VARp[i]
	print(VAR[i])
	}

# recompute y to get within-subject CIs (cf. Cousineau, 2005)  

bo <- merge(bo, data.frame(m_id=tapply(bo[,y], bo$id, mean)), by.x="id", by.y=0, sort=FALSE)
bp <- merge(bp, data.frame(m_id=tapply(bp[,y], bp$id, mean)), by.x="id", by.y=0, sort=FALSE)

mo_overall <- mean(bo[,y])
mp_overall <- mean(bp[,y])
bo[,y] <- bo[,y] - bo$m_id + mo_overall
bp[,y] <- bp[,y] - bp$m_id + mp_overall

# recombine date from two conditions
b <- rbind(bo, bp)

# prepare dataframe
p <- data.frame(y=b[,y])
p$Condition <- b[,CND]
levels(p$Condition) <- c("easy-question", "hard-question")
p[,VARp] <- b[,VARp]
p[,VAR] <- b[,VAR]
	
i <- 1 # plot current frequency effect (left panel)

plot1 <- ggplot(aes(y=y, group=Condition), data=p) + 
	stat_summary(aes(x=f.qm, group=Condition, shape=Condition), geom="point", fun.y=mean, size=4) + 
	stat_summary(aes(x=f.qm, group=Condition), geom="errorbar", fun.data="mean_cl_normal", conf.int=0.95, width= (max(b[VARp[i]])-min(b[VARp[i]])) /(NB[i]*3)) + 
	geom_smooth(aes(x=f, group=Condition, linetype=Condition), method="loess", se=FALSE, color="grey30") + 
	scale_y_continuous(name="Single fixation duration [ms]", breaks=seq(200,280,length=5)) +
	scale_x_continuous(name="Current Word Frequency [log]", limits=c(0,4.51)) + 
	coord_cartesian(ylim=c(195,285)) +
	theme_bw(base_size=15) + 
	theme(legend.position="none", 
		panel.border = element_rect(colour = "black", size=0.8), 
		panel.grid.major = element_line(colour = "grey80")) 

quartz(width=11, height=6)
pushViewport(viewport(layout = grid.layout(1, 2, widths=unit(c(0.47,0.63),"null") )))
print(plot1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))


i <- 2 # plot next word frequency effect (right panel)
	
plot2 <- ggplot(aes(y=y, group=p$Condition), data=p) + 
	stat_summary(aes(x=f2.qm, group=Condition, shape=Condition), geom="point", fun.y=mean, size=4) + 
	stat_summary(aes(x=f2.qm, group=Condition), geom="errorbar", fun.data="mean_cl_normal", conf.int=0.95, width= (max(b[VARp[i]])-min(b[VARp[i]])) /(NB[i]*3)) + 
	geom_smooth(aes(x=f2, group=Condition, linetype=Condition), method="loess", se=FALSE, color="grey30") + 
	scale_y_continuous(name="", breaks=seq(200,280,length=5)) + 
	scale_x_continuous(name="Next Word Frequency [log]", limits=c(0, 4.51)) +  
	coord_cartesian(ylim=c(195,285)) + 
	theme_bw(base_size=15) + 
	theme(legend.position="right", 
		panel.border = element_rect(colour = "black", size=0.8), 
		panel.grid.major = element_line(colour = "grey80")) 


print(plot2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))

#--------------------------------------------
# FIGURE 2: FREQUENCY EFFECTS (OLD ADULTS)
#--------------------------------------------
# 09/2009

rm(list=ls())
load("old_bd3.rda")

# VARIABLES
CND <- "cnd" # experimental condition (0 = easy-question, 1 = hard-question)
y    <- "dur" # single fixation duration
VAR   <- c("f","f2") # frequency word n, frequency word n+1
NB    <- c( 9, 9) 	#Number of Bins per VAR
			 
					
# for each condition seperately
bo <- b[which(b$cnd==0),]
bp <- b[which(b$cnd==1),]

# create categorized varibales
VARq <- paste(VAR,".q" ,sep="") 
VARp <- paste(VAR,".qm",sep="")

for (i in 1:length(VAR)){   
	bo[,VARq[i]] <- quantcut(bo[,VAR[i]], q=seq(0,1, length.out=NB[i]+1) ) 
	bp[,VARq[i]] <- quantcut(bp[,VAR[i]], q=seq(0,1, length.out=NB[i]+1) )
	bo <- merge(bo, data.frame(tmp =tapply(bo[,VAR[i]], bo[,VARq[i]], mean)), by.x=VARq[i], by.y=0, sort=FALSE)
	bp <- merge(bp, data.frame(tmp =tapply(bp[,VAR[i]], bp[,VARq[i]], mean)), by.x=VARq[i], by.y=0, sort=FALSE)
	colnames(bo)[colnames(bo)=="tmp" ] <- VARp[i]
	colnames(bp)[colnames(bp)=="tmp" ] <- VARp[i]
	print(VAR[i])
	}

# recompute y to get within-subject CIs (cf. Cousineau, 2005)  

bo <- merge(bo, data.frame(m_id=tapply(bo[,y], bo$id, mean)), by.x="id", by.y=0, sort=FALSE)
bp <- merge(bp, data.frame(m_id=tapply(bp[,y], bp$id, mean)), by.x="id", by.y=0, sort=FALSE)

mo_overall <- mean(bo[,y])
mp_overall <- mean(bp[,y])
bo[,y] <- bo[,y] - bo$m_id + mo_overall
bp[,y] <- bp[,y] - bp$m_id + mp_overall

# recombine date from two conditions
b <- rbind(bo, bp)

# prepare dataframe
p <- data.frame(y=b[,y])
p$Condition <- b[,CND]
levels(p$Condition) <- c("easy-question", "hard-question")
p[,VARp] <- b[,VARp]
p[,VAR] <- b[,VAR]

	
i <- 1 # plot current frequency effect (left panel)

plot1 <- ggplot(aes(y=y, group=Condition), data=p) + 
	stat_summary(aes(x=f.qm, group=Condition, shape=Condition), geom="point", fun.y=mean, size=4) + 
	stat_summary(aes(x=f.qm, group=Condition), geom="errorbar", fun.data="mean_cl_normal", conf.int=0.95, width= (max(b[VARp[i]])-min(b[VARp[i]])) /(NB[i]*3)) + 
	geom_smooth(aes(x=f, group=Condition, linetype=Condition), method="loess", se=FALSE, color="grey30") + 
	scale_y_continuous(name="Single fixation duration [ms]", breaks=seq(200,280,length=5)) +
	scale_x_continuous(name="Current Word Frequency [log]", limits=c(0,4.51)) + 
	coord_cartesian(ylim=c(195,285)) +
	theme_bw(base_size=15) + 
	theme(legend.position="none", 
		panel.border = element_rect(colour = "black", size=0.8), 
		panel.grid.major = element_line(colour = "grey80")) 

quartz(width=11, height=6)
pushViewport(viewport(layout = grid.layout(1, 2, widths=unit(c(0.47,0.63),"null") )))
print(plot1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))


i <- 2 # plot next word frequency effect (right panel)
	
plot2 <- ggplot(aes(y=y, group=p$Condition), data=p) + 
	stat_summary(aes(x=f2.qm, group=Condition, shape=Condition), geom="point", fun.y=mean, size=4) + 
	stat_summary(aes(x=f2.qm, group=Condition), geom="errorbar", fun.data="mean_cl_normal", conf.int=0.95, width= (max(b[VARp[i]])-min(b[VARp[i]])) /(NB[i]*3)) + 
	geom_smooth(aes(x=f2, group=Condition, linetype=Condition), method="loess", se=FALSE, color="grey30") + 
	scale_y_continuous(name="", breaks=seq(200,280,length=5)) + 
	scale_x_continuous(name="Next Word Frequency [log]", limits=c(0, 4.51)) +  
	coord_cartesian(ylim=c(195,285)) + 
	theme_bw(base_size=15) + 
	theme(legend.position="right", 
		panel.border = element_rect(colour = "black", size=0.8), 
		panel.grid.major = element_line(colour = "grey80")) 

print(plot2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))