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)

if (1==0){
# 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)){
		Lowscores<- rnorm(n, mean=mean(LowPerformingGroup), sd=sd(LowPerformingGroup))
		Highscores<- rnorm(n, mean=mean(HighPerformingGroup), sd=sd(HighPerformingGroup))
		
		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 ', EvidenceThreshold,' is ', powerEstimate, '\n')
	}
	
	if(powerEstimate>= DesiredPower){
		cat('Bayes Factor, smallest sample size with desired power is ', n, ' with power of ', powerEstimate, '\n')
		break
	}
}

}

if(1==1){

 
AweightThreshold = 0.9

# Stan/Ulam with model comparison 
library(rethinking)
numReps=20

for(n in seq(20, 1000)){

	CountSufficientEvidence=0
	# Generate simulated data sets
	for(rep in seq(1, numReps)){

		Lowscores<- rnorm(n, mean=mean(LowPerformingGroup), sd=sd(LowPerformingGroup))
		Highscores<- rnorm(n, mean=mean(HighPerformingGroup), sd=sd(HighPerformingGroup))
		
		# set up data frame
		data <- data.frame(Group=c(rep(0, n), rep(1,n)),  Score =c(Lowscores, Highscores) )		
		
		# Model comparison
		# 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  ) 
			
		
		WAICtable <-compare(NullModelWelch, AltModelWelch)
	
		# Identify which model wins the competition
		RowNames <- rownames(WAICtable)
		if(RowNames[1]=="AltModelWelch"){
			# Check Akaike weight
			thisWeight = WAICtable$weight[1]
			if(thisWeight > AweightThreshold){
				CountSufficientEvidence = CountSufficientEvidence +1		
			}
		}
	}
	powerEstimate = CountSufficientEvidence/numReps
	if(verbose){
		cat('With n=', n, ' WAIC power for threshold ', AweightThreshold,' is ', powerEstimate, '\n')

	}
	
	if(powerEstimate>= DesiredPower){
		cat('For model comparison, smallest sample size with desired  power is ', n, ' with power of ', powerEstimate, '\n')
		break
	}	
}



}