rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Serial Position
# Greg Francis
# PSY 646
# 15 February 2024

# load full data file
SPdata<-read.csv(file="SerialPosition.csv",header=TRUE,stringsAsFactors=FALSE)

## Re-organize with one score for each row
Participant<-c()
Position<-c()
Score<-c()
# Loop through each participant
for( i in c(1:109)){
	# loop through each position
	for(j in c(1:10)){
		Participant<-c(Participant, SPdata$Participant[i])
		Position<-c(Position, j)
		Score<-c(Score, SPdata[i, (j+1)])
	}
}

SPdata2 <- data.frame(Participant<-Participant, Position<-Position, Score<-Score)
colnames(SPdata2) <-c("Participant", "Position", "Score")

# load the rethinking library
library(rethinking)
		
# Positions: 			
SPsimple <- quap(
			alist( Score ~ dbinom(15, p), 
			logit(p) <- a[Position],
			a[Position] ~  dnorm(0, 10)
			), data= SPdata2  ) 
					
	
modelCoefficients <- coef(SPsimple)		
aValuesSimple<-c()
for(p in unique(SPdata2$Position)){
	code<-sprintf("a[%d]", p)
	aValuesSimple <- c(aValuesSimple, modelCoefficients[code])	
}
cat("Position probability of correct response: ", logistic(aValuesSimple))

dev.new()
sp.sequence <- c(1:10)
plot(sp.sequence, logistic(aValuesSimple))

logitPEst<-link(SPsimple, data=data.frame(Position=sp.sequence))
logitPEst.mean <- apply(logitPEst, 2, mean)
logitPEst.HPDI <-apply(logitPEst, 2, HPDI, prob=0.89)
lines(sp.sequence, logitPEst.mean)
shade(logitPEst.HPDI, sp.sequence)

# Curve for positions: different means and sds for each condition			
SPcurve <- quap(
			alist( Score ~ dbinom(15, p), 
			logit(p) <- a + b*(Position-c)*(Position-c),
			a ~  dnorm(0, 10),
			b ~ dnorm(0, 10),
			c ~ dunif(2, 9)
			), data= SPdata2 ) 
			
save(SPcurve, file="SPcurve.Rpd")		

dev.new()
# Get the posterior
post <- extract.samples(SPcurve)	
dens(post$c)
cat("Bottom of curve posterior: Mean=", mean(post$c), " sd=", sd(post$c))			
			

dev.new()
# Plot points for data means
dataMeans<-c()
for(j in c(1:10)){
	subData <- subset(SPdata2, SPdata2$Position==j)
	dataMeans<-c(dataMeans, mean(subData$Score))
}
plot(sp.sequence, dataMeans/15)

# Add model curve fit and HPDI
logitPEst<-link(SPcurve, data=data.frame(Position=sp.sequence))
logitPEst.mean <- apply(logitPEst, 2, mean)
logitPEst.HPDI <-apply(logitPEst, 2, HPDI, prob=0.89)
lines(sp.sequence, logitPEst.mean)
shade(logitPEst.HPDI, sp.sequence)



modelCoefficients <- coef(SPcurve)			
probValueCurve<-c()
for(p in unique(SPdata2$Position)){
	code<-sprintf("a[%d]", p)
	probValueCurve <- c(probValueCurve, (modelCoefficients["a"]+ modelCoefficients["b"]* (p-modelCoefficients["c"])*(p-modelCoefficients["c"])))	
}

cat("Curve probability of correct response: ", logistic(probValueCurve))