Merge branch 'release-0.1'

This is the first working release of bigsimr with the main function, rvec, working to generate data from specified marginals and parameters.

.Rbuildignore

@@ -1,2 +1,3 @@
 ^bigsimr\.Rproj$
 ^\.Rproj\.user$
+\tests

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: bigsimr
 Title: Simulate Multivariate Negative Binomial Data
-Version: 0.0.0.9000
+Version: 0.1.0
 Authors@R: 
     person(given = "Alexander",
            family = "Knudson",

R/crude_helper_funcs.R

@@ -1,15 +1,15 @@
 #' Generate a random Pearson correlation matrix
-#' 
+#'
 #' @param d A positive integer.
 #' @param constant_rho A boolean value.
 #' @return A \code{dxd} symmetric matrix with ones on the diagonal.
 #' @examples
-#' rcorr(4)
-#' rcorr(10, TRUE)
-rcorr <- function(d, constant_rho = FALSE) {
+#' rcor(4)
+#' rcor(10, TRUE)
+rcor <- function(d, constant_rho = FALSE) {
   if (constant_rho) {
     tmp_rho <- runif(n = 1, min = 0.25, max = 0.9)
-    sigma <- matrix(data = tmp_rho, nrow = d, ncol = d) 
+    sigma <- matrix(data = tmp_rho, nrow = d, ncol = d)
     diag(sigma) <- 1
   } else {
     tmp <- matrix(rnorm(d^2), d, d)
@@ -22,16 +22,16 @@ rcorr <- function(d, constant_rho = FALSE) {
 
 
 #' Generate negative binomial marginal parameters
-#' 
+#'
 #' @param d A positive integer.
 #' @param shape The shape parameter of a negative binomial distribution.
 #' @param id_margins A boolean value specifying if the margins should be identical.
 #' @return A matrix with columns containing (in order) lambda, alpha, mean, and variance.
 #' @examples
-#' rnegbin_params(4)
-#' rnegbin_params(10, shape = 40)
-#' rnegbin_params(7, 44, TRUE)
-rnegbin_params <- function(d, shape = 100, id_margins = FALSE) {
+#' rnbinom_params(4)
+#' rnbinom_params(10, shape = 40)
+#' rnbinom_params(7, 44, TRUE)
+rnbinom_params <- function(d, shape = 100, id_margins = FALSE) {
   if (!id_margins) {
     p <- runif(d, min = 0, max = 0.1)
     lambda <- (1 - p) / p
@@ -49,3 +49,14 @@ rnegbin_params <- function(d, shape = 100, id_margins = FALSE) {
   colnames(to_return) <- paste0("Var", 1:d)
   return(to_return)
 }
+
+
+#' Transforms a [multivariate]normal vector to a different marginal via a
+#' uniform intermediate transformation.
+#'
+#' @param x A normal random vector.
+#' @param param A list ccontaining the marginal and its parameters.
+normal2marginal <- function(x, param) {
+  do.call(what = paste0("q", param[[1]]),
+          args = c(list(p = pnorm(x)), param[-1]))
+}

R/helper_funcs.R

@@ -1,5 +1,5 @@
 #' Converts Pearson correlation to Spearman. Used for normal random variables.
-#' 
+#'
 #' @param R A a square symmetric Pearson correlation matrix.
 #' @return A Spearman correlation matrix.
 convertPearsonSpearman <- function(R) {
@@ -11,7 +11,7 @@ convertPearsonSpearman <- function(R) {
 
 
 #' Converts Spearman correlation to Pearson. Used for normal random variables.
-#' 
+#'
 #' @param R A a square symmetric Spearman correlation matrix.
 #' @return A Pearson correlation matrix.
 convertSpearmanPearson <- function(Rho) {
@@ -23,7 +23,7 @@ convertSpearmanPearson <- function(Rho) {
 
 
 #' Converts Pearson correlation to Kendall. Used for normal random variables.
-#' 
+#'
 #' @param R A a square symmetric Pearson correlation matrix.
 #' @return A Kendall correlation matrix.
 convertPearsonKendall <- function(R) {
@@ -35,7 +35,7 @@ convertPearsonKendall <- function(R) {
 
 
 #' Converts Kendall correlation to Pearson. Used for normal random variables.
-#' 
+#'
 #' @param R A a square symmetric Kendall correlation matrix.
 #' @return A Pearson correlation matrix.
 convertKendallPearson <- function(Rho) {
@@ -47,16 +47,16 @@ convertKendallPearson <- function(Rho) {
 
 
 #' Adjust an input correlation matrix for our method to match
-#' 
+#'
 #' @param R A a square symmetric correlation matrix.
 #' @param params The parameters of the marginals.
 #' @param cores The number of cores to utilize.
 adjustInputR <- function(R, params, cores = 1){
-  
+
   ## 1. find the pairs
   d <- NROW(R)
   index_mat <- combn(x = 1:d, 2)
-  
+
   ## 2. compute the first order approximation
   if (cores == 1) {
     input_rho_vector <- apply(index_mat, 2, function(tmp_index){
@@ -78,9 +78,10 @@ adjustInputR <- function(R, params, cores = 1){
     })
   } else {
     cl <- parallel::makeCluster(cores, type = "FORK")
-    
+
     ## parallel version of the above code
     doParallel::registerDoParallel(cl)
+    `%dopar%` <- foreach::`%dopar%`
     input_rho_vector <- foreach::foreach(i = 1:ncol(index_mat), .combine = 'c') %dopar% {
       ## extract and compute simple quantities
       tmp_index <- index_mat[,i]
@@ -97,25 +98,25 @@ adjustInputR <- function(R, params, cores = 1){
       tmp_nb_rho <- R[tmp_index[1], tmp_index[2]]
       ## convert to rho for gamma (exact)
       target_gamma_rho <- tmp_nb_rho * ((sd_nb1 * sd_nb2) / (scale_factor * lambda_nb1 * lambda_nb2))
-      
+
       ## now approximate rho for input corr for MVN
       return(min(1, target_gamma_rho))
     }
     ## close cluster
     parallel::stopCluster(cl)
   }
-  
+
   ## Put input_rho_vector into a matrix
   to_return <- diag(d)
   to_return[lower.tri(to_return)] <- input_rho_vector
   to_return <- to_return + t(to_return) - diag(diag(to_return))
-  
+
   ## Ensure positive semi-definiteness
   eigen_obj <- eigen(to_return)
   if (any(eigen_obj$values < 0)) {
     to_return <- as.matrix(Matrix::nearPD(to_return)$mat)
   }
-  
+
   ## Label and return
   rownames(to_return) <- colnames(to_return) <- rownames(R)
   return(to_return)

R/rvec.R

@@ -0,0 +1,66 @@
+#' Creates \code{n} observations from a multivariate distribution with the
+#' given marginals and correlation matrix.
+#'
+#' @param n The number random vectors to generate.
+#' @param R The input correlation matrix.
+#' @param params The parameters of the marginals.
+#' @param cores The number of cores to utilize.
+#' @param type The type of correlation matrix that is being passed.
+#' @return A matrix of random vectors generated from the specified marginals and parameters.
+#' @examples
+#' d <- 5
+#' rho <- rcor(d)
+#' margins <- list(
+#'   list("nbinom", size = 5, prob = 0.3),
+#'   list("exp", rate = 4),
+#'   list("binom", size = 5, prob = 0.7),
+#'   list("norm", mean = 10, sd = 3),
+#'   list("pois", lambda = 10)
+#' )
+#'
+#' margins2 <- list(
+#'   list("nbinom", 5, 0.3),
+#'   list("exp", 4),
+#'   list("binom", 5, 0.7),
+#'   list("norm", 10, 3),
+#'   list("pois", 10)
+#' )
+#'
+#' rvec(10, rho, margins, cores = 1, type = "pearson")
+#' rvec(10, rho, margins2, cores = 1, type = "pearson")
+rvec <- function(n, R, params, cores = 1,
+                 type = c("pearson", "kendall", "spearman")){
+  # Handle different types of dependencies
+  if (type == "spearman") {
+    R <- convertSpearmanPearson(R)
+  }
+  if (type == "kendall") {
+    R <- convertKendallPearson(R)
+  }
+
+  # determine the dimension d
+  d <- NROW(R)
+
+  # generate MVN sample
+  mvn_sim <- mvnfast::rmvn(n = n,
+                           mu = rep(0, d),
+                           sigma = R,
+                           ncores = cores,
+                           isChol = FALSE)
+
+  if (cores == 1) {
+    mv_sim <- sapply(1:d, function(i){
+      normal2marginal(mvn_sim[,i], params[[i]])
+    })
+  } else {
+    `%dopar%` <- foreach::`%dopar%`
+    cl <- parallel::makeCluster(cores, type = "FORK")
+    doParallel::registerDoParallel(cl)
+    mv_sim <- foreach::foreach(i = 1:d, .combine = 'cbind') %dopar% {
+      normal2marginal(mvn_sim[,i], params[[i]])
+    }
+    parallel::stopCluster(cl)
+  }
+  colnames(mv_sim) <- rownames(R)
+  return(mv_sim)
+}