rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics

# Greg Francis
# PSY 626
# 04 November 2025

# Stan/Ulam with rethinking 
library(rethinking)	

# Set up data frame to work with ulam
Participant <-c(seq(1, 80), seq(1, 80))  # two scores for each participant

# Correct or not for item 5, then correct or not for item 7
Correct<- c(rep(1, 12), rep(0, 30), rep(1, 31), rep(0, 7), rep(0, 12), rep(0, 30), rep(1, 31), rep(1, 7))

# Condition 0 for item 5, 1 for item 7
Condition <- c(rep(0, 80), rep(1, 80))

data<- data.frame(Participant=Participant, Correct=Correct, Condition=Condition)

# Note logit(0.5)=0, and logistic(0)=0.5
# Model comparison of logit models (converts proportions to regression)
# Null model (items have the same proportion correct)
NullModel <- ulam(
			alist( Correct ~ dbinom(1, p), 
			logit(p) <- a[Participant],
			a[Participant] ~ dnorm(logit(0.5), 10)
			), data= data, log_lik=TRUE, chains=4  ) 
			

# Alt model (possibly different proportions for each item)
AltModel <- ulam(
			alist( Correct ~ dbinom(1, p), 
			logit(p) <- a[Participant] + b*Condition,
			a[Participant] ~ dnorm(logit(0.5), 10),
			b ~ dnorm(logit(0.5), 10)
			), data= data, log_lik=TRUE, chains=4  ) 
			
print(compare(NullModel, AltModel) )


# Analyzing posterior of alt model	

#plot posterior of b (change due to item being 7 rather than 5)
numSamples=2000
post<-extract.samples(AltModel, n= numSamples)

dev.new()
dens(post$b, col=rangi2, lwd=2, xlab="posterior for difference of logit proportion ")

# Probability that a > 0.5 ?
cat("Probability that b greater than 0")
length(post$b[post$b >0])/length(post$b)
cat("Probability that b less than 0")
length(post$b[post$b <0])/length(post$b)