rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Proactive Interference
# Greg Francis
# PSY 626
# 12 September 2025

# Bayesian dependent t-test


# load data file
PIdata<-read.csv(file="ProactiveInterference.csv",sep="\t", header=TRUE,stringsAsFactors=TRUE)

# Dummy variable to indicate condition
PIdata$ConditionIndex <-0* PIdata$NumCorrect
PIdata$ConditionIndex[PIdata$Trial =="1"] =1
PIdata$ConditionIndex[PIdata$Trial =="5"] =2



#----------------------
# Run traditional dependent t-test

FirstTrialData = subset(PIdata, Trial=="1")
FifthTrialData = subset(PIdata, Trial=="5")
traditional <- t.test(FirstTrialData$NumCorrect, FifthTrialData$NumCorrect,  paired=TRUE)

print(traditional)


#------------------------
# Bayesian models

# load the rethinking library
library(rethinking)

# Null model: one mean for both conditions
# Different baseline (intercept) for each participant
PImodel0 <- quap(
			alist( NumCorrect ~ dnorm(mu, sigma), 
			mu <- a[Participant],
			a[Participant] ~ dnorm(5, 5),
			sigma ~ dunif(0, 7)
			), data= PIdata , control=list(maxit=10000) ) 
cat("Finished PImodel0\n")

# Alternative model 1: intercept is mean for trial 1, slope is deviation for trial 5	

PIdata$IsFifth <-0* PIdata$Trial
PIdata$IsFifth[PIdata$Trial =="5"] =1

PImodel1 <- quap(
			alist( NumCorrect ~ dnorm(mu, sigma), 
			mu <- a[Participant] +b*IsFifth,
			a[Participant] ~ dnorm(5, 5),
			b ~ dnorm(5, 5),
			sigma ~ dunif(0, 7)
			), data= PIdata , control=list(maxit=10000) ) 

cat("Finished PImodel1\n")


# Compare models

print(compare(PImodel0, PImodel1) )