rm(list=ls(all=TRUE))  # clear all variables
# Programmed by Greg Francis (gfrancis@purdue.edu) October, 2019

# Exploring excess linearity in Forster & Denzler (2012)

# define properties of experiments 
# Taken from Table 3 of investigative report
# Sample size (per cell)
n<-c(20, 20, 20, 20, 20, 20, 20, 20, 20, 15, 20, 15)
MeanLocal <-c(2.47, 2.51, 2.40, 2.41, 2.14, 3.19, 2.63, 2.87, 2.35, 2.55, 2.66, 2.42)
MeanControl<-c(3.04, 2.95, 2.90, 2.98, 2.82, 4.01, 3.73, 3.83, 3.66, 3.72, 3.69, 3.73)
MeanGlobal <-c(3.68, 3.35, 3.45, 3.64, 3.41, 4.79, 4.73, 4.79, 4.76, 4.78, 4.81, 5.02)

SDLocal<-c(1.21, 0.71, 0.86, 1.07, 1.20, 1.07, 1.49, 1.24, 1.01, 1.16, 1.21, 0.82)
SDControl <-c(0.72, 0.49, 0.51, 0.51, 0.78, 1.21, 1.21, 1.09, 1.19, 1.00, 1.30, 1.28)
SDGlobal<-c(0.68, 0.64, 0.80, 0.95, 0.71, 0.82, 1.55, 1.53, 1.71, 1.47, 1.54, 1.45)

numExperiments = 12

# Linearity of original data
# How different is control from the mean of the local and global?
# |((Global+Local)/2  – Control)/(Global-Local)| See http://datacolada.org/21#identifier_0_559
sum=0
cat("Linearity of original means\n")
for(studyIndex in c(1:numExperiments)){	
		diffPercent =  abs( ( (MeanGlobal[studyIndex] +MeanLocal[studyIndex])/2 - MeanControl[studyIndex] )/(MeanGlobal[studyIndex] - MeanLocal[studyIndex]) )
		cat(studyIndex, diffPercent, "\n")
		sum = sum + diffPercent
}
OriginalMeanDiff = sum/numExperiments

cat("Mean across all 12 experiments = ", OriginalMeanDiff)
cat("-------\n\n")


numRepeats = 5000
simulatedMeanDiffs <-c()

cat("Generating ", numRepeats, "simulated experiments...\n")


for(rep in c(1:numRepeats)){

	sum=0
	for(studyIndex in c(1:numExperiments)){
		
		# Generate random sample of data for each condition using the sample mean, standard deviation, and sample size
		xLocal <- rnorm(n=n[studyIndex], mean=MeanLocal[studyIndex], sd=SDLocal[studyIndex])
		xControl <- rnorm(n=n[studyIndex], mean=MeanControl[studyIndex], sd=SDControl[studyIndex])
		xGlobal <- rnorm(n=n[studyIndex], mean=MeanGlobal[studyIndex], sd=SDGlobal[studyIndex])


		diffPercent =  abs( ( (mean(xGlobal) +mean(xLocal))/2 - mean(xControl) )/(mean(xGlobal) - mean(xLocal)) )
#		cat(studyIndex, diffPercent, "\n")	
		sum = sum + diffPercent
	}
	simMeanDiff = sum/numExperiments
#	cat("Mean across all 12 experiments = ", simMeanDiff)
	
	simulatedMeanDiffs<- c(simulatedMeanDiffs, simMeanDiff)

	if(rep%%10000==0){cat("rep ", rep,"\n")}		
}


# Plot of simulated MeanDiffs
hist(simulatedMeanDiffs, xlim=c(0, 0.5), breaks= 2000)

# Compute how many simulated experiments have MeanDiffs less than observed
cat("There were ", length(simulatedMeanDiffs[simulatedMeanDiffs<OriginalMeanDiff]), "simulations with MeanDiff less than the observed value of ", OriginalMeanDiff, ".\n")
cat("The original data MeanDiff is at the ", 100*length(simulatedMeanDiffs[simulatedMeanDiffs<OriginalMeanDiff])/length(simulatedMeanDiffs), "percentile.\n")