data<-read.table("Exp3 data.txt",header=T)
str(data)


data$cb<-as.factor(data$cb)
data$subject<-as.factor(data$subject)
data$code<-as.factor(data$code)
data$pos1<-as.factor(data$pos1)
data$trial<-as.factor(data$trial)
data$occurrence<-as.factor(data$occurrence)
str(data)


table(data$subject)
with(data, aggregate(acc=="cor",by=list(subject=subject),mean))->subacc
plot(subacc$x)
	
subacc$subject[subacc$x<0.75]->under75


#high subject accuracy

subset(data,!subject%in%under75)->dh


round(dim(dh[!dh$acc=="cor",])[1]/dim(dh)[1],1)

subset(dh, acc=="cor"&rt>150&rt<1500&trial_type=="same")->trimmed	#same trials trimmed

round(dim(subset(dh,trial_type=="same"&(rt<150|rt>1500)))[1]/dim(subset(dh,trial_type=="same"))[1]*100,2)#0.47% outliers 

round(dim(subset(dh,trial_type=="same"&acc!="cor"&rt>150&rt<1500))[1]/dim(subset(dh,trial_type=="same"))[1]*100,2) #7.43% incorrect responses


subset(dh, acc=="cor"&rt>150&rt<1500&trial_type=="different")->trimmedDif #different trials trimmed

round(dim(subset(dh,trial_type=="different"&(rt<150|rt>1500)))[1]/dim(subset(dh,trial_type=="different"))[1]*100,2)#0.7% outliers 

round(dim(subset(dh,trial_type=="different"&acc!="cor"&rt>150&rt<1500))[1]/dim(subset(dh,trial_type=="different"))[1]*100,2) #11.21% incorrect responses 

with(trimmed,xtabs(rt~prime_type+trial_type)/xtabs(~prime_type+trial_type))
round(with(trimmed,xtabs(rt~prime_type+trial_type)/xtabs(~prime_type+trial_type)))

round(with(trimmedDif,xtabs(rt~prime_type+trial_type)/xtabs(~prime_type+trial_type)))


library(lme4)
library(phia)



###taken from https://github.com/lme4/lme4/issues/98
###for improved optimization

library(nloptr)
defaultControl <- list(algorithm="NLOPT_LN_BOBYQA",xtol_rel=1e-6,maxeval=1e5)
nloptwrap2 <- function(fn,par,lower,upper,control=list(),...) {
    for (n in names(defaultControl)) 
      if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]]
    res <- nloptr(x0=par,eval_f=fn,lb=lower,ub=upper,opts=control,...)
    with(res,list(par=solution,
                  fval=objective,
                  feval=iterations,
                  conv=if (status>0) 0 else status,
                  message=message))
}			


#RT same trials

lmer(rt~prime_type+(prime_type|subject)+(prime_type|trial),data=trimmed,control=lmerControl(optimizer="nloptwrap2"))->m1_nloptr #converged
Anova(m1_nloptr)

testInteractions(m1_nloptr,pairwise=c("prime_type"),adjustment="none")


#Accuracy same trials

100-round(with(dh,xtabs(acc=="cor"~prime_type+trial_type)/xtabs(~prime_type+trial_type))*100,1)

glmer(acc=="cor"~prime_type+(prime_type|subject)+(prime_type|trial),data=subset(dh,trial_type=="same"),family="binomial",control=glmerControl(optimizer="nloptwrap2"))->m1_nloptrAcc #converged

Anova(m1_nloptrAcc)

testInteractions(m1_nloptrAcc,pairwise=c("prime_type"),adjustment="none")

 

#RT different trials

#subset(dh, acc=="cor"&rt>150&rt<1500&trial_type=="different")->trimmedDif

with(trimmedDif,xtabs(rt~prime_type+trial_type)/xtabs(~prime_type+trial_type))


lmer(rt~prime_type+(1|subject)+(prime_type|trial),data=trimmedDif,control=lmerControl(optimizer="nloptwrap2"))->m1_nloptrDiff #Converged, the full model failed to converge.

Anova(m1_nloptrDiff)


#Accuracy different trials

glmer(acc=="cor"~prime_type+(prime_type|subject)+(prime_type|trial),data=subset(dh,trial_type=="different"),family="binomial",control=glmerControl(optimizer="bobyqa"))->m1_AccDif #Converged

Anova(m1_AccDif)
testInteractions(m1_AccDif,pairwise=c("prime_type"),adjustment="none")

