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

# Greg Francis
# PSY 626
# 03 November 2025


# Bayes Factor 

library(BayesFactor)

bf = proportionBF(y=135, N=255, p=.5, nullInterval = c(0, .5))

print(bf[2]/bf[1])

# Stan/Ulam with rethinking 
library(rethinking)

# set up data frame

# 1-indicates Republication, 0 - indicates Democrat
scores<- c(rep(1, 135), rep(0, 90))
data <- data.frame(Party =scores )		

# Note logit(0.5)=0, and logistic(0)=0.5
# Model comparison of logit models (converts proportions to regression)
NullModel <- ulam(
			alist( Party ~ dbinom(1, 0.5)

			), data= data, fixed_param=TRUE, log_lik=TRUE, chains=4  ) 			
			

# Alt model (broad prior bigger than null hypothesis of 0.5)
AltModel <- ulam(
			alist( Party ~ dbinom(1, p), 
			logit(p) <- a,
			a ~ dnorm(logit(0.5), 10)
			), data= data, constraints=list(a="lower=0"), log_lik=TRUE, chains=4  ) 
			
print(compare(NullModel, AltModel) )


# Analyzing posterior of full model

# Full model (broad prior but not restricted to be bigger than 0.5)
FullModel <- ulam(
			alist( Party ~ dbinom(1, p), 
			logit(p) <- a,
			a ~ dnorm(logit(0.5), 10)
			), data= data, log_lik=TRUE, chains=4  ) 	
	
#plot posterior of logistic(a)
numSamples=2000
post<-extract.samples(FullModel, n= numSamples)

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

# Probability that a > 0.5 ?
cat("Probability that a greater than 0 (Same as proportion greater than 0.5)")
length(post$a[post$a >0])/length(post$a)