# Model Residuals and Powertransformation Stuff
# for Bodner & Masson (1997, Exp1 + Exp2a) + Kinoshita (2006, Exp. 2) data

#getwd()
#setwd( ...) 

rm(list=ls())

library(ggplot2)
library(MASS)

load("MaskedPriming.rda")

d.rs <- melt(d, id=c("e", "subjects", "words", "f", "p"), measure=c("srt", "lrt", "rrt"), variable_name="RT")

# WITHIN-CELL DEVIATIONS
source("functions/parseq.R")

# ... FOR SUBJECTS
d$srt.id <- NA
d$lrt.id <- NA
d$rrt.id <- NA

d <- d[order(d$e, d$subjects, d$f, d$p),]  # e is between stimulus factor
ix <- parseq(as.integer(d$p))
nc <- nrow(ix)

for (i in 1:nc ) {
		d$srt.id[ix[i,1]:ix[i,2]] <- scale(d$srt[ix[i,1]:ix[i,2]], scale=FALSE)
		d$lrt.id[ix[i,1]:ix[i,2]] <- scale(d$lrt[ix[i,1]:ix[i,2]], scale=FALSE)
		d$rrt.id[ix[i,1]:ix[i,2]] <- scale(d$rrt[ix[i,1]:ix[i,2]], scale=FALSE)
}

# durations
Msrt.id <- tapply(d$srt,    list(d$p, d$f, d$subjects), mean)
Vsrt.id <- tapply(d$srt.id, list(d$p, d$f, d$subjects), var)

# log durations
Mlrt.id <- tapply(d$lrt,    list(d$p, d$f, d$subjects), mean)
Vlrt.id <- tapply(d$lrt.id, list(d$p, d$f, d$subjects), var)

# reciprocal transformation of time to speed, i.e., -1/s
Mrrt.id <- tapply(d$rrt,    list(d$p, d$f, d$subjects), mean)
Vrrt.id <- tapply(d$rrt.id, list(d$p, d$f, d$subjects), var)

# ... FOR STIMULI
d$srt.st <- NA
d$lrt.st <- NA
d$rrt.st <- NA

d <- d[order(d$f, d$words, d$e, d$p),]  # f is between stimulus factor
ix <- parseq(as.integer(d$p))
nc <- nrow(ix)

for (i in 1:nc ) {
		d$srt.st[ix[i,1]:ix[i,2]] <- scale(d$srt[ix[i,1]:ix[i,2]], scale=FALSE)
		d$lrt.st[ix[i,1]:ix[i,2]] <- scale(d$lrt[ix[i,1]:ix[i,2]], scale=FALSE)
		d$rrt.st[ix[i,1]:ix[i,2]] <- scale(d$rrt[ix[i,1]:ix[i,2]], scale=FALSE)
}

# durations
Msrt.st <- tapply(d$srt,    list(d$p, d$e, d$words), mean)
Vsrt.st <- tapply(d$srt.st, list(d$p, d$e, d$words), var)

# log durations
Mlrt.st <- tapply(d$lrt,    list(d$p, d$e, d$words), mean)
Vlrt.st <- tapply(d$lrt.st, list(d$p, d$e, d$words), var)

# reciprocal transformation of time to speed, i.e., -1/s
Mrrt.st <- tapply(d$rrt,    list(d$p, d$e, d$words), mean)
Vrrt.st <- tapply(d$rrt.st, list(d$p, d$e, d$words), var)

#pdf(file="~/documents/projects/masson/kmr08/FigA1.08.pdf",paper="letter")
par(mfrow = c(3,2))
plot(Msrt.id, Vsrt.id,
  xlab="Cell means",
  ylab = "Cell variances",
  main="Subject Cells: RT")
plot(Msrt.st, Vsrt.st,
  xlab="Cell means",
  ylab = "Cell variances",
  main="Word Cells: RT")

plot(Mlrt.id, Vlrt.id,
  xlab="Cell means",
  ylab = "Cell variances",
  main="Subject Cells: log(RT)")
plot(Mlrt.st, Vlrt.st,
  xlab="Cell means",
  ylab = "Cell variances",
  main="Word Cells: log(RT)")

plot(Mrrt.id, Vrrt.id,
  xlab="Cell means",
  ylab = "Cell variances",
  main="Subject Cells: -1/RT")
plot(Mrrt.st, Vrrt.st,
  xlab="Cell Means",
  ylab = "Cell variances",
  main="Word Cells: -1/RT")

#dev.off()

#############################
# Figure 2 - Model residuals
#load("models.rda")
library(lme4)
print(m2.srt <- lmer(srt ~ 1 + p*f*e + (p + f | subjects) + (1 | words), data=d), cor=FALSE)
print(m2.lrt <- lmer(lrt ~ 1 + p*f*e + (p + f | subjects) + (1 | words), data=d), cor=FALSE)
print(m2 <- lmer(rrt ~ 1 + p*f*e + (p + f | subjects) + (1 | words), data=d), cor=FALSE)

srt <- data.frame(res = resid(m2.srt), fit=scale(fitted(m2.srt)))
srt$DV <- 1

lrt <- data.frame(res = resid(m2.lrt), fit=scale(fitted(m2.lrt)))
lrt$DV <- 2

rrt <- data.frame(res = resid(m2), fit=scale(fitted(m2)))
rrt$DV <- 3

xx <- rbind(srt,lrt,rrt)
xx$DV <- factor(xx$DV)
levels(xx$DV) <- c("RT", "log(RT)", "-1/RT")

p <-qplot(x=fit, y=res, data=xx, geom="point",  shape=I("."), facets= . ~  DV, 
      xlab="Fitted values", ylab="Standardized residuals") 
p  + geom_hline(yintercept=0) + theme_bw() 
#ggsave(file="figures/Fig2.08.pdf")

##################################
# Figure 3 - Box-Cox power transformation
#pdf(file="~/documents/projects/masson/kmr08/Fig3.08.pdf",paper="letter")
par(mfrow = c(1,1))
boxcox( srt ~ f*p*subjects, data=d, lambda=seq(-0.8, -1.2, length=31 ) )   # best lambda for subjects: -1.01