# Example for estimating time-accuracy functions with nlsList and nlme with data from
# Oberauer & Kliegl (2006, Journal of Memory and Language, 55,601-626)
#
rm(list=ls())

# 
library(nlme)
library(lattice)
set.seed(1234321)

load("exampleTAF.rda")
str(tafx)
# tafx:	92 obs. of  9 variables:
# $ ID    : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...  # Subject
# $ STOR  : Factor w/ 2 levels "2","3": 1 1 1 1 1 1 1 1 1 1 ...          # Memory load
# $ PT    : num  0.272 0.308 0.596 0.668 0.776 ...                       # Presentation time
# $ PC    : num  0.237 0.289 0.321 0.607 0.792 ...                       # Probability correct
# $ NITEMS: int  19 19 14 14 12 14 16 17 17 17 ...                       # N of trials

# Treatment contrasts for storage
tafx$s2 <- ifelse(tafx$STOR=="2", 1, 0)
tafx$s3 <- ifelse(tafx$STOR=="3", 1, 0)

# Base estimation
tafx.lis  <- nlsList(PC ~ SSasympOff(PT, c, b, a) | ID, data=tafx)
tafx.0 <- nlme(tafx.lis)
        
# Subject-specific fits
tafx.1 <- update(tafx.0,
                fixed = list(a ~ s3, b ~ s3, c ~ s3),
                random = list(a ~ 1, b+c ~ s3-1),
                start = c(0.03, 0.24, 0.22, 1.27, 1.00, -0.19),
                verbose = TRUE )
        
# Predictions using conditonal modes of random effects (formerly known as BLUPs)
# ... generate data.frame to hold predictions
k <-  4     # n of ids
l <-  2     # n of storage demands
n <- 121    # n obs per id x storage demand = 0 : 6 s pres time
N = k*l*n

tafpred <- expand.grid(ID=1:k, STOR=2:3, PT=seq(0,6,length.out=n))
tafpred$s2   <- ifelse(tafpred$STOR==2, 1, 0)
tafpred$s3   <- ifelse(tafpred$STOR==3, 1, 0)
tafpred$ID   <- as.factor(tafpred$ID)
tafpred$STOR <- as.factor(tafpred$STOR)

# ... generate predictions
tafpred$PC <- predict(tafx.1,tafpred,level=1)
tafpred <- tafpred[tafpred$PC>.05,]

# Draw graph
tafobs <- tafx[,c(1:3,6:7,4)]

tafx2 <- rbind(tafobs,tafpred)
tafx2$op <- c(rep(10,nrow(tafobs)),rep(20,nrow(tafpred) ))
tafx2$op <- as.factor(tafx2$op+as.numeric(tafx2$STOR))

xyplot(PC ~ PT | ID, data=tafx2, groups = op, layout=c(2,2),
        panel = panel.superpose.2,
        type = c("p","p","l","l"), col=c(1:2), lty = 1)