R/simulations.R

In StatsIMUWr/slobeC: What the package does (short line)

#' Run a single simulation
#'
#' @param num_observations Number of observation in a data matrix
#' @param correlation_matrix Object of type "cor_desc" specifying correlation
#' structure
#' @param signal Object of class "signal" specifying the beta vector.
#' @param FDR FDR level for SLOPE
#' @param tolerance Tolerance to be used in admm algorithm
#' @param metric named list of metrics to be calculated
#'
#' @return List of simulation parameters and simulation results.
#'
#' @export
#'

single_simulation <- function(num_observations, correlation_matrix, signal,
                              FDR = 0.1, tolerance = 1e-4,
                              metrics = list("mse" = mse, "tp" = tp, "fp" = fp,
                                             "pred_err" = pred_error)) {

    correlation_matrix <- make_correlation_matrix(
        correlation_matrix[["type"]],
        correlation_matrix[["correlation"]],
        signal[["dimension"]]
    )
    num_variables <- signal[["dimension"]]
    beta <- signal[["beta"]]
    X <- anMC::mvrnormArma(num_observations, mu = rep(0, num_variables),
                       sigma = correlation_matrix[["matrix"]], 0)
    X <- X / sqrt(num_observations)
    Y <- X %*% beta + rnorm(num_observations)
    lasso_fit <- glmnet::cv.glmnet(X,Y)
    initial_beta <- as.numeric(coef(lasso_fit))[-1]
    initial_lambda <- qnorm(1 - FDR*(1:num_variables)/(2 * num_variables))

    slob_fit <- SLOB(initial_beta, Y, X, initial_lambda,
                     a = .01*num_observations,
                     b = .1*num_observations)[["beta"]]
    slope_fit <- SLOPE::SLOPE(X, Y, FDR = FDR, normalize = FALSE)[["beta"]]
    slob_cpp_fit <- slobeC::SLOBE_ADMM_approx(initial_beta, X, Y,
                                              a_prior = 0.01*num_observations,
                                              b_prior = 0.1*num_observations,
                                              FDR = FDR, verbose = F,
                                              tol = tolerance)[["beta"]]
    fits <- list("slob" = slob_fit, "slope" = slope_fit,
                 "slob_cpp" = slob_cpp_fit)

    simulation_results <- lapply(names(fits), function(fitted) {
        lapply(names(metrics), function(metric) {
            list(model = fitted,
                 statistic = metric,
                 value = metrics[[metric]](beta, fits[[fitted]], X))
        })
    })

    list(simulation_params = list(n_variables = num_variables,
                                  n_observations = num_observations,
                                  cor_type = correlation_matrix[["type"]],
                                  cor_val = correlation_matrix[["correlation"]],
                                  signal = signal[["signal"]],
                                  n_nonzero = sum(!dplyr::near(beta, 0)),
                                  FDR = FDR),
         simulation_results = unlist(simulation_results,
                                     recursive = FALSE, use.names = FALSE))
}


#' Multiple runs of the simulation for a single set of parameters
#'
#' @param num_iterations Number of replications
#' @inheritParams single_simulation
#'
#' @return A data frame of results
#'
#' @export
#'

slob_simulation <- function(num_iterations, num_observations,
                            correlation_matrix, signal,
                            FDR = 0.1, tolerance = 1e-4,
                            metrics = list("mse" = mse, "tp" = tp, "fp" = fp,
                                           "pred_err" = pred_error)) {
    raw_results <- lapply(1:num_iterations, function(iteration) {
        c(list(iteration = iteration),
          single_simulation(num_observations, correlation_matrix, signal,
                            FDR, tolerance, metrics))
    })
    results <- lapply(raw_results, function(result) {
        df <- dplyr::bind_rows(lapply(result[["simulation_results"]],
                                      as.data.frame))
        dplyr::mutate(df, iteration = result[["iteration"]])
    })
    results <- dplyr::bind_rows(results)
    params <- as.data.frame(raw_results[[1]][["simulation_params"]])

    dplyr::bind_cols(results, params[rep(1, nrow(results)), ])
}

#' Compare execution times
#'
#' @param n_rep number of replication for microbenchmark function
#' @param observation_numbers A list of observation numbers
#' @param correlation_matrices A list of correlation matrices
#' @param signals A list of signals
#' @param FDRs A list of FDR values
#'
#' @return data.frame
#'
#' @export
#'

execution_time <- function(n_rep, observation_numbers,
                           correlation_matrices, signals, FDRs) {
    all_results <- lapply(observation_numbers, function(n_obs) {
        lapply(signals, function(signal) {
            lapply(correlation_matrices, function(cor_m) {
                lapply(FDRs, function(fdr) {
                    cor_m <- make_correlation_matrix(
                        cor_m[["type"]],
                        cor_m[["correlation"]],
                        signal[["dimension"]]
                    )

                    num_variables <- signal[["dimension"]]
                    beta <- signal[["beta"]]
                    X <- anMC::mvrnormArma(n_obs, mu = rep(0, num_variables),
                                       sigma = cor_m[["matrix"]], 0)
                    X <- X / sqrt(n_obs)
                    Y <- X %*% beta + rnorm(n_obs)
                    lasso_fit <- glmnet::cv.glmnet(X,Y)
                    initial_beta <- as.numeric(coef(lasso_fit))[-1]
                    initial_lambda <- qnorm(1 - fdr*(1:num_variables)/(2 * num_variables))

                    time_results <- microbenchmark::microbenchmark(
                        slope = SLOPE::SLOPE(X, Y, FDR = fdr, normalize = FALSE),
                        slobe = SLOB(initial_beta, Y, X, initial_lambda,
                                     a = .01*n_obs,
                                     b = .1*n_obs),
                        slobe_cpp = slobeC::SLOBE_ADMM_approx(initial_beta, X, Y,
                                                              a_prior = 0.01*n_obs,
                                                              b_prior = 0.1*n_obs,
                                                              FDR = fdr, verbose = F),
                        times = n_rep
                    )
                    dplyr::bind_cols(time_results,
                                     data.frame(n_variables = num_variables,
                                                n_observations = n_obs,
                                                cor_type = cor_m[["type"]],
                                                cor_val = cor_m[["correlation"]],
                                                signal = signal[["signal"]],
                                                n_nonzero = sum(!dplyr::near(beta, 0)),
                                                FDR = fdr)[rep(1, nrow(time_results)), ])
                })
            })
        })
    })
    time_result_df <- dplyr::bind_rows(
        lapply(
            unlist(
                unlist(
                    unlist(
                        all_results,
                        FALSE, FALSE),
                    FALSE, FALSE),
                FALSE, FALSE),
            as.data.frame))
    time_result_df <- dplyr::mutate(time_result_df, statistic = "time")
    colnames(time_result_df)[1:2] <- c("model", "value")
    time_result_df
}


#' Run simulations for many sets of parameters
#'
#' @param num_iterations Number of runs for a single set of parameters
#' @param observation_numbers A list of observation numbers
#' @param correlation_matrices A list of correlation matrices
#' @param signals A list of signals
#' @param FDRs A list of FDR values
#'
#' @return data frame
#'
#' @export
#'

simulations <- function(num_iterations, observation_numbers,
                        correlation_matrices, signals, FDRs) {
    all_results <- lapply(observation_numbers, function(n_obs) {
        lapply(signals, function(signal) {
            lapply(correlation_matrices, function(cor_m) {
                lapply(FDRs, function(fdr) {
                    suppressMessages(slob_simulation(num_iterations, n_obs,
                                                     cor_m, signal, fdr))
                })
            })
        })
    })
    dplyr::bind_rows(unlist(
        unlist(
            unlist(all_results, FALSE, FALSE),
            FALSE, FALSE),
        FALSE, FALSE))
}