############################################################################################################
# Loading the first data set for testing.  
############################################################################################################
myd<-read.csv("Simulated Choice data revised simplified.dat")

# Wileyto analysis to derive k for second stage
myd$ip1<-1-myd$LLmagn/myd$SSmagn
myrWileyto<-by(myd, myd$Subject, function(x) coef(glm(Waited~ip1+LLdelay-1, data=x, family=binomial)))
summary(myrWileyto)

# The following commands perform a logistic regression on each subject and extracts the parameters.
# This looks complicated because I'm extracting information for a bunch of individual regressions.  
# The multilevel approach is much more straightforward.
library(plyr)
fitList <- dlply(myd, .(Subject), glm, formula=Waited~ip1+LLdelay-1, family=binomial) 
jack<-as.data.frame(t(sapply(fitList,coef)))
jack$k<-jack$LLdelay/jack$ip1
jack$Subject<-as.numeric(row.names(jack))
jack$logk<-log(jack$k+1) # problem with some negative k values; this captures all but one of them

# Setting it up and testing for a group comparison
jack$Group<-ifelse(jack$Subject < 26, 1, 0)
jackr<-lm(logk~Group, data=jack)
summary(jackr)
densityplot(~logk|Group, data=jack)

wilcox.test(k~Group, data=jack) # nonparametric test of full set of 50 k values

######  Multilevel regression of the same data - the recommended analysis
library(lme4)
myr<-glmer(Waited~log(LLmagn/SSmagn)+log(LLdelay+1)-1+(log(LLmagn/SSmagn)+log(LLdelay+1)-1|Subject), data=myd, family=binomial)
summary(myr)

# Generating plots like those that appear in the paper. 
plot(betaM~I(-betaD), data=myd, xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]), xlim=c(-4, .4), ylim=c(0,5))
plot(coef(myr)$Subject[[1]]~coef(myr)$Subject[[2]], cex=jack$logk+.5, xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]), xlim=c(-4, .4), ylim=c(0,5))
# four points aren't being plotted due to NaN or strong negative k
coefSub <- coef(myr)$Subject[c(1,4,9),]
points(coefSub[[1]]~coefSub[[2]], pch = 3)


jill<- data.frame(aggregate(myd$betaM, by=list(myd$Subject), FUN=mean),aggregate(myd$betaD, by=list(myd$Subject), FUN=mean))
colnames(jill)<-c("Subject", "betaM", "S", "betaD")
jill<-jill[c(1,2,4)]
jill$fit.betaM <- coef(myr)$Subject[[1]]
jill$fit.betaD <- coef(myr)$Subject[[2]]

# Plot estimated by actual regression weights
par(mfrow=c(1,2))
plot(jill$fit.betaM~jill$betaM, xlim=c(0,5), ylim=c(0,5), xlab=expression(beta["Magnitude"]), ylab=expression(paste("Fitted ", beta["Magnitude"])))
plot(function(x) x, 0,5, add=T, lty=2)
plot(jill$fit.betaD~ I(-jill$betaD), xlim=c(-4, .25), ylim=c(-4,.25), xlab=expression(beta["Delay"]), ylab=expression(paste("Fitted ", beta["Delay"])))
plot(function(x) x, -4,.25, add=T, lty=2)

# Same plots but for fits of each individual rather than a multilevel model. 
myrindivid<-by(myd, myd$Subject, function(x) coef(glm(Waited~log(LLmagn/SSmagn)+log(LLdelay+1)-1, data=x, family=binomial)))
summary(myrindivid)
fitList <- dlply(myd, .(Subject), glm, formula=Waited~log(LLmagn/SSmagn)+log(LLdelay+1)-1, family=binomial) 
julie<-as.data.frame(t(sapply(fitList,coef)))
julie$Subject<-as.numeric(row.names(julie))
summary(julie)

julie$betaM<-jill$betaM
julie$betaD<-jill$betaD
plot(julie[[1]]~jill$betaM, xlim=c(0,5), ylim=c(0,6), xlab=expression(beta["Magnitude"]), ylab=expression(paste("Fitted ", beta["Magnitude"])))
plot(function(x) x, 0,6, add=T, lty=2)
plot(julie[[2]]~ I(-jill$betaD), xlim=c(-4, .25), xlab=expression(beta["Delay"]), ylab=expression(paste("Fitted ", beta["Delay"])))
plot(function(x) x, -4,.25, add=T, lty=2)

plot(coef(myr)$Subject[[1]]~row.names(coef(myr)$Subject))
plot(coef(myr)$Subject[[2]]~as.numeric(row.names(coef(myr)$Subject))%%5)

# Testing ability to identify strength of an effect that was held constant
myr1.1<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 1), family=binomial)
myr1.2<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 2), family=binomial)
myr1.3<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 4), family=binomial)
myr1.4<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 8), family=binomial)
myr1.5<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 16), family=binomial)

summary(myr1.1)
summary(myr1.2)
summary(myr1.3)
summary(myr1.4)
summary(myr1.5)


############################################################################################################
# Testing extension of above to make a group comparison (Subjects 1-25 vs. 26-50)
############################################################################################################
myd$Group<-ifelse(myd$Subject < 26, 1, 0)

# Multilevel Fit
myr2<-glmer(Waited~Group*log(LLmagn/SSmagn)+Group*log(LLdelay+1)-1-Group+(log(LLmagn/SSmagn)+log(LLdelay+1)-1|Subject), data=myd, family=binomial)
summary(myr2)
xyplot(fitted(myr2)~LLdelay|Group, data=myd, type="a")

# Group comparison for individually derived weights generated in earlier analysis.
julie$Group<-ifelse(julie$Subject < 26, 1, 0)
colnames(julie)[1]<-"fit.betaM"
colnames(julie)[2]<-"fit.betaD"
julierM<-lm(fit.betaM~Group, data=julie)
julierD<-lm(fit.betaD~Group, data=julie)
summary(julierD)
summary(julierM)
densityplot(~fit.betaM|Group, data=julie)

# Try again by omitting subject 50
julierM<-lm(fit.betaM~Group, data=subset(julie, Subject!=50))
summary(julierM)
densityplot(~fit.betaM|Group, data=julie)

############################################################################################################
# Trying a new data set with nonzero SS delays.
############################################################################################################
myd2<-read.csv("Simulated Choice data revised simplified nonzero.dat")

myr2<-glmer(Waited~log(LLmagn/SSmagn)+log((LLdelay+1)/(SSdelay+1))-1+(log(LLmagn/SSmagn)+log((LLdelay+1)/(SSdelay+1))-1|Subject), data=myd2, family=binomial)
summary(myr2)

xyplot(fitted(myr2)~SSdelay, group=Subject, data=myd2, type="a", ylab="P(LL)")

############################################################################################################
# Analysis of Peterson et al. (2015) data set
############################################################################################################
myd3<-read.csv("ESchoiceGE.txt")  # I do not have permission to release this data, but the syntax is below

myr3<-glmer(Choice~log(LL.magnitude/SS.magnitude)+log((LL_Delay+1)/3)-1+(log(LL.magnitude/SS.magnitude)+log((LL_Delay+1)/3)-1|Subject), data=myd3, family=binomial)

# The following is functionally equivalent but uses intercept to estimate log(LLmagn/SSmagn) effect.
# Thus, to obtain the same magnitude regression weight as that used in previous analysis, you just 
# multiply the intercept times 1/log(2/1).  
myr3.1<-glmer(Choice~log((LL_Delay+1)/3)+(log((LL_Delay+1)/3)|Subject), data=myd3, family=binomial)
summary(myr3)
summary(myr3.1)

plot(coef(myr3)$Subject[[1]]~coef(myr3)$Subject[[2]])
plot(coef(myr3.1)$Subject[[1]]*1/log(2)~coef(myr3.1)$Subject[[2]], xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]))

# Trying a linear relationship.  The +1 and /3 were dropped because these linear transformations do not alter the fit metrics of a linear model.
myr3.2<-glmer(Choice~LL_Delay+(LL_Delay|Subject), data=myd3, family=binomial)
AIC(myr3.1, myr3.2)
# AIC much worse for the linear model (9985 for log vs. 10044 for linear, a difference of 59)

# Adding Trial to assess learning effects. Create new centered variable for trial.
myd3$logdelay <- log((myd3$LL_Delay+1)/3)
myd3$c.trial <- myd3$Trial - mean(myd3$Trial)

# The following is used if I want to center trial. I'm not centering delay because I need log(1) to predict 50% P(LL) for both magnitude and delay
# ratios of 1.
myr3.31<-glmer(Choice~c.trial*logdelay+(c.trial+logdelay|Subject), data=myd3, family=binomial)

# The next command is a test of a variation of the Wileyto model
# This model assumes a linear relationship between delay and outcome.  It also assumes a linear relationship
# for magnitudes, but these were held constant and thus cannot be distinguished.
myr3.31W<-glmer(Choice~c.trial*LL_Delay+(c.trial+LL_Delay|Subject), data=myd3, family=binomial)
# This model is worse.  

# plotting individuals - the first shows the model fit, the second shows raw data.
xyplot(fitted(myr3.31)~Trial|as.factor(Subject), group=LL_Delay, data=myd3, ylab="P(LL)",
       type = "a", auto.key = list(points = FALSE, lines = TRUE), 
       par.settings = simpleTheme(lty = c(1, 2, 3, 5), lwd=c(2,2,2,2)))
xyplot(Choice~Trial|as.factor(Subject), group=LL_Delay, data=myd3, ylab="P(LL)",
       type = "a", auto.key = list(points = FALSE, lines = TRUE), 
       par.settings = simpleTheme(lty = c(1, 2, 3, 5), lwd=c(2,2,2,2)))


# Use lstrends to estimate slopes (aka simple slopes in some areas of psychology).  Could also use multcomp.
library(lsmeans)
lstrends(myr3.31, ~ c.trial, var="logdelay", at = list(c.trial = c(-26.5, 26.5)))  # to examine slope differences
lsmeans(myr3.31, ~ c.trial, at = list(c.trial = c(-26.5, 26.5), logdelay=1.9265), transform="response")  # to examine slope differences
lsmeans(myr3.31, ~ c.trial, at = list(c.trial = c(-26.5, 26.5), logdelay=1.9265))  # to examine slope differences
lsmeans(myr3.31, ~ c.trial, at = list(c.trial = c(-26.5, 26.5), logdelay=0), transform="response")  # to examine slope differences
lsmeans(myr3.31, ~ c.trial, at = list(c.trial = c(-26.5, 26.5), logdelay=0))  # to examine slope differences

# If I wanted to, I could evaluate whether the trial effect was logarithmic rather than linear.
myr3.33<-glmer(Choice~log(Trial)*logdelay+(log(Trial)+logdelay|Subject), data=myd3, family=binomial)
AIC(myr3.32, myr3.33)
# It wasn't - the linear model was better for this data set. 

# Finally, here's a model testing changes across sessions within a block of similar LL delays.
myd3$c.sessionwithin <- myd3$SessionWithin -3
myr3.4<-glmer(Choice~c.sessionwithin*c.logdelay+(c.sessionwithin+c.logdelay|Subject), data=myd3, family=binomial)
xyplot(fitted(myr3.4)~SessionWithin|as.factor(Subject), group=LL_Delay, data=myd3, type="a", auto.key=T, ylab="P(LL)")
lsmeans(myr3.4, ~ c.sessionwithin, at = list(c.sessionwithin = c(-2, 2)), transform="response")  # to examine slope differences


############################################################################################################
# Another example using random delays to the LL and a SS delay of 0.25 s: Video game data.
############################################################################################################
mydvg <- read.csv("VideoGame DDDG data.dat")
mydvg = subset(mydvg, Pred.Rand=="Predictable") # only focusing on the predictable delays
mydvg$MagRatio <- 100/5
contrasts(mydvg$DD.DG)<-contr.sum(2)            # using effect coding
mydvg$logdelay <- log((mydvg$Delay+1)/1.25)

# Using intercept to estimate constant effect of magnitude, but could include MagRatio and drop intercept. 
myrvg<-glmer(Chose.Later~DD.DG*logdelay+(logdelay|Subject), data=mydvg, family=binomial)
summary(myrvg)

# Examine specific outcomes (simple slopes and means derived from the model)
lstrends(myrvg, ~ DD.DG, var="logdelay")  # to examine slope differences across conditions
lsmeans(myrvg, ~ DD.DG, transform="response")  # to examine intercept differences across conditions
lsmeans(myrvg, ~ DD.DG)  # to examine intercept differences across conditions

plot(coef(myrvg)$Subject[[1]]*1/log(20)~coef(myrvg)$Subject[[3]], xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]))
plot(function(x) -x, -5,5, add=T, lty=2)

# Examining impact of the three subjects who always waited for the LL by omitting them.  Note, cannot compare AIC/BIC with full data
# set because the sample size is now different. 
myrvg2<-glmer(Chose.Later~DD.DG*c.logdelay+(c.logdelay|Subject), data=subset(mydvg, Subject!="07b" & Subject!="07c" & Subject!="35a"), family=binomial)
summary(myrvg)
summary(myrvg2)
plot(coef(myrvg2)$Subject[[1]]*1/log(20)~coef(myrvg2)$Subject[[3]], xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]))
plot(function(x) -x, -5,5, add=T, lty=2)


# Individual fits of video game data
fitList <- dlply(mydvg, .(Subject), glm, formula=Chose.Later~log((Delay+1)/1.25), family=binomial) 
# For individual summaries, use a command like: 
summary(fitList$`07b`)

############################################################################################################
# Jarmolowicz data - I do not have permission to release this data, but the syntax is below
############################################################################################################
mydj <- read.csv("../MoneyDisc/TotalChoice.dat")
summary(mydj)

# Multilevel analysis
myrj<-glmer(Waited~log((LLmagn+1)/(SSmagn+1))+log(LLdelay+1)-1+(log((LLmagn+1)/(SSmagn+1))+log(LLdelay+1)-1|Subject), data=mydj, family=binomial)
summary(myrj)

plot(coef(myrj)$Subject[[1]]~coef(myrj)$Subject[[2]], xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]), xlim=c(-1.3, 0), ylim=c(0,6))

# Indifference point analysis (see Young JEAB 2017? for method)
mydjI <- read.csv("../MoneyDisc/TotalIndifference.dat")
library(nlme)
myrjI<-nlme(IndifferencePoint~10/(1+exp(logk)*Delay), fixed=logk~1, random=logk~1|Subject, data=mydjI, start=c(logk=-2))
summary(myrjI)
plot(coef(myrj)$Subject[[1]]~coef(myrjI)[[1]], xlab=expression(logk), ylab=expression(beta["Magnitude"]), ylim=c(0,6))
plot(coef(myrj)$Subject[[2]]~coef(myrjI)[[1]], xlab=expression(logk), ylab=expression(beta["Delay"]), ylim=c(0,-1.3))
plot(coef(myrj)$Subject[[1]]~coef(myrj)$Subject[[2]], cex=(coef(myrjI)[[1]]+9)*.3, xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]))

mydj$c.trial <- mydj$Trial - mean(mydj$Trial)

# Checking for changes across trials with a new analysis.
myrj3<-glmer(Waited~c.trial*logmagn+c.trial*log(LLdelay+1)-1-c.trial+(c.trial*logmagn+c.trial*log(LLdelay+1)-1-c.trial|Subject), data=mydj, family=binomial)

# I'm having problems wth lstrends when var name is so long, so creating a summary variable
mydj$logmagn <- log((mydj$LLmagn+1)/(mydj$SSmagn+1))

lstrends(myrj3, ~ c.trial, var="log(LLdelay+1)", at = list(c.trial = c(-36, 13)))  # to examine slope differences
lstrends(myrj3, ~ c.trial, var="logmagn", at = list(c.trial = c(-36, 13)))  # to examine slope differences

plot(coef(myrj3)$Subject[[3]]~coef(myrj3)$Subject[[1]], ylab=expression(beta["Magnitude x Trial"]), xlab=expression(beta["Magnitude"]))
plot(coef(myrj3)$Subject[[4]]~coef(myrj3)$Subject[[2]], ylab=expression(beta["Delay x Trial"]), xlab=expression(beta["Delay"]), ylim=c(.008, -.008))