#this is GCA code for analysis of Frequency x Neighborhood x Cohort visual world
#paradigm data from Magnuson et al., 2007

#make nlme library available
library(nlme)

#read in experiment data
data<-read.table("FNCdata.txt",col.names=c("subjID", "FreqCond", "NebCond", "CohCond", "time", "fixProp"))
#make a subject-by-condition variable
data$subcond<-factor(data$subjID*1000+data$FreqCond*100+data$NebCond*10+data$CohCond)
#treat condition and subject variables as factors, not ints
data$FreqCond<-factor(data$FreqCond)
data$NebCond<-factor(data$NebCond)
data$CohCond<-factor(data$CohCond)
data$subjID<-factor(data$subjID)

#scale fixProp up by 1000
#this can help avoid floating point problems and convergence issues without changing the substantive results
data$fixProp<-data$fixProp*1000

#create 3rd order orthogonal polynomial time variable
#there are 33 time steps in the data
t<-poly(1:33,3) 
#make a data-time sized array of orthogonal times
ot<-t[data$time,1:3]
#put the time values into separate variables in the experiment data frame
data$ot1<-ot[,1]
data$ot2<-ot[,2]
data$ot3<-ot[,3]

#base model for examining change in -2LL 
m.base<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3, 
data=data, random= ~ot1 | subcond, method="ML")
m.base$logLik*-2 #-2LL

#add frequency effect on intercept
m.f0<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ FreqCond, 
data=data, random= ~ot1 | subcond, method="ML")
m.f0$logLik*-2 #-2LL
#add frequency effect on slope
m.f1<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ FreqCond + FreqCond:ot1, 
data=data, random= ~ot1 | subcond, method="ML")
m.f1$logLik*-2 #-2LL
#frequency effect on slope, no intercept term
m.f1b<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ FreqCond:ot1, 
data=data, random= ~ot1 | subcond, method="ML")
m.f1b$logLik*-2 #-2LL

#add cohort effect on intercept
m.c0<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ CohCond, 
data=data, random= ~ot1 | subcond, method="ML")
m.c0$logLik*-2 #-2LL
#add cohort effect on slope
m.c1<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ CohCond + CohCond:ot1, 
data=data, random= ~ot1 | subcond, method="ML")
m.c1$logLik*-2 #-2LL
#cohort effect on slope, no intercept term
m.c1b<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ CohCond:ot1, 
data=data, random= ~ot1 | subcond, method="ML")
m.c1b$logLik*-2 #-2LL

#add neighborhood effect on intercept
m.n0<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ NebCond, 
data=data, random= ~ot1 | subcond, method="ML")
m.n0$logLik*-2 #-2LL
#add neighborhood effect on slope
m.n1<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ NebCond + NebCond:ot1, 
data=data, random= ~ot1 | subcond, method="ML")
m.n1$logLik*-2 #-2LL
#neighborhood effect on slope, no intercept term
m.n1b<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ NebCond:ot1, 
data=data, random= ~ot1 | subcond, method="ML")
m.n1b$logLik*-2 #-2LL