#this is GCA code for individual differences 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 last subject first to match SAS analysis
data$subjID<-replace(data$subjID,data$subjID==115,100)
#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]

#model for examining individual difference effects 
m.id<-lme(fixProp~ot1+ot2+ot3
+ subjID+subjID:ot1+subjID:ot2+subjID:ot3
+ CohCond + CohCond:ot1 + FreqCond + FreqCond:ot1 + NebCond + NebCond:ot1, 
data=data, random= ~ot1 | subcond, method="ML")
m.id$logLik*-2 #-2LL

#compute individual subject fixed effect correlations
#where values will be stored
r<-matrix(-1,nrow=4,ncol=4) #correlation matrix
p<-matrix(-1,nrow=4,ncol=4) #correlation significance matrix
#indices for individual parameter estimates
startInd<-c(5,22,36,50)
stopInd<-c(18,35,49,63)
#build correlation table
for (i in 1:4){
	for (j in 1:4){
		x<-cor.test(c(0,summary(m.id)$tTable[startInd[i]:stopInd[i]]),c(0,summary(m.id)$tTable[startInd[j]:stopInd[j]]))
		r[i,j]<-x$estimate
		p[i,j]<-x$p.value
	}
}
#NOTE: estimates in correlation are padded with 0 for the reference participant
#	individual fixed effects are computed relative to a reference, which gets a 0 estimate, which is not output in the t-Table

#compute random (residual) effects: averaging as described in appendix 
x<-m.id$coefficients
y<-x$random
int<-y$subcond[,1]
freqResInt<-(int[(1:60)*2]+int[((1:60)*2)-1])/2
cohResInt<-c((int[1]+int[5])/2,(int[3]+int[7])/2,(int[2]+int[6])/2,(int[4]+int[8])/2)
for (i in 1:14){
	temp<-c((int[(i*8)+1]+int[(i*8)+5])/2,(int[(i*8)+3]+int[(i*8)+7])/2,(int[(i*8)+2]+int[(i*8)+6])/2,(int[(i*8)+4]+int[(i*8)+8])/2)
	cohResInt<-c(cohResInt,temp)
}
linear<-y$subcond[,2]
freqResLin<-(linear[(1:60)*2]+linear[((1:60)*2)-1])/2
cohResLin<-c((linear[1]+linear[5])/2,(linear[3]+linear[7])/2,(linear[2]+linear[6])/2,(linear[4]+linear[8])/2)
for (i in 1:14){
	temp<-c((linear[(i*8)+1]+linear[(i*8)+5])/2,(linear[(i*8)+3]+linear[(i*8)+7])/2,(linear[(i*8)+2]+linear[(i*8)+6])/2,(linear[(i*8)+4]+linear[(i*8)+8])/2)
	cohResLin<-c(cohResLin,temp)
}

#compute correlations and store in correlation matrix, only top triangle is filled in
rres<-matrix(-1,nrow=4,ncol=4) #correlation matrix
pres<-matrix(-1,nrow=4,ncol=4) #correlation significance matrix
x<-cor.test(cohResInt,cohResLin)
rres[1,2]<-x$estimate
pres[1,2]<-x$p.value
x<-cor.test(cohResInt,freqResInt)
rres[1,3]<-x$estimate
pres[1,3]<-x$p.value
x<-cor.test(cohResInt,freqResLin)
rres[1,4]<-x$estimate
pres[1,4]<-x$p.value
x<-cor.test(cohResLin,freqResInt)
rres[2,3]<-x$estimate
pres[2,3]<-x$p.value
x<-cor.test(cohResLin,freqResLin)
rres[2,4]<-x$estimate
pres[2,4]<-x$p.value
x<-cor.test(freqResInt,freqResLin)
rres[3,4]<-x$estimate
pres[3,4]<-x$p.value