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


## @knitr LinearityExampleVisualization, fig.show='hold'
library(languageR)
data(lexdec)

g = glm(RT ~ I(Frequency^2) + I(Frequency^3), data=lexdec)
summary(g)$coefficients

library(ggplot2)
ggplot(lexdec, aes(x=I(Frequency^2), y=RT)) +
  geom_point() +
  geom_abline(
    intercept = coef(g)['(Intercept)'], 
    slope = coef(g)['I(Frequency^2)'],
    color = "blue")

ggplot(lexdec, aes(x=I(Frequency^3), y=RT)) +
  geom_point() +
  geom_abline(
    intercept = coef(g)['(Intercept)'], 
    slope = coef(g)['I(Frequency^3)'],
    color = "red")


## @knitr LinearityExampleVisualization3, echo=FALSE, out.width='8cm', out.height='8cm'
p = function(x) {
  y = coef(g)['(Intercept)'] + 
    coef(g)['I(Frequency^2)'] * x^2 +
    coef(g)['I(Frequency^3)'] * x^3 
  return(y)
}
ggplot(lexdec, aes(x=Frequency, y=RT)) +  
  geom_point() +
  geom_smooth(method='lm', fill=NA, color="black", linetype="dashed") +
  stat_function(fun= p, geom="line", color = "green", size=2) +
  stat_function(fun= function(x) {coef(g)['(Intercept)'] + 
    coef(g)['I(Frequency^2)'] * x^2}, geom="line", color = "blue") +
  stat_function(fun= function(x) {coef(g)['(Intercept)'] + 
    coef(g)['I(Frequency^3)'] * x^3}, geom="line", color = "red") +
  annotate("text", label = "x^2", x=7, y=5.9, color = "blue", size=10) +
  annotate("text", label = "x^3", x=7, y=7, color = "red", size=10) 


## @knitr FrequencyTrialModel, echo=T, results='hide', size='scriptsize'
summary(glm(RT ~ Frequency + Trial, data= lexdec))


## @knitr Centering, echo=T, size='scriptsize'
lexdec$cFrequency = lexdec$Frequency - mean(lexdec$Frequency)
summary(lexdec[,c('Frequency','cFrequency')])


## @knitr FrequencyTrialModel2, echo=T, size='scriptsize'
summary(glm(RT ~ Frequency * Trial, data= lexdec))$coefficients


## @knitr Predict, echo=T, eval=F, results='hide', size='scriptsize'
# Beyond the coefficients, another interesting thing to compare is 
# the models predictions (the predicted RTs)
g1 =  glm(RT ~ Frequency * Trial, data= lexdec)
g2 =  glm(RT ~ cFrequency * cTrial, data= lexdec)

# if the models make the same prediction the difference in their 
# predictions should be zero
predict(g1) - predict(g2)


## @knitr NativeSexModel, echo=T, results='hide', size='scriptsize'
summary(glm(RT ~ NativeLanguage * Sex, data= lexdec))$coefficients


## @knitr ReorderingFactorLevels, echo=T, size='scriptsize'
# by default R assumes alphanumeric ordering of factor levels
# but we can always change that by explicitly stating what
# order we would like:
lexdec$NativeLanguage = 
  factor(lexdec$NativeLanguage, levels=c('Other','English'))

# Btw, c() is simply referring to a vector in R. They can 
# contain numbers, strings, or any (combination) of other 
# objects, e.g. c(4, "seven", "Hans the Wurst", c(log(2), "hai"))


## @knitr AnovaCoding, echo=T, size='scriptsize'
# confirm that the default contrast was treatment-coding, e.g.
contrasts(lexdec$Sex)

# change contrasts
contrasts(lexdec$Sex) = contr.sum(2)
contrasts(lexdec$Sex)
# notice that default sum-coding in R assigns 1 to the 
# alphanumerically first level and -1 to the second level

# let's label the contrast so that we don't get confused 
# (cbind() makes a matrix out of multiple vectors; in this 
# case 1 vector, but it allows us to label the columns/vectors)
contrasts(lexdec$Sex) = cbind("Male" = c(-1, 1))


## @knitr AnovaCoding2, echo=T, size='scriptsize'
contrasts(lexdec$Sex) = contr.sum(2) / 2
contrasts(lexdec$Sex)

contrasts(lexdec$NativeLanguage) = contr.sum(2) / 2
contrasts(lexdec$NativeLanguage)