PosteriorPredictive_PValues.Rev

################################################################################
#
# RevBayes Example: Bayesian model testing using posterior predictive simulation
# This file calculates the Tree Summary Statistics
#
# authors: Lyndon M. Coghill, Sebastian Hoehna and Jeremy M. Brown
#
################################################################################



################### start of effect size calculations ########################

empData = readDataDelimitedFile(emp_pps_file,header=TRUE,delimiter=",")
simData = readDataDelimitedFile(sim_pps_file,header=TRUE,delimiter=",")

## setup the pvalue output file
write(file=outfileName, "Statistic", "Lower 1-tailed", "Upper 1-tailed", "Two-tailed", "Midpoint", "Effect Size", sep=",", append=FALSE)
write(file=outfileName, "\n", sep=",", append=TRUE)  

for (x in 2:simData.ncolumns()) {

    ## transform the data for easily manipulation
    ## we need to retrieve a column
    numbers <- simData.column(x)

    ## calculate median value of PPD here
    m = median(numbers)

    ## calculate effect size here
    empValue = empData[1][x-1]
    effect_size = abs((m - empValue) / stdev(numbers))
    
    ## Calculate and return a vector of lower, equal, and upper pvalues for a given test statistic
    p_values <- posteriorPredictiveProbability(numbers, empValue)
    
    ## 1-tailed
    lower_p_value <- p_values[1] 
	equal_p_value <- p_values[2] 
    upper_p_value <- p_values[3] 
	
    ## mid-point
    midpoint_p_value = lower_p_value + 0.5*equal_p_value

    ## 2-tailed
    two_tail_p_value = 2 * (min(v(lower_p_value+equal_p_value, upper_p_value+equal_p_value)))

    write(file=outfileName, statID[x], lower_p_value+equal_p_value, upper_p_value+equal_p_value, two_tail_p_value, midpoint_p_value, effect_size, sep=",", append=TRUE)
    write(file=outfileName, "\n", sep=",", append=TRUE)  

}
################### end of effect size calculations ##########################