#This code reproduces the behavioural results reported in #Kievit, R.A., Scholte, H.S., Waldorp, L. J. & Borsboom, D. (submitted).Inter- and intra-individual differences in fluid reasoning show distinct cortical responses #Date uploaded: 10-02-2016. The four associated raw datafiles are also in the figshare folder. #setwd("Your WD") #Load packages # install.packages("ltm") # install.packages("eRm") # install.packages("lattice") # install.packages("ggplot2") library(ltm) library(eRm) library(lattice) library(ggplot2) library(BayesFactor) #####DATA MANAGEMENT #Import behavioural data data=read.csv ("Accuracydata.csv", header=T) #Accuracy raschdata=data[,3:74] #Import RT data rtdatafull=read.csv('RTdata.csv',header=T) rtdatafull[rtdatafull==0] <- NA #Replace 0 rt's with NA (0 means no response within timeframe) rtdata=rtdatafull[,3:74] #Import IRT constraints constraints=read.csv("IRTconstraints.csv", header=F) constraints[,1]<-1:72 constraints<-as.matrix(constraints) betamanual=constraints[,2] #Beta parameters based on Raven's Matrix manual (Raven, Court, & Raven, 1996) #Import demographics Demos=read.csv('Demographics.csv') sum(Demos$Sex=='f') #Number of females range(Demos$Age) mean(Demos$Age) #Mean age sd(Demos$Age) #SD age #Descriptives correct/incorrect responses mean(rowSums(raschdata)) #Mean items correct range(rowSums(raschdata)) #Range items correct sd(rowSums(raschdata)) #SD items correct range(betamanual) #Range of difficulties based on manual #Descriptives Reaction time mean(rowMeans(rtdata,na.rm=T)) # mean average RT sd(rowMeans(rtdata,na.rm=T)) #SD average RT #Mean RT range(rtdata,na.rm=T) #RT range at the item level sd(as.matrix(rtdata),na.rm=T) #RT sd ###### Rasch analysis with constrained beta parameters Raschfit=rasch(raschdata, constraint=constraints, IRT.param = TRUE) Ability=factor.scores(Raschfit,raschdata)$score.dat$z1 #Fit rasch model without constraints Fitnoconstraints<-rasch(raschdata) #Andersen likelihood ratio test. Excludes items iteratively with full 0 or 1 responses w=RM(raschdata) LRtest(w) #Correlate estimated betas with betas from manual estimatedbeta=-1*Fitnoconstraints$coefficients[,1] #Times -1 because item easiness is estimated estimatedbeta[c(3,72)]=5 #Fixing two extreme outliers ((-25 and +25) items to 5 cor.test(estimatedbeta,betamanual,method='spearman') #Correlate theta estimated with or without constraints: virtually identical cor.test(Ability,factor.scores(Fitnoconstraints,raschdata)$score.dat$z1) #Create matrix with correct responses, difficulty, rt and nullresponses correctvector=c(t(raschdata)) difficultyvector=c(t(matrix(rep(betamanual,times=34),nrow=34,ncol=72,byrow=T))) rtvector=c(t(rtdata)) nullresponse=as.numeric(is.na(rtvector)) cors=data.frame(correctvector,difficultyvector,rtvector,nullresponse) #relate RT, correct and nullresponses to difficulty cor.test(cors$difficultyvector,cors$rtvector,method='spearman') #Difficulty and RT cor.test(cors$difficultyvector,cors$correctvector,method='spearman') #Difficulty and correct responses cor.test(cors$difficultyvector,cors$nullresponse,method='spearman') #Difficulty and probability of null response #Plotting #Figure 2 plot(Raschfit, annot=F, col=rainbow(172)[1:100],lwd=2,xlim=c(-5,5)) #plot ICC qplot(Ability,binwidth = .45)+theme(text=element_text(size=22))+xlim(-4,4)+theme_bw() ggsave('theta.eps') ggsave('theta.tiff') ################ COMPUTE CORRECTED ITEMS ################ rep.row<-function(x,n){ matrix(rep(x,each=n),nrow=n) } rep.col<-function(x,n){ matrix(rep(x,each=n), ncol=n, byrow=TRUE) } allsubstemp=read.csv('Accuracydata.csv',header=T) #Grab all accuracy data allsubs=allsubstemp[,1:74] # strip header ordered=allsubs[order(allsubs[,2]),1:74] #Order people by ability #calculating relative difficulties betas=constraints[,2] thetasmat=rep.col(ordered[,2],72) betasmat=rep.row(unlist(betas),34) adjustedbetas=betasmat-thetasmat #This matrix gives the betas conditional on theta, i.e. the relative difficulty of each item for each person interval=30 #Define number of items to include newmat=matrix(NA,nrow=34,ncol=interval) #create newmat[1,]=adjustedbetas[1,5:(4+interval)] newmat[2,]=adjustedbetas[2,7:(6+interval)] newmat[3,]=adjustedbetas[3,8:(7+interval)] newmat[4,]=adjustedbetas[4,15:(14+interval)] newmat[5,]=adjustedbetas[5,18:(17+interval)] newmat[6,]=adjustedbetas[6,19:(18+interval)] newmat[7,]=adjustedbetas[7,20:(19+interval)] newmat[8,]=adjustedbetas[8,21:(20+interval)] newmat[9,]=adjustedbetas[9,23:(22+interval)] newmat[10,]=adjustedbetas[10,24:(23+interval)] newmat[11,]=adjustedbetas[11,24:(23+interval)] newmat[12,]=adjustedbetas[12,25:(24+interval)] newmat[13,]=adjustedbetas[13,25:(24+interval)] newmat[14,]=adjustedbetas[14,27:(26+interval)] newmat[15,]=adjustedbetas[15,27:(26+interval)] newmat[16,]=adjustedbetas[16,28:(27+interval)] newmat[17,]=adjustedbetas[17,29:(28+interval)] newmat[18,]=adjustedbetas[18,29:(28+interval)] newmat[19,]=adjustedbetas[19,31:(30+interval)] newmat[20,]=adjustedbetas[20,32:(31+interval)] newmat[21,]=adjustedbetas[21,32:(31+interval)] newmat[22,]=adjustedbetas[22,33:(32+interval)] newmat[23,]=adjustedbetas[23,33:(32+interval)] newmat[24,]=adjustedbetas[24,34:(33+interval)] newmat[25,]=adjustedbetas[25,35:(34+interval)] newmat[26,]=adjustedbetas[26,35:(34+interval)] newmat[27,]=adjustedbetas[27,35:(34+interval)] newmat[28,]=adjustedbetas[28,35:(34+interval)] newmat[29,]=adjustedbetas[29,36:(35+interval)] newmat[30,]=adjustedbetas[30,39:(38+interval)] newmat[31,]=adjustedbetas[31,40:(39+interval)] newmat[32,]=adjustedbetas[32,40:(39+interval)] newmat[33,]=adjustedbetas[33,41:(40+interval)] newmat[34,]=adjustedbetas[34,42:(41+interval)] meansdiff=rowMeans(newmat) #Mean corrected difficulty sddiff=apply(newmat,1,sd) #SD corrected difficulty #Computing comparability of corrected difficulties: Bayes factors below 1 indicates evidence in favour of 'no difference in corrected item difficulty' bayesfactorsmatrix=matrix(NA,34,33) freqmattrix=matrix(NA,34,33) for (i in 1:34) { compdat = newmat[i,] for (j in 1:33) { compmat = newmat[-i,] bayesttest = ttestBF(compdat,compmat[j,]) freqtesttemp =t.test(compdat,compmat[j,]) freqtestp=freqtesttemp$p.value bayesfactemp = exp(bayesttest@bayesFactor$bf) bayesfactorsmatrix[i,j] = bayesfactemp freqmattrix[i,j] = freqtestp } } mean(bayesfactorsmatrix) max(bayesfactorsmatrix)