# ANALYSES for BoundaryN2-Experiments (Manuscript):
# EXP.1: n+2 frequency preview manipulation
# EXP.2: n+2 alternating case manipulation


# Sarah Risse, October 2011

rm(list=ls())

# SET WORKING DIRECTORY:
#setwd("")

# INPUT:
ifile_em   <- c("n2AC_em_filtered.rda")

# OUTPUT:


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

# FUNCTIONS:
source("remef.R")

# LOAD DATA:
load(ifile_em)
ls()			# em_n0, em_n1, em_n2, em_n3, em_nm1



# (1) COMPUTE TABLE


###############################
# WORD N
###############################

xdat <- em_n0

#---------------------------------------
# ANALYSES WITHOUT n2ab:
#----------------------------------------


# (1) get residuals of data without random effects (use remef.R):
lm1 <- lmer(log(gzd) ~ (sn1+lxn1+n2bb)^3 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n0gzd <- lm1@frame
n0gzd$gzd.adj <- exp(remef(lm1, fix = c(), ran = c("id", "wid", "sn")))


lm2 <- lmer(log(ffd) ~ (sn1+lxn1+n2bb)^3 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n0ffd <- lm2@frame
n0ffd$ffd.adj <- exp(remef(lm2, fix = c(), ran = c("id", "wid", "sn")))


lm3 <- lmer(log(sfd) ~ (sn1+lxn1+n2bb)^3 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n0sfd <- lm3@frame
n0sfd$sfd.adj <- exp(remef(lm3, fix = c(), ran = c("id", "wid", "sn")))


lm4 <- lmer(psk ~ (sn1+lxn1+n2bb)^3 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE, family=binomial) 

n0psk <- lm4@frame


# AGGREGATION OF DATA:
# GZD ... descriptives/ effect sizes:
n0gzd.rs <- melt(n0gzd, id=c("id","lxn1","n2bb","sn1")
		, measure=c("gzd.adj")
		, na.rm=TRUE)

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

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

n0gzd.rs$n2bb <- as.factor(n0gzd.rs$n2bb)
levels(n0gzd.rs$n2bb) <- c("easy","diff.")


# FFD .... descriptives/ effect sizes:
n0ffd.rs <- melt(n0ffd, id=c("id","lxn1","n2bb","sn1")
		, measure=c("ffd.adj")
		, na.rm=TRUE)

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

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

n0ffd.rs$n2bb <- as.factor(n0ffd.rs$n2bb)
levels(n0ffd.rs$n2bb) <- c("easy","diff.")


# SFD ... descriptives/ effect sizes:
n0sfd.rs <- melt(n0sfd, id=c("id","lxn1","n2bb","sn1")
		, measure=c("sfd.adj")
		, na.rm=TRUE)

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

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

n0sfd.rs$n2bb <- as.factor(n0sfd.rs$n2bb)
levels(n0sfd.rs$n2bb) <- c("easy","diff.")


# PSK ... descriptives/ effect sizes:
n0psk.rs <- melt(n0psk, id=c("id","lxn1","n2bb","sn1")
		, measure=c("psk")
		, na.rm=TRUE)

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

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

n0psk.rs$n2bb <- as.factor(n0psk.rs$n2bb)
levels(n0psk.rs$n2bb) <- c("easy","diff.")


# TABLE 5
mw <- data.frame(cast(n0gzd.rs, n2bb ~ sn1+lxn1
	, subset=variable==c("gzd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0gzd.rs, n2bb ~ sn1+lxn1
	, subset=variable==c("gzd.adj")
	, function(x) c(SD=signif(sd(x),2))))
mw


mw <- data.frame(cast(n0ffd.rs, n2bb ~ sn1+lxn1
	, subset=variable==c("ffd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0ffd.rs, n2bb ~ sn1+lxn1
	, subset=variable==c("ffd.adj")
	, function(x) c(SD=signif(sd(x),2))))
mw


mw <- data.frame(cast(n0sfd.rs, n2bb ~ sn1+lxn1
	, subset=variable==c("sfd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0sfd.rs, n2bb ~ sn1+lxn1
	, subset=variable==c("sfd.adj")
	, function(x) c(SD=signif(sd(x),2))))
mw


mw <- data.frame(cast(n0psk.rs, n2bb ~ sn1+lxn1
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0psk.rs, n2bb ~ sn1+lxn1
	, function(x) c(SD=signif(sd(x),2))))
mw

# INDIVIDUAL EFFECTS:
mw <- data.frame(cast(n0gzd.rs, lxn1 ~ variable
	, subset=variable==c("gzd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0ffd.rs, lxn1 ~ variable
	, subset=variable==c("ffd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0sfd.rs, lxn1 ~ variable
	, subset=variable==c("sfd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0gzd.rs, sn1 ~ variable
	, subset=variable==c("gzd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0ffd.rs, sn1 ~ variable
	, subset=variable==c("ffd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n0sfd.rs, sn1 ~ variable
	, subset=variable==c("sfd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw





###############################
# WORD N+1
###############################

xdat <- em_n1


# (1) get residuals of data without random effects (use remef.R):
lm1 <- lmer(log(gzd) ~ (lxn1+n2bb+n2ab)^3 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n1gzd <- lm1@frame
n1gzd$gzd.adj <- exp(remef(lm1, fix = c(), ran = c("id", "wid", "sn")))


lm2 <- lmer(log(ffd) ~ (lxn1+n2bb+n2ab)^3 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n1ffd <- lm2@frame
n1ffd$ffd.adj <- exp(remef(lm2, fix = c(), ran = c("id", "wid", "sn")))


lm3 <- lmer(log(sfd) ~ (lxn1+n2bb+n2ab)^3 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n1sfd <- lm3@frame
n1sfd$sfd.adj <- exp(remef(lm3, fix = c(), ran = c("id", "wid", "sn")))


lm4 <- lmer(psk ~ (lxn1+n2bb+n2ab)^3 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE, family=binomial) 

n1psk <- lm4@frame




# AGGREGATION OF DATA:
# GZD ... descriptives/ effect sizes:
n1gzd.rs <- melt(n1gzd, id=c("id","lxn1","n2bb","n2ab")
		, measure=c("gzd.adj")
		, na.rm=TRUE)

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

n1gzd.rs$n2bb <- as.factor(n1gzd.rs$n2bb)
levels(n1gzd.rs$n2bb) <- c("easy","diff.")

n1gzd.rs$n2ab <- as.factor(n1gzd.rs$n2ab)
levels(n1gzd.rs$n2ab) <- c("easy","diff.")


# FFD .... descriptives/ effect sizes:
n1ffd.rs <- melt(n1ffd, id=c("id","lxn1","n2bb","n2ab")
		, measure=c("ffd.adj")
		, na.rm=TRUE)

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

n1ffd.rs$n2bb <- as.factor(n1ffd.rs$n2bb)
levels(n1ffd.rs$n2bb) <- c("easy","diff.")

n1ffd.rs$n2ab <- as.factor(n1ffd.rs$n2ab)
levels(n1ffd.rs$n2ab) <- c("easy","diff.")


# SFD ... descriptives/ effect sizes:
n1sfd.rs <- melt(n1sfd, id=c("id","lxn1","n2bb","n2ab")
		, measure=c("sfd.adj")
		, na.rm=TRUE)

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

n1sfd.rs$n2bb <- as.factor(n1sfd.rs$n2bb)
levels(n1sfd.rs$n2bb) <- c("easy","diff.")

n1sfd.rs$n2ab <- as.factor(n1sfd.rs$n2ab)
levels(n1sfd.rs$n2ab) <- c("easy","diff.")


# PSK ... descriptives/ effect sizes:
n1psk.rs <- melt(n1psk, id=c("id","lxn1","n2bb","n2ab")
		, measure=c("psk")
		, na.rm=TRUE)

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

n1psk.rs$n2bb <- as.factor(n1psk.rs$n2bb)
levels(n1psk.rs$n2bb) <- c("easy","diff.")

n1psk.rs$n2ab <- as.factor(n1psk.rs$n2ab)
levels(n1psk.rs$n2ab) <- c("easy","diff.")



# TABLE 6
mw <- data.frame(cast(n1gzd.rs, n2bb+n2ab ~ lxn1
	, subset=variable==c("gzd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n1gzd.rs, n2bb+n2ab ~ lxn1
	, subset=variable==c("gzd.adj")
	, function(x) c(SD=signif(sd(x),2))))
mw


mw <- data.frame(cast(n1ffd.rs, n2bb+n2ab ~ lxn1
	, subset=variable==c("ffd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n1ffd.rs, n2bb+n2ab ~ lxn1
	, subset=variable==c("ffd.adj")
	, function(x) c(SD=signif(sd(x),2))))
mw


mw <- data.frame(cast(n1sfd.rs, n2bb+n2ab ~ lxn1
	, subset=variable==c("sfd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n1sfd.rs, n2bb+n2ab ~ lxn1
	, subset=variable==c("sfd.adj")
	, function(x) c(SD=signif(sd(x),2))))
mw

mw <- data.frame(cast(n1psk.rs, n2bb+n2ab ~ lxn1
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n1psk.rs, n2bb+n2ab ~ lxn1
	, function(x) c(SD=signif(sd(x),2))))
mw


# INDIVIDUAL EFFECTS:
mw <- data.frame(cast(n1gzd.rs, n2bb ~ variable
	, subset=variable==c("gzd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n1ffd.rs, n2bb ~ variable
	, subset=variable==c("ffd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n1sfd.rs, n2bb ~ variable
	, subset=variable==c("sfd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

###############################
# WORD N+2
###############################

xdat <- em_n2

# (1) get residuals of data without random effects (use remef.R):
lm1 <- lmer(log(gzd) ~ (sn1+lxn1+n2bb+n2ab)^4 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n2gzd <- lm1@frame
n2gzd$gzd.adj <- exp(remef(lm1, fix = c(), ran = c("id", "wid", "sn")))


lm2 <- lmer(log(ffd) ~ (sn1+lxn1+n2bb+n2ab)^4 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n2ffd <- lm2@frame
n2ffd$ffd.adj <- exp(remef(lm2, fix = c(), ran = c("id", "wid", "sn")))


lm3 <- lmer(log(sfd) ~ (sn1+lxn1+n2bb+n2ab)^4 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE) 

n2sfd <- lm3@frame
n2sfd$sfd.adj <- exp(remef(lm3, fix = c(), ran = c("id", "wid", "sn")))


lm4 <- lmer(psk ~ (sn1+lxn1+n2bb+n2ab)^4 + (1|id) + (1|wid) + (1|sn)
, data=xdat, REML=TRUE, family=binomial) 

n2psk <- lm4@frame





# AGGREGATION OF DATA:
# GZD ... descriptives/ effect sizes:
n2gzd.rs <- melt(n2gzd, id=c("id","lxn1","n2bb","n2ab","sn1")
		, measure=c("gzd.adj")
		, na.rm=TRUE)

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

n2gzd.rs$n2bb <- as.factor(n2gzd.rs$n2bb)
levels(n2gzd.rs$n2bb) <- c("easy","diff.")

n2gzd.rs$n2ab <- as.factor(n2gzd.rs$n2ab)
levels(n2gzd.rs$n2ab) <- c("easy","diff.")

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


# FFD .... descriptives/ effect sizes:
n2ffd.rs <- melt(n2ffd, id=c("id","lxn1","n2bb","n2ab","sn1")
		, measure=c("ffd.adj")
		, na.rm=TRUE)

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

n2ffd.rs$n2bb <- as.factor(n2ffd.rs$n2bb)
levels(n2ffd.rs$n2bb) <- c("easy","diff.")

n2ffd.rs$n2ab <- as.factor(n2ffd.rs$n2ab)
levels(n2ffd.rs$n2ab) <- c("easy","diff.")

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


# SFD ... descriptives/ effect sizes:
n2sfd.rs <- melt(n2sfd, id=c("id","lxn1","n2bb","n2ab","sn1")
		, measure=c("sfd.adj")
		, na.rm=TRUE)

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

n2sfd.rs$n2bb <- as.factor(n2sfd.rs$n2bb)
levels(n2sfd.rs$n2bb) <- c("easy","diff.")

n2sfd.rs$n2ab <- as.factor(n2sfd.rs$n2ab)
levels(n2sfd.rs$n2ab) <- c("easy","diff.")

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



# PSK ... descriptives/ effect sizes:
n2psk.rs <- melt(n2psk, id=c("id","lxn1","n2bb","n2ab","sn1")
		, measure=c("psk")
		, na.rm=TRUE)

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

n2psk.rs$n2bb <- as.factor(n2psk.rs$n2bb)
levels(n2psk.rs$n2bb) <- c("easy","diff.")

n2psk.rs$n2ab <- as.factor(n2psk.rs$n2ab)
levels(n2psk.rs$n2ab) <- c("easy","diff.")

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


# TABLE 7
mw <- data.frame(cast(n2gzd.rs, n2bb+n2ab ~ sn1+lxn1
	, subset=variable==c("gzd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n2gzd.rs, n2bb+n2ab ~ sn1+lxn1
	, subset=variable==c("gzd.adj")
	, function(x) c(SD=signif(sd(x),3))))
mw


mw <- data.frame(cast(n2ffd.rs, n2bb+n2ab ~ sn1+lxn1
	, subset=variable==c("ffd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n2ffd.rs, n2bb+n2ab ~ sn1+lxn1
	, subset=variable==c("ffd.adj")
	, function(x) c(SD=signif(sd(x),3))))
mw


mw <- data.frame(cast(n2sfd.rs, n2bb+n2ab ~ sn1+lxn1
	, subset=variable==c("sfd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n2sfd.rs, n2bb+n2ab ~ sn1+lxn1
	, subset=variable==c("sfd.adj")
	, function(x) c(SD=signif(sd(x),3))))
mw


mw <- data.frame(cast(n2psk.rs, n2bb+n2ab ~ sn1+lxn1
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n2psk.rs, n2bb+n2ab ~ sn1+lxn1
	, function(x) c(SD=signif(sd(x),3))))
mw


# INDIVIDUAL EFFECTS:
mw <- data.frame(cast(n2gzd.rs, n2ab ~ variable
	, subset=variable==c("gzd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n2ffd.rs, n2ab ~ variable
	, subset=variable==c("ffd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw

mw <- data.frame(cast(n2sfd.rs, n2ab ~ variable
	, subset=variable==c("sfd.adj")
	, function(x) c(M=signif(mean(x),3))))
mw