R/internal.PLSRejection.R

In bomeara/treevo: Using ABC to Understand Trait Evolution

#  Partial Least Squares Rejection
#
#  This function performs the ABC-rejection analysis using an input simulation
#  data. Particles are accepted is they fall sufficiently close to the target
#  data (within the tolerance). Distances are calculated using \code{abcDistance}.
#

#  @param summaryValuesMatrix Matrix of summary statistics from simulations.

#  @param trueFreeValuesMatrix Matrix of true free values from simulations.

#  @inheritParams doSimulation
#  @inheritParams doRun

#  @param abcTolerance Proportion of accepted simulations.

#  @param verbose option to message progress to screen.

#  @param validation Cross Validation procedure for ABC.

#  @param scale Scale for \code{pls.model.list}.

#  @param variance.cutoff Variance cutoff for \code{pls.model.list}.

#  @return Returns a list of the particle data frame and ABC distances.

#  @author Brian O'Meara and Barb Banbury

# @references O'Meara and Banbury, unpublished
# @keywords PLSRejection doRun doRun_rej abc

# @examples
#
# PLSRejection(
#   summaryValuesMatrix, 
#   trueFreeValuesMatrix, 
#   phy = simPhyExample, 
#   traits = simCharExample, 
#   abcTolerance = resultsBMExample[[1]]$abcTolerance
#   )
#

# This is a difficult function to write an example for, as mainly an algorithm for doRun_rej; unclear if even needs to be exported
# Will internalize and see if it makes any difference


#  @name PLSRejection
#  @rdname PLSRejection
#  @export
PLSRejection <- function(
		summaryValuesMatrix, 
		trueFreeValuesMatrix, 
		phy, traits, abcTolerance,
		verbose = TRUE, 
		validation = "CV", 
		scale = TRUE, 
		variance.cutoff = 95
		) {
	##############################################################################
	originalSummaryValues <- summaryStatsLong(phy = phy, traits = traits)
	if (verbose) {
		message("Done getting originalSummaryValues")
		}
	abcDistancesRaw <- sapply(
		sequence(dim(trueFreeValuesMatrix)[2]), 
		SingleParameterPLSDistanceSquared,
		summaryValuesMatrix = summaryValuesMatrix, 
		trueFreeValuesMatrix = trueFreeValuesMatrix, 
		originalSummaryValues = originalSummaryValues, 
		validation = validation, 
		scale = scale, 
		variance.cutoff = variance.cutoff
		)
	abcDistancesRawTotal <- apply(abcDistancesRaw, 1, sum)
	abcDistances <- sqrt(abcDistancesRawTotal) #Euclid rules.
	#
	# here's where we do abc - calculate quantile on distances
		# drop what distances are greater than that cutoff
	abcQuantile <- quantile(abcDistances, 
		prob = abcTolerance
		)
	acceptedParticles <- trueFreeValuesMatrix[
		which(abcDistances <= abcQuantile), , drop = FALSE] 
	#
	acceptedDistances <- abcDistances[
		which(abcDistances <= quantile(abcDistances, prob = abcTolerance))
		]
	#
	particleDataFrame <- data.frame(cbind(
		rep(1, dim(acceptedParticles)[1]), 
		as.vector(which(abcDistances <= quantile(abcDistances, prob = abcTolerance))), 
		seq(1:dim(acceptedParticles)[1]), 
		rep(0, dim(acceptedParticles)[1]), 
		acceptedDistances, rep(1, 
		dim(acceptedParticles)[1]), 
		acceptedParticles
		))
	colnames(particleDataFrame) <- c(
		"generation", "attempt", "id", 
		"parentid", "distance", "weight", 
		paste("param", 
			seq(dim(trueFreeValuesMatrix)[2]), 
			sep = ""
			)
		)
	res <- list(particleDataFrame = particleDataFrame, 
		abcDistances = abcDistances)
	return(res)
	}