######################################################################
#    DISCLAIMER
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
######################################################################
#
#    The R code below was tested on R version 2.13.0
#    September 2011
#
#    Analyses for the study:
#
#    Yan, M., Zhou, W., Shu, H. & Kliegl, R. (in press).
#	 Lexical and Sub-lexical Semantic Preview Benefits in Chinese Reading.
#	 Journal of Experimental Psychology: Learning, Memory, and Cognition.

rm(list=ls())
load("cpd4pmr2.rda")
library(lme4)

target$sfd = ifelse(target$n_fix==1, target$ffd, NA)

# fixation filters
# 1, duration filter
# 2, no regression from N or N+1
# 3, trials in which display change triggers were sent during first 80% of saccades
target1 = target[which(target$ffd > 60 & target$ffd < 600 
			& target$out_sac_len_n0 > 0 & target$out_sac_len > 0
			& target$dc_info!=-1 & target$dc_info<4/5),]
 
target1$id   <- factor(target1$id)
target1$sn   <- factor(target1$sn)
target1$cond <- factor(target1$cond, levels=c( "unr", "tpt", "opq", "idt"))


## basic models
# log-transformed durations
print(lmer(log(ffd) ~ cond + (1|sn) + (1|id), data=target1), cor=FALSE)
print(lmer(log(sfd) ~ cond + (1|sn) + (1|id), data=target1), cor=FALSE)
print(lmer(log(gd) ~ cond + (1|sn) + (1|id), data=target1), cor=FALSE)

# original durations
print(lmer(ffd ~ cond + (1|sn) + (1|id), data=target1), cor=FALSE)
print(lmer(sfd ~ cond + (1|sn) + (1|id), data=target1), cor=FALSE)
print(lmer(gd ~ cond + (1|sn) + (1|id), data=target1), cor=FALSE)


## additional analyses not included in the paper, excluding trials target words were foveally processed 
# log-transformed durations
print(lmer(log(ffd) ~ cond + (1|sn) + (1|id), data=target1, subset=fst_lp<1.2), cor=FALSE)
print(lmer(log(sfd) ~ cond + (1|sn) + (1|id), data=target1, subset=fst_lp<1.2), cor=FALSE)
print(lmer(log(gd) ~ cond + (1|sn) + (1|id), data=target1, subset=fst_lp<1.2), cor=FALSE)

# original durations
print(lmer(ffd ~ cond + (1|sn) + (1|id), data=target1, subset=fst_lp<1.2), cor=FALSE)
print(lmer(sfd ~ cond + (1|sn) + (1|id), data=target1, subset=fst_lp<1.2), cor=FALSE)
print(lmer(gd ~ cond + (1|sn) + (1|id), data=target1, subset=fst_lp<1.2), cor=FALSE)


target1$opq <- ifelse(target1$cond=="opq", 1, 0)
target1$tpt <- ifelse(target1$cond=="tpt", 1, 0)
target1$unr <- ifelse(target1$cond=="unr", 1, 0)
target1$idt <- ifelse(target1$cond=="idt", 1, 0)
sub.target1 <- subset(target1, dc_info!=-1 & dc_info<4/5)

## using likelihood ratio tests for the individual planned comparison, drop1-tests, also use ML instead of REML
# FFD, full model
model.ML.all <- lmer(log(ffd) ~ tpt + opq + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)

# check each of the three terms
model.ML.tpt <- lmer(log(ffd) ~ opq + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.tpt, model.ML.all)   # significant drop in goodness of fit when leaving out tpt = 0.0001112 *** 

model.ML.opq <- lmer(log(ffd) ~ tpt + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.opq, model.ML.all)   # significant drop in goodness of fit when leaving out opq = 0.0008796 *** 

model.ML.idt <- lmer(log(ffd) ~ tpt + opq + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.idt, model.ML.all)   # significant drop in goodness of fit when leaving out idt < 2.2e-16 *** 


# SFD, full model
model.ML.all <- lmer(log(sfd) ~ tpt + opq + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)

# check each of the three terms
model.ML.tpt <- lmer(log(sfd) ~ opq + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.tpt, model.ML.all)   # significant drop in goodness of fit when leaving out tpt = 0.0001175 ***

model.ML.opq <- lmer(log(sfd) ~ tpt + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.opq, model.ML.all)   # significant drop in goodness of fit when leaving out opq = 0.001688 **

model.ML.idt <- lmer(log(sfd) ~ tpt + opq + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.idt, model.ML.all)   # significant drop in goodness of fit when leaving out idt < 2.2e-16 *** 


# GD, full model
model.ML.all <- lmer(log(gd) ~ tpt + opq + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)

# check each of the three terms
model.ML.tpt <- lmer(log(gd) ~ opq + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.tpt, model.ML.all)   # significant drop in goodness of fit when leaving out tpt = 3.034e-06 ***

model.ML.opq <- lmer(log(gd) ~ tpt + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.opq, model.ML.all)   # significant drop in goodness of fit when leaving out opq = 0.01021 *

model.ML.idt <- lmer(log(gd) ~ tpt + opq + (1|sn) + (1|id), REML=FALSE, data=sub.target1)
anova(model.ML.idt, model.ML.all)   # significant drop in goodness of fit when leaving out idt < 2.2e-16 *** 


## use opaque as ref category
target1$cond1 <- factor(target1$cond, levels=c( "opq", "tpt", "unr", "idt"))
sub.target1 <- subset(target1, dc_info!=-1 & dc_info<4/5)

print(model <- lmer(log(gd) ~ cond1 + (1|sn) + (1|id), data=target1), cor=FALSE)
print(model <- lmer(gd ~ cond1 + (1|sn) + (1|id), data=target1), cor=FALSE)

# additional analyses not included in the paper, excluding trials target words were foveally processed 
print(model <- lmer(log(gd) ~ cond1 + (1|sn) + (1|id), data=target1, subset=fst_lp<1.2), cor=FALSE)
print(model <- lmer(gd ~ cond1 + (1|sn) + (1|id), data=target1, subset=fst_lp<1.2), cor=FALSE)

## using likelihood ratio tests
print(model.ML.all <- lmer(log(gd) ~ tpt + unr + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1), cor=FALSE)
print(model.ML.tpt <- lmer(log(gd) ~ unr + idt + (1|sn) + (1|id), REML=FALSE, data=sub.target1), cor=FALSE)
anova(model.ML.tpt, model.ML.all)   # significant drop in goodness of fit when leaving out tpt = 0.03378 *