# We require the function "sim_run" which is in the previous script
# If you don't run that script the function won't work
# We could add the following command to the end of the function script
# but I think it's more important to think about generating this output

# sapply takes a vector and for each value in the vector runs the function
# You add additional argument in the apply command, with the names for required arguments

# rep .2, 100 times, means that we have a vector that looks like this:
# [.2 .2 .2 .2 .... .2] and there are 50 elements
# The number of elements here is the number of simulations we will conduct

sapply(rep(.2,50),sim_run,asq=.3,esq=.5,nmz=500,ndz=500)->chi_from_c

# Now we have the object chi_from_c which are the -2*log-likelihood differences
# We can then manipulate the chi_from_c object in a couple of ways

# Let's plot the results

split.screen(c(1,1),erase=T)

plot.density(density(chi_from_c),xlab="Chi Square",ylab="Density",main="Testing for C")

# We can figure out our power from the chi_from_c object as well
# We do so by seeing how often our chi square is greater than our significance threshold

# Using R we can actually figure out what our significance threshold is by
# the qchisq command, which gives you the quantile of the Chi Square distribution for a given 
# probability. Our probability is the alpha level, which we will set at 0.05
# R by default considers the probability from left to right (i.e. it starts at minus infinity)
# We want to know what the threshold is to leave just the top 5% of the distribution

qchisq(0.05, 1, lower.tail=F)

# This gives us back a value of 3.84......

# We can add a horizontal line to our plot of chis to visually assess our power
# abline adds a line to the plot, either defined by a slope and intercept
# or by adding a horizontal or vertical line at a given value (h= for horizontal) (v = for vertical)

abline(v=qchisq(0.05,1,lower.tail=F))

# Now we need to define the subset of our simulated results which are more significant than
# our Chi Square as defined by the alpha level
# To do so, we will specify that the subset of the chis can be defined using the '[]' on the object

length(chi_from_c[chi_from_c>qchisq(0.05,1,lower.tail=F)])/length(chi_from_c)

# The length commmand counts up the number of observations of the object (subsetted or not)
# The above command should give a single value between 0 and 1