library(languageR)

## ---------------------------------------------------
# Lexical decision time data
#
# Q: To what extent are lexical decisions influenced by
#    (a) properties of the participants, (b) words (items) 
#    used in the experiment, (c) design properties of the 
#    experiment, or (d) more general factors?
## ---------------------------------------------------
data(lexdec)
library(lme4, keep.source = FALSE)
library(arm)

lm.simp <- lm(RT ~ 1, data=lexdec)
display(lm.simp)

lm.1 <- lm(RT ~ meanWeight + Frequency, data=lexdec)
summary(lm.1)


lm.2 <- lm(RT ~ 1 + Length + meanWeight + Frequency + FamilySize + SynsetCount
		   , data=lexdec)
summary(lm.2)

# comparison of models
anova(lm.2, lm.1, lm.simp)

## ---------------------------------------------------
# investigate the factors in the data set and build a good model 
# ONLY USING LINGUISTIC FACTORS
## ---------------------------------------------------

## ---------------------------------------------------
# Do subjects differ?
## ---------------------------------------------------

ll.1 <- lmList(RT ~ 1 + meanWeight + Frequency | Subject, data=lexdec)
str(ll.1)

# a little function to calculate the standard error of a 
# sample (in a vector)
se <- function(x) {
	var(x) / sqrt(length(x))
}

# coef and standard error of coef from each model
paste("[", sapply(ll.1, coef), ", ", sapply(ll.1, se.coef), "]", sep="")

# coef and standard error across model
paste("[", mean(sapply(ll.1, coef)[2,]), ", ", se(sapply(ll.1, se.coef)[,2]), "]", sep="")


# Trellis plot of individual participant coefficients, xyplot
trellis.device(color=F)
xyplot(RT ~ Frequency | Subject,
	data=lexdec, 
	main="By-subject LMs",
	ylab="log reaction times",
	xlab="",
#	ylim= c(-200,200),
	panel=function(x, y){
		panel.xyplot(x, y, col=1)
		panel.lmline(x, y, lty=4, col="blue", lwd=3)
	}	
)

old.prompt = grid::grid.prompt(T)
qqmath( ~ RT | Subject, data= lexdec, layout=c(4,4),
       prepanel = prepanel.qqmathline,
       panel = function(x, ...) {
          panel.qqmathline(x, col= 2, lty=1, ...)
          panel.qqmath(x, ...)
       }
)
grid::grid.prompt(old.prompt)

xylowess.fnc(RT ~ meanWeight | Subject, data=lexdec)


## ---------------------------------------------------
# multilevel (aka mixed) model
## ---------------------------------------------------

lmer.1 <- lmer(RT ~ 1 + meanWeight + Frequency +
		    (1 | Subject)
		   , data=lexdec)
summary(lmer.1)

# the random effects
ranef(lmer.1)
var(ranef(lmer.1)[[1]])


## ---------------------------------------------------
# R2
## ---------------------------------------------------
summary(lm.1)
cor(fitted(lm.1), lexdec$RT)^ 2

cor(fitted(lmer.1), lexdec$RT) ^ 2

## ---------------------------------------------------
# use R2 comparisons to investigate
# item effects
# design effects
# effects of subject differences
## ---------------------------------------------------




lexdec.lmer = lmer(RT ~ 1 + Correct + Trial + PrevType * meanWeight + 
   Frequency + NativeLanguage * Length + (1|Subject) + (1|Word), 
   data = lexdec)







lexdec.lmer = lmer(RT ~ 1 + Correct + Trial + PrevType * meanWeight + 
   Frequency + NativeLanguage * Length + (1|Subject) + (1|Word), 
   data = lexdec)


lexdec.lmer = lmer(RT ~ 1 + Correct + Trial + PrevType * meanWeight + 
   Frequency + NativeLanguage * Length + (1|Subject) + (1|Word), 
   data = lexdec)
pvals.fnc(lexdec.lmer)$summary