# Note that the files used in the following analyses are not provided.  All are read into the
# variable 'myd'.

library(nlme)

# Throughout, you can use "summary(myr)" and "coef(myr)" to see results where "myr" is the 
# variable holding the results of the analysis being summarized.


################################################################################
# Testing of SIMPLE fits of Kaplan & Reed data
################################################################################
# Individual fit
myrList<-nlsList(Value~500/(1+exp(logk)*Delay)|Subject, data=myd, start=c(logk=-2))

# Multilevel fit
myr<-nlme(Value~500/(1+exp(logk)*Delay), fixed=logk~1, random=logk~1|Subject, data=myd, start=c(logk=-2))

################################################################################
# Testing of SIMULATED data set.  Original dataset without noise is myd.  Noisy dataset is mynd.
################################################################################
mynd<-myd  
mynd$ran<-runif(dim(mynd)[[1]]) # add a column to be used for random selection

ProportionToOmit <- .0  # Note, no missing data in these analyses
AmountOfNoise <- 25 # sd of indifference points; set to 200 for the bottom of Figure 4.

mynd$Indifference<- mynd$Indifference + rnorm(dim(mynd)[[1]], sd=AmountOfNoise) # used to add noise
mynd$Indifference<- ifelse(mynd$Indifference>1000, 1000, floor(mynd$Indifference)) # ensure that values are not greater than max
mynd$Indifference<- ifelse(mynd$Indifference<0, 0, mynd$Indifference) # ensure that values are not greater than max
mynd$Indifference<-ifelse(mynd$ran<ProportionToOmit, NA, mynd$Indifference) # used to add missing values

myr<-nlme(Indifference~1000/(1+exp(logk)*Delay), fixed=logk~1, random=logk~1|Subject, data=mynd, start=c(logk=-2))
myrList<-nlsList(Indifference~1000/(1+exp(logk)*Delay)|Subject, data=mynd, start=c(logk=-2))


################################################################################
# Testing of SIMULATED data set with between-subject CONDITION
################################################################################
# Same code as above for preparing data except with ProportionToOmit = .30 and
# AmountOfNoise = 200.  Second set of simulations set ProportionToOmit = .50.
myr<-nlme(Indifference2~1000/(1+exp(logk)*Delay), fixed=logk~Condition, random=logk~1|Subject, data=mynd, na.action=na.omit, start=c(logk=-3, 0))

################################################################################
# Testing of HYPERBOLOID fits of Kaplan & Reed data
################################################################################
#### Individual fits
myrList<-nlsList(Value~500/(1+exp(logk)*Delay) | Subject, data=myd, start=c(logk=-2))
coef(myrList)

# Myerson/Green
myrListMG<-nlsList(Value~500/(1+exp(logk)*Delay)^s | Subject, data=myd, start=c(logk=-2, s=1))
coef(myrListMG)

# Rachlin
myrListR<-nlsList(Value~500/(1+exp(logk)*Delay^s) | Subject, data=myd, start=c(logk=-2, s=1))
coef(myrListR)

#### Multilevel fits
myr<-nlme(Value~500/(1+exp(logk)*Delay), fixed=logk~1, random=logk~1|Subject, data=myd, start=c(k=-2))

# Needed to increase error tolerance to get MG to produce a fit
myrMG<-nlme(Value~500/(1+exp(logk)*Delay)^s, fixed=logk+s~1, random=logk+s~1|Subject, data=myd, start=c(logk=-2, s=1), control=c(nlmeControl(tol=.5)))

myrRach<-nlme(Value~500/(1+exp(logk)*Delay^s), fixed=logk+s~1, random=logk+s~1|Subject, data=myd, start=c(logk=-2, s=1))

AIC(myr, myrMG, myrRach)