## @knitr Options, include=F, cache=F, message=F
opts_chunk$set(size='footnotesize', eval=T, cache=T, prompt=F, message=F, tidy=F, fig.width=4, fig.height=4)


## @knitr LinearMixedModelSimulation1, echo=T, results='markup', cache=T
## simulate some data
sigma.b <- 125       # inter-subject variation larger than
sigma.e <- 40        # intra-subject, inter-trial variation
alpha <- 500
beta <- 12
M <- 6                            # number of participants
n <- 50                           # trials per participant
b <- rnorm(M, 0, sigma.b)         # individual differences
nneighbors <- rpois(M*n,3) + 1    # generate num. neighbors
subj <- rep(1:M,n)
RT <- alpha + beta * nneighbors + # simulate RTs!
  b[subj] + rnorm(M*n,0,sigma.e)  #


## @knitr LinearMixedModel1, echo=T, results='markup', cache=T
library(languageR)
library(lme4)
data(lexdec)
l1 = lmer(RT ~ Frequency + (1 | Subject), data=lexdec)


## @knitr LinearMixedModel2, echo=F, results='markup', cache=T, dependson='LinearMixedModel1'
l1


## @knitr LinearMixedModelMCMC, echo=T, results='markup', cache=T, dependson='LinearMixedModel1'
library(languageR)
p = pvals.fnc(l1, nsim = 10000, addPlot=F)
print(p)


## @knitr LinearMixedModelR2, echo=T, results='markup', dependson='LinearMixedModel1'
cor(fitted(l1), lexdec$RT)^2


## @knitr LinearMixedModelMCMC2, echo=F, results='hide', cache=T, dependson='LinearMixedModel1', fig.width=3, fig.height=3
library(languageR)
pvals.fnc(l1, nsim = 10000)


## @knitr Ranef, echo=T, results='markup', dependson='LinearMixedModel1'
mean(ranef(l1)$Subject[[1]])
head(ranef(l1)$Subject[1], 5)


## @knitr RanefPredictedVsObserved, echo=F, results='markup', fig.width=3, fig.height=3, dependson='LinearMixedModel1'
# un-log RTs
par(cex=.7)
o = exp(as.vector(tapply(lexdec$RT, lexdec$Subject, mean)))
plot(o, type = "n", xlab = "Subject Index", 
  ylab = "Response latency in msecs", 
	main = "Subject as random effect")
# place by-subject mean on plot
text(o, as.character(unique(lexdec$Subject)), col = "blue")
# add overall mean
abline(mean(exp(lexdec$RT)),0, col = "gray", lwd=2)

# un-log model's predictions
p = as.numeric(unlist(exp(coef(l1)$Subject[1] + mean(lexdec$Frequency)*coef(l1)$Subject[2])))
# place model's predictions on plot
text(p, as.character(unique(lexdec$Subject)), col = "red")
legend(x=2, y=850, 
       legend=c("Predicted", "Observed"), 
       col=c("red", "blue"), pch=1
)


## @knitr LinearMixedModel3, echo=T, results='markup'
l2 = lmer(RT ~ 1 + (1 | Subject) + (1 | Word), 
          data = lexdec)
head(ranef(l2)$Subject, 5)
head(ranef(l2)$Word, 5)


## @knitr RanefPredictedVsObserved2, echo=F, results='markup', fig.width=5, fig.height=3.8
# un-log RTs
par(cex=.7)
s = exp(as.vector(tapply(lexdec$RT, lexdec$Subject, mean)))
o = exp(as.vector(tapply(lexdec$RT, lexdec$Word, mean)))
plot(o, ylim = c(min(s,o),max(s,o)), type = "n", xlab = "Word Index", 
  ylab = "Response latency in msecs", 
  main = "Word as random effect")
# place by-subject mean on plot
text(o, as.character(unique(lexdec$Word)), col = "blue")
# add overall mean
abline(mean(exp(lexdec$RT)),0, col = "gray", lwd=2)

# un-log model's predictions
p = as.numeric(unlist(exp(coef(l2)$Word[1])))
# place model's predictions on plot
segments(x0=1:length(o), y0=o, x1=1:length(p), 
  y1=p, col="red", lwd=2)
points(p, pch=ifelse(p > o, "^", ifelse(p < o, "v", ".")), col="red", lwd=2)
legend(x=2, y=850, 
       legend=c("Predicted", "Observed"), 
       col=c("red", "blue"), pch=1
)


## @knitr RandomSlopes1, echo=T, results='markup'
l3 = lmer(RT ~ 1 + Frequency + (1 + Frequency | Subject) + (1 | Word), 
            data = lexdec)


## @knitr RandomSlopes1b, echo=T, results='markup'
print(l3, corr=F)


## @knitr RandomSlopes2, echo=T, results='markup'
head(ranef(l3)$Subject, 5)
head(ranef(l3)$Word, 5)


## @knitr RandomSlopes4, echo=T, eval=FALSE
## dotplot(ranef(l3))


## @knitr RandomSlopes5, echo=T, eval=FALSE
## qqmath(ranef(l3))


## @knitr RandomCorrelations1, echo=T, results='markup', fig.width=2, fig.height=2
par(cex=.5)
plot(ranef(l3)$Subject, main="Random Effect Correlation")


## @knitr RandomCorrelations1b, echo=T, results='markup', fig.width=2, fig.height=2
lexdec$cFrequency = lexdec$Frequency - mean(lexdec$Frequency)
l3b = lmer(RT ~ 1 + cFrequency + (1 + cFrequency | Subject),
          data= lexdec)
print(l3b, corr=F)


## @knitr RandomCorrelations2, echo=T, results='markup', fig.width=2, fig.height=2
lexdec$HighFrequency = 
    ifelse(lexdec$Frequency > median(lexdec$Frequency), 1, 0)
l4 = lmer(RT ~ 1 + HighFrequency + (1 + HighFrequency | Subject),
          data= lexdec)
print(l4, corr=F)


## @knitr RandomCorrelations2b, echo=F, results='markup', fig.width=2, fig.height=2
par(cex=.5)
plot(ranef(l4)$Subject, main="Random Effect Correlation")


## @knitr RandomCorrelations3, echo=T, results='markup', fig.width=2, fig.height=2
lexdec$HighFrequency = ifelse(lexdec$Frequency > median(lexdec$Frequency), 
          .5, -.5)
l4b = lmer(RT ~ 1 + HighFrequency + (1 + HighFrequency | Subject),
          data= lexdec)
print(l4b, corr=F)
par(cex=.5)
plot(ranef(l4b)$Subject, main="Random Effect Correlation")


## @knitr RandomCorrelations3b, echo=F, results='markup', fig.width=2, fig.height=2
par(cex=.5)
plot(ranef(l4b)$Subject, main="Random Effect Correlation")


## @knitr RandomEffectExtraction, echo=T, results='markup'
# use vcov(model) for variance-covariance matrix of fixed effects
# these determiner the standard errors of the fixed effects.
vcov(l4b)

# get standard errors of fixed effects:
sqrt(diag(vcov(l4b)))


## @knitr RandomEffectExtraction2, echo=T, results='markup'
# use VarCorr(model) for variance-covariance matrix of random effects
# unlike the fixed effect variance-covariance matrix, only the *structure*
# of the VarCorr is known, the values are parameters that are estimated
# when we fit the model. 
VarCorr(l4b)

# for further details, see 
# https://stat.ethz.ch/pipermail/r-help/2006-July/109308.html


## @knitr RandomEffectExtraction3, echo=T, results='hide', dependson='LinearMixedModel3'
s = mcmcsamp(l2, 50000)

# a package for output analysis of MCMC
library(coda)
HPDinterval(s)


## @knitr L2repeated, echo=F, results='markup'
l2


## @knitr RandomEffectExtraction4, echo=F, results='markup', dependson='RandomEffectExtraction3'
HPDinterval(s)

# The ST slot contains estimates that determine the variance-covariance
# matrix of the random effects (these are the parameters of the Cholesky 
# decomposition; see next slide)
# (to understand what this all means, see # http://sciencemeanderthal.
# wordpress.com/2012/06/28/cholesky-decomposition-of-variance-covariance-
# matrices-in-the-classic-twin-study/)


## @knitr RandomEffectExtraction5, echo=T, results='markup', dependson='RandomEffectExtraction3'
HPDinterval(VarCorr(s, type = "varcov"))

# (given in the order they appear in l2: variance for Word intercept, 
# Subject intercept, Residual)