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

verbose=TRUE
DesiredPower=0.9


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)){
		scores<- rnorm(n, mean=6.2, sd=1.2)
		
		# Subtract the value in the null hypothesis (create difference scores)
		diffScores = scores - 7.2
		
		test = ttestBF(x = diffScores)
		
		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)){

		scores<- rnorm(n, mean=6.2, sd=1.2)	
		# set up data frame
		data <- data.frame(Participant=seq(1,n),  Score =scores )		
		
		# Model comparison
		# Null model
		NullModel <- ulam(
			alist( Score ~ dnorm(7.2, sigma), 
			sigma ~ dnorm(0, 10)
			), data= data, constraints=list(sigma="lower=0"), log_lik=TRUE, chains=4  ) 
			
		
		# Alt model
		AltModel <- ulam(
			alist( Score ~ dnorm(mu, sigma), 
			mu ~ dnorm(0, 10),
			sigma ~ dnorm(0, 10)
			), data= data,  constraints=list(sigma="lower=0"), log_lik=TRUE, chains=4  ) 
			
		
		WAICtable <-compare(NullModel, AltModel)
	
		# Identify which model wins the competition
		RowNames <- rownames(WAICtable)
		if(RowNames[1]=="AltModel"){
			# 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){
#		sink()
		cat('For model comparison, smallest sample size with desired  power is ', n, ' with power of ', powerEstimate, '\n')
		break
	}	
}



}