Clear work space

rm(list = ls(all = TRUE))

#SetDirectory and have a look at the files

#Go to session->set working directory->to source file location

#Load relevant libraries and functions - these might need to be installed individually

    library(languageR)
    library(lattice)
    library(plotrix)
    library(ggplot2)    #for graphs
    library(psych)
    library(plyr)
    library(doBy)   #to reshape
  library(lme4)
  library(knitr)

#Load Helper Functions

source("loadFunctionsGithub.R")

urlFolder <- 'https://api.github.com/repos/n400peanuts/languagelearninglab/git/trees/master?recursive=1'
urlRaw <- 'https://raw.githubusercontent.com/n400peanuts/languagelearninglab/master/tools/'

listFunctions <- c("summarySEwithin.R", "summarySE.R", "normDataWithin.R", "myCenter.R", "lizCenter.R", "bayesFactorUpdated.R", "Bf.R", "Bf_powercalc.R", "Bf_range.R", "logodds.R")

loadFunctionsGithub(urlFolder, urlRaw, listFunctions)
## [1] "----loading. Please wait----"

#IGL body-rime (2015-2016)

##Exposure Reaction Times analyses

###Load database

#Create Dataframes that we will be working on
#IGL data children

    IGL_expoRT <- read.csv("RT_trimimed_expo.csv", header=TRUE)
    IGL_expoRT$Subject = factor(IGL_expoRT$Subject)
    IGL_expoRT$Condition = factor(IGL_expoRT$Condition)
    IGL_expoRT$Block = factor(IGL_expoRT$Block)

###Reaction Times at exposure data analyses: Onset vs. Rime; by Blocks

#aggregate
    aggregated.IGLexpoRT= aggregate(TrimmedRTs ~ Subject + Condition + Block, IGL_expoRT, FUN=mean)

    ##calculate means and CORRECTED SEs using summarySEwithin (within; btw factors) to add error bars to my graphs
    Means_IGL_expo <- summarySEwithin(aggregated.IGLexpoRT, measurevar="TrimmedRTs", withinvars=c("Block"), betweenvar=c("Condition"), idvar="Subject", na.rm=FALSE, conf.interval=.95)
    Means_IGL_expo
##    Condition Block  N TrimmedRTs TrimmedRTs_norm        sd       se       ci
## 1       body     1 25   989.4043        979.3267  61.14178 12.22836 25.23808
## 2       body     2 25   989.9359        979.8583  56.55866 11.31173 23.34627
## 3       body     3 25  1009.0612        998.9836  63.80373 12.76075 26.33688
## 4       body     4 24   948.6935        908.5187  50.23401 10.25397 21.21196
## 5       body     5 24  1003.1552        962.9804  66.34339 13.54229 28.01435
## 6       body     6 24  1012.0110        971.8362  66.53290 13.58097 28.09438
## 7       rime     1 33   929.7668        948.1931  63.97587 11.13677 22.68486
## 8       rime     2 33   960.1193        978.5456  64.44946 11.21921 22.85279
## 9       rime     3 33   980.0030        998.4293  58.26744 10.14306 20.66074
## 10      rime     4 33   914.7249        933.1512  58.75902 10.22863 20.83504
## 11      rime     5 33   948.1677        966.5940  67.78292 11.79949 24.03478
## 12      rime     6 33   960.5079        978.9342 100.46654 17.48898 35.62389
    p = ggplot(Means_IGL_expo, aes(x = Block, y = TrimmedRTs, group=Condition, colour= Condition))
    p = p + geom_line(stat="identity", size=1)
    p = p + geom_point(stat="identity", size=2.5)
    p = p + geom_errorbar(aes(ymin=TrimmedRTs-ci, ymax=TrimmedRTs+ci), width=.2)
    p = p + xlab("Block") + ylab("Mean Reaction Time in milliseconds")
    p = p + geom_hline(yintercept=1000, linetype="dashed", color="black")
    p = p + theme_bw()
    p = p + coord_cartesian(ylim = c(0, 1500))
    p = p + ggtitle('IGL Reaction Times at Exposure')
    p = p+ theme(panel.background = element_rect(fill = 'white'),
    panel.grid.major = element_line(colour = "white"))
    print(p)

  ggsave(file="IGL_expoRT.jpeg")
  ggsave(file="IGL_expoRT.pdf")

###LMEs for expo RTs

#make sure that accuracy is always a factor
IGL_expoRT$Condition = factor(IGL_expoRT$Condition)
IGL_expoRT$Block = factor(IGL_expoRT$Block)
IGL_expoRT$TrimmedRTs = factor(IGL_expoRT$TrimmedRTs)

#centres the factor
IGL_expoRT <- lizCenter(IGL_expoRT, list("Condition"))
IGL_expoRT <- lizCenter(IGL_expoRT, list("Block"))
IGL_expoRT <- lizCenter(IGL_expoRT, list("TrimmedRTs"))

##does glmer (because accuracy data is binary)
#.ct are the centred factors
modelExpoRT <- glmer(TrimmedRTs ~  Condition.ct *Block.ct + (Condition.ct*Block.ct|Subject), control=glmerControl(optimizer = "bobyqa"), family = binomial, data = IGL_expoRT)
summary(modelExpoRT)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: TrimmedRTs ~ Condition.ct * Block.ct + (Condition.ct * Block.ct |  
##     Subject)
##    Data: IGL_expoRT
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     44.6    152.3     -8.3     16.6    16166 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -38.884   0.000   0.000   0.000   0.026 
## 
## Random effects:
##  Groups  Name                  Variance Std.Dev. Corr             
##  Subject (Intercept)           1.8093   1.3451                    
##          Condition.ct          1.6416   1.2813    0.02            
##          Block.ct              0.4315   0.6569    0.75  0.61      
##          Condition.ct:Block.ct 0.3768   0.6139    0.19 -0.64 -0.45
## Number of obs: 16180, groups:  Subject, 60
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)
## (Intercept)              48.84      83.47   0.585    0.558
## Condition.ct             28.36     153.48   0.185    0.853
## Block.ct                 12.33      51.87   0.238    0.812
## Condition.ct:Block.ct    21.49     125.80   0.171    0.864
## 
## Correlation of Fixed Effects:
##             (Intr) Cndtn. Blck.c
## Conditin.ct -0.048              
## Block.ct     0.471 -0.036       
## Cndtn.ct:B.  0.145  0.494 -0.690
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# provides a nice table overview
kable(summary(modelExpoRT)$coefficients, 
      digits = 3)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 48.836 83.467 0.585 0.558
Condition.ct 28.358 153.483 0.185 0.853
Block.ct 12.334 51.869 0.238 0.812
Condition.ct:Block.ct 21.493 125.803 0.171 0.864

##Test Reaction Times analyses

###Load database

#Create Dataframes that we will be working on
#IGL data children

    IGL_testRT <- read.csv("RT_trimimed_test_accurate.csv", header=TRUE)
    IGL_testRT$Subject = factor(IGL_testRT$Subject)
    IGL_testRT$Condition = factor(IGL_testRT$Condition)
    IGL_testRT$Legality = factor(IGL_testRT$Legality)

###Reaction Times at test data analyses: Onset vs. Rime and by legality

#aggregate
    aggregated.IGLtestRT= aggregate(TrimmedRTs ~ Subject + Condition + Legality, IGL_testRT, FUN=mean)

    ##calculate means and CORRECTED SEs using summarySEwithin (within; btw factors) to add error bars to my graphs
    Means_IGL_test <- summarySEwithin(aggregated.IGLtestRT, measurevar="TrimmedRTs", withinvars=c("Legality"), betweenvar=c("Condition"), idvar="Subject", na.rm=FALSE, conf.interval=.95)
    Means_IGL_test
##   Condition Legality  N TrimmedRTs TrimmedRTs_norm       sd        se       ci
## 1      body  illegal 21   2150.327        2205.097 554.6356 121.03140 252.4671
## 2      body    legal 23   2320.983        2271.305 529.0122 110.30668 228.7620
## 3      rime  illegal 29   2361.145        2353.791 444.8543  82.60736 169.2135
## 4      rime    legal 30   2122.567        2129.424 438.1085  79.98730 163.5924
    p = ggplot(Means_IGL_test, aes(x = Legality, y = TrimmedRTs, group=Condition, colour= Condition))
    p = p + geom_line(stat="identity", size=1)
    p = p + geom_point(stat="identity", size=3.5)
    p = p + geom_errorbar(aes(ymin=TrimmedRTs-ci, ymax=TrimmedRTs+ci), width=0.2)
    p = p + xlab("Legality") + ylab("Mean Reaction Time in milliseconds")
    p = p + geom_hline(yintercept=1000, linetype="dashed", color="black")
    p = p + theme_bw()
    p = p + coord_cartesian(ylim = c(0, 3000))
    p = p + ggtitle('IGL Reaction Times at Test')
    p = p+ theme(panel.background = element_rect(fill = 'white'),
    panel.grid.major = element_line(colour = "white"))
    print(p)

  ggsave(file="IGL_testRT.jpeg")
  ggsave(file="IGL_testRT.pdf")

###LMEs for test RTs

#make sure that accuracy is always a factor
IGL_testRT$Condition = factor(IGL_testRT$Condition)
IGL_testRT$Legality = factor(IGL_testRT$Legality)
IGL_testRT$TrimmedRTs = factor(IGL_testRT$TrimmedRTs)

#centres the factor
IGL_testRT <- lizCenter(IGL_testRT, list("Condition"))
IGL_testRT <- lizCenter(IGL_testRT, list("Legality"))
IGL_testRT <- lizCenter(IGL_testRT, list("TrimmedRTs"))

##does glmer (because accuracy data is binary)
#.ct are the centred factors
modelTestRT <- glmer(TrimmedRTs ~  Condition.ct *Legality.ct + (Condition.ct*Legality.ct|Subject), control=glmerControl(optimizer = "bobyqa"), family = binomial, data = IGL_testRT)
summary(modelTestRT)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## TrimmedRTs ~ Condition.ct * Legality.ct + (Condition.ct * Legality.ct |  
##     Subject)
##    Data: IGL_testRT
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##     38.6     96.2     -5.3     10.6      438 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -8.544  0.000  0.000  0.000  0.117 
## 
## Random effects:
##  Groups  Name                     Variance Std.Dev. Corr             
##  Subject (Intercept)              1.4882   1.2199                    
##          Condition.ct             0.5928   0.7699    0.64            
##          Legality.ct              2.1398   1.4628    0.93  0.40      
##          Condition.ct:Legality.ct 0.8671   0.9312   -0.54  0.02 -0.44
## Number of obs: 452, groups:  Subject, 56
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                21.684    322.871   0.067    0.946
## Condition.ct                9.334    428.122   0.022    0.983
## Legality.ct                 8.839    377.004   0.023    0.981
## Condition.ct:Legality.ct  -21.251    270.375  -0.079    0.937
## 
## Correlation of Fixed Effects:
##             (Intr) Cndtn. Lglty.
## Conditin.ct  0.661              
## Legality.ct  0.609 -0.131       
## Cndtn.ct:L. -0.324  0.102 -0.196
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# provides a nice table overview
kable(summary(modelExpoRT)$coefficients, 
      digits = 3)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 48.836 83.467 0.585 0.558
Condition.ct 28.358 153.483 0.185 0.853
Block.ct 12.334 51.869 0.238 0.812
Condition.ct:Block.ct 21.493 125.803 0.171 0.864

##Accuracy at test analyses

###Load database

  IGL_TEST_accuracy <- read.csv("IGL_body_rime_master_2015_2016_test.csv", header=TRUE)
  IGL_TEST_accuracy$Condition = factor(IGL_TEST_accuracy$Condition)
  IGL_TEST_accuracy$Subject = factor(IGL_TEST_accuracy$Subject)

#aggregate
    aggregated.IGL_TEST_accuracy= aggregate(accuracy ~ Subject + Condition, IGL_TEST_accuracy, FUN=mean)

###Calculate means and CORRECTED SEs using summarySEwithin (within; btw factors) to add error bars to my graphs

    Means_IGL_TEST_accuracy <- summarySE(aggregated.IGL_TEST_accuracy, measurevar="accuracy", groupvars ="Condition", na.rm=FALSE, conf.interval=.95)
    Means_IGL_TEST_accuracy
##   Condition  N  accuracy        sd         se         ci
## 1      body 27 0.5324074 0.1254444 0.02414178 0.04962414
## 2      rime 33 0.5435606 0.1366494 0.02378761 0.04845377
    p = ggplot(Means_IGL_TEST_accuracy, aes(x = Condition, y = accuracy, group=Condition))
    p = p + geom_bar(stat="identity", size=.5)
    p = p + geom_errorbar(aes(ymin=accuracy-ci, ymax=accuracy+ci), width=0.2)
    p = p + xlab("Condition") + ylab("Mean Proportion Correct")
    p = p + geom_hline(yintercept=0.5, linetype="dashed", color="black")
    p = p + theme_bw()
    p = p + coord_cartesian(ylim = c(0, 1))
    p = p + ggtitle('IGL accuracy at Test')
    p = p+ theme(panel.background = element_rect(fill = 'white'),
    panel.grid.major = element_line(colour = "white"))
    print(p)

  ggsave(file="IGL_test_accuracy.jpeg")
  ggsave(file="IGL_test_accuracy.pdf")

###T-Tests

#subsetting to get Body

  IGL_TEST_accuracy_body = subset(IGL_TEST_accuracy, Condition !="rime") 
  accuracy.body <- IGL_TEST_accuracy_body$accuracy

#subsetting to get Rime

  IGL_TEST_accuracy_rime = subset(IGL_TEST_accuracy, Condition !="body") 
  accuracy.rime <- IGL_TEST_accuracy_rime$accuracy
  
  accuracy.both = IGL_TEST_accuracy$accuracy
  
#t-test
  t.test(accuracy.body, mu=0.5, alternative = "greater")
## 
##  One Sample t-test
## 
## data:  accuracy.body
## t = 1.3484, df = 431, p-value = 0.08911
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.4927906       Inf
## sample estimates:
## mean of x 
## 0.5324074
  t.test(accuracy.rime, mu=0.5, alternative = "greater")
## 
##  One Sample t-test
## 
## data:  accuracy.rime
## t = 2.0076, df = 527, p-value = 0.0226
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.5078085       Inf
## sample estimates:
## mean of x 
## 0.5435606
  t.test(accuracy.both, mu=0.5, alternative = "greater")
## 
##  One Sample t-test
## 
## data:  accuracy.both
## t = 2.3942, df = 959, p-value = 0.008423
## alternative hypothesis: true mean is greater than 0.5
## 95 percent confidence interval:
##  0.5120375       Inf
## sample estimates:
## mean of x 
## 0.5385417

###LMEs

#make sure that accuracy is always a factor
IGL_TEST_accuracy$Condition = factor(IGL_TEST_accuracy$Condition)

#centres the factor
IGL_TEST_accuracy <- lizCenter(IGL_TEST_accuracy, list("Condition"))

##does glmer (because accuracy data is binary)
#.ct are the centred factors
modelTestAccuracy <- glmer(accuracy ~  Condition.ct* + (Condition.ct|Subject), control=glmerControl(optimizer = "bobyqa"), family = binomial, data = IGL_TEST_accuracy)
summary(modelTestAccuracy)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: accuracy ~ Condition.ct * +(Condition.ct | Subject)
##    Data: IGL_TEST_accuracy
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1334.9   1359.2   -662.5   1324.9      955 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1340 -1.0804  0.8980  0.9247  0.9750 
## 
## Random effects:
##  Groups  Name         Variance Std.Dev. Corr
##  Subject (Intercept)  0.00000  0.0000       
##          Condition.ct 0.06282  0.2506    NaN
## Number of obs: 960, groups:  Subject, 60
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.15508    0.06686   2.319   0.0204 *
## Condition.ct  0.04482    0.13449   0.333   0.7390  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Conditin.ct -0.011
## convergence code: 0
## boundary (singular) fit: see ?isSingular
# provides a nice table overview
kable(summary(modelTestAccuracy)$coefficients, 
      digits = 3)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.155 0.067 2.319 0.020
Condition.ct 0.045 0.134 0.333 0.739
#Looks at the mean for each group and condition
round(with(IGL_TEST_accuracy, tapply(accuracy, list(Condition), mean, na.rm=T)),3)
##  body  rime 
## 0.532 0.544

##D-Prime analysis

###Load database

  IGL_TEST_dprime <- read.csv("D_prime_all.csv", header=TRUE)
  IGL_TEST_dprime$Condition = factor(IGL_TEST_dprime$Condition)
  IGL_TEST_dprime$Subject = factor(IGL_TEST_dprime$Subject)

##calculate means and CORRECTED SEs using summarySEwithin (within; btw factors) to add error bars to my graphs

    Means_IGL_TEST_dprime <- summarySE(IGL_TEST_dprime, measurevar="DPRIME", groupvars ="Condition", na.rm=FALSE, conf.interval=.95)
    Means_IGL_TEST_dprime
##   Condition  N    DPRIME        sd        se        ci
## 1      body 27 0.1883091 0.7917437 0.1523712 0.3132034
## 2      rime 33 0.2775233 0.8164994 0.1421343 0.2895181
    p = ggplot(Means_IGL_TEST_dprime, aes(x = Condition, y = DPRIME, group=Condition))
    p = p + geom_bar(stat="identity", size=.5)
    p = p + geom_errorbar(aes(ymin=DPRIME-ci, ymax=DPRIME+ci), width=0.22)
    p = p + xlab("Condition") + ylab("D-prime")
    p = p + geom_hline(yintercept=0, linetype="dashed", color="black")
    p = p + theme_bw()
    p = p + theme(text = element_text(size=22 ,color="black"))
    p = p + ggtitle('IGL Signal Detection/Sensitivity at Test')
    p = p + theme(panel.background = element_rect(fill = 'white'))
    panel.grid.major = element_line(colour = "white")
    print(p)

  ggsave(file="IGL_test_dprime.jpeg")
  ggsave(file="IGL_test_dprime.pdf")

###T-Tests on D-prime

#subsetting to get Body

  IGL_TEST_dprime_body = subset(IGL_TEST_dprime, Condition !="rime") 
  dprime.body <- IGL_TEST_dprime_body$DPRIME

#subsetting to get Rime

  IGL_TEST_dprime_rime = subset(IGL_TEST_dprime, Condition !="body") 
  dprime.rime <- IGL_TEST_dprime_rime$DPRIME
  
  dprime.both = IGL_TEST_dprime$DPRIME
  
#t-test
  t.test(dprime.body, mu=0, alternative = "greater")
## 
##  One Sample t-test
## 
## data:  dprime.body
## t = 1.2359, df = 26, p-value = 0.1138
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.07157785         Inf
## sample estimates:
## mean of x 
## 0.1883091
  t.test(dprime.rime, mu=0, alternative = "greater")
## 
##  One Sample t-test
## 
## data:  dprime.rime
## t = 1.9525, df = 32, p-value = 0.02983
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.03676358        Inf
## sample estimates:
## mean of x 
## 0.2775233
  t.test(dprime.both, mu=0, alternative = "greater")
## 
##  One Sample t-test
## 
## data:  dprime.both
## t = 2.2987, df = 59, p-value = 0.01254
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.06481004        Inf
## sample estimates:
## mean of x 
## 0.2373769