rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Greg Francis
# PSY 626
# 19 November 2025

# load data file
LWdata<-read.csv(file="LinkWord.csv",header=TRUE,stringsAsFactors=FALSE)

# Stan/Ulam with rethinking 

library(rethinking)

# Model comparison
# Null model (same mean for both conditions), may vary across participants
NullModel <- ulam(
	alist( Recalled ~ dnorm(mu, sigma), 
	mu <- a[Participant],
	a[Participant] ~ dnorm(GrandMean, GrandSigma),
	GrandMean ~ dnorm(12.5, 10),
	GrandSigma ~ dnorm(0, 10),
	sigma ~ dnorm(0, 10)
	), data= LWdata,  constraints=list(GrandSigma="lower=0", sigma="lower=0"), log_lik=TRUE, chains=4  ) 


# Alt model (allows for variation across conditions), may vary across participants
AltModel <- ulam(
	alist( Recalled ~ dnorm(mu, sigma), 
	mu <- a[Participant] + b*Condition,
	a[Participant] ~ dnorm(GrandMean, GrandSigma),
	GrandMean ~ dnorm(12.5, 10),
	GrandSigma ~ dnorm(0, 10),
	b ~ dnorm(0, 5),
	sigma ~ dnorm(0, 10)
	), data= LWdata,  constraints=list(GrandSigma="lower=0", sigma="lower=0"), log_lik=TRUE, chains=4  ) 

compare(NullModel, AltModel)


# Analyzing posterior

#plot posterior of a
numSamples=2000
post<-extract.samples(AltModel, n= numSamples)

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

# Convert b values to score improvements for a 100 word quiz
post$QuizImprove = 4*post$b

length(post$QuizImprove[post$QuizImprove >5])/length(post$QuizImprove)


length(post$QuizImprove[post$QuizImprove >2])/length(post$QuizImprove)


# Utilities

UtilityRecall = (post$b)^5

dens(UtilityRecall, col=rangi2, lwd=2, xlab="posterior for UtilityRecall")

UtilityDuration = -(0.5)^2

TotalUtility = UtilityRecall + UtilityDuration

mean(TotalUtility)

dev.new()
dens(TotalUtility, col=rangi2, lwd=2, xlab="posterior for TotalUtility")

# Probability of positive utility
length(TotalUtility[TotalUtility >0])/length(TotalUtility)

length(TotalUtility[TotalUtility >0])/length(TotalUtility)