#' Conduct an ABC on a range of effect size distribution hyperparameters
#'
#' Runs Approximate Bayesian Computation across a range of marker effect size
#' distribution hyperparameters by comparing the distributions of phenotypes produced to those in the real data.
#'
#' Please note that no missing phenotypes or missing genotype data is permitted. Missing genotypes should be imputed
#' (such as with \code{\link{impute_and_phase_beagle}}) before this function is run.
#'
#' @param x matrix. Input genotypes, SNPs as rows, columns as individuals. Genotypes formatted as 0,1,2 for the major homozygote, heterozygote, and minor homozygote, respectively.
#' @param phenotypes numeric vector. Observed phenotypes, one per individual.
#' @param iters numeric. Number of ABC permutations to run.
#' @param effect_distribution function, default rbayesB. A function for a distribution of effect sizes. The first argument, n,
#' should be the number of effects to draw. Other arguments must be named to match the names of the parameter distribution
#' functions given in the parameter_distributions argument.
#' @param parameter_distributions List containing named functions, defualt
#' list(pi = function(x) rbeta(x, 25, 1), d.f = function(x) runif(x, 1, 100), scale = function(x) rbeta(x, 1, 3)*100).
#' A list containing functions for each parameter in the effect size distribution. The given functions should be named to
#' match the parameter for which they produce distributions, and should take a single argument, x, which holds the number of
#' parameter values to draw.
#' @param h_dist function, default function(x) rep(.5, x). A function for the distribution of heritability from which to draw
#' h^2 values for each iteration. Must take a single argument, x, which holds the number of h values to draw. The default
#' produces only heritabilities of .5.
#' @param center logical, default T. Determines if the phenotypes provided and generated during each iteration
#' should be centered (have their means set to 0).
#' @param par numeric or FALSE, default FALSE. If numeric, the number of cores on which to run the ABC.
#'
#' @export
<- (x, phenotypes, iters,
effect_distribution = ,
parameter_distributions = ( = (x) (x, 25, 1),
d.f = (x) (x, 1, 100),
= (x) (x, 1, 3)*100),
h_dist = (x) (.5, x),
center = ,
joint_res = ,
joint_acceptance = ,
joint_res_dist = ,
= , phased = ,
save_effects = , = 2000,
save_effects_nkeep = (iters * .1),
save_effects_keep_var = "ks"){
(!(save_effects)){
look_for_files <- (save_effects, (".bk", ".rds"))
files_here <- (look_for_files)
((files_here)){
(("File(s) ", (look_for_files[files_here], collapse = ", "), " already exist(s).\n"))
}
}
#============ABC_scheme function for one rep=============
scheme_D <- (x, phenotypes, effect_distribution, parameters, h, center = center, phased = F, save_effects = ){
pseudo <- (x, effect_distribution, parameters, h, center = center, phased = phased)
<- (phenotypes, pseudo$p)
(save_effects){
(( = , = pseudo$e))
}
()
}
#============initialization======================================
# center the distribution if requested
(center){
phenotypes <- phenotypes - (phenotypes)
}
# get the random values to run
run_parameters <- sample_parameters_from_distributions(parameter_distributions = parameter_distributions,
joint_res = joint_res, iters = iters,
joint_acceptance = joint_acceptance,
joint_res_dist = joint_res_dist,
reg_res = , = )
nparms <- (run_parameters)
# run_parameters <- vector("list", length = length(parameter_distributions))
# names(run_parameters) <- names(parameter_distributions)
# for(i in 1:length(run_parameters)){
# run_parameters[[i]] <- parameter_distributions[[i]](iters)
# }
# run_parameters <- as.data.frame(run_parameters)
h <- h_dist(iters)
# initialize storage
dist_output <- (0, iters, = ())
(dist_output) <-
#============ABC loop===================
# run the ABC
## serial
( == F | == 1){
(i 1:iters){
("Iter:", i, "\n")
d <- scheme_D(x = x,
phenotypes = phenotypes,
effect_distribution = effect_distribution,
parameters = (run_parameters[i,, = F]),
h = h[i],
center = center,
phased = phased, save_effects = !(save_effects))
(!(save_effects)){
data.table::fwrite((run_parameters[i,, = F], h = h[i], data.table::as.data.table((d$, = 1))),
(save_effects, ".tmp"), sep = "\t", col.names = , = , = )
dist_output[i,] <- d$
}
{
dist_output[i,] <- d
}
}
dist_output <- (dist_output)
# filter saved effects
# read effects into an FBM and save if requested (only the keepers!)
(!(save_effects)){
("Collating written effects into an FBM.\n")
dist_output$keep <-
dist_output$keep[order(dist_output,save_effects_keep_var])[1:save_effects_nkeep]] <- # keep only the nkeep best runs
effect_fbm <- bigstatsr::big_read((save_effects, ".tmp"),
= (dist_output$keep),
select = 1:((x) + nparms + 1))
((save_effects, ".tmp"))
# save fbm
effect_fbm <- effect_fbm$()
}
((as.data.table(run_parameters), h = h, as.data.table(dist_output)))
}
# parallel
{
<- (, iters)
cl <- snow::makeSOCKcluster()
doSNOW::registerDoSNOW(cl)
# divide up into ncore chunks
it_par <- (1:iters)%%
chunks <- (parms = .smart_split(run_parameters, it_par),
dist_output = .smart_split((dist_output), it_par),
h = .smart_split((h), it_par))
(run_parameters, dist_output)
# prepare reporting function
progress <- (n) (("Chunk %d out of", n), , "is complete.\n")
opts <- (progress=progress)
need_removal <- (save_effects, ".", 1:, ".tmp.part")
need_removal <- need_removal[which((need_removal))]
((need_removal) > 0){
(need_removal)
}
output <- foreach::foreach( = 1:, .inorder = , .errorhandling = "pass",
.options.snow = opts, = ("data.table", "GeneArchEst", "bigstatsr")
) %dopar% {
parm_chunk <- chunks$parm[q]]
h_chunk <- (chunks$h[q]])
dist_chunk <- chunks$dist_output[q]]
is.err <- (0)
errs <- (0)
# run once per iter in this chunk
(i 1:(parm_chunk)){
b <- (scheme_D(x = x,
phenotypes = phenotypes,
effect_distribution = effect_distribution,
parameters = (parm_chunk[i,, = F]),
h = h_chunk[i],
center = center,
phased = phased,
save_effects = !(save_effects)), silent = T)
((b) == "try-error"){
is.err <- (is.err, i)
errs <- (errs, b)
}
{
(!(save_effects)){
data.table::fwrite((parm_chunk[i,, = F], h = h_chunk[i], data.table::as.data.table((b$, = 1))),
(save_effects, ".", , ".tmp.part"), sep = "\t", col.names = , = , = )
dist_chunk[i,] <- b$
}
{
dist_chunk[i,] <- b
}
}
}
(!(save_effects)){
dist_chunk$q_id <-
dist_chunk$iter_id <- 1:(parm_chunk)
}
out <- ("list", 2)
(out) <- ("successes", "fails")
((is.err) > 0){
out[1]] <- (parm_chunk-is.err,], h = h_chunk-is.err], dist_chunk-is.err,])
out[2]] <- (error_msg = errs,
error_parms = (parm_chunk[is.err,, = F], h = h_chunk[is.err]))
}
{
out[1]] <- (parm_chunk, h = h_chunk, dist_chunk)
out[2]] <- (error_msg = (0),
error_parms = ((, = (parm_chunk) + 1, = 1)))
(out[2]]$error_parms) <- ((parm_chunk), "h")
}
out
}
# release cores and clean up
::(cl)
doSNOW::registerDoSNOW()
();()
("Run complete. Number of good runs per job:\n\t")
dims <- ((purrr::map(output, 1), ))
((1:(dims), ": ", dims, "\n\t"))
("Colnames:\n\t")
cns <- (purrr::map(output, 1), )
cns <- (cns, , collapse = "\t")
((1:(dims), ": ", cns, "\n\t"))
dist_output <- dplyr::bind_rows(purrr::map(output, 1))
errs <- purrr::map(output, 2)
err_parms <- dplyr::bind_rows(purrr::map(errs, 2))
err_msgs <- (purrr::map(errs, 1))
err_parms <- (err_parms)
errs <- (parms = err_parms, msgs = err_msgs)
# read effects into an FBM and save if requested (only the keepers!)
(!(save_effects)){
("Collating written effects into an FBM.\n")
dist_output$keep <-
dist_output$keep[order(dist_output,save_effects_keep_var])[1:save_effects_nkeep]] <- # keep only the nkeep best runs
effect_fbm <- bigstatsr::FBM(save_effects_nkeep, (x) + nparms + 1, backingfile = save_effects)
progress <- 0
(i 1:){
correct_iter <- (dist_output$q_id == i)
write_end <- (dist_output$keep[correct_iter])
(write_end == 0){
}
tmpfile_bk <- stringi::stri_rand_strings(1, 20)
((tmpfile_bk)){
tmpfile_bk <- stringi::stri_rand_strings(1, 20)
}
tmp_fbm <- bigstatsr::big_read((save_effects, ".", i, ".tmp.part"),
select = 1:((x) + nparms + 1),
= (dist_output$keep[correct_iter]),
backingfile = tmpfile_bk)
effect_fbm(progress + 1):(progress + write_end),] <- tmp_fbm]
(tmp_fbm)
();()
((tmpfile_bk, ".bk"))
((tmpfile_bk, ".rds"))
((save_effects, ".", i, ".tmp.part"))
progress <- progress + write_end
}
# clean dist_output
dist_output <- dplyr::arrange(dist_output, q_id, iter_id) # sort by q and i to make sure our orders all match up.
dist_output$q_id <-
dist_output$iter_id <-
# save fbm
effect_fbm <- effect_fbm$()
}
((errs$msgs) > 0){
("Errors occured on some ABC iterations. See named element 'errors' in returned list for details.\n")
((ABC_res = dist_output, errors = errs))
}
{
("Complete, no errors on any iterations.\n")
(dist_output)
}
}
}
#' Draw parameter values from the posteriour distribution of an ABC using a density kernal
#' @param x numeric, the number of values to draw.
#' @param res data.frame, the results of an ABC run, must include a "dist" column.
#' @param num_accepted numeric, either the proportion of runs to accept or the number of runs to accept.
#' @param parameters character, the parameter names for which to draw values.
#'
#' @export
<- (x, res, num_accepted, parameters, = 2000, dist.var = "ks"){
((parameters) != 2){
("Only two parameters accepted at the moment.\n")
}
# assign accepted runs
(num_accepted < 1){
qcut <- num_accepted
}
{
qcut <- num_accepted/(res)
}
# remove results with NA values for the dist var
bads <- ((res,dist.var]))
((bads) > 0){
res <- res-bads,]
}
res$hits <- (res,dist.var] <= (res,dist.var], qcut), 1, 0)
# generate kernal
op <- GenKern::KernSur(res[res$hits == 1,((res) == parameters[1])],
res[res$hits == 1,((res) == parameters[2])],
range.x = ((res,((res) == parameters[1])]),
(res,((res) == parameters[1])])),
range.y = ((res,((res) == parameters[2])]),
(res,((res) == parameters[2])])),
ygridsize = , xgridsize = )
# sample from kernal and readjust to rows/columns
cells <- (1:(op$zden), x, = T, prob = op$zden)
rows <- cells %% (op$zden)
cols <- (cells/(op$zden)) + 1
((rows == 0)){
cols[rows == 0] <- cols[rows == 0] - 1
rows[rows == 0] <- (op$zden)
}
# grab parameter values and return
vals <- (p1 = op$xords[rows], p2 = op$yords[cols])
(vals) <- parameters
(vals)
}
