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

# Greg Francis
# PSY 626
# 03 November 2025


# Bayes Factor 

library(BayesFactor)

data <- data.frame(Education=c(4, 6, 8, 11, 12, 14, 16, 17, 20), Income=c(6, 12, 14, 10, 17, 16, 13, 16, 19))

bf = correlationBF(data$Education, data$Income, nullInterval = c(-1, 0))

bf

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

# Stan/Ulam with rethinking 

# Linear regression models

library(rethinking)
	
# Null model (education has no effect on income)
NullModel <- ulam(
			alist( Income ~ dnorm(mu, sigma), 
			mu<- a,
			a ~ dnorm(15, 5),
			sigma~dnorm(0, 10)
			), data= data, constraints=list(sigma="lower=0"), log_lik=TRUE, chains=4  ) 
			

# Alt model (education has effect on income, and effect is positive)
AltModel <- ulam(
			alist( Income ~ dnorm(mu, sigma), 
			mu<- a + b*Education,
			a ~ dnorm(15, 5),
			b ~ dnorm(0, 5),
			sigma~dnorm(0, 10)
			), data= data, constraints=list(sigma="lower=0", b="lower=0"), log_lik=TRUE, chains=4  ) 
			
print(compare(NullModel, AltModel) )


# Analyzing posterior of full model
# Full model (education has effect on income, does not have to be positive)
FullModel <- ulam(
			alist( Income ~ dnorm(mu, sigma), 
			mu<- a + b*Education,
			a ~ dnorm(15, 5),
			b ~ dnorm(0, 5),
			sigma~dnorm(0, 10)
			), data= data, constraints=list(sigma="lower=0"), log_lik=TRUE, chains=4  ) 
	
numSamples=2000
post<-extract.samples(FullModel, n= numSamples)

# Probability that slope is positive ?
cat("Probability that slope is postive:")
length(post$b[post$b >0])/length(post$b)

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

# Can convert slope to correlation
Sx = sd(data$Education)
Sy= sd(data$Income)

# posterior for correlation calculated from posterior of slope
correlation <- post$b * Sy/Sx

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

# Probability that correlation is positve ?
cat("Probability that correlation is postive:")

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