

```{r config, echo = FALSE, warning = FALSE, results = 'hide', message=FALSE}
library(plyr)
library(lmSupport)
library(lme4)
library(reshape2)
library(ggplot2)
library(cowplot)
library(gridExtra)

# set working directory to your current folder that this Rmd file and all the logs are saved in. 
setwd("~/Downloads/osfstorage-archive")

# this analysis has the following sections:
# data preprocessing
# FC AC V: Forced Choice Vocabulary Accuracy incl covariate
# FC RT V: Forced Choice Vocabulary Reaction Time
# FC AC A: Forced Choice Suffix Understanding Accuracy
# FC RT A: Forced Choice Suffix Understanding Reaction Time
# Error monitoring task d prime
# EM AC W: Error Monitoring Word Order Accuracy
# EM RT W: Error Monitoring Word Order Reaction Time
# EM AC A: Error Monitoring Agreement Accuracy
# EM RT A: Error Monitoring Agreement Reaction Time
# results figures
# probabilistic co-occurrence (SOM-U)
```
## data processing

```{r, echo = FALSE, warning = FALSE, results = 'hide'}

# alternatively to running through all of the R code below, I have saved my workspace, and you can just load it in. 
load(file = "WorkspaceHopmanMacDonaldPsychScience.Rdata")

rawDataFC = read.csv('FClogsusable.csv', header=TRUE)
rawDataEM = read.csv('EMlogsusable.csv', header=TRUE)

d <-rbind.fill(rawDataFC, rawDataEM)

# recode learning condition to be centered around zero so that Comprehension = -0.5 and Production = 0.5
d$learningCondition <- varRecode(d$conditionnr,c(1,2), c(-0.5,0.5))

# calculate the reaction time as measured from the onset of the critical word
d$ART = 0
d$ART[d$critword == 1] = d$RT[d$critword == 1] 
d$ART[d$critword == 2] = d$RT[d$critword == 2]  - d$cumtimew1[d$critword == 2] 
d$ART[d$critword == 3] = d$RT[d$critword == 3]  - d$cumtimew2[d$critword == 3] 
d$ART[d$critword == 4] = d$RT[d$critword == 4]  - d$cumtimew3[d$critword == 4] 
d$ART[d$critword == 5] = d$RT[d$critword == 5]  - d$cumtimew4[d$critword == 5] 
d$ART[d$critword == 6] = d$RT[d$critword == 6]  - d$cumtimew5[d$critword == 6] 
d$ART[d$critword == 7] = d$RT[d$critword == 7]  - d$cumtimew6[d$critword == 7] 
d <- rename(d, c("ART"="reactionTime"))

relevantColumnsFC = c("subjectnr", "learningCondition", "itemtype", "correctanswer", "reactionTime")
relevantColumnsEM = c("subjectnr", "learningCondition", "itemtype", "correctanswer", "reactionTime", "wassoundcorrect", "trialnr")

# write the relevant columns away as the full dataset for EM and FC trials for accuracy analysis
allDataFC <- subset(d, d$trialtype == 'FC', relevantColumnsFC)
allDataEM <- subset(d, d$trialtype == 'EM', relevantColumnsEM)

# for the RT analysis we only want trials that were answered correctly and trials with a positive RT (negative RTs happen when participants press before the critical word in a trial)
dforRTFC = subset(allDataFC, allDataFC$reactionTime > 0 & allDataFC$correctanswer == 1)
dforRTEM = subset(allDataEM, allDataEM$reactionTime > 0 & allDataEM$correctanswer == 1)

# for each of the two tasks, get a participants mean RT
mRTFC <- ddply(dforRTFC,~subjectnr,summarise,meanRT=mean(reactionTime),sdRT=sd(reactionTime))
mRTEM <- ddply(dforRTEM, ~subjectnr,summarise,meanRT=mean(reactionTime), sdRT=sd(reactionTime))

dforRTFC = merge(dforRTFC, mRTFC)
dforRTEM = merge(dforRTEM, mRTEM)

# remove trials where a participant is slower than their own M + 3*SD 
RTDataFC = subset(dforRTFC, dforRTFC$reactionTime < dforRTFC$meanRT + 3*dforRTFC$sdRT, relevantColumnsFC)
RTDataEM = subset(dforRTEM, dforRTEM$reactionTime < dforRTEM$meanRT + 3*dforRTEM$sdRT, relevantColumnsEM)

# calculate the % of original FC+EM trials that are retained for RT analysis, reported under data processing
percentageRTData = (nrow(RTDataFC)+ nrow(RTDataEM))/(nrow(allDataFC)+nrow(allDataEM))
```

```{r, echo = TRUE}
percentageRTData # reported in Data Processing
```
## FC AC V 
```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# the first data to analyze is the covariate, the accuracy in the vocabulary trials in the FC task. 
dataFCACV = subset(allDataFC, allDataFC$itemtype == "FCvocabcolr" | allDataFC$itemtype == "FCvocabmons" | allDataFC$itemtype == "FCvocabpatt" | allDataFC$itemtype == "FCvocabvvid" | allDataFC$itemtype == "FCvocablvid")

# see how accuracy on these vocabulary trials differs by condition
mFCACV <- glmer(correctanswer ~ learningCondition + (1 | subjectnr), data = dataFCACV, family = "binomial")
```

```{r, echo = TRUE}
summary(mFCACV) # table 3 in SOM-R
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pLearningCondition <- data.frame(learningCondition = c(-0.5,0.5))
pFCACV <- modelPredictions(mFCACV, pLearningCondition)
pFCACV$Itemtype = "a) FC Vocabulary"

# get each participants' accuracy score on this test as a covariate for all other analyses. 
Covariate <- ddply(dataFCACV,~subjectnr + learningCondition,summarise,vocabCovariate=mean(correctanswer))
meanCovariate = mean(Covariate$vocabCovariate)
Covariate$vocabCovariateC <- Covariate$vocabCovariate - meanCovariate

# and add the covariate to all dataframes (both FC and EM tasks and both the data used for accuracy and RT analyses)
allDataFC = merge(allDataFC, Covariate)
allDataEM = merge(allDataEM, Covariate)
RTDataFC = merge(RTDataFC, Covariate)
RTDataEM = merge(RTDataEM, Covariate)

## FC RT V

# next, let's look at RT on these same vocabulary trials from the FC test
dataFCRTV = subset(RTDataFC, RTDataFC$itemtype == "FCvocabcolr" | RTDataFC$itemtype == "FCvocabmons" | RTDataFC$itemtype == "FCvocabpatt" | RTDataFC$itemtype == "FCvocabvvid" | RTDataFC$itemtype == "FCvocablvid")

# see how RT on these vocabulary trials differs by condition
mFCRTV = lmer(reactionTime ~  vocabCovariateC + learningCondition+ (1|subjectnr), data = dataFCRTV)

```

```{r, echo = TRUE}
modelSummary(mFCRTV) # table 4 in SOM-R
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pLearningConditionCovariate <- data.frame(learningCondition = c(-0.5, 0.5), vocabCovariateC = 0)
pFCRTV <- modelPredictions(mFCRTV, pLearningConditionCovariate)
pFCRTV$Itemtype = "b) FC Vocabulary"

## FC AC A 

# next, let's look at accuracy on the suffix understanding trials in the FC task
dataFCACA = subset(allDataFC, allDataFC$itemtype == "FCnumbritem" | allDataFC$itemtype == "FCsemanitem")

# recode itemtype to code for whether an item is testing understanding of the semantic or number aspect of the suffixes
dataFCACA$snItemtype = .5 # semantic items coded as .5
dataFCACA$snItemtype[dataFCACA$itemtype == 'FCnumbritem'] = -.5 # number items coded as -.5

# see how accuracy on these suffix understanding trials differs by learning condition and semantic/number itemtype
mFCACA = glmer(correctanswer ~ vocabCovariateC + learningCondition*snItemtype  + (1 + snItemtype|subjectnr), data=dataFCACA, family = "binomial")

```

```{r, echo = TRUE}
summary(mFCACA) # table 5 in SOM-R

```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}

# for making figures later on, save model predictions
pLearningConditionCovariateSNItemtype <- data.frame(learningCondition = c(-0.5, 0.5), snItemtype = c(0,0), vocabCovariateC = 0) # note that for the picture we are only interested in learningCondition, so we take the average of all other model predictors
pFCACA <- modelPredictions(mFCACA, pLearningConditionCovariateSNItemtype)
pFCACA$Itemtype = "c) FC Suffix"

## FC RT A

# next, let's look at RT on these same suffix understanding trials from the FC test
dataFCRTA = subset(RTDataFC, RTDataFC$itemtype == "FCnumbritem" | RTDataFC$itemtype == "FCsemanitem")

# recode itemtype to code for whether an item is testing understanding of the semantic or number aspect of the suffixes
dataFCRTA$snItemtype = .5 # semantic items coded as -.5
dataFCRTA$snItemtype[dataFCRTA$itemtype == 'FCnumbritem'] = -.5 # number items coded as .5

# see how RT on these suffix understanding trials differs by condition
mFCRTA = lmer(reactionTime ~ vocabCovariateC + learningCondition*snItemtype + (1 + snItemtype|subjectnr), data = dataFCRTA)

```

```{r, echo = TRUE}
modelSummary(mFCRTA) # table 6 in SOM-R
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}

# for making figures later on, save model predictions
pFCRTA <- modelPredictions(mFCRTA, pLearningConditionCovariateSNItemtype)
pFCRTA$Itemtype = "d) FC Suffix"

## d prime analysis on the EM data as a whole

# code for hits and false alarms
allDataEM[,"Hits"] <- 0
allDataEM[,"FA"] <- 0

# set Hits to 1 when the sentence was grammatically correct and the participant answered it correctly saying it was grammatical
allDataEM$Hits[allDataEM$wassoundcorrect==1 & allDataEM$correct == 1] <- 1  
# set FA to 1 when the sentence was grammatically incorrect and the participant answered it incorrectly saying it was grammatical
allDataEM$FA[allDataEM$wassoundcorrect== 0 & allDataEM$correct == 0] <- 1 

# in order to create a wide dataframe for Hits
vars=c("subjectnr","learningCondition","Hits", "trialnr")
dataH=allDataEM[vars] 
# cast into a wide dataframe
dH<-dcast(dataH,subjectnr + learningCondition~ trialnr,value.var="Hits") 
# sum the hits
dH$Hits=rowSums(dH[,3:length(dH)]) 
vars=c("subjectnr","Hits", "learningCondition") 
dH=dH[vars]

# do the same for FAs
vars=c("subjectnr","learningCondition","FA", "trialnr")
dataFA=allDataEM[vars] 
dFA<-dcast(dataFA,subjectnr~trialnr,value.var="FA") 
dFA$FA=rowSums(dFA[,2:length(dFA)]) # sum false alarms
vars=c("subjectnr","FA") # select variables
dFA=dFA[vars] # create dataframe with just those vars

# merge and look at wide dataframe
dW = merge(dH,dFA,by="subjectnr")

# make hit rate into proportion
Ncorrect = length(allDataEM$wassoundcorrect[allDataEM$wassoundcorrect==1 & allDataEM$subjectnr==101]) # nr of trials where sound was correct, should be 44
dW$Hprop = dW$Hits/Ncorrect
# make FA rate into proportion
Nincorrect = length(allDataEM$wassoundcorrect[allDataEM$wassoundcorrect==0 & allDataEM$subjectnr==101]) # nr of trials where sound was incorrect, should be 80
dW$FAprop = dW$FA/Nincorrect

# check for values of 0 and 1
varDescribe(dW) 
# ok, so someone has a hit rate of 1, that needs to be replaced by 1-1/2N with N the nr of stimuli (in this case Ncorrect, 44)
correctionh = 1-1/(2*Ncorrect)
dW$Hprop[dW$Hprop ==1] <- correctionh

# calculate d' for each participant
dW$dprime <- qnorm(dW$Hprop)-qnorm(dW$FAprop)

# calculate response bias for each participant
dW$c <- -0.5*(qnorm(dW$Hprop)+qnorm(dW$FAprop))

# see whether dprime differs by condition, report this in the Results section for Error Monitoring Tests
mD1 <- lm(dprime ~ learningCondition, data=dW)

```

```{r, echo = TRUE}
modelSummary(mD1) # d prime analysis reported in the first paragraph of Error Monitoring Tests
varDescribeBy(dW$dprime, dW$learningCondition) # d prime analysis reported in the first paragraph of Error Monitoring Tests
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}

## EM AC W

# next, let's look at accuracy on the word order error trials in the EM task
dataEMACW = subset(allDataEM, allDataEM$itemtype == "EMsw1" | allDataEM$itemtype == "EMsw2" | allDataEM$itemtype == "EMsw3" | allDataEM$itemtype == "EMsw4")

# see how accuracy on these word order errors differs by learning condition
mEMACW = glmer(correctanswer ~ vocabCovariateC + learningCondition + (1|subjectnr), data=dataEMACW, family = "binomial")
```

```{r, echo = TRUE}
summary(mEMACW) # table 7 in SOM-R
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pEMACW <- modelPredictions(mEMACW, pLearningConditionCovariate)
pEMACW$Itemtype = "e) EM Word Order"

## EM RT W

# next, let's look at RT on these same word order error trials
dataEMRTW = subset(RTDataEM, RTDataEM$itemtype == "EMsw1" | RTDataEM$itemtype == "EMsw2" | RTDataEM$itemtype == "EMsw3" | RTDataEM$itemtype == "EMsw4")

# see how RT on these word order error trials differs by condition
mEMRTW = lmer(reactionTime ~ vocabCovariateC + learningCondition + (1|subjectnr), data = dataEMRTW)
```

```{r, echo = TRUE}
modelSummary(mEMRTW) # table 8 in SOM-R
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pEMRTW <- modelPredictions(mEMRTW, pLearningConditionCovariate)
pEMRTW$Itemtype = "f) EM Word Order"

## EM AC A

# next, let's look at accuracy on the suffix agreement error trials in the EM task
dataEMACA = subset(allDataEM, allDataEM$itemtype == "EMnadj" | allDataEM$itemtype == "EMnnadj" | allDataEM$itemtype == "EMsadj" | allDataEM$itemtype == "EMsnadj")

# recode itemtype to code for whether an item is testing understanding of the semantic or number aspect of the suffixes
dataEMACA$snItemtype = .5 # semantic items coded as .5
dataEMACA$snItemtype[dataEMACA$itemtype == 'EMnadj' | dataEMACA$itemtype == 'EMnnadj'] = -.5 # number items coded as -.5

# recode itemtype to code for whether an item contains an adjacent (error-suffix on the monster word) or a non-adjacent (error-suffix on the verb) error
dataEMACA$adjacencyItemtype = .5 # non-adjacent error items coded as .5
dataEMACA$adjacencyItemtype[dataEMACA$itemtype == 'EMnadj' | dataEMACA$itemtype == 'EMsadj'] = -.5 # adjacent error items coded as -.5

# see how accuracy on these word order errors differs by learning condition
mEMACA1 = glmer(correctanswer ~  vocabCovariateC + learningCondition*adjacencyItemtype*snItemtype  + (1 + adjacencyItemtype*snItemtype|subjectnr), data=dataEMACA, family = "binomial")
summary(mEMACA1) ## does not converge

# following the recommendations in Barr et al. 2013 I took the following steps:
# 1. leave out the random slopes for adjacent and sntype, but keep the random slope for their interaction (still does not converge)
# 2. leave out all covariates between random slopes (still does not converge)
# 3. remove the random intercept (model converges)
mEMACA2 = glmer(correctanswer ~  vocabCovariateC + learningCondition*adjacencyItemtype*snItemtype  + (0 + adjacencyItemtype:snItemtype|subjectnr), data=dataEMACA, family = "binomial")
```

```{r, echo = TRUE}
summary(mEMACA2) # table 9 in SOM-R
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pLearningConditionCovariateSNItemtypeAdjItemtype <- data.frame(learningCondition = c(-0.5, 0.5), snItemtype = c(0,0), adjacencyItemtype = c(0,0), vocabCovariateC = 0)
pEMACA <- modelPredictions(mEMACA2, pLearningConditionCovariateSNItemtypeAdjItemtype)
pEMACA$Itemtype = "g) EM Suffix"

## EM RT A

# next, let's look at reaction time on the suffix agreement error trials in the EM task
dataEMRTA = subset(RTDataEM, RTDataEM$itemtype == "EMnadj" | RTDataEM$itemtype == "EMnnadj" | RTDataEM$itemtype == "EMsadj" | RTDataEM$itemtype == "EMsnadj")

# recode itemtype to code for whether an item is testing understanding of the semantic or number aspect of the suffixes
dataEMRTA$snItemtype = .5 # semantic items coded as .5
dataEMRTA$snItemtype[dataEMRTA$itemtype == 'EMnadj' | dataEMRTA$itemtype == 'EMnnadj'] = -.5 # number items coded as -.5

# recode itemtype to code for whether an item contains an adjacent (error-suffix on the monster word) or a non-adjacent (error-suffix on the verb) error
dataEMRTA$adjacencyItemtype = .5 # non-adjacent error items coded as .5
dataEMRTA$adjacencyItemtype[dataEMRTA$itemtype == 'EMnadj' | dataEMRTA$itemtype == 'EMsadj'] = -.5 # adjacent error items coded as -.5

# see how accuracy on these word order errors differs by learning condition
mEMRTA = lmer(reactionTime ~ vocabCovariateC + learningCondition*adjacencyItemtype*snItemtype + (1 + adjacencyItemtype*snItemtype|subjectnr), data=dataEMRTA)
```

```{r, echo = TRUE}
modelSummary(mEMRTA) # table 10 in SOM-R
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pEMRTA <- modelPredictions(mEMRTA, pLearningConditionCovariateSNItemtypeAdjItemtype)
pEMRTA$Itemtype = "h) EM Suffix"

## results figures

# put all model predictions for accuracy data in the same dataframe
dataAccuracy <- rbind.fill(pFCACV, pFCACA, pEMACW, pEMACA)
dataAccuracy$conditionStr = varRecode(dataAccuracy$learningCondition,c(-0.5,0.5),c('Comprehension', 'Production'))

# get dataframes with covariate data to plot for each participant in both conditions
Covariate$conditionStr = varRecode(Covariate$learningCondition,c(-0.5,0.5),c('Comprehension', 'Production'))
CovariateComprehension <- Covariate[ which(Covariate$learningCondition==-0.5), ]
CovariateProduction <- Covariate[ which(Covariate$learningCondition==0.5), ]

# first figure is for accuracy
Accuracy <- ggplot(dataAccuracy, aes(x = Itemtype, y = Predicted, fill = conditionStr), width = 0.5) + 
  geom_bar(mapping = aes(fill = as.factor(conditionStr),x=Itemtype,y=Predicted), 
           stat = "identity",
           position=position_dodge(.9), color = "black") + 
  geom_point(data = CovariateComprehension, aes(y = vocabCovariate, x = 0.78),colour='grey14',
             position = position_jitter(w = 0.15, h = 0), size = 2)+
  geom_point(data = CovariateProduction, aes(y = vocabCovariate, x = 1.22),colour='grey14', 
             position = position_jitter(w = 0.15, h = 0), size = 2)+
  geom_errorbar(width=.25, 
                aes(fill = as.factor(conditionStr),y = Predicted, x = Itemtype, 
                    ymin = CILo, 
                    ymax = CIHi), stat="identity", width=0.75, position = position_dodge(.9)) + 
  labs(y = 'Proportion Correct') + 
  theme_bw(base_size = 15) +     
  theme(legend.position = c(0.62,0.86)) +
  guides(fill=guide_legend(title="Learning Condition"))+theme(legend.title = element_text(size= 11, face = "bold"), legend.text = element_text(size = 11, face = "bold"))+
  coord_cartesian(ylim = c(0.47,1.02))+  
  scale_fill_manual(values=c("darkgreen", "lightgreen"))+
  theme( axis.title.x=element_blank(), axis.ticks.x=element_blank())+  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  annotate("text", x=4, y=0.8505, label="*", size=10) +
  geom_label(aes(x = 1, y = 0.465, label = "Covariate"), fill = "white", size = 5, face = "bold")+
  annotate("text", x=2, y=1.005, label="*", size=10) + 
  geom_segment(aes(x=1.78, y=0.97, xend=1.78, yend=1)) +
  geom_segment(aes(x=1.78, y=1, xend=2.22, yend=1)) +
  geom_segment(aes(x=2.22, y=1, xend=2.22, yend=0.99)) +
  geom_segment(aes(x=3.78, y=0.66, xend=3.78, yend=0.85)) +
  geom_segment(aes(x=3.78, y=0.85, xend=4.22, yend=0.85)) +
  geom_segment(aes(x=4.22, y=0.85, xend=4.22, yend=0.81)) +
  theme(axis.text = element_text(colour = "black", face = "bold")) +
  theme(axis.title.y = element_text(colour = "black", face = "bold")) +
  theme(aspect.ratio=2.5/5)
```

```{r, echo = TRUE}
Accuracy # part of Figure 5 with results
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
ggsave("Accuracy.png") # my current window is 8.71 x 6.86 in RStudio

# put all model predictions for RT data in the same dataframe
dataRT <- rbind.fill(pFCRTV, pFCRTA, pEMRTW, pEMRTA)
dataRT$conditionStr = varRecode(dataRT$learningCondition,c(-0.5,0.5),c('Comprehension', 'Production'))

# make RT figure
RT <- ggplot(dataRT, aes(x = Itemtype, y = Predicted, fill = conditionStr), width = 0.5) +
  geom_bar(mapping = aes(fill = as.factor(conditionStr),x=Itemtype,y=Predicted),
           stat = "identity",
           position=position_dodge(.9), color = "black") + 
  geom_errorbar(width=.25, 
                aes(fill = as.factor(conditionStr),y = Predicted, x = Itemtype, 
                    ymin = CILo, 
                    ymax = CIHi), stat="identity", width=0.75, position = position_dodge(.9)) + 
  labs(y = 'RT (S)') + 
  theme_bw(base_size = 15) + 
  coord_cartesian(ylim = c(1.35,2.82))+  
  scale_fill_manual(values=c("darkgreen", "lightgreen"))+
  theme( axis.title.x=element_blank(), axis.ticks.x=element_blank())+  
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  geom_segment(aes(x=0.78, y=1.64, xend=0.78, yend=1.8)) +
  geom_segment(aes(x=0.78, y=1.8, xend=1.22, yend=1.8)) +
  geom_segment(aes(x=1.22, y=1.8, xend=1.22, yend=1.46)) +
  geom_segment(aes(x=1.78, y=2.59, xend=1.78, yend=2.7)) +
  geom_segment(aes(x=1.78, y=2.7, xend=2.22, yend=2.7)) +
  geom_segment(aes(x=2.22, y=2.7, xend=2.22, yend=1.925)) +
  geom_segment(aes(x=2.78, y=2.71, xend=2.78, yend=2.76)) +
  geom_segment(aes(x=2.78, y=2.76, xend=3.22, yend=2.76)) +
  geom_segment(aes(x=3.22, y=2.76, xend=3.22, yend=2.32)) +
  geom_segment(aes(x=3.78, y=2.21, xend=3.78, yend=2.4)) +
  geom_segment(aes(x=3.78, y=2.4, xend=4.22, yend=2.4)) +
  geom_segment(aes(x=4.22, y=2.4, xend=4.22, yend=1.87)) +
  annotate("text", x=1, y=1.81, label="*", size=10) +
  annotate("text", x=2, y=2.71, label="*", size=10) +
  annotate("text", x=3, y=2.77, label="*", size=10)  +
  annotate("text", x=4, y=2.41, label="*", size=10) +
  theme(axis.text = element_text(colour = "black", face = "bold")) +
  theme(axis.title.y = element_text(colour = "black", face = "bold")) +
  theme(aspect.ratio=2/5) +
  theme(legend.position = c(0.135,0.82)) +
  guides(fill=guide_legend(title="Learning Condition"))+theme(legend.title = element_text(size= 11, face = "bold"), legend.text = element_text(size = 11, face = "bold"))
```

```{r, echo = TRUE}
RT # part of Figure 5 with results 
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
ggsave("RT.png") # my current window is 8.71 x 6.86 in RStudio

# save an image of the workspace so that, alternative to running through all of these commands, you can just load in the workspace
save.image(file = "WorkspaceHopmanMacDonaldPsychScience.Rdata")

## probabilistic co-occurrence (SOM-U)

# let's first look at accuracy on the probabilistic co-occurence trials in the FC task
dataFCACP = subset(allDataFC, allDataFC$itemtype == "FCprobitems" | allDataFC$itemtype == "FCimprbitem")

# recode itemtype to code for whether an item is testing understanding of probable or improabable items
dataFCACP$probItemtype = .5 # probable items coded as .5
dataFCACP$probItemtype[dataFCACP$itemtype == 'FCimprbitem'] = -.5 # improbable items coded as -.5

# see how accuracy on these probabilitic co-occurrence trials differs by learning condition and probable/improbable itemtype
mFCACP1 = glmer(correctanswer ~ vocabCovariateC + learningCondition*probItemtype  + (1 + probItemtype|subjectnr), data=dataFCACP, family = "binomial")
summary(mFCACP1) # model doesn't initially converge, which is funny, cause it used to converge in my older version of R
# to solve convergence, I tried a solution found here (accessed on April 19th, 2017)
# https://rstudio-pubs-static.s3.amazonaws.com/33653_57fc7b8e5d484c909b615d8633c01d51.html
# i restarted from the previous fit - since it used to converge just fine in a previous version of R i figured this would work
ss <- getME(mFCACP1,c("theta","fixef"))
mFCACP2 <- update(mFCACP1,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
summary(mFCACP2)
ss2 <- getME(mFCACP2,c("theta","fixef"))
mFCACP3 <- update(mFCACP2,start=ss2,control=glmerControl(optCtrl=list(maxfun=2e4)))
```

```{r, echo = TRUE}
summary(mFCACP3) # table 1 in SOM-U
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pLearningConditionCovariateProbItemtype <- data.frame(learningCondition = c(-0.5,-0.5, 0.5,0.5), probItemtype = c(-0.5, 0.5, -0.5, 0.5), vocabCovariateC = 0) 
pFCACP <- modelPredictions(mFCACP3, pLearningConditionCovariateProbItemtype)

# next, let's look at RT on these same probabilistic co-occurrence trials from the FC test
dataFCRTP = subset(RTDataFC, RTDataFC$itemtype == "FCprobitems" | RTDataFC$itemtype == "FCimprbitem")

# recode itemtype to code for whether an item is testing understanding of probable or improabable items
dataFCRTP$probItemtype = .5 # probable items coded as .5
dataFCRTP$probItemtype[dataFCRTP$itemtype == 'FCimprbitem'] = -.5 # improbable items coded as -.5

# see how RT on these probabilistic co-occurrence trials differs by condition
mFCRTP = lmer(reactionTime ~ vocabCovariateC + learningCondition*probItemtype + (1 + probItemtype|subjectnr), data = dataFCRTP)
```

```{r, echo = TRUE}
modelSummary(mFCRTP) # table 2 in SOM-U
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pFCRTP <- modelPredictions(mFCRTP, pLearningConditionCovariateProbItemtype)

# now let's look at accuracy on error monitoring probabilistic co-occurrence trials
dataEMACP = subset(allDataEM, allDataEM$itemtype == "EMimpr" | allDataEM$itemtype == "EMprob")

# recode itemtype to code for whether an item is testing understanding of probable or improabable items
dataEMACP$probItemtype = .5 # probable items coded as .5
dataEMACP$probItemtype[dataEMACP$itemtype == 'EMimpr'] = -.5 # improbable items coded as -.5

# see how accuracy on these probabilistic co-occurrence trials differs by condition
mEMACP = glmer(correctanswer ~ vocabCovariateC + learningCondition*probItemtype  + (1 + probItemtype|subjectnr), data=dataEMACP, family = "binomial")
```

```{r, echo = TRUE}
summary(mEMACP)  # table 3 in SOM-U
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pEMACP <- modelPredictions(mEMACP, pLearningConditionCovariateProbItemtype)

# next, let's look at RT on these same probabilistic co-occurrence trials from the EM task
dataEMRTP = subset(RTDataEM, RTDataEM$itemtype == "EMimpr" | RTDataEM$itemtype == "EMprob") 

# recode itemtype to code for whether an item is testing understanding of probable or improabable items
dataEMRTP$probItemtype = .5 # probable items coded as .5
dataEMRTP$probItemtype[dataEMRTP$itemtype == 'EMimpr'] = -.5 # improbable items coded as -.5

# see how accuracy on these probabilistic co-occurrence trials differs by condition
mEMRTP = lmer(reactionTime ~ vocabCovariateC + learningCondition*probItemtype  + (1 + probItemtype|subjectnr), data=dataEMRTP)
```

```{r, echo = TRUE}
modelSummary(mEMRTP)  # table 4 in SOM-U
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
# for making figures later on, save model predictions
pEMRTP <- modelPredictions(mEMRTP, pLearningConditionCovariateProbItemtype)

pFCACP$Itemtype = varRecode(pFCACP$probItemtype,c(-0.5,0.5),c('Improbable', 'Probable'))

legendplot <- ggplot(pFCACP, aes(x = learningCondition, y = Predicted, fill = Itemtype), width=0) + # set data frame, x and y vars
  geom_bar(mapping = aes(fill = as.factor(Itemtype),x=learningCondition,y=Predicted), # set the color based on dichotomous var
           data = pFCACP, # set bars to the predicted values
           stat = "identity", # this means you are giving the data
           position=position_dodge(.9), color = "black") + # width of the bars
  theme_bw(base_size = 12) +  # remove the grey background     
  guides(fill=guide_legend(title="Itemtype"))+
  coord_cartesian(ylim = c(1,3))+  # set the scale so it's the same as the other plots
  scale_fill_manual(values=c("darkblue", "lightblue"))+
  theme(aspect.ratio=4/2)


get_legend<-function(myggplot){
  tmp <- ggplot_gtable(ggplot_build(myggplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)
}

legend <- get_legend(legendplot)

pFCACP$conditionStr = varRecode(pFCACP$learningCondition,c(-0.5,0.5),c('Comprehension', 'Production'))

FCACP <- ggplot(pFCACP, aes(x = conditionStr, y = Predicted, fill = Itemtype), width=0) +
  geom_bar(mapping = aes(fill = as.factor(Itemtype),x=conditionStr,y=Predicted), data = pFCACP, stat = "identity", position=position_dodge(.9), color = "black") + 
  geom_errorbar(data = pFCACP, width=.25, aes(y = Predicted, x = conditionStr, ymin = CILo,ymax = CIHi), stat="identity", width=0.75, position = position_dodge(.9)) + 
  labs(y = 'Proportion Correct') + 
  theme_bw(base_size = 15) +  
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0.8,1))+  
  theme(axis.title.x=element_blank(), axis.ticks.x=element_blank())+
  scale_fill_manual(values=c("darkblue", "lightblue"))+
  theme(aspect.ratio=1/1)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())


pEMACP$conditionStr = varRecode(pEMACP$learningCondition,c(-0.5,0.5),c('Comprehension', 'Production'))
pEMACP$Itemtype = varRecode(pEMACP$probItemtype,c(-0.5,0.5),c('Improbable', 'Probable'))

EMACP <- ggplot(pEMACP, aes(x = conditionStr, y = Predicted, fill = Itemtype), width=0) + 
  geom_bar(mapping = aes(fill = as.factor(Itemtype),x=conditionStr,y=Predicted), data = pEMACP, stat = "identity",position=position_dodge(.9),color = "black") + 
  geom_errorbar(data = pEMACP, width=.25, aes(y = Predicted, x = conditionStr,  ymin = CILo,  ymax = CIHi), stat="identity", width=0.75, position = position_dodge(.9)) + 
  labs(y = 'Proportion Correct') + 
  theme_bw(base_size = 15) +   
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0.8,1))+  
  theme(axis.title.x=element_blank(), axis.ticks.x=element_blank())+
  scale_fill_manual(values=c("darkblue", "lightblue"))+
  theme(aspect.ratio=1/1)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())


pFCRTP$conditionStr = varRecode(pFCRTP$learningCondition,c(-0.5,0.5),c('Comprehension', 'Production'))
pFCRTP$Itemtype = varRecode(pFCRTP$probItemtype,c(-0.5,0.5),c('Improbable', 'Probable'))

FCRTP <- ggplot(pFCRTP, aes(x = conditionStr, y = Predicted, fill = Itemtype), width=0) + 
  geom_bar(mapping = aes(fill = as.factor(Itemtype),x=conditionStr,y=Predicted), data = pFCRTP, stat = "identity", position=position_dodge(.9), color = "black") + 
  geom_errorbar(data = pFCRTP, width=.25, aes(y = Predicted, x = conditionStr, ymin = CILo, ymax = CIHi), stat="identity", width=0.75, position = position_dodge(.9)) + 
  labs(y = 'RT (s)') + 
  theme_bw(base_size = 15) +   
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(1,1.6))+ 
  theme(axis.title.x=element_blank(), axis.ticks.x=element_blank())+
  scale_fill_manual(values=c("darkblue", "light blue"))+
  geom_segment(aes(x=0.75, y=1.37, xend=1.25, yend=1.37)) +
  geom_segment(aes(x=1, y=1.37, xend=1, yend=1.5)) +
  geom_segment(aes(x=1, y=1.5, xend=2, yend=1.5)) +
  geom_segment(aes(x=2, y=1.5, xend=2, yend=1.25)) +
  geom_segment(aes(x=1.75, y=1.25, xend=2.25, yend=1.25)) +
  annotate("text", x=1.5, y=1.5, label="*", size=10) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(aspect.ratio=1/1)



pEMRTP$conditionStr = varRecode(pEMRTP$learningCondition,c(-0.5,0.5),c('Comprehension', 'Production'))
pEMRTP$Itemtype = varRecode(pEMRTP$probItemtype,c(-0.5,0.5),c('Improbable', 'Probable'))

EMRTP <- ggplot(pEMRTP, aes(x = conditionStr, y = Predicted, fill = Itemtype), width=0) + 
  geom_bar(mapping = aes(fill = as.factor(Itemtype),x=conditionStr,y=Predicted), data = pEMRTP, stat = "identity", position=position_dodge(.9), color = "black") + 
  geom_errorbar(data = pEMRTP, width=.25, aes(y = Predicted, x = conditionStr,  ymin = CILo,  ymax = CIHi), stat="identity", width=0.75, position = position_dodge(.9)) + 
  labs(y = 'RT (s)') + 
  theme_bw(base_size = 15) +      
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(1,1.6))+ 
  theme(axis.title.x=element_blank(), axis.ticks.x=element_blank())+
  scale_fill_manual(values=c("darkblue", "lightblue"))+
  geom_segment(aes(x=0.75, y=1.52, xend=1.25, yend=1.52)) +
  geom_segment(aes(x=1, y=1.52, xend=1, yend=1.58)) +
  geom_segment(aes(x=1, y=1.58, xend=2, yend=1.58)) +
  geom_segment(aes(x=2, y=1.58, xend=2, yend=1.29)) +
  geom_segment(aes(x=1.75, y=1.29, xend=2.25, yend=1.29)) +
  annotate("text", x=1.5, y=1.58, label="*", size=10) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(aspect.ratio=1/1)


# put all four plots together in one figure

ProbPlot <- ggdraw(xlim=c(0,2.1),ylim=c(0,2.4)) +
  draw_plot(FCACP, 0, 1, 1, 1) + draw_plot(EMACP, 0.85, 1, 1, 1)+
  draw_plot(FCRTP, 0, 0, 1, 1) + draw_plot(EMRTP, 0.85, 0, 1, 1)


PP <- ProbPlot + draw_plot(legend,1.1,1.63,1,1)


PPF <- PP + draw_label("Forced Choice",x=0.3, y=2.2, vjust=0, hjust = 0, size = 15, fontfamily = "Arial")+
  draw_label("Error Monitoring",x=1.0, y=2.2, vjust=0, hjust = 0, size = 15, fontfamily = "Arial")


PPFL <- PPF + draw_label("a)",x=0.11, y=1.05, vjust=0, hjust = 0, size = 15, fontfamily = "Arial")+
  draw_label("b)",x=0.11, y=0.05, vjust=0, hjust = 0, size = 15, fontfamily = "Arial")+
  draw_label("c)",x=0.96, y=1.05, vjust=0, hjust = 0, size = 15, fontfamily = "Arial")+
  draw_label("d)",x=0.96, y=0.05, vjust=0, hjust = 0, size = 15, fontfamily = "Arial")
```

```{r, echo = TRUE}
PPFL # figure 2 in SOM-U
```

```{r, echo = FALSE, warning = FALSE, results = 'hide'}
ggsave("ProbPlotFinal.png")
# 1139 by 655 works 

```


 
