## @knitr preamble, echo=FALSE, results='hide', cache=FALSE, message=F
library(plyr)
library(reshape2)

library(ggplot2)
library(stringr)

# set knitr chunk option defaults for this document
opts_chunk$set(results='markup', echo=TRUE, cache=T, autodep=T, warning=F, message=F, prompt=F,
               fig.width=4, fig.height=4, out.width='0.5\\textwidth', fig.pos='!h', fig.lp='fig:',
               tidy=F, size='scriptsize')

data(lexdec, package='languageR')

# set ggplot theme globally
theme_set(theme_bw())



## @knitr 
freqEffects <- ddply(lexdec, .(Subject), function(df) {coef(lm(RT~Frequency, data=df))})


## @knitr echo=F
head(freqEffects)


## @knitr lexdec-freq-effect-bysub, echo=F, out.width='0.9\\textwidth'
ggplot(lexdec, aes(x=Frequency, y=RT)) +
  geom_point() +
  facet_wrap(~Subject) + 
  geom_abline(data=freqEffects, aes(slope=Frequency, intercept=`(Intercept)`), color='red')


## @knitr 
mean(runif(100))


## @knitr 
f <- mean
f(runif(100))


## @knitr 
function(x) sum(x)/length(x)


## @knitr 
mymean <- function(x) sum(x)/length(x)
mymean(runif(100))
mymean(runif(100, 1, 2))


## @knitr 
(function(x) sum(x)/length(x)) (runif(100))


## @knitr 
make.power <- function(n) {
  return(function(x) x^n)
}
my.square <- make.power(2)
my.square(3)
(make.power(4)) (2)


## @knitr 
sd(rnorm(sd=5, mean=1, 100))


## @knitr 
rnorm


## @knitr 
list.o.nums <- list(runif(100), rnorm(100), rpois(100, lambda=1))
lapply(list.o.nums, mean)


## @knitr 
sapply(1:10, function(n) rpois(5, lambda=n))


## @knitr 
sapply(1:10, function(n) mean(rnorm(n=5, mean=0, sd=1)))


## @knitr 
data(lexdec, package='languageR')
RT.bysub <- with(lexdec, split(RT, Subject))
RT.means.bysub <- sapply(RT.bysub, mean)
head(data.frame(RT.mean=RT.means.bysub))


## @knitr 
library(plyr) 
head(ddply(lexdec, .(Subject), function(df) data.frame(RT.mean=mean(df$RT)))) 


## @knitr eval=F
ddply(lexdec, .(Subject), function(df) data.frame(RT.mean=mean(df$RT)))


## @knitr ddply-output-exercises, eval=F
library(plyr)
data(lexdec, package='languageR')       # load the dataset if it isn't already
ddply(lexdec, .(PrevType, Class), function(df) with(df, data.frame(meanRT=mean(RT))))


## @knitr 
slowestTrials <- ddply(lexdec, .(Subject), subset, RT==max(RT))
head(slowestTrials[, c('Subject', 'RT', 'Trial', 'Word')])


## @knitr eval=FALSE
## ddply(lexdec, .(Subject), function(df, ...) subset(df, ...), RT==max(RT))
## ddply(lexdec, .(Subject), function(df) subset(df, RT==max(RT)))


## @knitr eval=F
ddply(lexdec, .(Subject), function(df) {
  df$RT.s <- scale(df$RT)
  return(df)
})


## @knitr 
lexdecScaledRT <- ddply(lexdec, .(Subject), transform, RT.s=scale(RT))


## @knitr summarise-RT-dist-plots, out.width='0.45\\textwidth', echo=FALSE, results='hide'
ggplot(lexdecScaledRT, aes(x=RT, group=Subject)) + geom_density()
ggplot(lexdecScaledRT, aes(x=RT.s, group=Subject)) + geom_density()


## @knitr eval=FALSE
## summarise(.data, summVar1=expr1, summVar2=expr2, ...)


## @knitr 
lexdec.RTsumm <- ddply(lexdec, .(Subject), summarise, mean=mean(RT), var=var(RT))
head(lexdec.RTsumm)


## @knitr eval=F
ddply(lexdec, .(Subject), function(df) with(df, data.frame(meanRT=mean(RT))))


## @knitr 
lexdec <- ddply(lexdec, .(Subject), transform, RTrank=rank(RT))


word.rt.ranks <- ddply(lexdec, .(Word, Frequency), summarise,
                       RTrank=mean(RTrank),
                       RTrank25=quantile(RTrank, 0.25), RTrank75=quantile(RTrank, 0.75))


## @knitr out.width='0.9\\textwidth'
ggplot(word.rt.ranks, aes(x=Frequency, y=RTrank, ymin=RTrank25, ymax=RTrank75)) +
  geom_pointrange()


## @knitr 
subset(word.rt.ranks,
       RTrank %in% c(max(RTrank), min(RTrank)))


## @knitr 
ddply(lexdec, .(NativeLanguage), summarise, acc=mean(Correct=='correct'))


## @knitr 
daply(lexdec, .(NativeLanguage, Correct), nrow)


## @knitr 
correctCounts <- daply(lexdec, .(NativeLanguage, Correct), nrow)
chisq.test(correctCounts)


## @knitr 
subjectSlopes <- ddply(lexdec, .(Subject), function(df) {coef(lm(RT~Frequency, data=df))})


## @knitr 
summary(subjectSlopes$Frequency)


## @knitr subject-slopes-rt
ggplot(lexdec, aes(x=Frequency, y=RT)) +
  geom_point() +
  facet_wrap(~Subject) + 
  geom_abline(data=subjectSlopes, aes(slope=Frequency, intercept=`(Intercept)`), color='red')


## @knitr 
lexdec <- transform(lexdec,
                    FreqHiLo=factor(ifelse(Frequency>median(Frequency), 
                                           'high', 'low'),
                                    levels=c('low', 'high')))


## @knitr results='hide', echo=FALSE
# this is necessary because of something strange in how knitr/plyr are
# using environments
se <- function(x) sd(x)/sqrt(length(x))


## @knitr 
se <- function(x) sd(x)/sqrt(length(x))
langfreq.summ <- ddply(lexdec,
                       .(NativeLanguage, FreqHiLo),
                       summarise, mean=mean(RT), se=se(RT))


## @knitr fig.width=6, out.width='0.75\\textwidth'
ggplot(langfreq.summ, aes(x=FreqHiLo, color=NativeLanguage,
                          y=mean, ymin=mean-1.96*se, ymax=mean+1.96*se)) +
  geom_point() +
  geom_errorbar(width=0.1) + 
  geom_line(aes(group=NativeLanguage))


## @knitr 
langfreq.bysub <- ddply(lexdec, .(NativeLanguage, FreqHiLo, Subject),
                        summarise, RT=mean(RT))


## @knitr 
langfreq.bysub.summ <- ddply(langfreq.bysub,
                             .(NativeLanguage, FreqHiLo),
                             summarise, mean=mean(RT), se=se(RT))


## @knitr by-subject-errorbars, fig.width=6, out.width='0.75\\textwidth'
ggplot(langfreq.bysub.summ, aes(x=FreqHiLo, color=NativeLanguage,
                                y=mean, ymin=mean-1.96*se, ymax=mean+1.96*se)) +
  geom_point() +
  geom_errorbar(width=0.1) + 
  geom_line(aes(group=NativeLanguage))


## @knitr 
langfreq.byitem <- ddply(lexdec, .(NativeLanguage, FreqHiLo, Word), summarise, RT=mean(RT))
langfreq.byitem.summ <- ddply(langfreq.byitem,
                             .(NativeLanguage, FreqHiLo),
                             summarise, mean=mean(RT), se=se(RT))

## ggplot(langfreq.byitem.summ, aes(x=FreqHiLo, color=NativeLanguage,
##                           y=mean, ymin=mean-2*se, ymax=mean+2*se)) +
##   geom_point() +
##   geom_errorbar(width=0.1) + 
##   geom_line(aes(group=NativeLanguage))


## @knitr echo=F, results='hide', fig.show='hide'
library(plyr)
library(mvtnorm)
library(lme4)

make.data.generator <- function(true.effects=c(0,0),
                                resid.var=1,
                                ranef.var=diag(c(1,1)),
                                n.subj=24,
                                n.obs=24
                                ) 
{   
  # create design matrix for our made up experiment
  data.str <- data.frame(freq=factor(c(rep('high', n.obs/2), rep('low', n.obs/2))))
  contrasts(data.str$freq) <- contr.sum(2)
  model.mat <- model.matrix(~ 1 + freq, data.str)
  
  generate.data <- function() {
    # sample data set
    simulated.data <- rdply(n.subj, {
      beta <- t(rmvnorm(n=1, sigma=ranef.var)) + true.effects
      expected.RT <- model.mat %*% beta
      epsilon <- rnorm(n=length(expected.RT), mean=0, sd=sqrt(resid.var))
      data.frame(data.str,
                 RT=expected.RT + epsilon)
    })
    names(simulated.data)[1] <- 'subject'
    simulated.data
  }
}

fit.models <- function(simulated.data) {
  # fit models and extract coefs
  lm.coefs <- coefficients(summary(lm(RT ~ 1+freq, simulated.data)))[, 1:3]
  rand.int.coefs <- summary(lmer(RT ~ 1+freq + (1|subject), simulated.data))@coefs
  rand.slope.coefs <- summary(lmer(RT ~ 1+freq + (1+freq|subject), simulated.data))@coefs
  # format output all pretty
  rbind(data.frame(model='lm', predictor=rownames(lm.coefs), lm.coefs),
        data.frame(model='rand.int', predictor=rownames(rand.int.coefs), rand.int.coefs),
        data.frame(model='rand.slope', predictor=rownames(rand.slope.coefs), rand.slope.coefs))
}


gen.dat <- make.data.generator()

simulations <- rdply(.n=100,
                     fit.models(gen.dat()),
                     .progress='text')

head(simulations, n=20)
sim.m <- melt(simulations, measure.var=c('Estimate', 'Std..Error', 't.value'))
head(sim.m)

ddply(simulations, .(model, predictor), summarise, type1err=mean(abs(t.value)>1.96))



## @knitr data-simulation-model-fitting-preview, fig.width=5, fig.height=3.5, out.width='\\textwidth', echo=F, results='hide'

ggplot(simulations, aes(x=t.value, color=model)) + 
  geom_vline(xintercept=c(-1.96, 1.96), color='#888888', linetype=3)  +
  scale_x_continuous('t value') + 
  geom_density() +
  facet_grid(predictor~.)



## @knitr eval=F
library(plyr)
library(mvtnorm)
library(lme4)

make.data.generator <- function(true.effects=c(0,0),
                                resid.var=1,
                                ranef.var=diag(c(1,1)),
                                n.subj=24,
                                n.obs=24
                                ) 
{   
  # create design matrix for our made up experiment
  data.str <- data.frame(freq=factor(c(rep('high', n.obs/2), rep('low', n.obs/2))))
  contrasts(data.str$freq) <- contr.sum(2)
  model.mat <- model.matrix(~ 1 + freq, data.str)
  
  generate.data <- function() {
    # sample data set under mixed effects model with random slope/intercepts
    simulated.data <- rdply(n.subj, {
      beta <- t(rmvnorm(n=1, sigma=ranef.var)) + true.effects
      expected.RT <- model.mat %*% beta
      epsilon <- rnorm(n=length(expected.RT), mean=0, sd=sqrt(resid.var))
      data.frame(data.str,
                 RT=expected.RT + epsilon)
    })
    names(simulated.data)[1] <- 'subject'
    simulated.data
  }
}


## @knitr eval=F
fit.models <- function(simulated.data) {
  # fit models and extract coefs
  lm.coefs <- coefficients(summary(lm(RT ~ 1+freq, simulated.data)))[, 1:3]
  rand.int.coefs <- summary(lmer(RT ~ 1+freq + (1|subject), simulated.data))@coefs
  rand.slope.coefs <- summary(lmer(RT ~ 1+freq + (1+freq|subject), simulated.data))@coefs
  # format output all pretty
  rbind(data.frame(model='lm', predictor=rownames(lm.coefs), lm.coefs),
        data.frame(model='rand.int', predictor=rownames(rand.int.coefs), rand.int.coefs),
        data.frame(model='rand.slope', predictor=rownames(rand.slope.coefs), rand.slope.coefs))
}


## @knitr eval=F
gen.dat <- make.data.generator()
simulations <- rdply(.n=100,
                     fit.models(gen.dat()),
                     .progress='text')


## @knitr 
head(simulations)
daply(simulations, .(model, predictor), function(df) type1err=mean(abs(df$t.value)>1.96))


## @knitr fig.width=5, fig.height=3.5, out.width='0.75\\textwidth'
# use reshape2::melt to get the data into a more convenient format (see next section)
ggplot(simulations, aes(x=t.value, color=model)) +
  geom_vline(xintercept=c(-1.96, 1.96), color='#888888', linetype=3) + 
  scale_x_continuous('t value') + 
  geom_density() +
  facet_grid(predictor~.)


## @knitr 
head(params <- expand.grid(n.obs=c(4, 16, 64), n.subj=c(4, 16, 64)))


## @knitr eval=F
man.simulations <- mdply(params, function(...) {
                           make.data <- make.data.generator(...)
                           rdply(.n=100, fit.models(make.data()))
                         }, .progress='text')


## @knitr eval=F
mdply(params, function(...) {
  make.data <- make.data.generator(...)
  rdply(.n=100, fit.models(make.data()))
}, .progress='text')


## @knitr eval=FALSE
## ddply(params, .(n.obs, n.subj), function(df) {
##   make.data <- make.data.generator(n.obs=df$n.obs, n.subj=df$n.subj)
##   rdply(.n=100, .fun=fit.models(make.data()))
## }, .progress='text')


## @knitr echo=F
ld.wide[1:7, 1:7]


## @knitr echo=F
ld.long[1:5, ]


## @knitr echo=F
ld.wide[1:5, 1:4]


## @knitr 
ld.wide[1:5, 1:7]


## @knitr 
head(ld.m <- melt(ld.wide, id.var='Subject', na.rm=T))


## @knitr 
require(stringr)
trials.and.vars <- ldply(str_split(ld.m$variable, '_'))
names(trials.and.vars) <- c('Trial', 'measure')


## @knitr 
head(ld.m <- cbind(ld.m, trials.and.vars))


## @knitr 
head(dcast(ld.m, Subject+Trial ~ measure))


## @knitr 
head(dcast(ld.m, ... ~ measure))


## @knitr 
ld.m$variable <- NULL
head(dcast(ld.m, ... ~ measure))


## @knitr 
head(dcast(ld.m, Subject ~ measure))


## @knitr results='hide'
melt(ld.wide, id.var='Subject', na.rm=T)


## @knitr results='hide'
ddply(ld.wide, .(Subject), function(df) {
  vars <- names(df)
  vals <- t(df)
  dimnames(vals) <- NULL
  return(subset(data.frame(variable=vars, value=vals),
                variable != 'Subject' & !is.na(value)))
})


## @knitr results='hide'
head(dcast(ld.m, Subject+Trial ~ measure))


## @knitr 
head(ddply(ld.m, .(Subject, Trial), function(df) {
  res <- data.frame(t(df$value))
  names(res) <- df$measure
  return(res)
}))