# Graphics for conditional means (BLUPs), shrinkage effects on correlations

#setwd( ...) 

library(MASS)
library(lattice)
library(reshape)
library(ggplot2)
library(Matrix)
library(lme4)

rm(list=ls())
load("MaskedPriming.rda")
# Final model for reciprocal RTs (-1/srt); alternative variants in p1.R
# ... Helmert contrast for experiments
cmat <- matrix(c(-1,  1,  0,
	             -1, -1,  2),  3,  2, 
	             dimnames=list(c("BM1", "BM2", "SK"), 
	                              c(".BM1-2", ".BM-SK")))
# ... LMM
contrasts(d$e) <- cmat
print(m2 <- lmer(rrt ~ p*f*e + (p + f | subjects) + (1 | words), data=d), cor=FALSE)

# Coefficients (i.e., fixed effect + random effect)
m2.coef <- unlist(coef(m2)[[2]][1:3])  # only id CMs: intercept, priming, frequency
dim(m2.coef) <- c(72, 3)
m2.coef <- data.frame(m2.coef)
names(m2.coef) <- c("Mean", "Priming", "Frequency")

# ... histograms of effects; leave out mean
m2.coef.long <- melt(m2.coef, measure=names(m2.coef)[2:3])
histogram( ~ value | variable, data = m2.coef.long,  
          xlab="Conditional Means", type = "density",
          panel = function(x, ...) {
              panel.histogram(x, ...)
              panel.mathdensity(dmath = dnorm, col = "black",
                                args = list(mean=mean(x),sd=sd(x)))
          } )

# ... splom of coefficiencts 
splom(~ m2.coef[,1:3])

# ... random effects sorted by intercept
m2.ranef <- ranef(m2, postVar=TRUE)   # see ?ranef for explanation of postVar
dotplot(m2.ranef,layout=c(3,1), par.settings = simpleTheme(col = 1), scales = list(x = list(relation = 'free', rot=0, cex=0.7), y = list(cex=0) ), 
        strip = strip.custom(bg = "grey90") )[[2]][c(1,2,3)]

# ... demonstrate shrinkage, compared to within-subject regression (adopted from Bates PotsdamGLMM/Longitudinal.Rnw)
# ... ... within-subject OLS estimates
df <- coef(lmList(rrt ~ f+p | subjects, data=d))[, 1:3]

# ... ... merge within-subject and mixed-model estimates
m2.coef <- as.data.frame(m2.coef)
m2.coef$id <- factor(1:72)

op <- options(digits=2)
cor(m2.coef)
sd(m2.coef)
options(op)

df <- cbind(df, m2.coef)

# Saved as figures/Fig5a.pdf
# ... ... mean rrt and priming effects
with(df,
     print(xyplot(p ~ `(Intercept)`, aspect = 1,
            x1 = Mean, y1 = Priming, xlab="Mean -1/RT", ylab="Priming Effect",
            panel = function(x, y, x1, y1, subscripts, ...) {
                panel.grid(h = -1, v = -1)
                x1 <- x1[subscripts]
                y1 <- y1[subscripts]
                panel.arrows(x, y, x1, y1, type = "open", length = 0.1, col="gray50", lty=1,
                              angle = 15, ...)
                panel.points(x, y, col="black")
                panel.points(x1, y1, col="black", cex=0.5, pch=19)
            } #, key = list(space = "top", columns = 2, text = list(c("Mixed model", "Within-subject")),
#                points = list(col = trellis.par.get("superpose.symbol")$col[1],
#                pch = trellis.par.get("superpose.symbol")$pch[1]))
               ), split = c(1, 1, 3, 1)))

# ... ... mean rrt and frequency effects
with(df,
     plot(xyplot(f ~ `(Intercept)`, aspect = 1, 
            x1 = Mean, y1 = Frequency, xlab="Mean -1/RT", ylab="Frequency Effect",
            panel = function(x, y, x1, y1, subscripts, ...) {
                panel.grid(h = -1, v = -1)
                x1 <- x1[subscripts]
                y1 <- y1[subscripts]
                panel.arrows(x, y, x1, y1, type = "open", length = 0.1, col="gray50", lty=1,
                              angle = 15, ...)
                panel.points(x, y, col="black")
                panel.points(x1, y1, col="black", cex=0.5, pch=19)
            }
               ), split = c(2, 1, 3, 1), newpage=FALSE) )

# Saved as figures/Fig5c.pdf
# ... ... priming and frequency effects
with(df,
     print(xyplot(f ~ p, aspect = 1,
            x1 = Priming, y1 = Frequency, xlab="Priming Effect", ylab="Frequency Effect",
            panel = function(x, y, x1, y1, subscripts, ...) {
                panel.grid(h = -1, v = -1)
                x1 <- x1[subscripts]
                y1 <- y1[subscripts]
                panel.arrows(x, y, x1, y1, type = "open", length = 0.1, col="gray50", lty=1,
                              angle = 15, ...)
                panel.points(x, y, col="black")
                panel.points(x1, y1, col="black", cex=0.5, pch=19)
            }
               ), split = c(3, 1, 3, 1), newpage=FALSE) )