# R-scripts for manuscript:
# Göthe, Esser, Gendt, Kliegl, Dev. Psych. 2011
# Working memory in children: Tracing age differences and special 
# educational needs to parameters of a formal model


# C)  ANALYZING EFFECTS OF STORAGE DEMAND, MATERIAL
# IN WORKING MEMORY, 2 GROUPS OF CHILDREN (CG2 vs. GLD5-6)
# NLME- test difference vs. developmental retardation hypothesis 

# 0 Preliminaries
rm(list=ls())
set.seed(1234321)
library(nlme)

# 1 Input
load("GAAred.rda")
taf<- tafred

taf$group[(taf$ID<200)]=1
taf$group[(taf$ID>200 &taf$ID<300)]=2
taf$group[(taf$ID>300&taf$ID<400)]=3
taf$group[(taf$ID>400)]=4

taf <- groupedData(PC ~ PT | ID, inner = ~ STOR*TASK,
             data=taf,
		 order.groups = F )

# Treatment contrasts for condition
taf$cond[(taf$STOR==1 & taf$TASK==1) ] <-  1
taf$cond[(taf$STOR==2 & taf$TASK==1) ] <-  2
taf$cond[(taf$STOR==1 & taf$TASK==2) ] <-  3
taf$cond[(taf$STOR==2 & taf$TASK==2) ] <-  4


# Treatment contrast for storage (memory demand)
taf$s2 <- 0
taf$s2 [taf$STOR==2] <-  +1

# Treatment contrast for material
taf$m2<- 0			   # verbal
taf$m2[taf$TASK==2] <-  +1 # spatial

# Reduction of data frame, select group 1 & 4
taf<-taf[(taf$group==1 | taf$group==4),]

# Treatment contrasts for Group, testing developmental vs. difference hypothesis (101-120=CG2, 401-419=GLD5-6)
taf$diffdev[(taf$group==1)]<- 0	# CG2
taf$diffdev[(taf$group==4)]<- +1	# GLD5-6

# 4 Base function per ID (specified in taf via groupedData)
source("FIM2.R")#(t, condition, s2, m2, c,  rate, std,diffdev)
mutaf.lis <- nlsList(PC ~ FIM2 (PT, cond, s2, m2, C ,RATE, STD,diffdev) | ID,
              data=taf,start=c(C=0.3, RATE=1, STD=0.15) )

mutaf2.0 <- nlme(mutaf.lis,
            fixed = list(C+ RATE+ STD ~1),
		start = c(0.2,3,  2,),
		random = list(C~1), verbose=TRUE, control=list(maxIter=200))
		summary(mutaf2.0)
		mean(coef(mutaf2.0))
		rs <- FRsquared(mutaf2.0, taf$PC)
		rs
		
mutaf2.1 <- update(mutaf2.0,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2+diffdev),
		start = c(0.75,.03,
		1.5,3,-0.4,0.7, -1.3, 
		0.2,.001,.001),
		random = list(C~1), verbose=TRUE, control=list(maxIter=200))
		summary(mutaf2.1)
		coef(mutaf2.1)
		rs <- FRsquared(mutaf2.1, taf$PC)
		rs
anova(mutaf2.0, mutaf2.1)

# Test for interaction C & diffdev--> no!
mutaf2.2 <- update(mutaf2.0,
            fixed = list(C~m2* diffdev, RATE~m2*s2+diffdev, STD ~m2+diffdev),
		start = c(0.75,.03,0.1,0.1,
		1.5,3,-0.4,0.7, -1.3, 
		0.2,.001,.001),
		random = list(C~1), verbose=TRUE, control=list(maxIter=200))
		summary(mutaf2.2)
		coef(mutaf2.2)
		rs <- FRsquared(mutaf2.2, taf$PC)
		rs
anova(mutaf2.1, mutaf2.2)

# increase random effects
mutaf2.1b <- update(mutaf2.0,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2+diffdev),
		start = c(0.75,.03,
		1.5,3,-0.4,0.7, -1.3, 
		0.2,.001,.001),
		random = pdDiag(c(C~1, STD~1)), verbose=TRUE, control=list(maxIter=200))
		summary(mutaf2.1b)
		mean(coef(mutaf2.1b))
		rs <- FRsquared(mutaf2.1b, taf$PC)
		rs
anova(mutaf2.1, mutaf2.1b)

mutaf2.1c <- update(mutaf2.1b,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2+diffdev),
		start = c(0.557,.3,
		1.5,2.7,-0.37,0.288, -1.6, 
		0.2,.045,-.0227),
		random = pdDiag(c(C~1, STD~1+m2)), verbose=TRUE, control=list(maxIter=300))
		summary(mutaf2.1c)
		coef(mutaf2.1c)
		rs <- FRsquared(mutaf2.1c, taf$PC)
		rs
anova(mutaf2.1b, mutaf2.1c)

mutaf2.1c2 <- update(mutaf2.1b,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2),
		start = c(0.557,.3,
		1.5,2.7,-0.37,0.288, -1.6, 
		0.2,.045),
		random = pdDiag(c(C~1, STD~1+m2)), verbose=TRUE, control=list(maxIter=300))
		summary(mutaf2.1c2)
		coef(mutaf2.1c2)
		rs <- FRsquared(mutaf2.1c2, taf$PC)
		rs
anova(mutaf2.1c, mutaf2.1c2)

mutaf2.1d <- update(mutaf2.1b,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2+diffdev),
		start = c(0.557,.3,
		1.5,2.7,-0.37,0.288, -1.6, 
		0.2,.045,-.0227),
		random = pdDiag(c(C~1, RATE~1, STD~1+m2)), verbose=TRUE, control=list(maxIter=300))
		summary(mutaf2.1d)
		coef(mutaf2.1d)
		rs <- FRsquared(mutaf2.1d, taf$PC)
		rs
anova(mutaf2.1c, mutaf2.1d)


mutaf2.1d3 <- update(mutaf2.1b,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2+diffdev),
		start = c(0.557,.3,
		1.5,2.7,-0.37,0.288, -1.6, 
		0.2,.045,-.0227),
		random = pdDiag(c(C~1, RATE~1, STD~1+m2+diffdev)), verbose=TRUE, control=list(maxIter=300))
		summary(mutaf2.1d3)
		coef(mutaf2.1d3)
		rs <- FRsquared(mutaf2.1d3, taf$PC)
		rs
anova(mutaf2.1d, mutaf2.1d3)

mutaf2.1d2 <- update(mutaf2.1b,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2),
		start = c(0.557,.3,
		1.5,2.7,-0.37,0.288, -1.6, 
		0.2,.045),
		random = pdDiag(c(C~1, RATE~1, STD~1+m2)), verbose=TRUE, control=list(maxIter=300))
		summary(mutaf2.1d2)
		coef(mutaf2.1d2)
		rs <- FRsquared(mutaf2.1d2, taf$PC)
		rs
anova(mutaf2.1d, mutaf2.1d2)
anova(mutaf2.1c2, mutaf2.1d2)

mutaf2.1e2 <- update(mutaf2.1b,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2),
		start = c(0.557,.3,
		1.5,2.7,-0.37,0.288, -1.6, 
		0.2,.045),
		random = pdDiag(c(C~1, RATE~1+m2, STD~1+m2)), verbose=TRUE, control=list(maxIter=300))
		summary(mutaf2.1e2)
		coef(mutaf2.1e2)
		rs <- FRsquared(mutaf2.1e2, taf$PC)
		rs
anova(mutaf2.1d2, mutaf2.1e2)

# Test for correlations
mutaf2.1.correl <- update(mutaf2.0,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2+diffdev),
		start = c(0.75,.03,
		1.5,3,-0.4,0.7, -1.3, 
		0.2,.001,.001),
		random = list(C~1), verbose=TRUE, control=list(maxIter=200))
		summary(mutaf2.1)
		coef(mutaf2.1)
		rs <- FRsquared(mutaf2.1, taf$PC)
		rs

############################################################
# C--> fixed, when sign.--> corresp. random
############################################################
taf1 <- update(mutaf2.0,
            fixed = list(C~m2, RATE+ STD ~1),
		start = c(0.7, 0.4,
		1.5,  
		.2),
		random = pdDiag(c(C~1)), verbose=TRUE, control=list(maxIter=200))
		summary(taf1)
		mean(coef(taf1))
		rs <- FRsquared(taf1, taf$PC)
		rs
anova(mutaf2.0, taf1)

taf2<- update(mutaf2.0,
            fixed = list(C~m2, RATE+ STD ~1),
		start = c(0.559, 0.4,
		3.22,  
		.2),
		random = pdDiag(c(C~1+m2)), verbose=TRUE, control=list(maxIter=200))
		summary(taf1)
		mean(coef(taf2))
		rs <- FRsquared(taf2, taf$PC)
		rs
anova(taf1, taf2)

############################################################
# RATE
############################################################
taf2.0 <- update(mutaf2.0,
            fixed = list(C~1, RATE~1, STD ~1),
		start = c(0.7,
		3.5,  
		.2),
		random = pdDiag(c(RATE~1)), verbose=TRUE, control=list(maxIter=200))
		summary(taf2.0)
		mean(coef(taf2.0))
		rs <- FRsquared(taf2.0, taf$PC)
		rs

taf2.1 <- update(mutaf2.0,
            fixed = list(C~1, RATE~m2, STD ~1),
		start = c(0.7,
		3.5,1.4,  
		.2),
		random = pdDiag(c(RATE~1)), verbose=TRUE, control=list(maxIter=200))
		summary(taf2.1)
		mean(coef(taf2.1))
		rs <- FRsquared(taf2.1, taf$PC)
		rs
anova(taf2.0, taf2.1)

taf2.1b <- update(mutaf2.0,
            fixed = list(C~1, RATE~m2, STD ~1),
		start = c(0.7,
		3.5,1.4,  
		.2),
		random = pdDiag(c(RATE~1+m2)), verbose=TRUE, control=list(maxIter=200))
		summary(taf2.1b)
		mean(coef(taf2.1b))
		rs <- FRsquared(taf2.1b, taf$PC)
		rs
anova( taf2.1, taf2.1b)

taf2.2 <- update(mutaf2.0,
            fixed = list(C~1, RATE~m2+s2, STD ~1),
		start = c(0.7,
		3.5,1.4,-0.5  ,
		.2),
		random = pdDiag(c(RATE~1)), verbose=TRUE, control=list(maxIter=200))
		summary(taf2.2)
		mean(coef(taf2.2))
		rs <- FRsquared(taf2.2, taf$PC)
		rs
anova(taf2.1, taf2.2)

taf2.2b <- update(mutaf2.0,
            fixed = list(C~1, RATE~m2+s2, STD ~1),
		start = c(0.7,
		3.5,1.4,-0.5  ,
		.2),
		random = pdDiag(c(RATE~1+s2)), verbose=TRUE, control=list(maxIter=200))
		summary(taf2.2b)
		mean(coef(taf2.2b))
		rs <- FRsquared(taf2.2b, taf$PC)
		rs
anova(taf2.2, taf2.2b)

taf2.3 <- update(mutaf2.0,
            fixed = list(C~1, RATE~m2*s2, STD ~1),
		start = c(0.237,
		3.5,-.31,-2.95 ,.9 ,
		.2),
		random = pdDiag(c(RATE~1+s2)), verbose=TRUE, control=list(maxIter=200))
		summary(taf2.3)
		mean(coef(taf2.3))
		rs <- FRsquared(taf2.3, taf$PC)
		rs
anova(mutaf2.2b, taf2.3)

############################################################
# NOISE
############################################################
taf3.0 <- update(mutaf2.0,
            fixed = list(C~1, RATE~1, STD ~1),
		start = c(0.7,
		3.5,
		.2),
		random = pdDiag(c(STD~1)), verbose=TRUE, control=list(maxIter=200))
		summary(taf3.0)
		mean(coef(taf3.0))
		rs <- FRsquared(taf3.0, taf$PC)
		rs

taf3.1 <- update(mutaf2.0,
            fixed = list(C~1, RATE~1, STD ~m2),
		start = c(0.7,
		3.5,
		.2,0.1),
		random = pdDiag(c(STD~1)), verbose=TRUE, control=list(maxIter=200))
		summary(taf3.1)
		mean(coef(taf3.1))
		rs <- FRsquared(taf3.1, taf$PC)
		rs
anova(taf3.0, taf3.1)

taf3.2 <- update(mutaf2.0,
            fixed = list(C~1, RATE~1, STD ~m2),
		start = c(0.7,
		3.5,
		.2,0.1),
		random = pdDiag(c(STD~1+m2)), verbose=TRUE, control=list(maxIter=200))
		summary(taf3.2)
		mean(coef(taf3.2))
		rs <- FRsquared(taf3.2, taf$PC)
		rs
anova(taf3.1, taf3.2)

taf3.3 <- update(mutaf2.0,
            fixed = list(C~1, RATE~1, STD ~m2+s2),
		start = c(0.87,
		2.5,
		.189,0.021,0.4),
		random = pdDiag(c(STD~1+m2)), verbose=TRUE, control=list(maxIter=200))
		summary(taf3.3)
		mean(coef(taf3.3))
		rs <- FRsquared(taf3.3, taf$PC)
		rs
anova(taf3.2, taf3.3)

############################################################
# Integration
############################################################

#no convergence--> intger2 reduce randoms
intger1<- update(mutaf2.0,
            fixed = list(C~m2, RATE~m2*s2, STD ~m2),
		start = c(0.559, 0.4,
		3.5,-.31,-2.95 ,.9 ,
		.189,0.021),
		random = pdDiag(c(C~1+m2, RATE~1+s2,STD~1+m2)), 
		verbose=TRUE, control=list(maxIter=500))
		summary(intger1)
		mean(coef(intger1))
		rs <- FRsquared(intger1, taf$PC)
		rs

# no convergence.
intger3<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2, STD ~m2),
		start = c(0.559, 0.4,
		3.5,-.31,-2.95 ,.9 ,
		.189,0.021),
		random = pdDiag(c(C~m2,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=500))
		summary(intger3)
		mean(coef(intger3))
		rs <- FRsquared(intger3, taf$PC)
		rs

intger2<- update(mutaf2.0,
            fixed = list(C~m2, RATE~m2*s2, STD ~m2),
		start = c(0.559, 0.4,
		3.5,-.31,-2.95 ,.9 ,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=500))
		summary(intger2)
		mean(coef(intger2))
		rs <- FRsquared(intger2, taf$PC)
		rs
anova(taf3.2, intger2)

############################################################
# Gruppeneffekt
############################################################

group1<- update(intger2,
            fixed = list(C~m2+diffdev, RATE~m2*s2, STD ~m2),
		start = c(0.559, 0.4,-0.1,
		3.5,-.31,-2.95 ,.9 ,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=500))
		summary(group1)
		mean(coef(group1))
		rs <- FRsquared(group1, taf$PC)
		rs
anova(intger2, group1)

group2<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2),
		start = c(0.559, 0.4,
		3.5,-.31,-2.95 ,.9 ,-0.1,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=500))
		summary(group2)
		mean(coef(group2))
		rs <- FRsquared(group2, taf$PC)
		rs
anova(intger2, group2)

group3<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev+m2*diffdev, STD ~m2),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,0.1,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group3)
		mean(coef(group3))
		rs <- FRsquared(group3, taf$PC)
		rs
anova(group2, group3)

group4<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev, STD ~m2),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~1+diffdev,STD~m2)), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group4)
		mean(coef(group4))
		rs <- FRsquared(group4, taf$PC)
		rs
anova(group2, group4)

group5<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev+s2*diffdev, STD ~m2),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,0.1,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group5)
		mean(coef(group5))
		rs <- FRsquared(group5, taf$PC)
		rs
anova(group2, group5)

group6<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev+s2*diffdev+m2*diffdev, STD ~m2),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,-0.1,-0.1,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group6)
		mean(coef(group6))
		rs <- FRsquared(group6, taf$PC)
		rs
anova(group6, group5)

group7<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev+s2*diffdev+m2*diffdev+m2*s2*diffdev, STD ~m2),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,-0.1,-0.1,-0.1,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group7)
		mean(coef(group7))
		rs <- FRsquared(group7, taf$PC)
		rs
anova(group7, group5)

# no convergence
group8<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev+s2*diffdev, STD ~m2),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,0.1,
		.189,0.021),
		random = pdDiag(c(C~1,RATE~s2,STD~m2)), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group8)
		mean(coef(group8))
		rs <- FRsquared(group8, taf$PC)
		rs

group9<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev+s2*diffdev, STD ~m2+diffdev),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,0.1,
		.189,-0.021, -0.01),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group9)
		mean(coef(group9))
		rs <- FRsquared(group9, taf$PC)
		rs
anova(group5, group9)

group10<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev+s2*diffdev, STD ~m2+diffdev+m2*diffdev),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,0.1,
		.189,-0.021, -0.01, -0.01),
		random = pdDiag(c(C~1,RATE~1,STD~m2)), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group10)
		mean(coef(group10))
		rs <- FRsquared(group10, taf$PC)
		rs
anova(group5, group10)

# group 5= winner --> Correlation test
group5.c<- update(intger2,
            fixed = list(C~m2, RATE~m2*s2+diffdev+s2*diffdev, STD ~m2),
		start = c(0.46, 0.4,
		1.35,2.66,-0.57 ,0.52 ,-1.22,0.1,
		.189,0.021),
		random = list(C~1,RATE~1,STD~m2), 
		verbose=TRUE, control=list(maxIter=600))
		summary(group5.c)
		mean(coef(group5.c))
		rs <- FRsquared(group5.c, taf$PC)
		rs
anova(group5, group5.c)
###############################################

anova(mutaf2.0, taf1, taf2)
anova(mutaf2.0, taf2.1)
anova( taf2.1, taf2.1b)
anova(taf2.1, taf2.2,taf2.2b, taf2.3)
anova(mutaf2.2b, taf2.3)

# noise
anova(mutaf2.0, taf3.1)
anova(taf3.1, taf3.2, taf3.3)
anova(taf3.2, intger2)
anova( group2, group5)
anova(group1, group2, group3, group4, group5, group5.c)


###############################################