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

verbose=TRUE
DesiredPower=0.9

# From original study
LowPerformingGroup<- c(8, 6, 4, 12, 16, 17, 12, 10, 11, 13)
HighPerformingGroup<- c(23, 11, 17, 16, 6, 14, 15, 19)
n1=length(LowPerformingGroup)
n2=length(HighPerformingGroup)
pooledSD = sqrt(((n1-1)*var(LowPerformingGroup) + (n2-1)*var(HighPerformingGroup))/(n1+n2-2))

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

# Estimate posterior of mean difference (low - high)
chains = posterior(test, iterations=10000)
diffposterior = chains[,2]	
				
# Bayes Factor 

library(BayesFactor)

EvidenceThreshold=3

for(n in seq(8,1000)){
	CountSufficientEvidence=0
	numReps=1000
	# Generate simulated data sets
	for(rep in seq(1, numReps)){
		# Pick a random difference of means from posterior
		diffMean= sample(diffposterior, 1)
		Lowscores<- rnorm(n, mean=mean(LowPerformingGroup), sd=pooledSD)
		Highscores<- rnorm(n, mean=mean(LowPerformingGroup)-diffMean, sd=pooledSD)
		
		test = ttestBF(x = Lowscores, y = Highscores)
		
		bf = extractBF(test, logbf=FALSE, onlybf=TRUE)
		
		if(bf>EvidenceThreshold){
			CountSufficientEvidence = CountSufficientEvidence +1
		}
	
	}
	powerEstimate = CountSufficientEvidence/numReps
	if(verbose){
		cat('With n=', n, ' Bayesian power for threshold based on posterior', EvidenceThreshold,' is ', powerEstimate, '\n')
	}
	
	if(powerEstimate>= DesiredPower){
		cat('Bayes Factor, smallest sample size with desired power based on posterior is ', n, ' with power of ', powerEstimate, '\n')
		break
	}
}