#this is GCA code for analysis of Cohort-Rhyme-Unrelated visual world
#paradigm data from Magnuson et al., 2003

#make nlme library available
library(nlme)

#read in experiment data
data<-read.table("CRUdata.txt",col.names=c("subjID", "day", "time", "targFreq", "compFreq", "Cond", "targFix", "compFix", "unrelFix"))
#remove time steps after 1400ms
data<-subset(data,data$time<43)
#make a subject-by-condition variable
data$subcond<-factor(data$subjID*1000+data$targFreq*100+data$compFreq*10+data$Cond)

#make the "unrelated" condition the default comparison condition
data$Cond<-replace(data$Cond,data$Cond==3,0)
#treat condition and subject variables as factors, not ints
data$Cond<-factor(data$Cond)
data$subjID<-factor(data$subjID)

#create 4th order orthogonal polynomial time variable
t<-poly(0:42,4)
#make a data-time sized array of orthogonal times, time adjusted for excluded values
ot<-t[data$time+1,1:4] #+1 reflects a change in indexing from 0 to indexing from 1
#put the time values into separate variables in the experiment data frame
data$ot1<-ot[,1]
data$ot2<-ot[,2]
data$ot3<-ot[,3]
data$ot4<-ot[,4]

#DAY 1 ANALYSIS
#gradually build-up model and monitor change in -2LL
#base model
m.base<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4, #subject effects
data=subset(data,data$day==1), random= ~ot1 | subcond, method="ML")
m.base$logLik*-2 #-2LL
#add condition effect on intercept
m.0<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond,
data=subset(data,data$day==1), random= ~ot1 | subcond, method="ML")
m.0$logLik*-2 #-2LL
#add condition effect on slope
m.1<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1,
data=subset(data,data$day==1), random= ~ot1 | subcond, method="ML")
m.1$logLik*-2 #-2LL
#add condition effect on quad
m.2<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2,
data=subset(data,data$day==1), random= ~ot1 | subcond, method="ML")
m.2$logLik*-2 #-2LL
#add condition effect on cube
m.3<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2+Cond:ot3,
data=subset(data,data$day==1), random= ~ot1 | subcond, method="ML")
m.3$logLik*-2 #-2LL
#add condition effect on quart
m.4<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2+Cond:ot3+Cond:ot4,
data=subset(data,data$day==1), random= ~ot1 | subcond, method="ML")
m.4$logLik*-2 #-2LL
#get parameter estimates and t-test results
summary(m.4)$tTable

#DAY 2 ANALYSIS
#gradually build-up model and monitor change in -2LL
#base model
m.base2<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4, #subject effects
data=subset(data,data$day==2), random= ~ot1 | subcond, method="ML")
m.base2$logLik*-2 #-2LL
#add condition effect on intercept
m.0<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond,
data=subset(data,data$day==2), random= ~ot1 | subcond, method="ML")
m.0$logLik*-2 #-2LL
#add condition effect on slope
m.1<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1,
data=subset(data,data$day==2), random= ~ot1 | subcond, method="ML")
m.1$logLik*-2 #-2LL
#add condition effect on quad
m.2<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2,
data=subset(data,data$day==2), random= ~ot1 | subcond, method="ML")
m.2$logLik*-2 #-2LL
#add condition effect on cube
m.3<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2+Cond:ot3,
data=subset(data,data$day==2), random= ~ot1 | subcond, method="ML")
m.3$logLik*-2 #-2LL
#add condition effect on quart
m.4<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2+Cond:ot3+Cond:ot4,
data=subset(data,data$day==2), random= ~ot1 | subcond, method="ML")
m.4$logLik*-2 #-2LL
#get parameter estimates and t-test results
summary(m.4)$tTable

#DIRECT COMPARISON OF COHORT AND RHYME FIXATION
dataRC<-subset(data,as.integer(data$Cond)>1) #note: in factor-to-integer conversion unrelated goes from 0 to 1
#DAY 1 ANALYSIS
#gradually build-up model and monitor change in -2LL
#base model
m.base<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4, #subject effects
data=subset(dataRC,dataRC$day==1), random= ~ot1 | subcond, method="ML")
m.base$logLik*-2 #-2LL
#add condition effect on intercept
m.0<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond,
data=subset(dataRC,dataRC$day==1), random= ~ot1 | subcond, method="ML")
m.0$logLik*-2 #-2LL
#add condition effect on slope
m.1<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1,
data=subset(dataRC,dataRC$day==1), random= ~ot1 | subcond, method="ML")
m.1$logLik*-2 #-2LL
#add condition effect on quad
m.2<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2,
data=subset(dataRC,dataRC$day==1), random= ~ot1 | subcond, method="ML")
m.2$logLik*-2 #-2LL
#add condition effect on cube
m.3<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2+Cond:ot3,
data=subset(dataRC,dataRC$day==1), random= ~ot1 | subcond, method="ML")
m.3$logLik*-2 #-2LL
#add condition effect on quart
m.4<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2+Cond:ot3+Cond:ot4,
data=subset(dataRC,dataRC$day==1), random= ~ot1 | subcond, method="ML")
m.4$logLik*-2 #-2LL
#DAY 2 ANALYSIS
#gradually build-up model and monitor change in -2LL
#base model
m.base<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4, #subject effects
data=subset(dataRC,dataRC$day==2), random= ~ot1 | subcond, method="ML")
m.base$logLik*-2 #-2LL
#add condition effect on intercept
m.0<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond,
data=subset(dataRC,dataRC$day==2), random= ~ot1 | subcond, method="ML")
m.0$logLik*-2 #-2LL
#add condition effect on slope
m.1<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1,
data=subset(dataRC,dataRC$day==2), random= ~ot1 | subcond, method="ML")
m.1$logLik*-2 #-2LL
#add condition effect on quad
m.2<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2,
data=subset(dataRC,dataRC$day==2), random= ~ot1 | subcond, method="ML")
m.2$logLik*-2 #-2LL
#add condition effect on cube
m.3<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2+Cond:ot3,
data=subset(dataRC,dataRC$day==2), random= ~ot1 | subcond, method="ML")
m.3$logLik*-2 #-2LL
#add condition effect on quart
m.4<-lme(compFix~ot1+ot2+ot3+ot4 #time effects
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3+subjID:ot4 #subject effects
+ Cond+Cond:ot1+Cond:ot2+Cond:ot3+Cond:ot4,
data=subset(dataRC,dataRC$day==2), random= ~ot1 | subcond, method="ML")
m.4$logLik*-2 #-2LL