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

# Greg Francis
# PSY 626
# 31 October 2025


# Bayes Factor 

library(BayesFactor)

LowPerformingGroup<- c(8, 6, 4, 12, 16, 17, 12, 10, 11, 13)
HighPerformingGroup<- c(23, 11, 17, 16, 6, 14, 15, 19)

bf = ttestBF(x = LowPerformingGroup, y = HighPerformingGroup)

bf

# Stan/Ulam with rethinking 

# set up data frame

data <- data.frame(Group=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1), Score =c(LowPerformingGroup, HighPerformingGroup) )		

library(rethinking)

# Model comparison
# Null model (same mean for each group)
NullModel <- ulam(
	alist( Score ~ dnorm(mu, sigma), 
	mu ~ dnorm(10, 10),
	sigma ~ dnorm(0, 10)
	), data= data, constraints=list(sigma="lower=0"), log_lik=TRUE, chains=4  ) 

# Alt model (different mean for each group, same standard deviation)
AltModel <- ulam(
	alist( Score ~ dnorm(mu, sigma), 
	mu <- a + b*Group,
	a ~ dnorm(10, 10),
	b ~ dnorm(0, 10),
	sigma ~ dnorm(0, 10)
	), data= data, constraints=list(sigma="lower=0"), log_lik=TRUE, chains=4  ) 

compare(NullModel, AltModel)


# Models with different variances

# Null model (same mean for each group)
NullModelWelch <- ulam(
	alist( Score ~ dnorm(mu, sigma), 
	mu ~ dnorm(10, 10),
	sigma <- d + e*Group,
	d ~ dnorm(0, 10),
	e ~ dnorm(0, 5)
	), data= data, constraints=list(d="lower=0"), log_lik=TRUE, chains=4  ) 
	
# Alt model (different mean for each group)
AltModelWelch <- ulam(
	alist( Score ~ dnorm(mu, sigma), 
	mu <- a + b*Group,
	a ~ dnorm(10, 10),
	b ~ dnorm(0, 10),
	sigma <- d + e*Group,
	d ~ dnorm(0, 10),
	e ~ dnorm(0, 5)
	), data= data, constraints=list(d="lower=0"), log_lik=TRUE, chains=4  )