# ANALYSES AND PLOTS:
# ARTICLE "Adult Age Differences in the Perceptual Span" (2010)
# Psychology and Aging

# Sarah Risse, August 2010



rm(list=ls())

# LIBRARIES:
library(lme4)
library(reshape)
library(Hmisc)

# LOAD DATA:
load("RK2010_data.rda")
ls()			# em_n0, em_n1, em_n2


# STRUCTURE DATA FRAMES:

# em_n0:	eye-movement data on preboundary word N
# em_n1:	... on postboundary word N+1 (three letters length)
# em_n2:	... on target word N+2

# (1)   id:		subject id (nested within age)
# (2)   sn:		sentence id
# (3)   wn: 	word number in sentence
# (4)   nw: 	number of words in sentence
# (5)   pvn2:	preview condition of word N+2, -0.5 -> identical preview; 0.5 -> nonword preview
# (6)   lxn1:	lexical status of word N+1 (i.e., processing difficulty), -0.5 -> function word (easy); 0.5 -> content word (difficulty)
# (7)   f:		current word (n) log10-frequency, centered around zero 
# (8)   wid:	unique word id in sentence (same word counted as same across sentences)
# (9)   wl:		1/word length of current word (n), centered around zero 
# (10)  wl1:	wl of previous word to the left (n-1)
# (11)  f1: 	f of (n-1)
# (12)  wl2:  	wl of upcoming word to the right (n+1)
# (13)  f2:		f of (n+1)
# (14)  wl3: 	wl of penultimate word to the left (n-2)
# (15)  f3:		f of (n-2)
# (16)  wl4:	wl of the second next word to the right (n+2)
# (17)  f4: 	f of (n+2)
# (18)  ffd:	first fixation duration
# (19)  gzd:	gaze duration
# (20)  sfd:	single fixation duration
# (21)  tvt:	total viewing time
# (22)  prx:	probability of refixations (0: no refixation, 1: refixated)
# (23)  psk:	probability of skipping (0: fixated, 1: skipped) 
# (24)  prg:	probability of regression (out of this word)
# (25)  ilp:	initial landing position, centered
# (26)  age:	age group, -0.5 -> young adults, 0.5 -> old adults
# (27)  sn1:	only for em_n0 and em_n2: skipping of word N+1, -0.5 -> fixated; 0.5 -> skipped


# (1) EXCLUDE SUBJECTS:
# (1.1) YOUNGER ADULTS:
# 2 subjects had to be excluded from analysis
# Vp20: raw data was incomplete thus not even preprocessed
# Vp30: noticed display changes (id 29 in R-data)
# (1.2) OLDER ADULTS:
# 2 subjects had to be excluded from analysis
# Vp09: noticed about 10 display changes (id 109 in R-data)
# Vp29: noticed 2-3 display changes (id 129 in R-data)

exclude <- which(em_n0$id==29 | em_n0$id==109 | em_n0$id==129)
em_n0 <- em_n0[-exclude,] 
em_n1 <- em_n1[-exclude,] 
em_n2 <- em_n2[-exclude,] 

#-------------------------------------
# (2) LMM ANALYSIS:
# (2.1) TARGET WORD N+2:
#-------------------------------------

xdat <- em_n2

# run lmer-models:


lm1.gzd <- lmer(log(gzd) ~ (age*lxn1*pvn2*sn1) + (1|id) + (1|wid) + (1|sn), data=xdat) 
print(lm1.gzd, cor=F)

# ADDITIONAL LMM -- CONTRASTS:
# SPECIFICATION OF NESTED CONTRASTS IN LMM:

xdat$LexStat <- as.factor(xdat$lxn1)
levels(xdat$LexStat) <- c("CW","FW")
xdat$Preview <- as.factor(xdat$pvn2)
levels(xdat$Preview) <- c("identical","non-preview")
xdat$Skipping <- as.factor(xdat$sn1)
levels(xdat$Skipping) <- c("fixated","skipped")
xdat$AgeGroup <- as.factor(xdat$age)
levels(xdat$AgeGroup) <- c("young","old")

# "cnd" is some kind of help variable that codes the combination of all 
# factor levels in your design (for each row of data frame)

xdat$cnd <- as.factor(paste(xdat$AgeGroup,xdat$Skipping, xdat$Preview))
str(xdat)
levels(xdat$cnd)		# shows you what combinations of factor levels 
# are now considered in this help variable

mat <- matrix(c( -1, -1, -1, -1,  1,  1,  1,  1, # main effect age
		     -1,  1,  0,  0,  0,  0,  0,  0, # PBE old n1 fixated
		      0,  0, -1,  1,  0,  0,  0,  0, # PBE old n1 skipped
                  0,  0,  0,  0, -1,  1,  0,  0, # PBE young n1 fixated
			0,  0,  0,  0,  0,  0, -1,  1), 8, 5)/2   


colnames(mat) <- c(".Age",".PBE:old:n1.fix",".PBE:old:n1.skp",".PBE:yng:n1.fix",".PBE:yng:n1.skp")

xdat$c1 <- C(xdat$cnd, mat, 5)

m0 <- lmer(log(gzd) ~ (lxn1*c1) + (1|id) + (1|wid) + (1|sn), data=xdat)
print(m0, cor=F)


#-------------------------------------
# (2) FIGURE 2:
# (2.2) TARGET WORD N+2:
#-------------------------------------

# Figure 2. Difference in N+2 preview benefit conditional on 
# skipping word N+1. Plotted is the mean GD on the target word N+2 
# for the N+2 preview conditions by skipping the preceding word N+1. 
# Error bars represent the 95 % confidence interval.  


# descriptives/ effect sizes:
n2.rs <- melt(xdat, id=c("id","age","lxn1","pvn2","sn1")
		, measure=c("gzd")
		, na.rm=TRUE)

n2.rs$age <- as.factor(n2.rs$age)
levels(n2.rs$age) <- c("young","older")

n2.rs$pvn2 <- as.factor(n2.rs$pvn2)
levels(n2.rs$pvn2) <- c("identical","nonword")

n2.rs$lxn1 <- as.factor(n2.rs$lxn1)
levels(n2.rs$lxn1) <- c("FW","CW")

n2.rs$sn1 <- as.factor(n2.rs$sn1)
levels(n2.rs$sn1) <- c("fixated","skipped")


mw <- cast(n2.rs, pvn2 ~ sn1
	, function(x) c(M=signif(mean(x),3)
	, SE=1.96*(sd(x)/(length(x)^.5)) ))

mw$fix.lower <- mw$fixated_M - mw$fixated_SE
mw$fix.upper <- mw$fixated_M + mw$fixated_SE

mw$skp.lower <- mw$skipped_M - mw$skipped_SE
mw$skp.upper <- mw$skipped_M + mw$skipped_SE

data  <- as.matrix(mw[,c("fixated_M","skipped_M")])
lower <- as.matrix(mw[,c("fix.lower","skp.lower")])
upper <- as.matrix(mw[,c("fix.upper","skp.upper")])


windows(1,width=4, height=5, pointsize=11)
par(mfrow=c(1,1),lwd=2, cex.axis=1.2, cex.lab=1.2, mar=c(5,4.2,2.5,2.5))


plot(c(1,2),data[1,], type='b',lty='solid'
	, pch=21
	, bg='black'
	, ylim=c(200,300)
	, xlim=c(0.8,2.2)
	, ylab=c('Gaze duration on word N+2 [ms]')
	, xlab=c('Skipping of word N+1')
	, axes=F)
	axis(side=1, at=c(1,2), labels=c("fixated","skipped")
	, tick=T)
	axis(side=2)
	errbar(c(1,2),data[1,], yplus=upper[1,]
	, yminus=lower[1,], add=T, col="black", lty="solid"
	, lwd=1.5)
	errbar(c(1,2),data[2,], yplus=upper[2,]
	, yminus=lower[2,], add=T, col="black", lty="solid"
	, lwd=1.5)
	lines(c(1,2),data[2,], type='b', lty='dashed'
	, pch=21, bg='white')
	lines(c(1,2),data[1,], type='b', lty='solid'
	, pch=21, bg='black')	
	box(lty='solid')

	legend(0.9,300, legend=c('identical', 'nonword')
	, title=c("Preview of word N+2")
	, lty=c('solid','dashed')
	, pch=c(21,21)
	, pt.bg=c('black','white')
	, merge=F
	, ncol=1, #bty='n')
	, bg='white'
	, box.col='white')


#-------------------------------------
# (2) FIGURE 3:
# (2.3) TARGET WORD N+2:
#-------------------------------------

# Figure 3. Age difference in post-skipping costs. 
# Plotted is the mean GD on the target word N+2 conditional on 
# skipping the preceding word N+1, both for young and old adults. 
# Error bars represent the 95 % confidence interval.  

mPlot <- cast(n2.rs, age ~ sn1 
	, subset=variable==c("gzd")
	, function(x) c(M=signif(mean(x),3)
	, SE=1.96*(sd(x)/(length(x)^.5)) ))

mPlot
mPlot$fl1.lower <- mPlot[,2] - mPlot[,3]
mPlot$fl1.upper <- mPlot[,2] + mPlot[,3]

mPlot$fl2.lower <- mPlot[,4] - mPlot[,5]
mPlot$fl2.upper <- mPlot[,4] + mPlot[,5]

data  <- as.matrix(mPlot[,c("fixated_M","skipped_M")])
lower <- as.matrix(mPlot[,c("fl1.lower","fl2.lower")])
upper <- as.matrix(mPlot[,c("fl1.upper","fl2.upper")])

windows(1,width=4, height=5, pointsize=11)
par(mfrow=c(1,1),lwd=2, cex.axis=1.2, cex.lab=1.2, mar=c(5,4.2,2.5,2.5))


plot(c(1,2),data[1,], type='b',lty='solid'
	, pch=21
	, bg='black'
	, ylim=c(200,300)
	, xlim=c(0.8,2.2)
	, ylab=c('Gaze duration on word N+2 [ms]')
	, xlab=c('Skipping of word N+1')
	, axes=F)
	axis(side=1, at=c(1,2), labels=c("fixated","skipped")
	, tick=T)
	axis(side=2)
	errbar(c(1,2),data[1,], yplus=upper[1,]
	, yminus=lower[1,], add=T, col="black", lty="solid"
	, lwd=1.5)
	errbar(c(1,2),data[2,], yplus=upper[2,]
	, yminus=lower[2,], add=T, col="black", lty="solid"
	, lwd=1.5)
	lines(c(1,2),data[2,], type='b', lty='dashed'
	, pch=21, bg='white')
	lines(c(1,2),data[1,], type='b', lty='solid'
	, pch=21, bg='black')	
	box(lty='solid')

	legend(0.9,300, legend=c("young adults", "     old adults")
	, title=c("                 Age group")
	, lty=c('solid','dashed')
	, pch=c(21,21)
	, pt.bg=c('black','white')
	, merge=F
	, ncol=1, #bty='n')
	, bg='white'
	, box.col='white')


#-------------------------------------
# (3) LMM ANALYSIS:
# (3.1) PREBOUNDARY WORD N:
#-------------------------------------

xdat <- em_n0


# skipping probability:
length(which(xdat$psk==1))/nrow(xdat)*100


# run lmer-models:

lm.gzd <- lmer(log(gzd) ~ (age*lxn1*pvn2*sn1) + (1|id) + (1|wid) + (1|sn), data=xdat, REML=T) 
print(lm.gzd, cor=F)


#-------------------------------------
# (3) FIGURE 4:
# (3.2) PREBOUNDARY WORD N:
#-------------------------------------

# Figure 4. Age difference in the parafoveal-on-foveal effect of the 
# neighbouring word N+1. Plotted is the mean GD on the pre-boundary word N 
# depending on the processing difficulty of word N+1, both for young and 
# old adults. Error bars represent the 95 % confidence interval. 
# FW: function word; CW: content word.  

# descriptives/ effect sizes:
n2.rs <- melt(xdat, id=c("id","age","lxn1","pvn2","sn1")
		, measure=c("gzd")
		, na.rm=TRUE)

n2.rs$age <- as.factor(n2.rs$age)
levels(n2.rs$age) <- c("young","older")

n2.rs$pvn2 <- as.factor(n2.rs$pvn2)
levels(n2.rs$pvn2) <- c("identical","nonword")

n2.rs$lxn1 <- as.factor(n2.rs$lxn1)
levels(n2.rs$lxn1) <- c("FW","CW")

n2.rs$sn1 <- as.factor(n2.rs$sn1)
levels(n2.rs$sn1) <- c("fixated","skipped")


mPlot <- cast(n2.rs, lxn1 ~ age
	, subset=variable==c("gzd")
	, function(x) c(M=signif(mean(x),3)
	, SE=1.96*(sd(x)/(length(x)^.5)) ))

mPlot$fl1.lower <- mPlot[,2] - mPlot[,3]
mPlot$fl1.upper <- mPlot[,2] + mPlot[,3]

mPlot$fl2.lower <- mPlot[,4] - mPlot[,5]
mPlot$fl2.upper <- mPlot[,4] + mPlot[,5]

data  <- t(as.matrix(mPlot[1:2,c("young_M", "older_M")]))
lower <- t(as.matrix(mPlot[1:2,c("fl1.lower","fl2.lower")]))
upper <- t(as.matrix(mPlot[1:2,c("fl1.upper","fl2.upper")]))

y.limits <- c(200, 300)


windows(1,width=4, height=5, pointsize=11)
par(mfrow=c(1,1),lwd=2, cex.axis=1.2, cex.lab=1.2, mar=c(5,4.2,2.5,2.5))


plot(c(1,2),data[1,], type='b',lty='solid'
	, pch=21
	, bg='black'
	, ylim=y.limits
	, xlim=c(0.8,2.2)
	, ylab=c('Gaze duration on word N [ms]')
	, xlab=c('Processing difficulty of word N+1')
	, axes=F)
	axis(side=1, at=c(1,2), labels=c("easy (FW)","difficult (CW)")
	, tick=T)
	axis(side=2)
	errbar(c(1,2),data[1,], yplus=upper[1,]
	, yminus=lower[1,], add=T, col="black", lty="solid"
	, lwd=1.5)
	errbar(c(1,2),data[2,], yplus=upper[2,]
	, yminus=lower[2,], add=T, col="black", lty="solid"
	, lwd=1.5)
	lines(c(1,2),data[2,], type='b', lty='dashed'
	, pch=21, bg='white')
	lines(c(1,2),data[1,], type='b', lty='solid'
	, pch=21, bg='black')	
	box(lty='solid')

	legend(0.9,300, legend=c("young adults", "     old adults")
	, title=c("                 Age group")
	, lty=c('solid','dashed')
	, pch=c(21,21)
	, pt.bg=c('black','white')
	, merge=F
	, ncol=1, #bty='n')
	, bg='white'
	, box.col='white')



#-------------------------------------
# (4) LMM ANALYSIS:
# (4.1) POSTBOUNDARY WORD N+1:
#-------------------------------------

xdat <- em_n1

# probability of skipping word N+1:
(length(which(xdat$psk==1)))/(length(which(xdat$psk==0 | xdat$psk==1)))*100

# run lmer-models:

lm.gzd <- lmer(log(gzd) ~ (age*lxn1*pvn2) + (1|id) + (1|wid) + (1|sn), data=xdat) 
print(lm.gzd, cor=F)


# POST-HOC LMM's SEPARATE FOR AGE GROUPS:
lm.yng <- lmer(log(gzd) ~ (lxn1*pvn2) + (1|id) + (1|wid) + (1|sn)
		, data=xdat[which(xdat$age==-0.5),]) 
print(lm.yng, cor=F)

lm.old <- lmer(log(gzd) ~ (lxn1*pvn2) + (1|id) + (1|wid) + (1|sn)
		, data=xdat[which(xdat$age==0.5),]) 
print(lm.old, cor=F)

#-------------------------------------
# (4) FIGURE 5:
# (4.2) POSTBOUNDARY WORD N+1:
#-------------------------------------

# Figure 5. N+2 preview effect on the post-boundary word N+1. 
# Plotted is the mean GD on word N+1 conditional on the preview condition 
# of word N+2 for young and old adults. 
# Error bars represent the 95 % confidence interval.  

# descriptives/ effect sizes:
n2.rs <- melt(xdat, id=c("id","age","lxn1","pvn2")
		, measure=c("ffd","gzd")
		, na.rm=TRUE)

n2.rs$age <- as.factor(n2.rs$age)
levels(n2.rs$age) <- c("young","older")

n2.rs$pvn2 <- as.factor(n2.rs$pvn2)
levels(n2.rs$pvn2) <- c("identical","nonword")

n2.rs$lxn1 <- as.factor(n2.rs$lxn1)
levels(n2.rs$lxn1) <- c("FW","CW")


mPlot <- cast(n2.rs, pvn2 ~ age
	, subset=variable==c("gzd")
	, function(x) c(M=signif(mean(x),3)
	, SE=1.96*(sd(x)/(length(x)^.5)) ))

mPlot$fl1.lower <- mPlot[,2] - mPlot[,3]
mPlot$fl1.upper <- mPlot[,2] + mPlot[,3]

mPlot$fl2.lower <- mPlot[,4] - mPlot[,5]
mPlot$fl2.upper <- mPlot[,4] + mPlot[,5]

data  <- t(as.matrix(mPlot[1:2,c("young_M", "older_M")]))
lower <- t(as.matrix(mPlot[1:2,c("fl1.lower","fl2.lower")]))
upper <- t(as.matrix(mPlot[1:2,c("fl1.upper","fl2.upper")]))

y.limits <- c(200, 260)

windows(1,width=4, height=5, pointsize=11)
par(mfrow=c(1,1),lwd=2, cex.axis=1.2, cex.lab=1.2, mar=c(5,4.2,2.5,2.5))


plot(c(1,2),data[1,], type='b',lty='solid'
	, pch=21
	, bg='black'
	, ylim=y.limits
	, xlim=c(0.8,2.2)
	, ylab=c('Gaze duration on word N+1 [ms]')
	, xlab=c('Preview of word N+2')
	, axes=F)
	axis(side=1, at=c(1,2), labels=c("identical","nonword")
	, tick=T)
	axis(side=2)
	errbar(c(1,2),data[1,], yplus=upper[1,]
	, yminus=lower[1,], add=T, col="black", lty="solid"
	, lwd=1.5)
	errbar(c(1,2),data[2,], yplus=upper[2,]
	, yminus=lower[2,], add=T, col="black", lty="solid"
	, lwd=1.5)
	lines(c(1,2),data[2,], type='b', lty='dashed'
	, pch=21, bg='white')
	lines(c(1,2),data[1,], type='b', lty='solid'
	, pch=21, bg='black')	
	box(lty='solid')

	legend(0.9,260, legend=c("young adults", "     old adults")
	, title=c("                 Age group")
	, lty=c('solid','dashed')
	, pch=c(21,21)
	, pt.bg=c('black','white')
	, merge=F
	, ncol=1, #bty='n')
	, bg='white'
	, box.col='white')