# Example for estimating TAFs with nlmer (lme4) with data from
# Oberauer & Kliegl (2006, Journal of Memory and Language, 55,601-626)
#
rm(list=ls());

# 
library(lme4)
library(ggplot2)
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)

# Change from prob correct to percentage correct
tafx$PC <- tafx$PC*100

xyplot(PC ~ PT | ID, tafx, groups = STOR, 
       auto.key = list(lines = TRUE, columns = 2),
       layout = c(4, 1),
       type = "o")

xyplot(PC ~ PT | ID, tafx, groups = STOR,
       auto.key = list(lines = TRUE, columns = 2),
       layout = c(4, 1),
       type = c("g","p","smooth"))

# Base estimations
(tafx.0 <- nlmer(PC ~ SSasympOff(PT, c, b, a) ~ (c + b  | ID), tafx,
				start = c(c = 100, b = 0.50, a = 0.10), verb=1) )

(tafx.0.1 <- nlmer(PC ~ SSlogis(PT, Asym,  xmid, scal) ~ (Asym + scal | ID), tafx,
				start = c(Asym=80, xmid=0.50, scal = 0.50), verb=1) )

# Base estimation via function Ftaf.V2.R, built after SSasympOff(), incl.gradient attribute
source("Ftaf.V2.R")  
(tafx.0.2 <- nlmer(PC ~ Ftaf(PT, c, b, a) ~  (c + b | ID), tafx,
				start = c(c = 100, b = 0.50, a = 0.10), verb=1) )
							
# Subject-specific effects: HOW WOULD I SPECIFY THESE in nlmer?
#  ... delta c: dc , if s3
source("Ftaf.V3.R")
(tafx.0.3 <- nlmer(PC ~ Ftaf(PT, s3, c, b, a, dc) ~  (c + b + dc | ID), tafx,
				start = c(c = 100, b = 0.50, a = 0.10, dc= -15), verb=FALSE) )

# ... delta b: db, if s3
source("Ftaf.V4.R")
(tafx.0.4 <- nlmer(PC ~ Ftaf(PT, s3, c, b, a, dc, db) ~  (c + b | ID) + (dc + db | ID), tafx,
				start = c(c = 100, b = 0.50, a = 0.10, dc= -15, db= -.20), verb=1) )

# Generate simple graph; predict() not available yet, but fitted values are
tafx$fit0 <- fitted(tafx.0.4,level=1)
tafL <- reshape(tafx,timevar="DV",drop=c("s2","s3"),
        varying=list(c("PC","fit0")),direction="long")  

tafL$DVS <- as.factor(tafL$DV*10+as.numeric(tafL$STOR))

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