#####################################################
# 			Bayesian A/B test on small sample sizes	#
#			Functions								#
#####################################################

#----------------
# Functions used in the main script

# Function to compute the analytical probability of a proportion a greater than another
# P(theta2 > theta1)
# Code from https://osf.io/anvg2
# this function is based on a formula published in
# Schmidt, M. N., & Morup, M. (2019). Efficient computation for Bayesian 
# comparison of two proportions. Statistics & Probability Letters, 145, 57-62.
# a1			[integer] : 			Count of successes in Group 1
# b1			[integer] : 			Count of failures in Group 1
# a2			[integer] : 			Count of successes in Group 2
# b2			[integer] : 			Count of failures in Group 2
prob.ab = function(a1, b1, a2 ,b2){
  
  if (a1*b2 > a2*b1){ 
    
    # normal series
    logz <- lgamma(a1 + a2) + lgamma(b1 + b2) - lgamma(a1 + b1 + a2 + b2 - 1) +
      log(genhypergeo(U = c(1, 1-a1, 1-b2), L = c(b1 + 1, a2 + 1), z = 1, 
                      tol = 1e-15,  series = T)) - log(b1*a2)
    
  } else {
    
    # flipped series
    logz <- lgamma(a1 + a2) + lgamma(b1 + b2) - lgamma(a1 + b1 + a2 + b2 - 1) +
      log(genhypergeo(U = c(1, 1-b1, 1-a2), L = c(a1 + 1, b2 + 1), z = 1, 
                      tol = 1e-15,  series = T)) - log(a1*b2)
  }
  
  # posterior probability of the event theta1 > theta2 
  return(exp(logz - lbeta(a1, b1) - lbeta(a2, b2)))
}



# Function using JAGS with the R2jags package to estimate two proportions and their difference with the IBE (Independent beta estimation) approach implemented by the MCMC algorithm
# data1				[list]: 				Information to estimate the first proportion
# data2				[list]: 				Information to estimate the second proportion
# sampleSize		[integer]:				Integer corresponding to the final chain sizes
# Inputs from data1 and data2
##	observations	[vector] : 				Binary vector with the observations (0 = success et 1 = failure)
##	alphaPrior		[numeric] : 			Prior value for the alpha parameter of the beta distribution
##	betaPrior		[numeric] : 			Prior value for the beta parameter of the beta distribution
# The output is a list containing 
# model1_jags								Output from the jags function applied to the first sample
# model2_jags								Output from the jags function applied to the second sample
# model1_mcmc								mcmc object version of model1_jags
# model2_mcmc								mcmc object version of model2_jags
# difference_mcmc							Differences of te estimates
bayesMCMC_2proportionsIBE = function(data1, data2, sampleSize = 1e5){
	# Parameters for the jags function
	nBurnin = 1000
	nThin = 10
	nIteration = sampleSize * nThin + nBurnin
	
	# JAGS Model
	bayes.mod = function() {
		for(i in 1:N){
		y[i] ~ dbern(pi1)
		}
		pi1 ~ dbeta(alpha, beta)
	}
	
	# Parameter of post distributions
	bayes.mod.params = c("pi1")
	
	# Initializing the processus by assigning random values to start with
	bayes.mod.inits = function(){
		list("pi1" = runif(1, min = 0, max = 1))
	}
	
	# Creation of the parameters to train the model
	y1 = data1$observations
	N1 = length(data1$observations)
	alpha1 = data1$alphaPrior
	beta1 = data1$betaPrior
	y2 = data2$observations
	N2 = length(data2$observations)
	alpha2 = data2$alphaPrior
	beta2 = data2$betaPrior
	
	sim.dat.jags1 = list(
		y = y1,
		N = N1,
		alpha = alpha1,
		beta = beta1
	)
	sim.dat.jags2 = list(
		y = y2,
		N = N2,
		alpha = alpha2,
		beta = beta2
	)
	
	# Adjust the models
	bayes.mod1.fit = jags(data = sim.dat.jags1, inits = bayes.mod.inits,
		parameters.to.save = bayes.mod.params, n.chains = 3, n.iter = nIteration,
		n.burnin = nBurnin, n.thin = nThin, model.file = bayes.mod)
	bayes.mod2.fit = jags(data = sim.dat.jags2, inits = bayes.mod.inits,
		parameters.to.save = bayes.mod.params, n.chains = 3, n.iter = nIteration,
		n.burnin = nBurnin, n.thin = nThin, model.file = bayes.mod)
	
	# Conversion into mcmc objects for a further checking of the output
	bayes.mod1.fit.mcmc = as.mcmc(bayes.mod1.fit)
	bayes.mod2.fit.mcmc = as.mcmc(bayes.mod2.fit)
	
	# Calculation of the differences
	deltas = list(
		bayes.mod2.fit.mcmc[[1]][,2] - bayes.mod1.fit.mcmc[[1]][,2],
		bayes.mod2.fit.mcmc[[2]][,2] - bayes.mod1.fit.mcmc[[2]][,2],
		bayes.mod2.fit.mcmc[[3]][,2] - bayes.mod1.fit.mcmc[[3]][,2]
	)

	bayes.diff.mcmc = bayes.mod1.fit.mcmc
	for(j in 1:length(bayes.diff.mcmc)){
		bayes.diff.mcmc[[j]][,2] = deltas[[j]]
	}
	
	return(list(
		model1_jags = bayes.mod1.fit, model2_jags = bayes.mod2.fit,
		model1_mcmc = bayes.mod1.fit.mcmc, model2_mcmc = bayes.mod2.fit.mcmc,
		difference_mcmc = bayes.diff.mcmc))
}



# Function to get the estimates of the prior densities
# donneesPrior		[list] : 				List containing the prior parameters
#												mu_gammaPrior : mean (mathematical expectation) of the gamma parameter (sometimes denoted by mu)
#												sigma_gammaPrior : standard deviation of the gamma parameter (sometimes denoted by mu)
#												mu_psiPrior : mean (mathematical expectation) of the psi parameter
#												sigma_psiPrior : standard deviation of the psi parameter
# The output is a list containing
# gamma										the 1e5 draws for the distribution of gamma
# psi										the 1e5 draws for the distribution of psi
# p1										the 1e5 draws for the distribution of pi1 computed from gamma and psi
# p2										the 1e5 draws for the distribution of pi2 computed from gamma and psi
# delta										the 1e5 draws for the distribution of delta = pi2 - pi1
priorEstimatesLTT = function(donneesPrior){
	parGamma = rnorm(1e5, mean = donneesPrior$mu_gammaPrior, sd = donneesPrior$sigma_gammaPrior)
	parPsi = rnorm(1e5, mean = donneesPrior$mu_psiPrior, sd = donneesPrior$sigma_psiPrior)
	
	aux1 = exp(parGamma - parPsi/2)
	aux2 = exp(parGamma + parPsi/2)
	resP1 = aux1 /(1 + aux1)
	resP2 = aux2 /(1 + aux2)
	resDelta = resP2 - resP1
	
	return(list(gamma = parGamma, psi = parPsi, p1 = resP1, p2 = resP2, delta = resDelta))
}



# Function that uses JAGS through the R2jags package fitted to estimate two proportions and their difference with the LTT (Log Transformation Testing) approach implemented with the MCMC algorithm
# dataInput			[list] : 				Information about the estimates of the two proportions
# dataInput content
##	data			[data.frame] : 			Sequential observations of the data
##												n1 : sequential count in Group 1
##												y1 : sequential count of successes in Group 1
##												n2 : sequential count in Group 2
##												y2 : sequential count of successes in Group 2
##	mu_gammaPrior		[numeric] : 		Prior value of the mu parameter of the gamma (or mu) distribution
##	sigma_gammaPrior	[numeric] : 		Prior value of the sigma parameter of the gamma (or mu) distribution
##	mu_psiPrior			[numeric] : 		Prior value of the mu parameter of the psi distribution
##	sigma_psiPrior		[numeric] : 		Prior value of the sigma parameter of the psi distribution
##	pH0					[numeric] : 		Prior probability value of H0 (between 0 and 1)
##	sampleSize			[integer] :			Large integer to obtain the desired number of draws at the output
# Return a list containing
# model_jags								Output of the jags function applied on the input data
# model_mcmc								mcmc object conversion from model_jags
bayesMCMC_2proportionsLTT = function(dataInput){
	# Parameters for the jags function
	nBurnin = 1000
	nThin = 10
	nIteration = dataInput$sampleSize * nThin + nBurnin

	# Model creation
	#-- Model with no constraints
	bayes.mod = function() {
		y1 ~ dbin(pi1, n1)
		y2 ~ dbin(pi2, n2)
		
		logit(pi1) <- mu - psi/2
		logit(pi2) <- mu + psi/2
		
		# In JAGS, the parameter of dnorm is the precision denoted by tau (not the standard deviance or variance)
		# sd = 1 / tau²
		tau_mu <- pow(sigma_mu, -2)
		tau_psi <- pow(sigma_psi, -2)
		mu ~ dnorm(mu_mu,tau_mu)
		psi ~ dnorm(mu_psi, tau_psi)
		
		# Parameter of interest
		delta <- pi2 - pi1
	}
	
	# Parameters of the post distributions
	bayes.mod.params = c("mu", "psi", "pi1", "pi2", "delta")

	# Initializing the processus by assigning random values to start with
	bayes.mod.inits = function(){
		list("mu" = runif(1, min = -1, max = 1), "psi" = runif(1, min = 0.1, max = 1.5))
	}
	
	# Creation of the parameters to train the model
	y1 = dataInput$data$y1
	y2 = dataInput$data$y2
	n1 = dataInput$data$n1
	n2 = dataInput$data$n2
	mu_mu = dataInput$mu_gammaPrior
	sigma_mu = dataInput$sigma_gammaPrior
	mu_psi = dataInput$mu_psiPrior
	sigma_psi = dataInput$sigma_psiPrior
		
	sim.dat.jags1 = list(
		y1 = y1,
		y2 = y2,
		n1 = n1,
		n2 = n2,
		mu_mu = mu_mu,
		sigma_mu = sigma_mu,
		mu_psi = mu_psi,
		sigma_psi = sigma_psi
	)
	
	# Adjust the models
	bayes.mod1.fit = jags(data = sim.dat.jags1, inits = bayes.mod.inits,
		parameters.to.save = bayes.mod.params, n.chains = 3, n.iter = nIteration,
		n.burnin = nBurnin, n.thin = nThin, model.file = bayes.mod)
	
	# Conversion into an mcmc object for further checking of the output
	bayes.mod1.fit.mcmc = as.mcmc(bayes.mod1.fit)
	
	return(list(model_jags = bayes.mod1.fit, model_mcmc = bayes.mod1.fit.mcmc))
}



# Function that computes the Bayes factors and the posterior probabilities of the H+ and H0 hypotheses
# dataInput			[list] : 				Information for the estimates of the two proportions
# dataInput content
##	samples_psi		[vector] : 				Réalizations of the psi parameter samples by the MCMC algorithm
##											It is the "psi" column of the model_jags$BUGSoutput$sims.matrix matrix from the JAGS process in the bayesMCMC_2proportionsLTT function
##	priorParameters	[list] : 				Prior parameters for the process
##												muPsi : 	Mean parameter for the psi density
##												sigmaPsi : 	Standard deviation for the psi density
##												p_H0 : 		Prior probability of H0 P(H0)
##												p_Hplus : 	Prior probability of H+ P(H+)
# Return a list containing
# BF_10										Bayes factor BF_10 (in favor of the bilateral alternative hypothesis over H0)
# BF_plus0									Bayes factor BF_+0 (in favor of the unilateral alternative hypothesis + over H0)
# BF_minus0									Bayes factor BF_-0 (in favor of the unilateral alternative hypothesis - over H0)
# p_Hplus_post								Post probability of H+ P(H+ | data)
# p_H0_post									Post probability of H0 P(H0 | data)
bayesMCMC_SD_post_BF = function(samples_psi, priorParameters = list(muPsi = 0, sigmaPsi = 1, p_H0 = 0.5, p_Hplus = 0.5)){
	# Density estimated by kernel
	dens_psi = density(samples_psi, adjust = 1)
	post_density_psi_at0_kde = approx(dens_psi$x, dens_psi$y, xout = 0)$y
	
	# Value of the prior density at 0
	prior_density_psi_at0 = dnorm(0, mean = priorParameters$muPsi, sd = priorParameters$sigmaPsi)

	# Bayes factor BF_01 (H0 vs H1) by Savage-Dickey
	BF_01_SD = post_density_psi_at0_kde / prior_density_psi_at0
	BF_10_SD = 1 / BF_01_SD

	# Prior probability of positivity: P(psi > 0) (depends mu_psi, sigma_psi)
	p_prior_psi_pos = pnorm(0, mean = priorParameters$muPsi, sd = priorParameters$sigmaPsi, lower.tail = FALSE)
	
	# Post probability of positivity: P(psi > 0 | data)
	p_post_psi_pos = length(which(samples_psi > 0)) / length(samples_psi)
	
	# Formula: BF_{+0} = BF_10_SD * (P(psi > 0 | data) / P(psi > 0))
	BF_plus0 = BF_10_SD * p_post_psi_pos / p_prior_psi_pos
	
	# Formula: BF_{-0} = BF_10_SD * (P(psi < 0 | data) / P(psi < 0))
	BF_minus0 = BF_10_SD * (1 - p_post_psi_pos) / (1 - p_prior_psi_pos)

	# Posterior probabilities for H+ and H0 given prior model probs (here 0.5/0.5)
	# Posterior odds H+ / H0 = prior odds * BF_{+0}; prior odds = p_Hplus_prior / p_H0_prior
	prior_odds_plus0 = priorParameters$p_Hplus / priorParameters$p_H0
	post_odds_plus0  = prior_odds_plus0 * BF_plus0

	p_Hplus_post = post_odds_plus0 / (1 + post_odds_plus0)
	p_H0_post = 1 / (1 + post_odds_plus0)
	names(p_H0_post) = "H0"

	result = list(BF_10 = BF_10_SD, BF_plus0 = BF_plus0, BF_minus0 = BF_minus0, 
		p_Hplus_post = p_Hplus_post, p_H0_post = p_H0_post)
	
	return(result)
}



# String conversion from a credible interval (CI)
# intervalle	[object] : 				R Object R from a CI calculation by the bayestestR::ci function
# precision		[integer] : 			Desired precision (number of digits)
# Return a string format of the CI: [CI_low; CI_high]
CI2string = function(intervalle, precision = 2) {
	result = paste0("[", round(intervalle$CI_low, digits = precision), "; ", round(intervalle$CI_high, digits = precision), "]")
	return(result)
}



# Function implementing our own simulation (i.e. the rbeta sampling) of the IBE approach, based on random samples of beta distributions (dbeta)
# contingencyData	[vector]: 			R vector containing the contingency table (i.e., count of successes and failures in each group)
# priorParameters	[vector]: 			R vector containing the parameters for the prior beta distributions
#										_ G1_alpha: alpha parameter of the beta distribution for Group 1
#										_ G1_beta: beta parameter of the beta distribution for Group 1
#										_ G2_alpha: alpha parameter of the beta distribution for Group 2
#										_ G2_beta: beta parameter of the beta distribution for Group 2
# sampleSize		[large integer]: 	Sample size during random draws
# ciLevel			[numeric]:			Numeric between 0 and 1 excluded to customize the credible interval level
# Outputs: The function provides the random draws and IBE indicators in a list
simulationIBE = function(contingencyData, priorParameters = c(1, 1, 1, 1), sampleSize = 1e5, ciLevel = 0.9) {
	# Count of successes and failures
	nbS_1 = contingencyData[2,1]
	nbS_2 = contingencyData[2,2]
	nbC_1 = contingencyData[1,1]
	nbC_2 = contingencyData[1,2]
	
	# Distribution parameters
	alphaG1 = priorParameters[1]
	betaG1 = priorParameters[2]
	alphaG2 = priorParameters[3]
	betaG2 = priorParameters[4]
	
	# Simulation of prior distributions
	simPrior = data.frame(
		probP1 = rbeta(sampleSize, alphaG1, betaG1),
		probP2 = rbeta(sampleSize, alphaG2, betaG2)
	)

	simPrior = simPrior %>%
		mutate(
			difference = probP2 - probP1,
			differenceR = (probP2 - probP1) / probP1,
		)
	
	# Simulation of posterior distributions
	simPosterior = data.frame(
		probP1 = rbeta(sampleSize, alphaG1 + nbS_1, betaG1 + nbC_1),
		probP2 = rbeta(sampleSize, alphaG2 + nbS_2, betaG2 + nbC_2)
	)

	simPosterior = simPosterior %>%
		mutate(
			difference = probP2 - probP1,
			differenceR = (probP2 - probP1) / probP1,
		)

	# Quantiles
	priorQuantilesG1 = quantile(simPrior$probP1)
	priorQuantilesG2 = quantile(simPrior$probP2)
	priorQuantilesDelta = quantile(simPrior$difference)
	postQuantilesG1 = quantile(simPosterior$probP1)
	postQuantilesG2 = quantile(simPosterior$probP2)
	postQuantilesDelta = quantile(simPosterior$difference)
	
	# Credible interval by HDI and ETI
	SimHDI_Pi1 = CI2string(bayestestR::ci(simPosterior$probP1, ci = ciLevel, method = "HDI"))
	SimHDI_Pi2 = CI2string(bayestestR::ci(simPosterior$probP2, ci = ciLevel, method = "HDI"))
	SimHDI_delta = CI2string(bayestestR::ci(simPosterior$difference, ci = ciLevel, method = "HDI"))

	SimETI_Pi1 = CI2string(bayestestR::ci(simPosterior$probP1, ci = ciLevel, method = "ETI"))
	SimETI_Pi2 = CI2string(bayestestR::ci(simPosterior$probP2, ci = ciLevel, method = "ETI"))
	SimETI_delta = CI2string(bayestestR::ci(simPosterior$difference, ci = ciLevel, method = "ETI"))
	
	# Prior probabilities, posterior probabilities and Bayes factor computation
	SimPriorProb_12 = length(which(simPrior$difference > 0)) / length(simPrior$difference)
	SimPosteriorProb_12 = length(which(simPosterior$difference > 0)) / length(simPosterior$difference)
	SimBF = SimPosteriorProb_12 / (1 - SimPosteriorProb_12) * (1 - SimPriorProb_12) / SimPriorProb_12
	
	# Return the results
	result = list(
		simulationPrior = data.frame(
			simulation_Pi1 = simPrior$probP1,
			simulation_Pi2 = simPrior$probP2,
			simulation_delta = simPrior$difference
		),
		simulationPosterior = data.frame(
			simulation_Pi1 = simPosterior$probP1,
			simulation_Pi2 = simPosterior$probP2,
			simulation_delta = simPosterior$difference
		),
		priorQuantiles_Pi1 = priorQuantilesG1,
		priorQuantiles_Pi2 = priorQuantilesG2,
		priorQuantiles_delta = priorQuantilesDelta,
		postQuantiles_Pi1 = postQuantilesG1,
		postQuantiles_Pi2 = postQuantilesG2,
		postQuantiles_delta = postQuantilesDelta,
		HDI_Pi1 = SimHDI_Pi1,
		HDI_Pi2 = SimHDI_Pi2,
		HDI_delta = SimHDI_delta,
		ETI_Pi1 = SimETI_Pi1,
		ETI_Pi2 = SimETI_Pi2,
		ETI_delta = SimETI_delta,
		priorProbabilities = SimPriorProb_12,
		posteriorProbabilities = SimPosteriorProb_12,
		bayesFactor = SimBF
	)
	
	return(result)
}



# Function implementing the bayesAB function and then summarizing the results into indicators (posterior probabilities, HDI, ETI, Bayes factors, etc.)
# dataG1			[vector]: 			R vector containing the binary outputs (0 for failures and 1 for sucesses) of Group 1
# dataG2			[vector]: 			R vector containing the binary outputs (0 for failures and 1 for sucesses) of Group 2
# priorParameters	[vector]: 			R vector containing the parameters for the prior beta distributions
#										_ alpha: alpha parameter of the beta distributions
#										_ beta: beta parameter of the beta distributions
# sampleSize		[large integer]: 	sample size during random draws
# ciLevel			[numeric]:			Numeric between 0 and 1 excluded to customize the credible interval level
# Outputs: The function provides the random draws and IBE indicators in a list
bayesABSummarize = function(dataG1, dataG2, priorParameters = c(1, 1), sampleSize = 1e5, ciLevel = 0.9) {
	# Distribution parameters
	alphaPar = priorParameters[1]
	betaPar = priorParameters[2]

	simulation = bayesTest(dataG2, dataG1,
		priors = c('alpha' = alphaPar, 'beta' = betaPar), 
		n_samples = sampleSize, distribution = 'bernoulli')

	# Quantiles
	resSimulation = summary(simulation)
	differenceDelta = simulation$posteriors$Probability$A - simulation$posteriors$Probability$B
	quantilesG1 = quantile(simulation$posteriors$Probability$B)
	quantilesG2 = quantile(simulation$posteriors$Probability$A)
	quantilesDelta = quantile(differenceDelta)
	
	# Posterior probabilities and Bayes factor computation
	posteriorProb_12 = resSimulation$probability$Probability
	BF = resSimulation$probability$Probability / (1 - resSimulation$probability$Probability)
	# In this case, the prior distributions for Pi1 and Pi2 are equal => P(Pi1 > Pi2) = P(Pi1 < Pi2) = 0.5
	
	# Credible interval by HDI and ETI
	HDI_Pi1 = CI2string(bayestestR::ci(simulation$posteriors$Probability$B, ci = ciLevel, method = "HDI"))
	HDI_Pi2 = CI2string(bayestestR::ci(simulation$posteriors$Probability$A, ci = ciLevel, method = "HDI"))
	HDI_delta = CI2string(bayestestR::ci(differenceDelta, ci = ciLevel, method = "HDI"))

	ETI_Pi1 = CI2string(bayestestR::ci(simulation$posteriors$Probability$B, ci = ciLevel, method = "ETI"))
	ETI_Pi2 = CI2string(bayestestR::ci(simulation$posteriors$Probability$A, ci = ciLevel, method = "ETI"))
	ETI_delta = CI2string(bayestestR::ci(differenceDelta, ci = ciLevel, method = "ETI"))
	
	# Return the results
	result = list(
		simulationPosterior = data.frame(
			simulation_Pi1 = simulation$posteriors$Probability$B,
			simulation_Pi2 = simulation$posteriors$Probability$A,
			simulation_delta = differenceDelta
		),
		postQuantiles_Pi1 = quantilesG1,
		postQuantiles_Pi2 = quantilesG2,
		postQuantiles_delta = quantilesDelta,
		HDI_Pi1 = HDI_Pi1,
		HDI_Pi2 = HDI_Pi2,
		HDI_delta = HDI_delta,
		ETI_Pi1 = ETI_Pi1,
		ETI_Pi2 = ETI_Pi2,
		ETI_delta = ETI_delta,
		posteriorProbabilities = posteriorProb_12,
		bayesFactor = BF
	)
	
	return(result)
}



# Function that puts the indicators from objects implementing the IBE approach into a table
# simulation		[object]: 			Output from the simulationIBE function containing indicators of the IBE approach implemented by our own simulation (i.e. the rbeta sampling)
# bayesAB			[object]: 			Output from the bayesABSummarize function containing indicators of the IBE approach implemented by the bayesTest function
# mcmc				[Object]:			Output from the mcmcSummarize function containing indicators of the OBE approach implemented into the jags function by the MCMC algorithm
# analytic			[numeric]: 			Output from the prob.ab computing the post probability of Pi2 > Pi1
# Return a talbe with the indicators
summarize_IBE = function(simulation, bayesAB, mcmc, analytic){
	nrowsRes = 5 + 1 + (5 + 4) * 3 + 1			# 5 prior quantiles + 1 prior probabilities + (5 posterior quantiles + posterior probability + BF + 90% HDI + 90% ETI) pour simulation, bayesAB et MCMC + 1 posterior probability par analytique
	ncolsRes = 3								# G1, G2 and delta
	rowNamesRes = c(paste0("Prior quantiles ", 0:4), "Sim prior Prob",
		paste0("Sim Posterior quantiles ", 0:4),
		"Sim HDI", "Sim ETI", "Sim posterior Prob", "Sim BF",
		paste0("bayesAB Posterior quantiles ", 0:4), 
		"bayesAB HDI", "bayesAB ETI", "bayesAB posterior Prob", "bayesAB BF",
		paste0("MCMC Posterior quantiles ", 0:4), 
		"MCMC HDI", "MCMC ETI", "MCMC posterior Prob", "MCMC BF", "Analytic posterior probability")
	colNamesRes = c("G1", "G2", "Delta")
	
	tableResIBE = array(c(NA), dim = c(nrowsRes, ncolsRes))
	rownames(tableResIBE) = rowNamesRes
	colnames(tableResIBE) = colNamesRes
	
	tableResIBE[1:5,1] = simulation$priorQuantiles_Pi1
	tableResIBE[1:5,2] = simulation$priorQuantiles_Pi2
	tableResIBE[1:5,3] = simulation$priorQuantiles_delta
	tableResIBE[6,1] = simulation$priorProbabilities
	tableResIBE[7:11,1] = simulation$postQuantiles_Pi1
	tableResIBE[7:11,2] = simulation$postQuantiles_Pi2
	tableResIBE[7:11,3] = simulation$postQuantiles_delta
	tableResIBE[12,1:3] = c(simulation$HDI_Pi1, simulation$HDI_Pi2, simulation$HDI_delta)
	tableResIBE[13,1:3] = c(simulation$ETI_Pi1, simulation$ETI_Pi2, simulation$ETI_delta)
	tableResIBE[14,1] = simulation$posteriorProbabilities
	tableResIBE[15,1] = simulation$bayesFactor
	
	if(!is.null(bayesAB)){
		tableResIBE[16:20,1] = bayesAB$postQuantiles_Pi1
		tableResIBE[16:20,2] = bayesAB$postQuantiles_Pi2
		tableResIBE[16:20,3] = bayesAB$postQuantiles_delta
		tableResIBE[21,1:3] = c(bayesAB$HDI_Pi1, bayesAB$HDI_Pi2, bayesAB$HDI_delta)
		tableResIBE[22,1:3] = c(bayesAB$ETI_Pi1, bayesAB$ETI_Pi2, bayesAB$ETI_delta)
		tableResIBE[23,1] = bayesAB$posteriorProbabilities
		tableResIBE[24,1] = bayesAB$bayesFactor
	}
	
	tableResIBE[25:29,1] = mcmc$postQuantiles_Pi1
	tableResIBE[25:29,2] = mcmc$postQuantiles_Pi2
	tableResIBE[25:29,3] = mcmc$postQuantiles_delta
	tableResIBE[30,1:3] = c(mcmc$HDI_Pi1, mcmc$HDI_Pi2, mcmc$HDI_delta)
	tableResIBE[31,1:3] = c(mcmc$ETI_Pi1, mcmc$ETI_Pi2, mcmc$ETI_delta)
	tableResIBE[32,1] = mcmc$posteriorProbabilities
	tableResIBE[33,1] = mcmc$bayesFactor
	tableResIBE[34,1] = analytic
	
	return(tableResIBE)
}



# Function implementing the bayesAB function and then summarizing the results into indicators (posterior probabilities, HDI, ETI, Bayes factors, etc.)
# dataG1			[vector]: 			R vector containing the binary outputs (0 for failures and 1 for sucesses) of Group 1
# dataG2			[vector]: 			R vector containing the binary outputs (0 for failures and 1 for sucesses) of Group 2
# priorParameters	[vector]: 			R vector containing the parameters for the prior beta distributions
#										_ G1_alpha: alpha parameter of the beta distribution for Group 1
#										_ G1_beta: beta parameter of the beta distribution for Group 1
#										_ G2_alpha: alpha parameter of the beta distribution for Group 2
#										_ G2_beta: beta parameter of the beta distribution for Group 2
# sampleSize		[large integer]: 	sample size during random draws
# ciLevel			[numeric]:			Numeric between 0 and 1 excluded to customize the credible interval level
# Outputs: The function provides the random draws and IBE indicators in a list
mcmcSummarize = function(dataG1, dataG2, priorParameters = c(1, 1, 1, 1), sampleSize = 1e5, ciLevel = 0.9){
	# Distribution parameters
	alphaG1 = priorParameters[1]
	betaG1 = priorParameters[2]
	alphaG2 = priorParameters[3]
	betaG2 = priorParameters[4]
	
	# Simulation of the prior distributions
	simPrior = data.frame(
		probP1 = rbeta(1e5, alphaG1, betaG1),
		probP2 = rbeta(1e5, alphaG2, betaG2)
	)

	simPrior = simPrior %>%
		mutate(
			difference = probP2 - probP1,
			differenceR = (probP2 - probP1) / probP1,
		)

	# Implementation of the MCMC algorithm
	donnees1 = list(observations = dataG1, alphaPrior = alphaG1, betaPrior = betaG1)
	donnees2 = list(observations = dataG2, alphaPrior = alphaG2, betaPrior = betaG2)
	AB_MCMC = bayesMCMC_2proportionsIBE(donnees1, donnees2, sampleSize)

	# Aggregation of the chains
	allChains_Pi1 = c(AB_MCMC$model1_mcmc[[1]][,2], AB_MCMC$model1_mcmc[[2]][,2], AB_MCMC$model1_mcmc[[3]][,2])
	allChains_Pi2 = c(AB_MCMC$model2_mcmc[[1]][,2], AB_MCMC$model2_mcmc[[2]][,2], AB_MCMC$model2_mcmc[[3]][,2])
	allChains_delta = c(AB_MCMC$difference_mcmc[[1]][,2], AB_MCMC$difference_mcmc[[2]][,2], AB_MCMC$difference_mcmc[[3]][,2])
	
	# Quantiles
	quantilesMCMCG1 = quantile(allChains_Pi1)
	quantilesMCMCG2 = quantile(allChains_Pi2)
	quantilesMCMCDelta = quantile(allChains_delta)

	# Computation of the posterior probabilities for each MCMC chain
	prob1MCMC_vect = c()
	for(j in 1:length(AB_MCMC$difference_mcmc)){
		prob1MCMC_vect[j] = length(which(AB_MCMC$difference_mcmc[[j]][,2] > 0)) / length(AB_MCMC$difference_mcmc[[j]][,2])
	}

	# Prior and posterior probabilities and Bayes factor computation
	SimPriorProb_12 = length(which(simPrior$difference > 0)) / length(simPrior$difference)
	prob_MCMC = mean(prob1MCMC_vect)
	BF_MCMC = prob_MCMC / (1 - prob_MCMC) * (1 - SimPriorProb_12) / SimPriorProb_12

	# Credible interval by HDI and ETI
	HDImcmc_Pi1 = CI2string(bayestestR::ci(allChains_Pi1, ci = ciLevel, method = "HDI"))
	HDImcmc_Pi2 = CI2string(bayestestR::ci(allChains_Pi2, ci = ciLevel, method = "HDI"))
	HDImcmc_delta = CI2string(bayestestR::ci(allChains_delta, ci = ciLevel, method = "HDI"))

	ETImcmc_Pi1 = CI2string(bayestestR::ci(allChains_Pi1, ci = ciLevel, method = "ETI"))
	ETImcmc_Pi2 = CI2string(bayestestR::ci(allChains_Pi2, ci = ciLevel, method = "ETI"))
	ETImcmc_delta = CI2string(bayestestR::ci(allChains_delta, ci = ciLevel, method = "ETI"))
	
	# Return the results
	result = list(
		simulationPosterior = data.frame(
			simulation_Pi1 = allChains_Pi1,
			simulation_Pi2 = allChains_Pi2,
			simulation_delta = allChains_delta
		),
		postQuantiles_Pi1 = quantilesMCMCG1,
		postQuantiles_Pi2 = quantilesMCMCG2,
		postQuantiles_delta = quantilesMCMCDelta,
		HDI_Pi1 = HDImcmc_Pi1,
		HDI_Pi2 = HDImcmc_Pi2,
		HDI_delta = HDImcmc_delta,
		ETI_Pi1 = ETImcmc_Pi1,
		ETI_Pi2 = ETImcmc_Pi2,
		ETI_delta = ETImcmc_delta,
		posteriorProbabilities = prob_MCMC,
		bayesFactor = BF_MCMC
	)
	
	return(result)
}



# Function that summarizes and puts the quantiles from priorEstimatesLTT function outputs
# res_prior			[object]					Object containing outputs from priorEstimatesLTT function
summarize_prior_LTT = function(res_prior){
	
	result = cbind(quantile(res_prior$gamma),
		quantile(res_prior$psi),
		quantile(res_prior$p1),
		quantile(res_prior$p2),
		quantile(res_prior$delta))
	
	rownames(result) = paste0("Prior quantiles ", 0:4)
	colnames(result) = c("Gamma", "Psi", "G1", "G2", "Delta")
	
	return(result)
}


# Function that summarizes and puts the indicators from ab_test function outputs
# res_ab_test			[object]				Object containing outputs from ab_test function
# ciLevel				[numeric]				Numeric between 0 and 1 excluded to customize the credible interval level
# Return a table containing the indicators
summarize_ab_test = function(res_ab_test, ciLevel = 0.9){
	# New object to sum up res_ab_test
	resSumm_ab_test = summary(res_ab_test)
	# Determine the hypotheses to compare
	hypotheses = names(which(resSumm_ab_test$input$prior_prob > 0))
	
	# Build the result table
	nrowsRes = 2 + 5 + 5		# 2 for prior probabilities, 5 quantiles and 5 indicators (90% HDI et ETI, post probabilities of H+ and H0, Bayes factor) for ab_test in order to compare H+ with H0, then H+ and H-
	
	ncolsRes = 5				# Gamma, Psi, G1, G2, delta
	rowNamesRes = c(paste0("Prior prob ", hypotheses[1]), paste0("Prior prob ", hypotheses[2]),
		paste0("ab_test Post quantiles ", 0:4), 
		"ab_test HDI", "ab_test ETI", paste0("ab_test posterior Prob ", hypotheses[1]), paste0("ab_test posterior Prob ", hypotheses[2]), "ab_test BF")
	colNamesRes = c("Gamma", "Psi", "G1", "G2", "Delta")
	result = array(c(NA), dim = c(nrowsRes, ncolsRes))
	rownames(result) = rowNamesRes
	colnames(result) = colNamesRes
	
	
	# Estimates of delta
	ab_test_difference = resSumm_ab_test$post$Hplus$p2 - resSumm_ab_test$post$Hplus$p1
	
	# Credible interval by HDI and ETI
	ab_testHDI_Gamma = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$beta, ci = ciLevel, method = "HDI"))
	ab_testHDI_Psi = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$psi, ci = ciLevel, method = "HDI"))
	ab_testHDI_Pi1 = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$p1, ci = ciLevel, method = "HDI"))
	ab_testHDI_Pi2 = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$p2, ci = ciLevel, method = "HDI"))
	ab_testHDI_delta = CI2string(bayestestR::ci(ab_test_difference, ci = ciLevel, method = "HDI"))
	ab_testHDI_OR = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$or, ci = ciLevel, method = "HDI"))
	ab_testHDI_Arisk = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$arisk, ci = ciLevel, method = "HDI"))
	ab_testHDI_Rrisk = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$rrisk, ci = ciLevel, method = "HDI"))

	ab_testETI_Gamma = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$beta, ci = ciLevel, method = "ETI"))
	ab_testETI_Psi = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$psi, ci = ciLevel, method = "ETI"))
	ab_testETI_Pi1 = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$p1, ci = ciLevel, method = "ETI"))
	ab_testETI_Pi2 = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$p2, ci = ciLevel, method = "ETI"))
	ab_testETI_delta = CI2string(bayestestR::ci(ab_test_difference, ci = ciLevel, method = "ETI"))
	ab_testETI_OR = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$or, ci = ciLevel, method = "ETI"))
	ab_testETI_Arisk = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$arisk, ci = ciLevel, method = "ETI"))
	ab_testETI_Rrisk = CI2string(bayestestR::ci(resSumm_ab_test$post$Hplus$rrisk, ci = ciLevel, method = "ETI"))

	# Save the results
	result[1:2,1] = c(resSumm_ab_test$input$prior_prob[hypotheses[1]], resSumm_ab_test$input$prior_prob[hypotheses[2]])
	result[3:7,1] = c(quantile(resSumm_ab_test$post$Hplus$beta))
	result[3:7,2] = c(quantile(resSumm_ab_test$post$Hplus$psi))
	result[3:7,3] = c(quantile(resSumm_ab_test$post$Hplus$p1))
	result[3:7,4] = c(quantile(resSumm_ab_test$post$Hplus$p2))
	result[3:7,5] = c(quantile(ab_test_difference))
	result[8,] = c(ab_testHDI_Gamma, ab_testHDI_Psi, ab_testHDI_Pi1, ab_testHDI_Pi2, ab_testHDI_delta)
	result[9,] = c(ab_testETI_Gamma, ab_testETI_Psi, ab_testETI_Pi1, ab_testETI_Pi2, ab_testETI_delta)
	result[10:11,1] = c(resSumm_ab_test$post_prob[hypotheses[1]], resSumm_ab_test$post_prob[hypotheses[2]])
	result[12,1] = resSumm_ab_test$bf$bfplus0
	
	return(result)
}



# Function that summarizes and puts the indicators from ab_test function outputs
# res_MCMC				[object]				Object containing outputs from bayesMCMC_2proportionsLTT function
# priorPsi				[vector]				Vector containing the prior settings for psi
#												_ psiMu for the density center
#												_ psiSigma for the standard deviation
# priorProbabilities	[vector]				Vector containing the prior probaility settings
#												_ H+
#												_ H0
# ciLevel				[numeric]				Numeric between 0 and 1 excluded to customize the credible interval level
# Return a table containing the indicators
summarize_MCMC_LTT = function(res_MCMC, priorPsi = c(psiMu = 0, psiSigma = 1), priorProbabilities = c("H+" = 0.5, "H0" = 0.5), ciLevel = 0.9){

	# Build the result table
	nrowsRes = 5 + 6			# 5 quantiles and 5 quantiles and 6 indicators (90% HDI et ETI, post probabilities of H+ and H0, Bayes factor) for MCMC in order to compare H+ with H0 and BF-0
	
	ncolsRes = 5				# Gamma, Psi, G1, G2, delta
	rowNamesRes = c(paste0("MCMC Post quantiles ", 0:4), 
		"MCMC HDI", "MCMC ETI", "MCMC posterior Prob H+", "MCMC posterior Prob H0", "MCMC BF+0", "MCMC BF-0")
	colNamesRes = c("Gamma", "Psi", "G1", "G2", "Delta")
	result = array(c(NA), dim = c(nrowsRes, ncolsRes))
	rownames(result) = rowNamesRes
	colnames(result) = colNamesRes
	
	
	# Quantiles and other indicators to save into the result table
	allChains_Gamma = c(res_MCMC$model_mcmc[[1]][,'mu'], res_MCMC$model_mcmc[[2]][,'mu'], res_MCMC$model_mcmc[[3]][,'mu'])
	allChains_Psi = c(res_MCMC$model_mcmc[[1]][,'psi'], res_MCMC$model_mcmc[[2]][,'psi'], res_MCMC$model_mcmc[[3]][,'psi'])
	allChains_Pi1 = c(res_MCMC$model_mcmc[[1]][,'pi1'], res_MCMC$model_mcmc[[2]][,'pi1'], res_MCMC$model_mcmc[[3]][,'pi1'])
	allChains_Pi2 = c(res_MCMC$model_mcmc[[1]][,'pi2'], res_MCMC$model_mcmc[[2]][,'pi2'], res_MCMC$model_mcmc[[3]][,'pi2'])
	allChains_delta = c(res_MCMC$model_mcmc[[1]][,'delta'], res_MCMC$model_mcmc[[2]][,'delta'], res_MCMC$model_mcmc[[3]][,'delta'])
	allChains_or = (allChains_Pi2/(1-allChains_Pi2)) / (allChains_Pi1/(1-allChains_Pi1))
	allChains_ar = allChains_Pi2 - allChains_Pi1
	allChains_rr = allChains_Pi2 / allChains_Pi1

	MCMC_HDI_Gamma = CI2string(bayestestR::ci(allChains_Gamma, ci = ciLevel, method = "HDI"))
	MCMC_HDI_Psi = CI2string(bayestestR::ci(allChains_Psi, ci = ciLevel, method = "HDI"))
	MCMC_HDI_Pi1 = CI2string(bayestestR::ci(allChains_Pi1, ci = ciLevel, method = "HDI"))
	MCMC_HDI_Pi2 = CI2string(bayestestR::ci(allChains_Pi2, ci = ciLevel, method = "HDI"))
	MCMC_HDI_delta = CI2string(bayestestR::ci(allChains_delta, ci = ciLevel, method = "HDI"))
	MCMC_HDI_OR = CI2string(bayestestR::ci(allChains_or, ci = ciLevel, method = "HDI"))
	MCMC_HDI_Arisk = CI2string(bayestestR::ci(allChains_ar, ci = ciLevel, method = "HDI"))
	MCMC_HDI_Rrisk = CI2string(bayestestR::ci(allChains_rr, ci = ciLevel, method = "HDI"))

	MCMC_ETI_Gamma = CI2string(bayestestR::ci(allChains_Gamma, ci = ciLevel, method = "ETI"))
	MCMC_ETI_Psi = CI2string(bayestestR::ci(allChains_Psi, ci = ciLevel, method = "ETI"))
	MCMC_ETI_Pi1 = CI2string(bayestestR::ci(allChains_Pi1, ci = ciLevel, method = "ETI"))
	MCMC_ETI_Pi2 = CI2string(bayestestR::ci(allChains_Pi2, ci = ciLevel, method = "ETI"))
	MCMC_ETI_delta = CI2string(bayestestR::ci(allChains_delta, ci = ciLevel, method = "ETI"))
	MCMC_ETI_OR = CI2string(bayestestR::ci(allChains_or, ci = ciLevel, method = "ETI"))
	MCMC_ETI_Arisk = CI2string(bayestestR::ci(allChains_ar, ci = ciLevel, method = "ETI"))
	MCMC_ETI_Rrisk = CI2string(bayestestR::ci(allChains_rr, ci = ciLevel, method = "ETI"))
	
	# Computation of post probabilities and Bayes factors with Savage-Dickey on MCMC outputs
	mcmc_samples = as.matrix(res_MCMC$model_jags$BUGSoutput$sims.matrix)

	postBySD = bayesMCMC_SD_post_BF(mcmc_samples[, "psi"], 
		priorParameters = list(muPsi = priorPsi["psiMu"], sigmaPsi = priorPsi["psiSigma"], p_H0 = priorProbabilities["H0"], p_Hplus = priorProbabilities["H+"]))

	
	# Save the results
	result[1:5,1] = c(quantile(allChains_Gamma))
	result[1:5,2] = c(quantile(allChains_Psi))
	result[1:5,3] = c(quantile(allChains_Pi1))
	result[1:5,4] = c(quantile(allChains_Pi2))
	result[1:5,5] = c(quantile(allChains_delta))
	result[6,] = c(MCMC_HDI_Gamma, MCMC_HDI_Psi, MCMC_HDI_Pi1, MCMC_HDI_Pi2, MCMC_HDI_delta)
	result[7,] = c(MCMC_ETI_Gamma, MCMC_ETI_Psi, MCMC_ETI_Pi1, MCMC_ETI_Pi2, MCMC_ETI_delta)
	
	result[8:9,1] = c(postBySD$p_Hplus_post, postBySD$p_H0_post)
	result[10,1] = postBySD$BF_plus0
	result[11,1] = postBySD$BF_minus0

	return(result)
}



# Function that summarizes and puts the indicators from objects implementing the LTT approach
# res_prior				[object]				Object containing outputs from priorEstimatesLTT function
# res_ab_test			[object]				Object containing outputs from ab_test function
# res_MCMC				[object]				Object containing outputs from bayesMCMC_2proportionsLTT function
# priorPsi				[vector]				Vector containing the parameters of the prior Psi
#												_ psiMu for the density center
#												_ psiSigma for the standard deviation
# priorProbabilities	[vector]				
# Return a table containing the indicators
summarize_LTT = function(res_prior, res_ab_test, res_MCMC, priorPsi = c(psiMu = 0, psiSigma = 1), priorProbabilities = c("H+" = 0.5, "H0" = 0.5), ciLevel = 0.9){
	
	# Build the result table
	nrowsRes = 2 + 5 + (5 + 5) * 2 + 1			# 5 for prior quantiles 2 for prior probabilities, 5 quantiles and 5 indicators (90% HDI et ETI, post probabilities of H+ and H0, Bayes factor) for ab_test et MCMC, all in order to compare H+ with H0 and one BF-0 for MCMC
												# we don't store the MCMC samples when comparing H- and H+ because they are similar to those from the comparison between H0 and H+ (so -7)
	ncolsRes = 5								# Gamma, Psi, G1, G2, delta
	rowNamesRes = c("Prior prob H+", "Prior prob H0",
		paste0("Prior quantiles ", 0:4), 
		paste0("ab_test Post quantiles ", 0:4), 
		"ab_test HDI", "ab_test ETI", "ab_test posterior Prob H+", "ab_test posterior Prob H0", "ab_test BF+0",
		paste0("MCMC Post quantiles ", 0:4), 
		"MCMC HDI", "MCMC ETI", "MCMC posterior Prob H+", "MCMC posterior Prob H0", "MCMC BF+0", "MCMC BF-0")
	colNamesRes = c("Gamma", "Psi", "G1", "G2", "Delta")
	result = array(c(NA), dim = c(nrowsRes, ncolsRes))
	rownames(result) = rowNamesRes
	colnames(result) = colNamesRes
	
	# Prior simulations
	prior_table = summarize_prior_LTT(res_prior)
	
	# ab_test outputs
	ab_test_table = summarize_ab_test(res_ab_test, ciLevel = 0.9)
	
	# MCMC outputs
	mcmc_table = summarize_MCMC_LTT(res_MCMC, priorPsi = priorPsi, priorProbabilities = priorProbabilities, ciLevel = 0.9)
	
	result[1:2,1] = ab_test_table[1:2,1]
	result[3:7,] = prior_table
	result[8:12,] = ab_test_table[paste0("ab_test Post quantiles ", 0:4),]
	result[13:14,] = ab_test_table[c("ab_test HDI", "ab_test ETI"),]
	result[15:17,1] = ab_test_table[10:12,1]
	rownames(result)[1:2] = rownames(ab_test_table)[1:2]
	rownames(result)[3:7] = rownames(prior_table)[1:5]
	rownames(result)[8:17] = rownames(ab_test_table)[3:12]
	
	result[18:22,] = mcmc_table[1:5,]
	result[23:24,] = mcmc_table[6:7,]
	result[25:28,1] = mcmc_table[8:11,1]
	rownames(result)[18:28] = rownames(mcmc_table)[1:11]
	
	return(result)
}



# Function that summarizes and puts the indicators from objects implementing the LTT approach
# res_prior				[object]				Object containing outputs from priorEstimatesLTT function
# res_ab_test			[object]				Object containing outputs from ab_test function
# res_MCMC				[object]				Object containing outputs from bayesMCMC_2proportionsLTT function
summarize_LTT_plus_minus = function(res_prior, res_ab_test, res_MCMC, priorPsi = c(psiMu = 0, psiSigma = 1), priorProbabilities = c("H+" = 0.5, "H-" = 0.5, "H0" = 0)){
	#-- Bayes factors and post probabilities for ab_test
	ab_test_BF_plus_minus = res_ab_test$post_prob["H+"]/res_ab_test$post_prob["H-"] * priorProbabilities["H-"] / priorProbabilities["H+"]
	
	#-- Indicators for MCMC
	# --- Computation of Bayes factors and post probabilities by Savage-Dickey
	mcmc_samples = as.matrix(res_MCMC$model_jags$BUGSoutput$sims.matrix)
	
	MCMC_postBySD = bayesMCMC_SD_post_BF(mcmc_samples[, "psi"], 
		priorParameters = list(muPsi = priorPsi["psiMu"], sigmaPsi = priorPsi["psiSigma"], p_H0 = priorProbabilities["H-"], p_Hplus = priorProbabilities["H+"]))

	
	# BF+- = BF+0 / BF-0
	MCMC_BF_plusMinus = MCMC_postBySD$BF_plus0 / MCMC_postBySD$BF_minus0
	# pPost_plus = (BF+0 * pPrior_plus) / (BF+0 * pPrior_plus + BF-0 * pPrior_minus + p_H0)
	MCMC_postProb_LTT_Hplus = MCMC_postBySD$BF_plus0 * priorProbabilities["H+"] / (MCMC_postBySD$BF_plus0 * priorProbabilities["H+"] + MCMC_postBySD$BF_minus0 * priorProbabilities["H-"] + priorProbabilities["H0"])
	MCMC_postProb_LTT_Hminus = 1 - MCMC_postProb_LTT_Hplus

	
	result = c(priorProbabilities["H+"], priorProbabilities["H-"],
		res_ab_test$post_prob["H+"], res_ab_test$post_prob["H-"], ab_test_BF_plus_minus,
		MCMC_postProb_LTT_Hplus, MCMC_postProb_LTT_Hminus, MCMC_BF_plusMinus)
		
	names(result) = c("Prior prob H+", "Prior prob H-",
		"ab_test posterior Prob H+", "ab_test posterior Prob H-", "ab_test BF+-",
		"MCMC posterior Prob H+", "MCMC posterior Prob H-", "MCMC BF+-")
	
	return(result)
}



# Function that simulates alterations in a dataset. The alteration consists in adding a unique observation in each group considered, and doubling the sample size
# dataset			[data.frame]: 		R dataframe containing the dataset with binary data (i.e., each group is identified by 0 or 1)
# priorParameters	[vector]: 			R vector containing the parameters for the prior beta distributions
#										_ G1_alpha: alpha parameter of the beta distribution for Group 1
#										_ G1_beta: beta parameter of the beta distribution for Group 1
#										_ G2_alpha: alpha parameter of the beta distribution for Group 2
#										_ G2_beta: beta parameter of the beta distribution for Group 2
# sampleSize		[large integer]: 	sample size during random draws
# ciLevel			[numeric]:			Numeric between 0 and 1 excluded to customize the credible interval level
# chosenSeed		[numeric]:			Numeric value used for the simulations reproducibility
# Outputs: The function provides the results of the simulationIBE function (IBE simulation) and prob.ab (analytic computation of post probabilities) in a list
#			Each element of the list corresponds to a case (G1S: Another success in Group 1, G1F: ANother failure in Group 1, G2S: Another success in Group 2, G2F: Another failure in Group 2, Double: Doubling the sample size)
#			Each element contain: 
#				_ prob_analytic: the post probabilities (H+ against H-) analytically computed
#				_ simulation: results of the simulationIBE function
simulationIBE_alterData_add = function(dataset, priorParameters = c(1, 1, 1, 1), sampleSize = 1e5, ciLevel = 0.9, chosenSeed = 202505){
	# Distribution parameters
	alphaG1 = priorParameters[1]
	betaG1 = priorParameters[2]
	alphaG2 = priorParameters[3]
	betaG2 = priorParameters[4]
	
	
	# 1st case: Another success in Group 1
	alteredData = rbind(dataset, c(0, 1))
	
	currentContingence = table(alteredData$adolescentWS, alteredData$parentsWS)
	# Count of successes and failures
	nbS_1 = currentContingence[2,1]
	nbS_2 = currentContingence[2,2]
	nbC_1 = currentContingence[1,1]
	nbC_2 = currentContingence[1,2]
	
	set.seed(chosenSeed)
	simulation_G1S = simulationIBE(currentContingence, priorParameters = c(alphaG1, betaG1, alphaG2, betaG2), sampleSize = 1e5, ciLevel = 0.9)
	prob_BsupA_G1S = prob.ab(alphaG1 + nbS_1, betaG1 + nbC_1, alphaG2 + nbS_2, betaG2 + nbC_2)

	
	# 2nd case: Another failure in Group 1
	alteredData = rbind(dataset, c(0, 0))
	
	currentContingence = table(alteredData$adolescentWS, alteredData$parentsWS)
	# Count of successes and failures
	nbS_1 = currentContingence[2,1]
	nbS_2 = currentContingence[2,2]
	nbC_1 = currentContingence[1,1]
	nbC_2 = currentContingence[1,2]
	
	set.seed(chosenSeed)
	simulation_G1F = simulationIBE(currentContingence, priorParameters = c(alphaG1, betaG1, alphaG2, betaG2), sampleSize = 1e5, ciLevel = 0.9)
	prob_BsupA_G1F = prob.ab(alphaG1 + nbS_1, betaG1 + nbC_1, alphaG2 + nbS_2, betaG2 + nbC_2)


	# 3rd case: Another success in Group 2
	alteredData = rbind(dataset, c(1, 1))
	
	currentContingence = table(alteredData$adolescentWS, alteredData$parentsWS)
	# Count of successes and failures
	nbS_1 = currentContingence[2,1]
	nbS_2 = currentContingence[2,2]
	nbC_1 = currentContingence[1,1]
	nbC_2 = currentContingence[1,2]

	set.seed(chosenSeed)
	simulation_G2S = simulationIBE(currentContingence, priorParameters = c(alphaG1, betaG1, alphaG2, betaG2), sampleSize = 1e5, ciLevel = 0.9)
	prob_BsupA_G2S = prob.ab(alphaG1 + nbS_1, betaG1 + nbC_1, alphaG2 + nbS_2, betaG2 + nbC_2)

	
	# 4th case: Another failure in Group 2
	alteredData = rbind(dataset, c(1, 0))
	
	currentContingence = table(alteredData$adolescentWS, alteredData$parentsWS)
	# Count of successes and failures
	nbS_1 = currentContingence[2,1]
	nbS_2 = currentContingence[2,2]
	nbC_1 = currentContingence[1,1]
	nbC_2 = currentContingence[1,2]

	set.seed(chosenSeed)
	simulation_G2F = simulationIBE(currentContingence, priorParameters = c(alphaG1, betaG1, alphaG2, betaG2), sampleSize = 1e5, ciLevel = 0.9)
	prob_BsupA_G2F = prob.ab(alphaG1 + nbS_1, betaG1 + nbC_1, alphaG2 + nbS_2, betaG2 + nbC_2)
	
	
	# 5th case: Doubling the sample size
	alteredData = rbind(dataset, dataset)
	
	currentContingence = table(alteredData$adolescentWS, alteredData$parentsWS)
	# Count of successes and failures
	nbS_1 = currentContingence[2,1]
	nbS_2 = currentContingence[2,2]
	nbC_1 = currentContingence[1,1]
	nbC_2 = currentContingence[1,2]

	set.seed(chosenSeed)
	simulation_Double = simulationIBE(currentContingence, priorParameters = c(alphaG1, betaG1, alphaG2, betaG2), sampleSize = 1e5, ciLevel = 0.9)
	prob_BsupA_Double = prob.ab(alphaG1 + nbS_1, betaG1 + nbC_1, alphaG2 + nbS_2, betaG2 + nbC_2)

	results = list(
		G1S = list(
			prob_analytic = prob_BsupA_G1S,
			simulation = simulation_G1S
		),
		G1F = list(
			prob_analytic = prob_BsupA_G1F,
			simulation = simulation_G1F
		),
		G2S = list(
			prob_analytic = prob_BsupA_G2S,
			simulation = simulation_G2S
		),
		G2F = list(
			prob_analytic = prob_BsupA_G2F,
			simulation = simulation_G2F
		),
		Double = list(
			prob_analytic = prob_BsupA_Double,
			simulation = simulation_Double
		)
	)

	return(results)
}



# Function to summarize the simulationIBE_alterData_add output in a table
# results_alterData		[list]				List containing the simulationIBE_alterData_add output
# Output:		An array with indicators for each parameter of interest (Pi1, Pi2, and Delta) and for each case of alteration (G1S, G1F, G2S, G2F). The indicators are:
#				_ Analytic post prob: The analytic posterior probability P(H+ | data)
#				_ Quantiles: Quantiles of the simulated parameter of interest
#				_ HDI: 90% Highest Density Interval
#				_ ETI: 90% Equal-Tailed Interval
#				_ Sim post prob: The simulated posterior probability P(H+ | data)
#				_ BF: Bayes factor
summarizeAlteration_IBE = function(results_alterData){
	# Result table with its characteristics
	nrowsRes = (1 + 5 + 4) * 5					# (analytic posterior probability + 5 posterior quantiles + posterior probability + BF + 90% HDI + 90% ETI) for Additional observation in Group 1, in Group 2 and doubled sample size
	ncolsRes = 3								# G1, G2 and delta
	rowNamesRes = c(paste0("G1S ", c("Analytic post prob", paste0("Quantiles ", 0:4), "HDI", "ETI", "Sim post prob", "BF")), 
		paste0("G1F ", c("Analytic post prob", paste0("Quantiles ", 0:4), "HDI", "ETI", "Sim post prob", "BF")),
		paste0("G2S ", c("Analytic post prob", paste0("Quantiles ", 0:4), "HDI", "ETI", "Sim post prob", "BF")),
		paste0("G2F ", c("Analytic post prob", paste0("Quantiles ", 0:4), "HDI", "ETI", "Sim post prob", "BF")),
		paste0("Doubling ", c("Analytic post prob", paste0("Quantiles ", 0:4), "HDI", "ETI", "Sim post prob", "BF")))

	colNamesRes = c("G1", "G2", "Delta")
	
	# Result table
	tableResIBE = array(c(NA), dim = c(nrowsRes, ncolsRes))
	rownames(tableResIBE) = rowNamesRes
	colnames(tableResIBE) = colNamesRes

	tableResIBE[1,1] = results_alterData$G1S$prob_analytic
	tableResIBE[2:6,1] = results_alterData$G1S$simulation$postQuantiles_Pi1
	tableResIBE[2:6,2] = results_alterData$G1S$simulation$postQuantiles_Pi2
	tableResIBE[2:6,3] = results_alterData$G1S$simulation$postQuantiles_delta
	tableResIBE[7,1:3] = c(results_alterData$G1S$simulation$HDI_Pi1, results_alterData$G1S$simulation$HDI_Pi2, results_alterData$G1S$simulation$HDI_delta)
	tableResIBE[8,1:3] = c(results_alterData$G1S$simulation$ETI_Pi1, results_alterData$G1S$simulation$ETI_Pi2, results_alterData$G1S$simulation$ETI_delta)
	tableResIBE[9,1] = results_alterData$G1S$simulation$posteriorProbabilities
	tableResIBE[10,1] = results_alterData$G1S$simulation$bayesFactor

	tableResIBE[11,1] = results_alterData$G1F$prob_analytic
	tableResIBE[12:16,1] = results_alterData$G1F$simulation$postQuantiles_Pi1
	tableResIBE[12:16,2] = results_alterData$G1F$simulation$postQuantiles_Pi2
	tableResIBE[12:16,3] = results_alterData$G1F$simulation$postQuantiles_delta
	tableResIBE[17,1:3] = c(results_alterData$G1F$simulation$HDI_Pi1, results_alterData$G1F$simulation$HDI_Pi2, results_alterData$G1F$simulation$HDI_delta)
	tableResIBE[18,1:3] = c(results_alterData$G1F$simulation$ETI_Pi1, results_alterData$G1F$simulation$ETI_Pi2, results_alterData$G1F$simulation$ETI_delta)
	tableResIBE[19,1] = results_alterData$G1F$simulation$posteriorProbabilities
	tableResIBE[20,1] = results_alterData$G1F$simulation$bayesFactor

	tableResIBE[21,1] = results_alterData$G2S$prob_analytic
	tableResIBE[22:26,1] = results_alterData$G2S$simulation$postQuantiles_Pi1
	tableResIBE[22:26,2] = results_alterData$G2S$simulation$postQuantiles_Pi2
	tableResIBE[22:26,3] = results_alterData$G2S$simulation$postQuantiles_delta
	tableResIBE[27,1:3] = c(results_alterData$G2S$simulation$HDI_Pi1, results_alterData$G2S$simulation$HDI_Pi2, results_alterData$G2S$simulation$HDI_delta)
	tableResIBE[28,1:3] = c(results_alterData$G2S$simulation$ETI_Pi1, results_alterData$G2S$simulation$ETI_Pi2, results_alterData$G2S$simulation$ETI_delta)
	tableResIBE[29,1] = results_alterData$G2S$simulation$posteriorProbabilities
	tableResIBE[30,1] = results_alterData$G2S$simulation$bayesFactor

	tableResIBE[31,1] = results_alterData$G2F$prob_analytic
	tableResIBE[32:36,1] = results_alterData$G2F$simulation$postQuantiles_Pi1
	tableResIBE[32:36,2] = results_alterData$G2F$simulation$postQuantiles_Pi2
	tableResIBE[32:36,3] = results_alterData$G2F$simulation$postQuantiles_delta
	tableResIBE[37,1:3] = c(results_alterData$G2F$simulation$HDI_Pi1, results_alterData$G2F$simulation$HDI_Pi2, results_alterData$G2F$simulation$HDI_delta)
	tableResIBE[38,1:3] = c(results_alterData$G2F$simulation$ETI_Pi1, results_alterData$G2F$simulation$ETI_Pi2, results_alterData$G2F$simulation$ETI_delta)
	tableResIBE[39,1] = results_alterData$G2F$simulation$posteriorProbabilities
	tableResIBE[40,1] = results_alterData$G2F$simulation$bayesFactor
	
	tableResIBE[41,1] = results_alterData$Double$prob_analytic
	tableResIBE[42:46,1] = results_alterData$Double$simulation$postQuantiles_Pi1
	tableResIBE[42:46,2] = results_alterData$Double$simulation$postQuantiles_Pi2
	tableResIBE[42:46,3] = results_alterData$Double$simulation$postQuantiles_delta
	tableResIBE[47,1:3] = c(results_alterData$Double$simulation$HDI_Pi1, results_alterData$Double$simulation$HDI_Pi2, results_alterData$Double$simulation$HDI_delta)
	tableResIBE[48,1:3] = c(results_alterData$Double$simulation$ETI_Pi1, results_alterData$Double$simulation$ETI_Pi2, results_alterData$Double$simulation$ETI_delta)
	tableResIBE[49,1] = results_alterData$Double$simulation$posteriorProbabilities
	tableResIBE[50,1] = results_alterData$Double$simulation$bayesFactor

	return(tableResIBE)
}



# Function that simulates alterations in a dataset. The alteration consists in adding a unique observation in each group considered and doubling the sample size
# dataset				[data.frame]: 		R dataframe containing the dataset compliant for the ab_test function (i.e., each group is identified by 1 or 2 and each outcome by 0 or 1)
# priorParameters		[vector]: 			R vector containing the parameters for the prior distributions
#											_ gammaMu: parameter for the center of the gamma distribution in the LTT approach
#											_ gammaSigma: parameter for the standard deviation of the gamma distribution in the LTT approach
#											_ psiMu: parameter for the center of the psi distribution in the LTT approach
#											_ psiSigma: parameter for the standard deviation of the psi distribution in the LTT approach
# priorProbabilities	[vector]: 			R vector containing the prior probabilities to compare H+ and H0
#											_ H+
#											_ H0
# priorProbabilities2	[vector]: 			R vector containing the prior probabilities to compare H+ and H-
#											_ H+
#											_ H-
# sampleSize			[large integer]: 	sample size during random draws
# ciLevel				[numeric]:			Numeric between 0 and 1 excluded to customize the credible interval level
# chosenSeed			[numeric]:			Numeric value used for the simulations reproducibility
# Outputs: The function provides the results of the ab_test function (LTT simulation) in a list
#			Each element of the list corresponds to results of post distributions related to a case (G1S: Another success in Group 1, G1F: ANother failure in Group 1, G2S: Another success in Group 2, G2F: Another failure in Group 2, Double: Doubling the sample size)
#			Each result element contains: 
#				_ a table with a content similar to the summarize_ab_test output and post probabilities and Bayes factors estimates to compare H+ and H- (P(H+|data), P(H-|data) and BF+-)
#			Each result distribution contains: 
#				_ a data.frame with 5 columns (gamma, psi, pi1, pi2 and delta) and sampleSize rows
simulationLTT_alterData_add = function(dataset,  priorParameters = c(gammaMu = 0, gammaSigma = 1, psiMu = 0, psiSigma = 1), 
	priorProbabilities = c("H+" = 0.5, "H0" = 0.5), priorProbabilities2 = c("H+" = 0.5, "H-" = 0.5), sampleSize = 1e5, ciLevel = 0.9, chosenSeed = 202505){
	# Distribution parameters
	gammaMu = priorParameters["gammaMu"]
	gammaSigma = priorParameters["gammaSigma"]
	psiMu = priorParameters["psiMu"]
	psiSigma = priorParameters["psiSigma"]
	
	
	#---
	# 1st case: Another success in Group 1
	alteredData = rbind(dataset, c(1, 1))
	
	#-- Comparison between H+ and H0
	set.seed(chosenSeed)
	ab_test_LTT_G1S = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities, nsamples = sampleSize)

	tableResLTT_G1S = summarize_ab_test(ab_test_LTT_G1S, ciLevel = ciLevel)

	#-- Comparison between H+ and H-
	set.seed(chosenSeed)
	ab_test_LTT2_G1S = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities2, nsamples = sampleSize)

	# Bayes factors and post probabilities for ab_test
	ab_test_BF_plus_minus_G1S = ab_test_LTT2_G1S$post_prob["H+"]/ab_test_LTT2_G1S$post_prob["H-"] * priorProbabilities2["H-"] / priorProbabilities2["H+"]
	
	tableFinalLTT_G1S = rbind(tableResLTT_G1S,
		c(ab_test_LTT2_G1S$post_prob["H+"], rep(NA, 4)),
		c(ab_test_LTT2_G1S$post_prob["H-"], rep(NA, 4)),
		c(ab_test_BF_plus_minus_G1S, rep(NA, 4)))
	rownames(tableFinalLTT_G1S)[13:15] = c("ab_test posterior Prob H+", "ab_test posterior Prob H-", "ab_test BF+-")
	
	# Save the distribution in order to draw the graphs
	res_G1S = summary(ab_test_LTT_G1S)
	distributions_G1S = data.frame(
		gamma = res_G1S$post$Hplus$beta,
		psi = res_G1S$post$Hplus$psi,
		pi1 = res_G1S$post$Hplus$p1,
		pi2 = res_G1S$post$Hplus$p2,
		delta = res_G1S$post$Hplus$p2 - res_G1S$post$Hplus$p1
	)
	
	#---
	# 2nd case: Another failure in Group 1
	alteredData = rbind(dataset, c(1, 0))
	
	#-- Comparison between H+ and H0
	set.seed(chosenSeed)
	ab_test_LTT_G1F = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities, nsamples = sampleSize)
	
	tableResLTT_G1F = summarize_ab_test(ab_test_LTT_G1F, ciLevel = ciLevel)

	#-- Comparison between H+ and H-
	set.seed(chosenSeed)
	ab_test_LTT2_G1F = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities2, nsamples = sampleSize)

	# Bayes factors and post probabilities for ab_test
	ab_test_BF_plus_minus_G1S = ab_test_LTT2_G1F$post_prob["H+"]/ab_test_LTT2_G1F$post_prob["H-"] * priorProbabilities2["H-"] / priorProbabilities2["H+"]
	
	tableFinalLTT_G1F = rbind(tableResLTT_G1F,
		c(ab_test_LTT2_G1F$post_prob["H+"], rep(NA, 4)),
		c(ab_test_LTT2_G1F$post_prob["H-"], rep(NA, 4)),
		c(ab_test_BF_plus_minus_G1S, rep(NA, 4)))
	rownames(tableFinalLTT_G1F)[13:15] = c("ab_test posterior Prob H+", "ab_test posterior Prob H-", "ab_test BF+-")
	
	# Save the distribution in order to draw the graphs
	res_G1F = summary(ab_test_LTT_G1F)
	distributions_G1F = data.frame(
		gamma = res_G1F$post$Hplus$beta,
		psi = res_G1F$post$Hplus$psi,
		pi1 = res_G1F$post$Hplus$p1,
		pi2 = res_G1F$post$Hplus$p2,
		delta = res_G1F$post$Hplus$p2 - res_G1F$post$Hplus$p1
	)
	
	#---
	# 3rd case: Another success in Group 2
	alteredData = rbind(dataset, c(2, 1))
	
	#-- Comparison between H+ and H0
	set.seed(chosenSeed)
	ab_test_LTT_G2S = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities, nsamples = sampleSize)
	
	tableResLTT_G2S = summarize_ab_test(ab_test_LTT_G2S, ciLevel = ciLevel)

	#-- Comparison between H+ and H-
	set.seed(chosenSeed)
	ab_test_LTT2_G2S = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities2, nsamples = sampleSize)

	# Bayes factors and post probabilities for ab_test
	ab_test_BF_plus_minus_G2S = ab_test_LTT2_G2S$post_prob["H+"]/ab_test_LTT2_G2S$post_prob["H-"] * priorProbabilities2["H-"] / priorProbabilities2["H+"]
	
	tableFinalLTT_G2S = rbind(tableResLTT_G2S,
		c(ab_test_LTT2_G2S$post_prob["H+"], rep(NA, 4)),
		c(ab_test_LTT2_G2S$post_prob["H-"], rep(NA, 4)),
		c(ab_test_BF_plus_minus_G2S, rep(NA, 4)))
	rownames(tableFinalLTT_G2S)[13:15] = c("ab_test posterior Prob H+", "ab_test posterior Prob H-", "ab_test BF+-")
	
	# Save the distribution in order to draw the graphs
	res_G2S = summary(ab_test_LTT_G2S)
	distributions_G2S = data.frame(
		gamma = res_G2S$post$Hplus$beta,
		psi = res_G2S$post$Hplus$psi,
		pi1 = res_G2S$post$Hplus$p1,
		pi2 = res_G2S$post$Hplus$p2,
		delta = res_G2S$post$Hplus$p2 - res_G2S$post$Hplus$p1
	)
	
	#---
	# 4th case: Another failure in Group 2
	
	#-- Comparison between H+ and H0
	alteredData = rbind(dataset, c(2, 0))
	set.seed(chosenSeed)
	ab_test_LTT_G2F = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities, nsamples = sampleSize)
	
	tableResLTT_G2F = summarize_ab_test(ab_test_LTT_G2F, ciLevel = ciLevel)

	#-- Comparison between H+ and H-
	set.seed(chosenSeed)
	ab_test_LTT2_G2F = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities2, nsamples = sampleSize)

	# Bayes factors and post probabilities for ab_test
	ab_test_BF_plus_minus_G2F = ab_test_LTT2_G2F$post_prob["H+"]/ab_test_LTT2_G2F$post_prob["H-"] * priorProbabilities2["H-"] / priorProbabilities2["H+"]
	
	tableFinalLTT_G2F = rbind(tableResLTT_G2F,
		c(ab_test_LTT2_G2F$post_prob["H+"], rep(NA, 4)),
		c(ab_test_LTT2_G2F$post_prob["H-"], rep(NA, 4)),
		c(ab_test_BF_plus_minus_G2F, rep(NA, 4)))
	rownames(tableFinalLTT_G2F)[13:15] = c("ab_test posterior Prob H+", "ab_test posterior Prob H-", "ab_test BF+-")
	
	# Save the distribution in order to draw the graphs
	res_G2F = summary(ab_test_LTT_G2F)
	distributions_G2F = data.frame(
		gamma = res_G2F$post$Hplus$beta,
		psi = res_G2F$post$Hplus$psi,
		pi1 = res_G2F$post$Hplus$p1,
		pi2 = res_G2F$post$Hplus$p2,
		delta = res_G2F$post$Hplus$p2 - res_G2F$post$Hplus$p1
	)
	
	#---
	# 5th case: Doubling the sample size
	alteredData = rbind(dataset, dataset)
	
	#-- Comparison between H+ and H0
	set.seed(chosenSeed)
	ab_test_LTT_Double = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities, nsamples = sampleSize)
	
	tableResLTT_Double = summarize_ab_test(ab_test_LTT_Double, ciLevel = ciLevel)

	#-- Comparison between H+ and H-
	set.seed(chosenSeed)
	ab_test_LTT2_Double = ab_test(data = alteredData,
		prior_par = list(mu_psi = psiMu, sigma_psi = psiSigma, mu_beta = gammaMu, sigma_beta = gammaSigma),
		prior_prob = priorProbabilities2, nsamples = sampleSize)

	# Bayes factors and post probabilities for ab_test
	ab_test_BF_plus_minus_Double = ab_test_LTT2_Double$post_prob["H+"]/ab_test_LTT2_Double$post_prob["H-"] * priorProbabilities2["H-"] / priorProbabilities2["H+"]
	
	tableFinalLTT_Double = rbind(tableResLTT_Double,
		c(ab_test_LTT2_Double$post_prob["H+"], rep(NA, 4)),
		c(ab_test_LTT2_Double$post_prob["H-"], rep(NA, 4)),
		c(ab_test_BF_plus_minus_Double, rep(NA, 4)))
	rownames(tableFinalLTT_Double)[13:15] = c("ab_test posterior Prob H+", "ab_test posterior Prob H-", "ab_test BF+-")
	
	# Save the distribution in order to draw the graphs
	res_Double = summary(ab_test_LTT_Double)
	distributions_Double = data.frame(
		gamma = res_Double$post$Hplus$beta,
		psi = res_Double$post$Hplus$psi,
		pi1 = res_Double$post$Hplus$p1,
		pi2 = res_Double$post$Hplus$p2,
		delta = res_Double$post$Hplus$p2 - res_Double$post$Hplus$p1
	)
	
	#---
	# Return the result
	return(list(
		G1S = tableFinalLTT_G1S,
		G1F = tableFinalLTT_G1F,
		G2S = tableFinalLTT_G2S,
		G2F = tableFinalLTT_G2F,
		Double = tableFinalLTT_Double,
		distributions_G1S = distributions_G1S,
		distributions_G1F = distributions_G1F,
		distributions_G2S = distributions_G2S,
		distributions_G2F = distributions_G2F,
		distributions_Double = distributions_Double
	))
}



#----------------
# Functions used to draw the figures

# Function used to draw Figure 2:
# Prior distributions used in the present study for each scenario investigated (1. Non-informative, 2. Informative, 3. Optimistic):
#   (a) βeta distributions used to estimate the π1 and π2 parameters in the IBE approach,
#   (b) Normal distributions used to estimate the γ and ψ parameters in the LTT approach.
# IBEparameters			[data.frame]:		R data.frame containing the parameters (alpha and beta) of the prior beta distributions with the following columns
#											alpha: the alpha parameter
#											beta: the beta parameter
# LTTparameters			[data.frame]:		R data.frame containing the LTT prior parameters
#											gamma_mu: the mean parameter for gamma
#											gamma_sigma: the sd parameter for gamma
#											psi_mu: the mean parameter for psi
#											psi_sigma: the sd parameter for psi
# Return the ggplot2 graph object
drawPriorDistributions = function(IBEparameters, LTTparameters){
	
	###---###
	# IBE approach	

	x = seq(0, 1, by = 0.01)									# x axis of the distribution on [0,1] interval
	
	alphaT1 = IBEparameters$alpha[1]
	betaT1 = IBEparameters$beta[1]
	alphaT2 = IBEparameters$alpha[2]
	betaT2 = IBEparameters$beta[2]
	alphaT3 = IBEparameters$alpha[3]
	betaT3 = IBEparameters$beta[3]
	alphaT4 = IBEparameters$alpha[4]
	betaT4 = IBEparameters$beta[4]
	
	# Object dedicated to the storage of x and y prior distributions
	dataDensity = data.frame(
		x = x,
		y1 = dbeta(x, alphaT1, betaT1),
		y2 = dbeta(x, alphaT2, betaT2),
		y3 = dbeta(x, alphaT3, betaT3),
		y4 = dbeta(x, alphaT4, betaT4)
	)
	
	dataGraph = dataDensity %>%
		pivot_longer(cols = c(y1:y4),
			names_to = "factorName",
			values_to = "values")

	# Texts in the legend
	legendLabels = c(paste0("Non-informative: \nBeta (", alphaT1, ", ", betaT1, ")"),
		paste0("Informative: \nBeta (", alphaT2, ", ", betaT2, ")"),
		paste0("Optimistic Group 1: \nBeta (", alphaT3, ", ", betaT3, ")"),
		paste0("Optimistic Group 2: \nBeta (", alphaT4, ", ", betaT4, ")"))
	
	# Draw the graph
	graphLine_IBE_priorDistributions = ggplot(dataGraph, 
			aes(x = x, y = values, color = factorName)) + 
		geom_line(linewidth = 1.2) + 
		scale_color_brewer(palette = "Set2", labels = legendLabels) + 
		labs(x = "", y = "", color = "IBE prior distribution", title = "") + 
		theme_minimal() +
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),
			title = element_text(size = titleSize))
	
	
	###---###
	# LTT approach
	
	x = seq(-4, 4, by = 0.01)								# x axis of the Gamma and Psi distributions on [-4,4] interval
	
	# Object dedicated to the storage of x and y prior distributions
	dataDensity = data.frame(
		x = x,
		gamma_1 = dnorm(x, LTTparameters$gamma_mu[1], LTTparameters$gamma_sigma[1]),				# Non-informative
		psi_1 = dnorm(x, LTTparameters$psi_mu[1], LTTparameters$psi_sigma[1]),
		gamma_2 = dnorm(x, LTTparameters$gamma_mu[2], LTTparameters$gamma_sigma[2]),				# Informative
		psi_2 = dnorm(x, LTTparameters$psi_mu[2], LTTparameters$psi_sigma[2]),
		gamma_3 = dnorm(x, LTTparameters$gamma_mu[3], LTTparameters$gamma_sigma[3]),				# Optimistic
		psi_3 = dnorm(x, LTTparameters$psi_mu[3], LTTparameters$psi_sigma[3])
	)
	
	dataGraph = dataDensity %>%
		pivot_longer(cols = c(gamma_1:psi_3),
			names_to = "factorName",
			values_to = "values") %>%
		mutate(
			parameter = case_when(
				str_detect(factorName, "gamma") ~ "Gamma",
				str_detect(factorName, "psi") ~ "Psi",
				TRUE ~ NA_character_
			),
			parameter2show = factor(parameter, levels = c("Gamma", "Psi"), 
				labels = c(bquote( ~ gamma), bquote( ~ psi))),
			scenario = case_when(
				str_detect(factorName, "_1") ~ "1. Non-informative",
				str_detect(factorName, "_2") ~ "2. Informative",
				str_detect(factorName, "_3") ~ "3. Optimistic",
				TRUE ~ NA_character_
			)
		)
	
	# Texts in the legend
	tildaString = "~"
	legendLabels = c(bquote(atop(atop(Non-informative: phantom(), ~ gamma ~ .(tildaString) ~ N(.(LTTparameters$gamma_mu[1]), .(LTTparameters$gamma_sigma[1])^2) ~ and ~ psi ~ .(tildaString) ~ N(.(LTTparameters$psi_mu[1]), .(LTTparameters$psi_sigma[1])^2)), phantom(0))),
		bquote(atop(atop(Informative: phantom(), ~ gamma ~ .(tildaString) ~ N(.(LTTparameters$gamma_mu[2]), .(LTTparameters$gamma_sigma[2])^2) ~ and ~ psi ~ .(tildaString) ~ N(.(LTTparameters$psi_mu[2]), .(LTTparameters$psi_sigma[2])^2)), phantom(0))),
		bquote(atop(atop(Optimistic: phantom(), ~ gamma ~ .(tildaString) ~ N(.(LTTparameters$gamma_mu[3]), .(LTTparameters$gamma_sigma[3])^2) ~ and ~ psi ~ .(tildaString) ~ N(.(LTTparameters$psi_mu[3]), .(LTTparameters$psi_sigma[3])^2)), phantom(0))))

	graphLine_LTT_priorDistributions = ggplot(dataGraph, aes(x = x, y = values, color = scenario)) + 
		facet_grid(. ~ parameter2show, labeller = labeller(parameter2show = label_parsed)) +
		geom_line(linewidth = 1.2) + 
		scale_color_brewer(palette = "Set2", labels = legendLabels) + 
		labs(x = "", y = "", color = "LTT prior distribution", title = "") + 
		theme_minimal() +
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize * 1.5),
			strip.text.x = element_text(size = titleSize),
			title = element_text(size = titleSize))
		
	###---###
	# Aggregation of the graphs
	
	graphArrange_graphLine_priorDistribution = grid.arrange(graphLine_IBE_priorDistributions, graphLine_LTT_priorDistributions,
		ncol=1, nrow=2)
	
	# Add the graph title and letters (a) and (b)
	graphArrange_graphLine_priorDistribution_Def = as_ggplot(graphArrange_graphLine_priorDistribution) +                                # transform to a ggplot
		draw_plot_label(label = c("Prior distributions"), size = titleSize,
			x = c(0), y = c(1)) +
		draw_plot_label(label = c("(a)", "(b)"), size = mean(c(contentSize, titleSize)),
			x = c(0.02, 0.02), y = c(0.95, 0.5))

	graphResult = graphLine_IBE_priorDistributions
	graphResult = graphLine_LTT_priorDistributions
	graphResult = graphArrange_graphLine_priorDistribution_Def
	
	# Return the graph
	return(graphResult)
}


# Function used to draw Figure 3:
# Prior and posterior βeta distributions of the IBE approach applied on Groups 1 (upper row) and 2 (lower row) for the Move Ô case study for each scenario investigated: 1. Non-informative, 2. Informative, 3. Optimistic.
# IBEparameters			[data.frame]:		R data.frame containing the parameters (alpha and beta) of the prior beta distributions with the following columns
#											scenario: the scenario considered
#											group: the group considered
#											alpha: the alpha parameter
#											beta: the beta parameter
# contingency			[table]:			contingency table: groups are in columns and successes/failures are in rows
# Return the ggplot2 graph object
drawPriorPosteriorDistributions_IBE = function(IBEparameters, contingency){
	
	# Parameters for the beta distributions
	alphaNIg1 = IBEparameters$alpha[which(IBEparameters$scenario == "1. Non-informative" & IBEparameters$group == "Group 1")]
	betaNIg1 = IBEparameters$beta[which(IBEparameters$scenario == "1. Non-informative" & IBEparameters$group == "Group 1")]
	alphaNIg2 = IBEparameters$alpha[which(IBEparameters$scenario == "1. Non-informative" & IBEparameters$group == "Group 2")]
	betaNIg2 = IBEparameters$beta[which(IBEparameters$scenario == "1. Non-informative" & IBEparameters$group == "Group 2")]
	
	alphaIg1 = IBEparameters$alpha[which(IBEparameters$scenario == "2. Informative" & IBEparameters$group == "Group 1")]
	betaIg1 = IBEparameters$beta[which(IBEparameters$scenario == "2. Informative" & IBEparameters$group == "Group 1")]
	alphaIg2 = IBEparameters$alpha[which(IBEparameters$scenario == "2. Informative" & IBEparameters$group == "Group 2")]
	betaIg2 = IBEparameters$beta[which(IBEparameters$scenario == "2. Informative" & IBEparameters$group == "Group 2")]
	
	alphaOg1 = IBEparameters$alpha[which(IBEparameters$scenario == "3. Optimistic" & IBEparameters$group == "Group 1")]
	betaOg1 = IBEparameters$beta[which(IBEparameters$scenario == "3. Optimistic" & IBEparameters$group == "Group 1")]
	alphaOg2 = IBEparameters$alpha[which(IBEparameters$scenario == "3. Optimistic" & IBEparameters$group == "Group 2")]
	betaOg2 = IBEparameters$beta[which(IBEparameters$scenario == "3. Optimistic" & IBEparameters$group == "Group 2")]
	
	# Number of successes and failures per group
	nbS_1 = currentContingence[2,1]
	nbS_2 = currentContingence[2,2]
	nbC_1 = currentContingence[1,1]
	nbC_2 = currentContingence[1,2]
	
	# Object dedicated to the storage of x and y prior and posterior distributions
	x = seq(0, 1, by = 0.01)										# x axis interval [0,1]
	dataDensity = data.frame(
		x = x,
		y1_1 = dbeta(x, alphaNIg1, betaNIg1),						# Non informative
		y1_2 = dbeta(x, alphaNIg2, betaNIg2),
		y2_1 = dbeta(x, alphaIg1, betaIg1),							# Informative
		y2_2 = dbeta(x, alphaIg2, betaIg2),
		y3_1 = dbeta(x, alphaOg1, betaOg1),							# Optimistic
		y3_2 = dbeta(x, alphaOg2, betaOg2),
		z1_1 = dbeta(x, alphaNIg1 + nbS_1, betaNIg1 + nbC_1),		# Non informative
		z1_2 = dbeta(x, alphaNIg2 + nbS_2, betaNIg2 + nbC_2),
		z2_1 = dbeta(x, alphaIg1 + nbS_1, betaIg1 + nbC_1),			# Informative
		z2_2 = dbeta(x, alphaIg2 + nbS_2, betaIg2 + nbC_2),
		z3_1 = dbeta(x, alphaOg1 + nbS_1, betaOg1 + nbC_1),			# Optimistic
		z3_2 = dbeta(x, alphaOg2 + nbS_2, betaOg2 + nbC_2)
	)
	
	dataGraph = dataDensity %>%
		pivot_longer(cols = c(y1_1:z3_2),
			names_to = "factorName",
			values_to = "values") %>%
		mutate(
			legende = case_when(
				factorName == "y1_1" ~ paste0("alpha = ", alphaNIg1, ", beta = ", betaNIg1),
				factorName == "y1_2" ~ paste0("alpha = ", alphaNIg2, ", beta = ", betaNIg2),
				factorName == "y2_1" ~ paste0("alpha = ", alphaIg1, ", beta = ", betaIg1),
				factorName == "y2_2" ~ paste0("alpha = ", alphaIg2, ", beta = ", betaIg2),
				factorName == "y3_1" ~ paste0("alpha = ", alphaOg1, ", beta = ", betaOg1),
				factorName == "y3_2" ~ paste0("alpha = ", alphaOg2, ", beta = ", betaOg2),
				factorName == "z1_1" ~ paste0("alpha = ", alphaNIg1 + nbS_1, ", beta = ", betaNIg1 + nbC_1),
				factorName == "z1_2" ~ paste0("alpha = ", alphaNIg2 + nbS_2, ", beta = ", betaNIg2 + nbC_2),
				factorName == "z2_1" ~ paste0("alpha = ", alphaIg1 + nbS_1, ", beta = ", betaIg1 + nbC_1),
				factorName == "z2_2" ~ paste0("alpha = ", alphaIg2 + nbS_2, ", beta = ", betaIg2 + nbC_2),
				factorName == "z3_1" ~ paste0("alpha = ", alphaOg1 + nbS_1, ", beta = ", betaOg1 + nbC_1),
				factorName == "z3_2" ~ paste0("alpha = ", alphaOg2 + nbS_2, ", beta = ", betaOg2 + nbC_2),
				TRUE ~ NA_character_
			),
			type = factor(case_when(
				str_detect(factorName, "y") ~ "Prior",
				str_detect(factorName, "z") ~ "Posterior",
				TRUE ~ NA_character_
			), c("Prior", "Posterior")),
			group = case_when(
				str_detect(factorName, "_1") ~ "Group 1",
				str_detect(factorName, "_2") ~ "Group 2",
				TRUE ~ NA_character_
			),
			scenario = case_when(
				str_detect(factorName, "1_") ~ "1. Non-informative",
				str_detect(factorName, "2_") ~ "2. Informative",
				str_detect(factorName, "3_") ~ "3. Optimistic",
				TRUE ~ NA_character_
			)
		)

	graphLine_IBE_priorAndPosteriorDistributions = ggplot(dataGraph, aes(x = x, y = values, color = type)) + 
		geom_line(linewidth = 1.2) + 
		facet_grid(group ~ scenario) +
		scale_color_brewer(palette = "Set2") +
		labs(x = "", y = "", color = "Distribution", title = "") + 
		theme_minimal() +
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),
			strip.text.x = element_text(size = contentSize, face="bold"), strip.text.y = element_text(size = contentSize, face="bold"),
			title = element_text(size = titleSize))

	return(graphLine_IBE_priorAndPosteriorDistributions)
}


# Function used to draw Figure 4:
# Distributions of the difference ( δ=π2−π1 ) simulated in the prior and posterior settings for the Move Ô case study by the three implementations (our own simulation coding, bayesAB package, and MCMC with R2jags) of the IBE approach for each scenario investigated: 1. Non-informative, 2. Informative, 3. Optimistic.
# data4graph			[data.frame]:		data.frame in a longer format design containing the following columns:
#											scenario: the scenario considered (non-informative, informative or optimistic)
#											method: the method used (prior, simulation, bayesAB or MCMC)
#											method2show: the method in a factor format instead of text
#											values: the delta values for each sample
# Return the ggplot2 graph object
drawDeltaDistributions_IBE = function(data4graph){
	
	# Text legend
	legendLabels = c("Prior simulation", "rbeta-based sampling", "bayesAB", "MCMC with R2jags")
	
	deltaDensities = ggplot(data4graph, aes(x = values, fill = method2show, color = method2show)) + 
		facet_grid(. ~ scenario) +
		geom_density(alpha = 0.4, linewidth = 1.2) + 
		labs(x = expression(paste("Difference", ~delta)), y = "", title = "", color = "Method", fill = "Method") +
		scale_color_manual(values = colors4delta_col, labels = legendLabels) + 
		scale_fill_manual(values = colors4delta_fill, labels = legendLabels) + 
		theme_minimal() + 
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),
			strip.text.x = element_text(size = contentSize, face="bold"), strip.text.y = element_text(size = contentSize, face="bold"),
			title = element_text(size = titleSize))
	
	return(deltaDensities)
}


# Function dedicated to the aggregation of results from LTT approach
# data					[list]:				list containing the results of the simulations, each element of the list basically corresponds to a scenario.
#											Each element of the list is constituted by:
#											_ nonInformative: the results returned by priorEstimatesLTT
#											_ informative: the results returned by ab_test
#											_ optimistic: the results returned by bayesMCMC_2proportionsLTT
# scenarios				[vector]:			Character vector describing the names of the scenarios. Its length must be the length of the data list
# Returns a list ready to be used in the drawDistributions_LTT function.
# The resulting list contains 3 elements:
# gamma_psi_data: the data.frame containing the data for the simulations of gamma and psi parameters
# pi1_pi2_data: the data.frame containing the data for the simulations of pi1 and pi2 parameters
# delta_data: the data.frame containing the data for the simulations of delta parameter
dataAggregation4graphLTT = function(data, scenarios = c("1. Non-informative", "2. Informative", "3. Optimistic")){
	nScenarios = length(scenarios)
	if(length(data) == nScenarios){
		for(j in 1:nScenarios){
			currentData = data[[j]]
			simPrior_res = currentData$prior
			ab_test_res = currentData$ab_test
			mcmc_res = as.matrix(currentData$mcmc$model_jags$BUGSoutput$sims.matrix)
			
			# Gamma and Psi distributions
			dataGrapheGammaPsiLTT = data.frame(
				scenario = scenarios[j],
				method = c(rep("prior", length(simPrior_res$gamma) * 2), 
					rep("ab_test", length(ab_test_res$post$Hplus$beta) * 2),
					rep("MCMC", length(mcmc_res[,"mu"]) * 2)),
				parameter = c(rep("gamma", length(simPrior_res$gamma)), rep("psi", length(simPrior_res$psi)), 
					rep("gamma", length(ab_test_res$post$Hplus$beta)), rep("psi", length(ab_test_res$post$Hplus$psi)),
					rep("gamma", length(mcmc_res[,"mu"])), rep("psi", length(mcmc_res[,"psi"]))),
				values = c(simPrior_res$gamma, simPrior_res$psi, 
					ab_test_res$post$Hplus$beta, ab_test_res$post$Hplus$psi,
					mcmc_res[,"mu"], mcmc_res[,"psi"])
			)
			
			# Pi1 and Pi2 distributions
			dataGraphePisLTT = data.frame(
				scenario = scenarios[j],
				method = c(rep("prior", length(simPrior_res$p1) * 2), 
					rep("ab_test", length(ab_test_res$post$Hplus$p1) * 2),
					rep("MCMC", length(mcmc_res[,"pi1"]) * 2)),
				parameter = c(rep("pi1", length(simPrior_res$p1)), rep("pi2", length(simPrior_res$p2)), 
					rep("pi1", length(ab_test_res$post$Hplus$p1)), rep("pi2", length(ab_test_res$post$Hplus$p2)),
					rep("pi1", length(mcmc_res[,"pi1"])), rep("pi2", length(mcmc_res[,"pi2"]))),
				values = c(simPrior_res$p1, simPrior_res$p2, 
					ab_test_res$post$Hplus$p1, ab_test_res$post$Hplus$p2,
					mcmc_res[,"pi1"], mcmc_res[,"pi2"])
			)
			
			# Delta distributions
			# Estimates of delta with ab_test
			ab_test_difference = ab_test_res$post$Hplus$p2 - ab_test_res$post$Hplus$p1
			
			dataGrapheDeltaLTT = data.frame(
				scenario = scenarios[j],
				method = c(rep("prior", length(simPrior_res$delta)), 
					rep("ab_test", length(ab_test_difference)),
					rep("MCMC", length(mcmc_res[,"delta"]))),
				parameter = "delta",
				values = c(simPrior_res$delta, 
					ab_test_difference,
					mcmc_res[,"delta"])
			)
			
			if(j == 1){
				res_dataGrapheGammaPsiLTT = dataGrapheGammaPsiLTT
				res_dataGraphePisLTT = dataGraphePisLTT
				res_dataGrapheDeltaLTT = dataGrapheDeltaLTT
			} else {
				res_dataGrapheGammaPsiLTT = rbind(res_dataGrapheGammaPsiLTT, dataGrapheGammaPsiLTT)
				res_dataGraphePisLTT = rbind(res_dataGraphePisLTT, dataGraphePisLTT)
				res_dataGrapheDeltaLTT = rbind(res_dataGrapheDeltaLTT, dataGrapheDeltaLTT)
			}
		}
		
		# Add columns as factors instead of characters
		res_dataGrapheGammaPsiLTT = res_dataGrapheGammaPsiLTT %>%
			mutate(method2show = factor(method, levels = c("prior", "ab_test", "MCMC"),
					labels = c("Prior simulation", "ab_test", "MCMC with R2jags")),
				parameter2show = factor(parameter, levels = c("gamma", "psi"), 
					labels = c(bquote( ~ gamma), bquote( ~ psi))))
		res_dataGraphePisLTT = res_dataGraphePisLTT %>%
			mutate(method2show = factor(method, levels = c("prior", "ab_test", "MCMC"),
					labels = c("Prior simulation", "ab_test", "MCMC with R2jags")),
				parameter2show = factor(parameter, levels = c("pi1", "pi2"), 
				labels = c(bquote( ~ pi[1]), bquote( ~ pi[2]))))
		res_dataGrapheDeltaLTT = res_dataGrapheDeltaLTT %>%
			mutate(method2show = factor(method, levels = c("prior", "ab_test", "MCMC"),
					labels = c("Prior simulation", "ab_test", "MCMC with R2jags")),
				parameter2show = factor(parameter, levels = c("delta"), 
				labels = c(bquote( ~ delta))))
		result = list(gamma_psi_data = res_dataGrapheGammaPsiLTT,
			pi1_pi2_data = res_dataGraphePisLTT,
			delta_data = res_dataGrapheDeltaLTT)
	} else {
		print("Problem: the number of scenarios should be equal the length of data")
		result = NULL
	}
	return(result)
}


# Function used to draw Figure 5:
# Distributions of the γ, ψ, π1 and π2 parameters and the difference δ=π2−π1 simulated under the prior and posterior settings for the Move Ô case study by the two implementations (ab_test from abtest package, and MCMC with R2jags) of the LTT approach for each scenario investigated: 1. Non-informative, 2. Informative, 3. Optimistic.
# data4graph			[list]:				list containing the results of the simulations. This is the result of the dataAggregation4graphLTT function.
drawDistributions_LTT = function(data4graph){
	# Graph with gamma and psi parameters
	dataGrapheGammaPsiLTT = data4graph$gamma_psi_data
	gammaPsiDensitiesLTT = ggplot(dataGrapheGammaPsiLTT, aes(x = values, fill = method2show, color = method2show)) + 
		facet_grid(parameter2show ~ scenario, labeller = labeller(parameter2show = label_parsed)) +
		geom_density(alpha = 0.2, linewidth = 1.2) + 
		labs(x = "", y = "", title = "", color = "Method", fill = "Method") +
		scale_color_manual(values = colors4delta_col[-2]) + 
		scale_fill_manual(values = colors4delta_fill[-2]) + 
		theme_minimal() + 
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),	#contentSize * 1.5
			strip.text.x = element_text(size = titleSize, face="bold"), strip.text.y = element_text(size = titleSize, face="bold"),
			title = element_text(size = titleSize))
	
	# Graph with pi1 and pi2 parameters
	dataGraphePisLTT = data4graph$pi1_pi2_data
	pisDensitiesLTT = ggplot(dataGraphePisLTT, aes(x = values, fill = method2show, color = method2show)) + 
		facet_grid(parameter2show ~ scenario, labeller = labeller(parameter2show = label_parsed)) +
		geom_density(alpha = 0.2, linewidth = 1.2) + 
		labs(x = "", y = "", title = "", color = "Method", fill = "Method") +
		scale_color_manual(values = colors4delta_col[-2]) + 
		scale_fill_manual(values = colors4delta_fill[-2]) + 
		theme_minimal() + 
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),	#contentSize * 1.5
			strip.text.x = element_blank(), strip.text.y = element_text(size = titleSize, face="bold"),
			title = element_text(size = titleSize))
	
	# Graph with delta parameter
	dataGrapheDeltaLTT = data4graph$delta_data
	deltaDensitiesLTT = ggplot(dataGrapheDeltaLTT, aes(x = values, fill = method2show, color = method2show)) + 
		facet_grid(parameter2show ~ scenario, labeller = labeller(parameter2show = label_parsed)) +
		geom_density(alpha = 0.2, linewidth = 1.2) + 
		labs(x = "", y = "", title = "", color = "Method", fill = "Method") +
		scale_color_manual(values = colors4delta_col[-2]) + 
		scale_fill_manual(values = colors4delta_fill[-2]) + 
		theme_minimal() + 
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),
			strip.text.x = element_blank(), strip.text.y = element_text(size = titleSize, face="bold"),
			title = element_text(size = titleSize))

	# Aggregation of the graphs with shared legend
	legend = g_legend(gammaPsiDensitiesLTT + theme(legend.position='right'))
	graphResult = grid_arrange_shared_legend(gammaPsiDensitiesLTT, pisDensitiesLTT, deltaDensitiesLTT,
		ncol = 1, nrow = 3, position='right', heights=c(2,2,1))
	
	return(graphResult)
}


# Function used to draw Figure 6:
# Posterior βeta distributions of the IBE approach applied to Groups 1 (upper row) and 2 (lower row) for the altered data from the Move Ô case study under each scenario investigated: 1. Non-informative, 2. Informative, 3. Optimistic.
# IBEparameters			[data.frame]:		R data.frame containing the parameters (alpha and beta) of the prior beta distributions with the following columns
#											scenario: the scenario considered
#											group: the group considered
#											alpha: the alpha parameter
#											beta: the beta parameter
# contingency			[table]:			contingency table: groups are in columns and successes/failures are in rows
# Return the ggplot2 graph object
drawDataAlterationDistributions_IBE = function(IBEparameters, contingency){
	
	# Parameters for the beta distributions
	alphaNIg1 = IBEparameters$alpha[which(IBEparameters$scenario == "1. Non-informative" & IBEparameters$group == "Group 1")]
	betaNIg1 = IBEparameters$beta[which(IBEparameters$scenario == "1. Non-informative" & IBEparameters$group == "Group 1")]
	alphaNIg2 = IBEparameters$alpha[which(IBEparameters$scenario == "1. Non-informative" & IBEparameters$group == "Group 2")]
	betaNIg2 = IBEparameters$beta[which(IBEparameters$scenario == "1. Non-informative" & IBEparameters$group == "Group 2")]
	
	alphaIg1 = IBEparameters$alpha[which(IBEparameters$scenario == "2. Informative" & IBEparameters$group == "Group 1")]
	betaIg1 = IBEparameters$beta[which(IBEparameters$scenario == "2. Informative" & IBEparameters$group == "Group 1")]
	alphaIg2 = IBEparameters$alpha[which(IBEparameters$scenario == "2. Informative" & IBEparameters$group == "Group 2")]
	betaIg2 = IBEparameters$beta[which(IBEparameters$scenario == "2. Informative" & IBEparameters$group == "Group 2")]
	
	alphaOg1 = IBEparameters$alpha[which(IBEparameters$scenario == "3. Optimistic" & IBEparameters$group == "Group 1")]
	betaOg1 = IBEparameters$beta[which(IBEparameters$scenario == "3. Optimistic" & IBEparameters$group == "Group 1")]
	alphaOg2 = IBEparameters$alpha[which(IBEparameters$scenario == "3. Optimistic" & IBEparameters$group == "Group 2")]
	betaOg2 = IBEparameters$beta[which(IBEparameters$scenario == "3. Optimistic" & IBEparameters$group == "Group 2")]
	
	# Number of successes and failures per group
	nbS_1 = currentContingence[2,1]
	nbS_2 = currentContingence[2,2]
	nbC_1 = currentContingence[1,1]
	nbC_2 = currentContingence[1,2]
	
	# Object dedicated to the storage of x and y prior and posterior distributions
	x = seq(0, 1, by = 0.01)												# x axis interval [0,1]
	dataDensity = data.frame(
		x = x,
		y1_1 = dbeta(x, alphaNIg1, betaNIg1),								# Non informative
		y1_2 = dbeta(x, alphaNIg2, betaNIg2),
		y2_1 = dbeta(x, alphaIg1, betaIg1),									# Informative
		y2_2 = dbeta(x, alphaIg2, betaIg2),
		y3_1 = dbeta(x, alphaOg1, betaOg1),									# Optimistic
		y3_2 = dbeta(x, alphaOg2, betaOg2),
		# Original data
		z1_1_OD = dbeta(x, alphaNIg1 + nbS_1, betaNIg1 + nbC_1),			# Non informative
		z1_2_OD = dbeta(x, alphaNIg2 + nbS_2, betaNIg2 + nbC_2),
		z2_1_OD = dbeta(x, alphaIg1 + nbS_1, betaIg1 + nbC_1),				# Informative
		z2_2_OD = dbeta(x, alphaIg2 + nbS_2, betaIg2 + nbC_2),
		z3_1_OD = dbeta(x, alphaOg1 + nbS_1, betaOg1 + nbC_1),				# Optimistic
		z3_2_OD = dbeta(x, alphaOg2 + nbS_2, betaOg2 + nbC_2),
		# G1S: Adding a success in Group 1
		z1_1_G1S = dbeta(x, alphaNIg1 + nbS_1 + 1, betaNIg1 + nbC_1),		# Non informative
		z2_1_G1S = dbeta(x, alphaIg1 + nbS_1 + 1, betaIg1 + nbC_1),			# Informative
		z3_1_G1S = dbeta(x, alphaOg1 + nbS_1 + 1, betaOg1 + nbC_1),			# Optimistic
		# G1F: Adding a failure in Group 1
		z1_1_G1F = dbeta(x, alphaNIg1 + nbS_1, betaNIg1 + nbC_1 + 1),		# Non informative
		z2_1_G1F = dbeta(x, alphaIg1 + nbS_1, betaIg1 + nbC_1 + 1),			# Informative
		z3_1_G1F = dbeta(x, alphaOg1 + nbS_1, betaOg1 + nbC_1 + 1),			# Optimistic
		# G2S: Adding a success in Group 2
		z1_2_G2S = dbeta(x, alphaNIg2 + nbS_2 + 1, betaNIg2 + nbC_2),		# Non informative
		z2_2_G2S = dbeta(x, alphaIg2 + nbS_2 + 1, betaIg2 + nbC_2),			# Informative
		z3_2_G2S = dbeta(x, alphaOg2 + nbS_2 + 1, betaOg2 + nbC_2),			# Optimistic
		# G2F: Adding a failure in Group 2
		z1_2_G2F = dbeta(x, alphaNIg2 + nbS_2, betaNIg2 + nbC_2 + 1),		# Non informative
		z2_2_G2F = dbeta(x, alphaIg2 + nbS_2, betaIg2 + nbC_2 + 1),			# Informative
		z3_2_G2F = dbeta(x, alphaOg2 + nbS_2, betaOg2 + nbC_2 + 1),			# Optimistic
		# Double
		z1_1_Double = dbeta(x, alphaNIg1 + nbS_1 * 2, betaNIg1 + nbC_1 * 2),	# Non informative
		z1_2_Double = dbeta(x, alphaNIg2 + nbS_2 * 2, betaNIg2 + nbC_2 * 2),
		z2_1_Double = dbeta(x, alphaIg1 + nbS_1 * 2, betaIg1 + nbC_1 * 2),		# Informative
		z2_2_Double = dbeta(x, alphaIg2 + nbS_2 * 2, betaIg2 + nbC_2 * 2),
		z3_1_Double = dbeta(x, alphaOg1 + nbS_1 * 2, betaOg1 + nbC_1 * 2),		# Optimistic
		z3_2_Double = dbeta(x, alphaOg2 + nbS_2 * 2, betaOg2 + nbC_2 * 2)
	)
	
	dataGraph = dataDensity %>%
		pivot_longer(cols = c(y1_1:z3_2_Double),
			names_to = "factorName",
			values_to = "values") %>%
		mutate(
			type = factor(case_when(
				str_detect(factorName, "y") ~ "Prior",
				str_detect(factorName, "z") ~ "Posterior",
				TRUE ~ NA_character_
			), c("Prior", "Posterior")),
			group = case_when(
				str_detect(factorName, "_1_") ~ "Group 1",
				str_detect(factorName, "_2_") ~ "Group 2",
				TRUE ~ NA_character_
			),
			scenario = case_when(
				(str_detect(factorName, "y1_") | str_detect(factorName, "z1_")) ~ "1. Non-informative",
				(str_detect(factorName, "y2_") | str_detect(factorName, "z2_")) ~ "2. Informative",
				(str_detect(factorName, "y3_") | str_detect(factorName, "z3_")) ~ "3. Optimistic",
				TRUE ~ NA_character_
			),
			alteration = factor(case_when(
				str_detect(factorName, "OD") ~ "Original data",
				str_detect(factorName, "G1S") ~ "Another success in Group 1",
				str_detect(factorName, "G1F") ~ "Another failure in Group 1",
				str_detect(factorName, "G2S") ~ "Another success in Group 2",
				str_detect(factorName, "G2F") ~ "Another failure in Group 2",
				str_detect(factorName, "Double") ~ "Doubled sample size",
				TRUE ~ NA_character_
			), c("Original data", "Another success in Group 1", "Another failure in Group 1", 
			"Another success in Group 2", "Another failure in Group 2", "Doubled sample size")),
		)

	graphLine_IBE_alteredDataDistributions = ggplot(dataGraph %>% filter(type == "Posterior"), aes(x = x, y = values, color = alteration)) + 
		geom_line(linewidth = 1.2) + 
		facet_grid(group ~ scenario) +
		scale_color_manual(values = colorsalteration_col) +
		labs(x = "", y = "", color = "Data alteration", title = "") + 
		theme_minimal() +
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),
			strip.text.x = element_text(size = contentSize, face="bold"), strip.text.y = element_text(size = contentSize, face="bold"),
			title = element_text(size = titleSize))

	return(graphLine_IBE_alteredDataDistributions)
}


# Function dedicated to the aggregation of results from LTT approach with data alteration
# data					[list]:				list containing the results of the simulations, each element of the list basically corresponds to an alteration.
#											Each element of the list is constituted by:
#											_ nonInformative: the results returned by priorEstimatesLTT
#											_ informative: the results returned by ab_test
#											_ optimistic: the results returned by bayesMCMC_2proportionsLTT
# scenarios				[vector]:			Character vector describing the names of the scenarios. Its length must be the length of the data list
# Returns a data.frame ready to be used in the drawAlterationDistributions_LTT function.
# The resulting list contains 6 columns:
# alteration: the alteration considered
# factorName: the name of the parameter considered (gamma, psi, pi1, pi2, delta)
# values: the observed values reported during the simulation
# scenario: the scenarion considered (Non-informative, Informative, Optimistic)
# alteration2show: the factor formated of alteration
# parameter2show: an alternative of the factorName, formated to be used shown on graphs
dataAggregation4graphDataAlterationLTT = function(data, scenarios = c("1. Non-informative", "2. Informative", "3. Optimistic")){
	nScenarios = length(scenarios)
	if(length(data) == nScenarios){
		for(j in 1:nScenarios){
			currentData = data[[j]]
			orignalData_res = currentData$originalData
			G1S_res = currentData$G1S
			G1F_res = currentData$G1F
			G2S_res = currentData$G2S
			G2F_res = currentData$G2F
			double_res = currentData$double			
			
			# Add the column of the alteration
			orignalData_res$alteration = "Original data"
			G1S_res$alteration = "Another success in Group 1"
			G1F_res$alteration = "Another failure in Group 1"
			G2S_res$alteration = "Another success in Group 2"
			G2F_res$alteration = "Another failure in Group 2"
			double_res$alteration = "Doubled sample size"
			
			# Aggregate and pivot the data.frame
			dataAlterationDistributions = rbind(orignalData_res, G1S_res, G1F_res, G2S_res, G2F_res, double_res) %>%
				pivot_longer(
					cols = c(gamma:delta),
					names_to = "factorName",
					values_to = "values"
				) %>%
				mutate(scenario = scenarios[j])

			if(j == 1){
				res_dataAlterationDistributions = dataAlterationDistributions
			} else {
				res_dataAlterationDistributions = rbind(res_dataAlterationDistributions, dataAlterationDistributions)
			}
		}
		
		# Add columns as factors instead of characters
		res_dataAlterationDistributions = res_dataAlterationDistributions %>%
			mutate(alteration2show = factor(alteration, c("Original data", "Another success in Group 1", "Another failure in Group 1",
				"Another success in Group 2", "Another failure in Group 2", "Doubled sample size")),
				parameter2show = factor(factorName, levels = c("gamma", "psi", "pi1", "pi2", "delta"), 
				labels = c(bquote( ~ gamma), bquote( ~ psi), bquote( ~ pi[1]), bquote( ~ pi[2]), bquote( ~ delta))))
		
		result = res_dataAlterationDistributions
	} else {
		print("Problem: the number of scenarios should be equal the length of data")
		result = NULL
	}
	return(result)
}


# Function used to draw Figure 7:
# Posterior distributions of the γ , ψ , π1 and π2 parameters and the difference δ=π2−π1 simulated for the altered data from the Move Ô case study by the abtest implementation of the LTT approach for each scenario investigated: 1. Non-informative, 2. Informative, 3. Optimistic.
# data4graph			[list]:				data.frame containing the distributions of the simulated parameters after a data alteration. This is the result of the dataAggregation4graphDataAlterationLTT function.
drawAlterationDistributions_LTT = function(data4graph){
	# Graph with gamma and psi parameters
	dataGrapheGammaPsiLTT = data4graph %>%
		filter(factorName %in% c("gamma", "psi"))
	gammaPsiDensitiesLTT = ggplot(dataGrapheGammaPsiLTT, aes(x = values, fill = alteration2show, color = alteration2show)) + 
		facet_grid(parameter2show ~ scenario, labeller = labeller(parameter2show = label_parsed)) +
		geom_density(alpha = 0.2, linewidth = 1.2) + 
		labs(x = "", y = "", title = "", color = "Data alteration", fill = "Data alteration") +
		scale_color_manual(values = colorsalteration_col) + 
		scale_fill_manual(values = colorsalteration_col) + 
		theme_minimal() + 
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),	#contentSize * 1.5
			strip.text.x = element_text(size = titleSize, face="bold"), strip.text.y = element_text(size = titleSize, face="bold"),
			title = element_text(size = titleSize))
	
	# Graph with pi1 and pi2 parameters
	dataGraphePisLTT = data4graph %>%
		filter(factorName %in% c("pi1", "pi2"))
	pisDensitiesLTT = ggplot(dataGraphePisLTT, aes(x = values, fill = alteration2show, color = alteration2show)) + 
		facet_grid(parameter2show ~ scenario, labeller = labeller(parameter2show = label_parsed)) +
		geom_density(alpha = 0.2, linewidth = 1.2) + 
		labs(x = "", y = "", title = "", color = "Data alteration", fill = "Data alteration") +
		scale_color_manual(values = colorsalteration_col) + 
		scale_fill_manual(values = colorsalteration_col) + 
		theme_minimal() + 
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),	#contentSize * 1.5
			strip.text.x = element_blank(), strip.text.y = element_text(size = titleSize, face="bold"),
			title = element_text(size = titleSize))
	
	# Graph with delta parameter
	dataGrapheDeltaLTT = data4graph %>%
		filter(factorName %in% c("delta"))
	deltaDensitiesLTT = ggplot(dataGrapheDeltaLTT, aes(x = values, fill = alteration2show, color = alteration2show)) + 
		facet_grid(parameter2show ~ scenario, labeller = labeller(parameter2show = label_parsed)) +
		geom_density(alpha = 0.2, linewidth = 1.2) + 
		labs(x = "", y = "", title = "", color = "Data alteration", fill = "Data alteration") +
		scale_color_manual(values = colorsalteration_col) + 
		scale_fill_manual(values = colorsalteration_col) + 
		theme_minimal() + 
		theme(axis.title.x = element_text(size = titleSize), axis.title.y = element_text(size = titleSize), legend.title = element_text(size = titleSize),
			axis.text.x = element_text(size = contentSize), axis.text.y = element_text(size = contentSize), legend.text = element_text(size = contentSize),
			strip.text.x = element_blank(), strip.text.y = element_text(size = titleSize, face="bold"),
			title = element_text(size = titleSize))

	# Aggregation of the graphs with shared legend
	legend = g_legend(gammaPsiDensitiesLTT + theme(legend.position='right'))
	graphResult = grid_arrange_shared_legend(gammaPsiDensitiesLTT, pisDensitiesLTT, deltaDensitiesLTT,
		ncol = 1, nrow = 3, position='right', heights=c(2,2,1))
	
	return(graphResult)
}


###########################################-- End --####################################################