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