/
interface_abc.R
94 lines (90 loc) · 3.36 KB
/
interface_abc.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#' Convert Simulation Results to abc's Parameter Format
#'
#' This function creates an object compatible with the \code{param}
#' argument of the \code{\link[abc]{abc}} function from coala's simulation
#' results.
#'
#' @param sim_results The simulation results as returned from
#' \code{\link[=simulate.coalmodel]{simulate}}.
#' @param model The model used for the simulations.
#'
#' @return A data.frame that can be used as \code{param}
#' argument of \code{\link[abc]{abc}}.
#' @export
#' @seealso For generating abc's summary statistics format:
#' \code{\link{create_abc_sumstat}}
#'
#' @examples
#' model <- coal_model(10, 1) +
#' feat_mutation(par_prior("theta", rnorm(1, 5, .5))) +
#' sumstat_sfs()
#' sim_results <- simulate(model, nsim = 2)
#' abc_param <- create_abc_param(sim_results, model)
#' print(abc_param)
create_abc_param <- function(sim_results, model) {
# First take care for special case when nsim was 1
if (!is.null(sim_results$pars)) {
sim_results <- list(sim_results)
}
# Now create parmeter data.frame
do.call(rbind, lapply(sim_results, function(x) {
data.frame(t(x$pars))
}))
}
#' Convert Simulation Results to abc's Summary Statistic Format
#'
#' This function creates an object compatible with the \code{sumstat}
#' argument of the \code{\link[abc]{abc}} function from coala's simulation
#' results. It converts all summary statistics that are in the simulation
#' results and expects that each of them is a numeric vector.
#' Use transformation functions to convert none vector-valued statistics
#' (e.g. \code{\link{sumstat_jsfs}}, \code{\link{sumstat_omega}} or
#' \code{\link{sumstat_trees}}) into a vector.
#'
#' @inheritParams create_abc_param
#'
#' @return A data.frame that can be used as \code{sumstat}
#' argument of \code{\link[abc]{abc}}.
#' @export
#' @seealso For generating abc's parameter format:
#' \code{\link{create_abc_param}}
#'
#' @examples
#' # Using the SFS:
#' model <- coal_model(10, 1) +
#' feat_mutation(par_prior("theta", rnorm(1, 5, .5))) +
#' sumstat_sfs()
#' sim_results <- simulate(model, nsim = 2)
#' abc_sumstat <- create_abc_sumstat(sim_results, model)
#' print(abc_sumstat)
#'
#' # Using the JSFS and converting it into a vector:
#' model <- coal_model(c(10, 10), 1) +
#' feat_mutation(par_prior("theta", rnorm(1, 5, .5))) +
#' feat_migration(par_prior("m", rnorm(1, .5, .1)), symmetri = TRUE) +
#' sumstat_jsfs(transformation = function(jsfs) {
#' c(sum(jsfs[1, ]), sum(jsfs[, 1]), sum(jsfs[-1, -1]))
#' })
#' sim_results <- simulate(model, nsim = 2)
#' abc_sumstat <- create_abc_sumstat(sim_results, model)
#' print(abc_sumstat)
create_abc_sumstat <- function(sim_results, model) {
# First take care for special case when nsim was 1
if (!is.null(sim_results$pars)) {
sim_results <- list(sim_results)
}
sumstats <- get_summary_statistics(model)
do.call(rbind, lapply(sim_results, function(x) {
stats_combined <- do.call(c, lapply(sumstats, function(stat) {
value <- x[[stat$get_name()]]
assert_that(is.numeric(value))
if (!is.vector(value)) {
warning("Value of the ", stat$get_name(), " statistic is not a vector. ",
"Converting it to vector, but this may not be correct. ",
"Better use transformation functions to convert the ",
"statistics to vectors.")
}
value
}))
}))
}