################################################
#Note:  First run the syntax in Power Paper Functions to load the functions.
# Required packages are MBESS and stats
# We are currently working on making the code general and robust and placing it as a package in R
# The PPOWER function currently requires wise choices for lowern and uppern.  The following examples work and replicate
# the analyses presented in 
# Biesanz, J. C., & Schrager, S. M. (2010).  Sample size planning with effect size estimates. Manuscript under review.
################################################

################################################
#Function for estimating average statistical power across the posterior distribution of the effect size estimate
# PPOWER(effect, n1, n2, DesiredPower, lowern, uppern, alpha, type, posterior)
#Variable definitions:
#  effect:  the observed effect size
#  n1			sample size for group 1 (or observed sample size for the correlation)
#  n2			sample size for group 2 (NA in the case of the correlation)
#  DesiredPower	value of desired average power
#  lowern		Lower sample size bounds to examine
#  uppern		Upper sample size bounds to examine
#  alpha		Type I error rate (e.g., .05)
#  type			effect size type.  Options are "d" and "r".
#  posterior	posterior distribution of the effect size.  May be user-supplied (i.e., from a HB analysis).
#				Otherwise specify as NA to generate a posterior from the noninformative prior. 
#
#Returns:
#  1.  Sample size required to achieve the specified level of DesiredPower
#  2.  Median 90% confidence interval for the effect size based on that sample size
#  3.  Quantiles for power.
#  4.  A graph of Average power and Confidence interval width across the range of lowern to uppern.
################################################
#Note:  Currently the EB and HB routines are written for standardize mean differences.  Correlations are an easy extension
################################################

#Noninformative prior example
PPOWER(effect=.5, n1=25,n2=25, DesiredPower=.80, lowern=100, uppern=240, alpha=.05, type="d", posterior=NA)




# Hierarchical Bayesian analysis using the data from Dunn et al. (2007)
# First we call the HB function to produce draws from the posterior distribution using the other 4 effect size estimates
# and then we input these posterior draws into the PPOWER function.

dunn<-read.table("lizcohend.txt", header=TRUE)
posHB <- HB(dunn$estimate, dunn$sd, dunn$df, 10000, 5, 100000)
PPOWER(effect=.769, n1=23,n2=23, DesiredPower=.80, lowern=25, uppern=200, alpha=.05, type="d", posterior= posHB)



#Empirical Bayesian analysis using the littering example and the attitude-behavior meta-analysis
attd<-read.table("attituded.txt", header=TRUE)
vardbeh <-((attd$n1+attd$n2)/(attd$n1*attd$n2)+(attd$dbeh*attd$dbeh)/
	      (2*(attd$n1+attd$n2-2)))*((attd$n1+attd$n2)/(attd$n1+attd$n2-2))
df <- (attd$n1+attd$n2-2)
dout <- reml(attd$dbeh, vardbeh)

#estimated standard error for the initial effect size estimate of d=.50.
sd <- sqrt(((25+25)/(25*25)+(.5* .5)/(2*(25+25-2)))*((25+25)/(25+25-2))) 

# First we call the EB function to produce draws from the posterior distribution using the REML estimates from the meta-analysis
# and then we input these posterior draws into the PPOWER function.
posEB <- EB(estimate=.5, sd=sd,  n1=25, n2=25, k=50000, mu=dout$mu, tau2=dout$vart, metan1 = attd$n1, metan2=attd$n2)
PPOWER(effect=.5, n1=25,n2=25, DesiredPower=.80, lowern=100, uppern=240, alpha=.05, type="d", posterior= posEB)