# Estimating a between-subject variance in ANOVA and LMM
# Compute shrinkage factor

#getwd()
#setwd( ...) 


library(reshape)
library(lme4)

rm(list=ls())
load("MaskedPriming.rda")
d.rs <- melt(d, id=c("subjects", "words", "e", "f", "p"),
                measure=c("srt", "lrt", "rrt"), variable="RT")

# ...subjects
d.id <- cast(d.rs, e+subjects+f+p ~ RT, mean)

# ##########################################
# ANOVA
op <- options(contrasts=c("contr.sum", "contr.poly"))
summary(m1 <- aov(rrt ~ 1 + subjects, data=d.id))
round(coef(m1))

# ... fixed-effects predictions
(feID <- unlist(model.tables(m1)[1]))

# ... estim for between-id variance: (MSid - MSerror)/Nt --from EMS stuff
# (SigmaID <- sqrt( (summary(m1)[[1]][1, 3] - summary(m1)[[1]][2, 3])/96 ))
(MSA <- summary(m1)[[1]][1, 3])
(MSE <- summary(m1)[[1]][2, 3])  # == SigmaError2
(SigmaID <- sqrt((MSA - MSE)/4))
(SigmaError <- sqrt(MSE))

############################################
# LMM
(m2 <- lmer(rrt ~ (1|subjects), data=d.id))
(reID.2 <- ranef(m2, postVar=TRUE))

# "Shrinkage"
# ... constant shrinkage factor for balanced design
feID/unlist(reID.2)
quartz()
plot(feID ~ unlist(reID.2))
abline(0,1)

# ... prediction intervals for random effects
dotplot(reID.2, scales = list(x = list(relation = 'free', rot=0, cex=1), y = list(cex=1) ) )

# ... combine various subject means in wide format and order by conditional means (reID.2)
xx <- data.frame(cbind(unlist(feID), unlist(reID.2)))
names(xx) <- c("feID", "reID.2")
xx <- xx[order(xx$reID.2), ]   # alternative: xx <- xx[order(xx$trueID.c), ]

# ... convert to long format
x <- data.frame(id=factor(rep(1:72, 2)))
x$coef=factor(rep(1:2, each=72))
levels(x$coef) <- c("fixed effect", "random effect")
x$Mean <- c(xx$feID, xx$reID.2)

# ... plot
dotplot(id ~ Mean, data=x, group=coef, ref=TRUE, auto.key = list(columns=2), #auto.key=TRUE,  
	    scales = list(x = list(relation = 'free', rot=0, cex=1), y = list(cex=0) ), pch=16 )

############################################
# Differential shrinkage due to different n
# ... ... unbalanced N of trials
ix <- which(as.numeric(d.id$subjects) <= 10)
b <- d.id[-sample(ix, length(ix)*3/4), ]

(m3 <- lmer(rrt ~ (1|subjects), data=b))
(reID.3 <- ranef(m3, postVar=TRUE))

# "Shrinkage"
# ... prediction intervals for random effects larger for subjects with fewer n 
dotplot(reID.3, scales = list(x = list(relation = 'free', rot=0, cex=1), y = list(cex=1) ) )