rm(list=ls())
#
# 0 Preliminaries
#
library(nlme)
#
library(lattice)
#
set.seed(1234321)
#
old.settings <- trellis.par.get()        
#
trellis.par.set(theme=col.whitebg())      
#
# trellis.par.set(theme = old.settings) 
#
grid::grid.prompt(TRUE)
ls()
taf <- groupedData(PC ~ PT | ID, inner = ~ STOR, 
#
             data=read.table("mutaf96.age.dat", header=TRUE, dec="."), 
#
             labels=list(x = "Time (s)", y = "p(Correct)") )
taf
taf$AGE <- taf$AGE - 1
#
#codes AGE as 0 = young, 1 = old 
#
taf$AGE <- as.factor(taf$AGE)
#
taf$ID <- as.factor(taf$ID)
#
taf$STOR <- as.factor(taf$STOR)
taf$s1 <- ifelse(taf$STOR==1, 1, 0)
#
taf$s2 <- ifelse(taf$STOR==2, 1, 0)#
taf$s3 <- ifelse(taf$STOR==3, 1, 0)#
taf$s4 <- ifelse(taf$STOR==4, 1, 0)#

#
# 3 Outlier removal; 7 data points
#
taf <- taf[taf$PC > .05,]
source("Fmutaf1rate.R")
mutaf.lis <- nlsList(PC ~ Fmutaf1rate(PT, STOR, s1, s2, s3, s4, C, RATE, STD) | ID,
#
              data=taf, start=c(C=0.2, RATE=1.2, STD=0.1519348))
mutaf.lis
mutaf0 <- nlme(mutaf.lis,
#
              fixed = list(C+RATE+STD ~ 1), random = list(C ~ 1), verbose=TRUE)
mutaf1 <- update(mutaf0, fixed = list(C ~ AGE, RATE+STD ~ 1), random = list(C ~ 1), 
#
              start = c(0.2, 0.05, 0.7, 0.13))
mutaf2 <- update(mutaf1, fixed = list(C ~ AGE, RATE+STD ~ 1), random = pdDiag(list(C+RATE+STD ~ 1)), 
#
              start = fixef(mutaf1))
mutaf3 <- update(mutaf2, fixed = list(C+RATE ~ AGE, STD ~ 1), random = pdDiag(list(C+RATE+STD ~ 1)), 
#
              start = c(0.2, 0.05, 0.7, 0, 0.13))
