library(gplots)
library(MASS)

#Preliminary analyses when we've ran 20 participants
read.table("exp03_words_prel.txt",header=T)->ld

#First, we check that the data was read ok
summary(ld)
str(ld)
head(ld)

#next we classify into right and left handed
ld$hnd_bin<-ifelse(ld$hnd>50,"R","L")

#next we classify into mono and bi handed
ld$lang_bin<-ifelse(ld$lng>50,"mono","bi")

#Check that all the participants are ok
boxplot(ld$rt~ld$ppnt) #ppnt 1 has some outlier really long response times -- on task? maybe exclude?
	#p4 has REALLY fast time, who is this guy??? I'm taking him out, he wasn't even thinking...
	
ld[ld$ppnt!="p_4",]->ld

hist(table(ld$resp,ld$ppnt)[1,]/240*100) #it looks like people make 4-8% errors, that's great

#now we take only correct trials
ld[ld$resp==1,]->ld

aggregate(ld$rt,by=list(ld$hnd_bin,ld$lang_bin,ld$job_type,ld$val,ld$ppnt,ld$key),mean)->meanrts
colnames(meanrts)<-c("hnd","lang","job","val","ppnt","key","rt")

#There are too many variables, so we start exploring the data taking 2 or 3 into account at a time

#Looking at the effects of valence, key, and handedness
plotmeans(rt~val,data=meanrts,pch="L",subset=c(meanrts$hnd=="L" & meanrts$key=="L"),main="left key")
plotmeans(rt~val,data=meanrts,pch="R",subset=c(meanrts$hnd=="R"& meanrts$key=="L"),add=T)

plotmeans(rt~val,data=meanrts,pch="L",subset=c(meanrts$hnd=="L" & meanrts$key=="R"),main="right key")
plotmeans(rt~val,data=meanrts,pch="R",subset=c(meanrts$hnd=="R"& meanrts$key=="R"),add=T)
#hmm looks like right handers are overall faster than left handers
#and looks like valence didn't have any effect?? 
#I also don't see any effect/interaction with key
# so let's collapse across valence and key since they didn't matter -- actually, let's clean up and have only the variables we are looking at in the table of averages
aggregate(ld$rt,by=list(ld$hnd_bin, ld$ppnt),mean)->temp
colnames(temp)<-c("hnd","ppnt","rt")

plotmeans(rt~hnd,data= temp,main="speed by handedness")


t.test(rt~hnd,data= temp) #not yet significant, perhaps with more subjects

#We keep exploring, now looking at language
plotmeans(rt~val,data=meanrts,pch="L",subset=c(meanrts$hnd=="L" & meanrts$lang=="mono"),main="monolinguals")
plotmeans(rt~val,data=meanrts,pch="R",subset=c(meanrts$hnd=="R"& meanrts$lang=="mono"),add=T)

plotmeans(rt~val,data=meanrts,pch="L",subset=c(meanrts$hnd=="L" & meanrts$lang=="bi"),main="bilinguals")
plotmeans(rt~val,data=meanrts,pch="R",subset=c(meanrts$hnd=="R"& meanrts$lang=="bi"),add=T)
#this analysis replicates the speed difference between left and right handers, and the fact that there are no effects on valence
#if the latter were true, this would be really interesting - it would mean that mono and bilinguals process valence the same!
#let's look at language collapsing across everything else
aggregate(ld$rt,by=list(ld$lang_bin, ld$ppnt),mean)->temp
colnames(temp)<-c("lang","ppnt","rt")

plotmeans(rt~lang,data=temp,main="language comparison")

#Actually no, there IS a language difference: monolinguals are overall faster than bilinguals -- makes perfect sense, but let's check with real statistics

t.test(rt~lang,data=temp) #not sig either, but getting there - cool!

#The only thing left to explore is the type of profession, but since valence and key*handedness didn't matter, I guess it won't play any role - but let's check
plotmeans(rt~val,data=meanrts,pch="L",subset=c(meanrts$hnd=="L" & meanrts$job==0),main="don't work with hands")
plotmeans(rt~val,data=meanrts,pch="R",subset=c(meanrts$hnd=="R"& meanrts$job==0),add=T)

plotmeans(rt~val,data=meanrts,pch="L",subset=c(meanrts$hnd=="L" & meanrts$job==1),main="work with hands")
plotmeans(rt~val,data=meanrts,pch="R",subset=c(meanrts$hnd=="R"& meanrts$job==1),add=T)
#How interesting, the handedness difference in RT only exists for people who don't work with their hands!
#lets check with proper stats

aggregate(ld$rt,by=list(ld$hnd_bin,ld$job_type,ld$ppnt),mean)->temp
colnames(temp)<-c("hnd","job","ppnt","rt")

summary(aov(rt~hnd*job,data= temp))
#weird -- there is no interaction... but perhaps it's lack of power, so I'll collect some more data

boxplot(rt~val*job,data=temp,subset=c(temp$hnd=="R")) #no, it seems not...


#We're done running subjects! I can do the final analyses now
read.table("exp03_words.txt",header=T)->ld

#As usual clean up and add cols first
ld[ld$ppnt!="p_4",]->ld
ld[ld$resp==1,]->ld
ld$hnd_bin<-ifelse(ld$hnd>50,"R","L")
ld$lang_bin<-ifelse(ld$lng>50,"mono","bi")

#then get participant-based averages
aggregate(ld$rt,by=list(ld$hnd_bin,ld$lang_bin,ld$job_type,ld$val,ld$ppnt,ld$key),mean)->meanrts
colnames(meanrts)<-c("hnd","lang","job","val","ppnt","key","rt")

#Following Peter's advice, I use a linear regression and stepwise deletion to find the right model
fit<-lm(rt~hnd*lang*job*val*key,data= meanrts)
stepAIC(fit)

best=lm(formula = rt ~ hnd + lang + job + hnd:lang + hnd:job + lang:job , data = meanrts)
summary(best)
#TERRIFIC this model is significant and accounts for about 11 percent of variance - not too shabby!!
#OK, so it looks like:
	#hndR negative estimate: Right-handers are, as I suspected above, FASTER than left-handers
	#langmono negative estimate: again, as I suspected, monolinguals are FASTER than bilinguals
	#also people who work with their hands are faster -- all of this is making perfect sense
	#there are also two interactions, I'm not sure how to explore them so I'll ask Peter
	
#Yikes Peter said I should have declared repeated measures in that analysis, or averaged over valence and key since they don't play a role
aggregate(ld$rt,by=list(ld$hnd_bin,ld$lang_bin,ld$job_type,ld$ppnt),mean)->meanrts
colnames(meanrts)<-c("hnd","lang","job","ppnt","rt")

fit<-lm(rt~hnd*lang*job,data= meanrts)
stepAIC(fit)
#what, the effect goes away?


#maybe I should have exclude those data, let's start over without exclusions
read.table("exp03_words.txt",header=T)->ld
ld$hnd_bin<-ifelse(ld$hnd>50,"R","L")
ld$lang_bin<-ifelse(ld$lng>50,"mono","bi")

aggregate(ld$rt,by=list(ld$hnd_bin,ld$lang_bin,ld$job_type,ld$ppnt),mean)->meanrts
colnames(meanrts)<-c("hnd","lang","job","ppnt","rt")

fit<-lm(rt~hnd*lang*job,data= meanrts)
stepAIC(fit)
best=lm(formula = rt ~ hnd + lang, data = meanrts)
summary(best)
#what??? still not significant?
#I don't get it -- let's inspect the data once more

#then get participant-based averages
aggregate(ld$rt,by=list(ld$hnd_bin,ld$lang_bin,ld$job_type,ld$val,ld$ppnt,ld$key),mean)->meanrts
colnames(meanrts)<-c("hnd","lang","job","val","ppnt","key","rt")

par(mar=c(10,4,1,0)+0.1)
boxplot(rt~hnd*lang*job*val*key,data=meanrts,las=2) #I definitely see a recurrent pattern... 
boxplot(rt~hnd*lang*val*key,data=meanrts,las=2) #yes, its there
boxplot(rt~hnd*lang*key,data=meanrts,las=2) #yes, its there
boxplot(rt~hnd*lang,data=meanrts,las=2) #yes, its DEFINITELY there, why doesn't it come out??

#oh maybe it's this repeated measures problem again??
aggregate(ld$rt,by=list(ld$hnd_bin,ld$lang_bin,ld$ppnt),mean)->meanrts
colnames(meanrts)<-c("hnd","lang","ppnt","rt")
boxplot(rt~hnd*lang,data=meanrts,las=2) #no, it's DEFINITELY there!!


t.test(meanrts$rt[meanrts$hnd=="L"& meanrts$lang=="bi"],meanrts$rt[meanrts$hnd=="R"& meanrts$lang=="mono"],alternative="greater")#still not significant
sum(meanrts$hnd=="L"& meanrts$lang=="bi")#oooh I see - there are only 4 participants in this group
#so we should collect more data