rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Sternberg Search
# Greg Francis
# PSY 626
# March 27, 2024


# load the rethinking library
library(rethinking)

# load full data file
SSdata<-read.csv(file="SternbergSearch.csv",header=TRUE,stringsAsFactors=FALSE)
SSdata$TargetPresent <- ifelse(SSdata$Condition =="Present", 1, 0)  


# Stan model for null (just intercept), with shrinkage
SSdataNull <- data.frame(RT= SSdata$RT, Participant=SSdata$Participant)		
SSmodelNull <- quap(
			alist( RT ~ dnorm(mu, sigma), 
			mu <- a[Participant],
			a[Participant] ~ dnorm(1000, 200),
			sigma ~ dunif(0, 1000)
			), data= SSdataNull  ) 
				
cat("Finished SSmodelNull\n")

# Stan model for main effect of MemorySetSize, with shrinkage
SSdataMSS <- data.frame(RT= SSdata$RT, MemorySetSize=SSdata$MemorySetSize, Participant=SSdata$Participant)		
SSmodelMSS <- quap(
			alist( RT ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]*MemorySetSize,
			a[Participant] ~ dnorm(1000, 200),
			b[Participant] ~ dnorm(0, 100),			
			sigma ~ dunif(0, 1000)
			), data= SSdataMSS ) 
				
cat("Finished SSmodelMSS\n")



# Stan model for main effect of Target condition
SSdataCondition <- data.frame(RT= SSdata$RT, Condition=SSdata$TargetPresent, Participant=SSdata$Participant)		
SSmodelCondition <- quap(
			alist( RT ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]*Condition,
			a[Participant] ~ dnorm(1000, 200),
			b[Participant] ~ dnorm(0, 100),
			sigma ~ dunif(0, 1000)
			), data= SSdataCondition ) 
				
cat("Finished SSmodelCondition\n")


# Stan model for additive effects 
SSdataAdditive <- data.frame(RT= SSdata$RT, Condition=SSdata$TargetPresent, MemorySetSize=SSdata$MemorySetSize, Participant=SSdata$Participant)		
SSmodelAdditive <- quap(
			alist( RT ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]*MemorySetSize + c[Participant]*Condition,
			a[Participant] ~ dnorm(1000, 200),
			b[Participant] ~ dnorm(0, 100),
			c[Participant] ~ dnorm(0, 100),
			sigma ~ dunif(0, 1000)
			), data= SSdataAdditive ) 
				
cat("Finished SSmodelAdditive\n")

# Stan model for additive effects and interaction
SSdataInteraction <- data.frame(RT= SSdata$RT, Condition=SSdata$TargetPresent, MemorySetSize=SSdata$MemorySetSize, Participant=SSdata$Participant)		
SSmodelInteraction <- quap(
			alist( RT ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]*MemorySetSize + c[Participant]*Condition+ d[Participant]*MemorySetSize*Condition,
			a[Participant] ~ dnorm(1000, 200),
			b[Participant] ~ dnorm(0, 100),
			c[Participant] ~ dnorm(0, 100),
			d[Participant] ~ dnorm(0, 100),
			sigma ~ dunif(0, 1000)
			), data= SSdataInteraction ) 
				
cat("Finished SSmodelInteraction\n")


# Stan model for interaction only (no additive effects)
SSdataInteraction <- data.frame(RT= SSdata$RT, Condition=SSdata$TargetPresent, MemorySetSize=SSdata$MemorySetSize, Participant=SSdata$Participant)		
SSmodelInteractionOnly <- quap(
			alist( RT ~ dnorm(mu, sigma), 
			mu <- a[Participant] + d[Participant]*MemorySetSize*Condition,
			a[Participant] ~ dnorm(1000, 200),
			d[Participant] ~ dnorm(0, 100),
			sigma ~ dunif(0, 1000)
			), data= SSdataInteraction ) 
	
cat("Finished SSmodelInteractionOnly\n")
				
	
# Stan model for Memory Set size and interaction
SSdataInteraction <- data.frame(RT= SSdata$RT, Condition=SSdata$TargetPresent, MemorySetSize=SSdata$MemorySetSize, Participant=SSdata$Participant)		
SSmodelInteractionMMS <- quap(
			alist( RT ~ dnorm(mu, sigma), 
			mu <- a[Participant] + b[Participant]*MemorySetSize + d[Participant]*MemorySetSize*Condition,
			a[Participant] ~ dnorm(1000, 200),
			b[Participant] ~ dnorm(0, 100),
			d[Participant] ~ dnorm(0, 100),
			sigma ~ dunif(0, 1000)
			), data= SSdataInteraction ) 
					
cat("Finished SSmodelInteractionMMS\n")

# Stan model for target condition and interaction
SSdataInteraction <- data.frame(RT= SSdata$RT, Condition=SSdata$TargetPresent, MemorySetSize=SSdata$MemorySetSize, Participant=SSdata$Participant)		
SSmodelInteractionCondition <- quap(
			alist( RT ~ dnorm(mu, sigma), 
			mu <- a[Participant] + c[Participant]*Condition+ d[Participant]*MemorySetSize*Condition,
			a[Participant] ~ dnorm(1000, 200),
			c[Participant] ~ dnorm(0, 100),
			d[Participant] ~ dnorm(0, 100),
			sigma ~ dunif(0, 1000)
			), data= SSdataInteraction ) 
	
cat("Finished SSmodelInteractionCondition\n")
				
		
print(compare(SSmodelInteractionCondition, SSmodelInteractionMMS, SSmodelInteractionOnly, SSmodelInteraction, SSmodelAdditive , SSmodelCondition, SSmodelMSS, SSmodelNull))