R/rstan-wrapper.R

In amarder/powerstan: Power Analysis using Stan

#' Create a random data frame.
#'
#' @param model The Stan model describing how the data is generated.
#' @param parameters The assumed parameter values to use in the data
#'     generation process.
#' @param n The number of samples to make.
#' @param data Any data that needs to be fed into the model - things
#'     like number of explanatory variables.
#'
#' @importFrom rstan sampling extract
#' @importFrom utils capture.output
#' @export
#'
simulate_data <- function(model, parameters, n, data = list()) {
    capture.output(
        fit <- sampling(
            object = model,
            data = data,
            chains = 1,
            iter = n,
            warmup = 0,
            init = function(i) parameters,
            algorithm = "Fixed_param"
        )
    )
    data <- as.data.frame(extract(fit))
    return(data)
}

#' Simulate a data set `n` times, return the proportion of simulations
#' where the null hypothesis is rejected.
#'
#' @param data_generator A function to simulate a data set.
#' @param test_null A function that tests the null hypothesis given an
#'     observed data set.
#' @param n The number of simulations to run.
#'
#' @export
#'
calculate_power <- function(data_generator, test_null, n = 100) {
    nrejected <- 0
    for (i in 1:n) {
        data <- data_generator()
        stopifnot(inherits(data, "data.frame"))
        rejected <- test_null(data)
        stopifnot(rejected %in% c(FALSE, TRUE))
        if (rejected) {
            nrejected <- nrejected + 1
        }
    }
    return(nrejected / n)
}

#' @importFrom rstan stan_model
#' @export
#'
rstan::stan_model