diff --git a/CRAN-RELEASE b/CRAN-RELEASE
deleted file mode 100644
index 4ba470f..0000000
--- a/CRAN-RELEASE
+++ /dev/null
@@ -1,2 +0,0 @@
-This package was submitted to CRAN on 2019-09-27.
-Once it is accepted, delete this file and tag the release (commit c1541e7383).
diff --git a/DESCRIPTION b/DESCRIPTION
index aa80393..c2e24c8 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: colorednoise
 Type: Package
 Title: Simulate Temporally Autocorrelated Populations
-Version: 1.0.5
-Date: 2019-09-26
+Version: 1.1.0
+Date: 2020-08-03
 Authors@R: c(
     person("Julia", "Pilowsky", email = "jap2178@caa.columbia.edu", 
     role = c("aut", "cre"), comment = c(ORCID = "0000-0002-6376-2585")
@@ -14,12 +14,11 @@ License: GPL-3
 Depends: R (>= 3.3.0)
 Imports:
   stats (>= 3.3.2),
-  dplyr (>= 0.7.3),
   purrr (>= 0.2.3),
-  tibble (>= 2.0.0),
-  tidyr (>= 1.0.0)
+  Rcpp (>= 1.0.5),
+  data.table (>= 1.12.8)
 LinkingTo: Rcpp, RcppArmadillo
-RoxygenNote: 6.1.1
+RoxygenNote: 7.1.1
 Encoding: UTF-8
 LazyData: true
 BugReports: http://github.com/japilo/colorednoise/issues
diff --git a/NAMESPACE b/NAMESPACE
index 0439521..ccd4145 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,14 +7,13 @@ export(colored_noise)
 export(cor2cov)
 export(matrix_model)
 export(multi_rnorm)
+export(stdev_transform)
 export(unstructured_pop)
-import(dplyr)
+import(data.table, except = transpose)
 import(purrr)
-import(tidyr)
+importFrom(Rcpp,evalCpp)
 importFrom(stats,acf)
 importFrom(stats,na.omit)
 importFrom(stats,plogis)
 importFrom(stats,sd)
-importFrom(tibble,as_tibble)
-importFrom(tibble,tibble)
 useDynLib(colorednoise, .registration = TRUE)
diff --git a/NEWS.md b/NEWS.md
index 2ccb6c7..b6835e4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,10 +1,15 @@
+## colorednoise 1.1.0
+
+* New function exported: `stdev_transform`, which adapts standard deviations to different probability distributions.
+* colorednoise now runs on `data.table` instead of `dplyr`, `tibble`, and `tidyr` for faster simulations.
+
 ## colorednoise 1.0.5
 
-* Updated to be compatible with tidyr v1.0.0
+* Updated to be compatible with `tidyr` v1.0.0
 
 ## colorednoise 1.0.4
 
-* Updated to be compatible with tibble v2.0.0
+* Updated to be compatible with `tibble` v2.0.0
 
 ## colorednoise 1.0.3
 
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 6229c60..abe163a 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -68,18 +68,41 @@ cor2cov <- function(sigma, corrMatrix) {
 #' @param covMatrix A valid covariance matrix. The number of rows/columns must match the length of the mu, sigma, and phi vectors.
 #' @return A matrix with as many rows as timesteps and as many columns as mu/sigma/phi values.
 #' @examples
-#' cov <- matrix(c(0.037, 0.044, -0.048, 0.044, 0.247, -0.008, -0.047, -0.008, 0.074), nrow = 3)
+#' cov <- matrix(c(1, 0.53, 0.73, 0.53, 1, 0.44, 0.73, 0.44, 1), nrow = 3)
 #' test <- colored_multi_rnorm(100, c(0, 3, 5), c(1, 0.5, 1), c(0.5, -0.3, 0), cov)
 #' var(test)
-#' library(dplyr)
-#' test %>% as.data.frame() %>% summarize_all(.funs = c("mean", "sd", "autocorrelation"))
+#' library(data.table)
+#' as.data.table(test)[, .(V1_mean = mean(V1), V2_mean = mean(V2), V3_mean = mean(V3),
+#' V1_sd = sd(V1), V2_sd = sd(V2), V3_sd = sd(V3),
+#' V1_autocorrelation = autocorrelation(V1), V2_autocorrelation = autocorrelation(V2),
+#' V3_autocorrelation = autocorrelation(V3))]
 #' @export
 colored_multi_rnorm <- function(timesteps, mean, sd, phi, covMatrix) {
     .Call(`_colorednoise_colored_multi_rnorm`, timesteps, mean, sd, phi, covMatrix)
 }
 
-variancefix <- function(mu, sigma, dist) {
-    .Call(`_colorednoise_variancefix`, mu, sigma, dist)
+#' Translate Standard Deviation from the Natural Scale to the Log or Logit Scale
+#'
+#' This function changes a given standard deviation so that when a vector of samples is drawn from the given distribution,
+#' the original standard deviation will be recovered once it is back-transformed from the log or logit scale. In effect,
+#' the function "translates" a standard deviation from the natural scale to the log or logit scale for the purposes of
+#' random draws from a probability distribution.
+#' @param mu The mean of the distribution on the natural scale.
+#' @param sigma The standard devation of the distribution on the natural scale.
+#' @param dist The distribution to which the standard deviation should be transformed.
+#' @return The standard deviation translated to the log or logit scale.
+#' @examples
+#' mean <- 10
+#' stdev <- 2
+#' mean_trans <- log(mean)
+#' stdev_trans <- stdev_transform(mean, stdev, "log")
+#' draws <- rnorm(50, mean_trans, stdev_trans)
+#' natural_scale <- exp(draws)
+#' mean(draws)
+#' sd(draws)
+#' @export
+stdev_transform <- function(mu, sigma, dist) {
+    .Call(`_colorednoise_stdev_transform`, mu, sigma, dist)
 }
 
 #' Simulated Time Series of an Unstructured Temporally Autocorrelated Population
diff --git a/R/colorednoise.R b/R/colorednoise.R
index bb5ded6..a9e1b46 100644
--- a/R/colorednoise.R
+++ b/R/colorednoise.R
@@ -8,13 +8,12 @@
 #' @name colorednoise
 #' @useDynLib colorednoise, .registration = TRUE
 #' @import purrr
-#' @import dplyr
+#' @importFrom Rcpp evalCpp
+#' @rawNamespace import(data.table, except = transpose)
 #' @importFrom stats sd acf na.omit plogis
-#' @importFrom tibble tibble as_tibble
-#' @import tidyr
 NULL
 
 ## quiets concerns of R CMD check re: the .'s that appear in
 ## pipelines
 if (getRversion() >= "2.15.1") utils::globalVariables(c(".", "mean.trans",
-    "sd.trans", "noise", "timestep", "dist", "zero"))
+    "sd.trans", "noise", "timestep", "dist", "zero", "ref"))
diff --git a/R/simulate_populations.R b/R/simulate_populations.R
index 18f6815..eaf9e2f 100644
--- a/R/simulate_populations.R
+++ b/R/simulate_populations.R
@@ -57,8 +57,8 @@ autocorr_sim <- function(timesteps, start, survPhi, fecundPhi, survMean,
         }
     }
     # Unnests the list and adds estimates of survival and fertility
-    sims <- labeled_sims %>% flatten() %>% map(~mutate(., est_surv = survivors/population,
-        est_fecund = newborns/survivors))
+    sims <- labeled_sims %>% flatten() %>% map(setDT) %>%
+      map(~.[, `:=`(est_surv = survivors/population, est_fecund = newborns/survivors)])
     return(sims)
 }
 
@@ -117,10 +117,10 @@ matrix_model <- function(data, initialPop, timesteps, covMatrix = NULL,
     colNames = NULL, matrixStructure = NULL, repeatElements = NULL,
     survivalOverflow = "scale") {
     stages <- length(initialPop)
-    # Regularize all valid data inputs to the same format
+    # Regularize all valid data input to the same format
     if (is.data.frame(data) == T) {
         if (is.null(colNames) == F) {
-            data <- data[, colNames] %>% rename(!(!(!colNames)))
+            data <- data %>% setDT() %>% setnames(old = colNames, new = names(colNames))
         }
         if (all(names(data) == c("mean", "sd", "autocorrelation")) ==
             F) {
@@ -135,7 +135,7 @@ matrix_model <- function(data, initialPop, timesteps, covMatrix = NULL,
         if (all(data$sd > 0) == F) {
             stop("Invalid values in SD column")
         }
-        dat <- data %>% as_tibble()
+        dat <- setDT(data)
     } else if (is.list(data) == T) {
         if (length(data) > 3) {
             stop("List data should only have 3 elements")
@@ -149,7 +149,7 @@ matrix_model <- function(data, initialPop, timesteps, covMatrix = NULL,
         if (length(unique(map_int(data, length)))>1) {
           stop("Matrices are not equal dimensions")
         }
-        dat <- tibble(mean = as.vector(t(data[[1]])), sd = as.vector(t(data[[2]])),
+        dat <- data.table(mean = as.vector(t(data[[1]])), sd = as.vector(t(data[[2]])),
             autocorrelation = as.vector(t(data[[3]])))
     } else {
         stop("Invalid data type. Must be a list of three matrices or a data frame with three columns.")
@@ -166,31 +166,35 @@ matrix_model <- function(data, initialPop, timesteps, covMatrix = NULL,
     if (is.null(repeatElements) == T) {
       repeatElements <- matrix(seq(1:stages^2), ncol = stages, byrow = T)
     }
-    dat <- dat %>% mutate(dist = dists, zero = mean==0&sd==0,
-                          ref = as.vector(t(repeatElements)))
+    dat <- dat[, `:=`(dist = dists, zero = mean == 0 & sd == 0, ref = as.vector(t(repeatElements)))]
     repeats <- repeatElements == matrix(seq(1:stages^2), ncol = stages, byrow = T)
     # Create version of data that can be used to generate colored noise
-    inputs <- dat %>% slice(which(t(repeats))) %>% filter(zero == F) %>%
-      rowwise() %>% mutate(
-        mean.trans = ifelse(mean == 0, 0, invoke(dist, list(mean))),
-        sd.trans = ifelse(sd == 0, 0, variancefix(mean, sd, dist))
-      )
+    unique <- dat[which(t(repeats))][zero == F]
+    unique$mean.trans <- map_dbl(c(1:nrow(unique)), function(x) {
+      if(dat[x,]$mean == 0) {0} else {
+        invoke(dat[x,]$dist, list(dat[x,]$mean))
+        }
+    })
+    unique$sd.trans <- map_dbl(c(1:nrow(unique)), function(x) {
+      if(dat[x,]$sd == 0) {0} else {
+        stdev_transform(dat[x,]$mean, dat[x,]$sd, dat[x,]$dist)
+      }
+    })
     if (is.null(covMatrix) == T) {
-      covMatrix <- cor2cov(inputs$sd, diag(nrow(inputs)))
+      covMatrix <- cor2cov(unique$sd, diag(nrow(unique)))
     }
     # Create colored noise, discard if invalid matrix
     if(survivalOverflow == "redraw") {
       repeat {
-        inputs$noise <- colored_multi_rnorm(timesteps, inputs$mean.trans, inputs$sd.trans,
-                                            inputs$autocorrelation, covMatrix) %>% split(col(.))
-        result <- left_join(dat, inputs, by = c("mean", "sd", "autocorrelation", "dist", "zero", "ref"))
+        unique$noise <- colored_multi_rnorm(timesteps, unique$mean.trans, unique$sd.trans,
+                                            unique$autocorrelation, covMatrix) %>% split(col(.))
+        result <- dat[unique, on = .(mean, sd, autocorrelation, dist, zero, ref), allow.cartesian = TRUE]
         result$noise[dat$zero==T] <- rep(list(rep.int(0, timesteps)), sum(dat$zero==T))
         # checking for >1 probability
-        result <- result %>% rowwise() %>% mutate(
-          natural.noise = ifelse(zero == T, list(noise),
-                                 ifelse(dist == "log", list(exp(noise)),
-                                        list(plogis(noise))))
-        )
+        result$natural.noise <- map(c(1:nrow(result)), function(x) {
+          if (result[x,]$zero) {result[x,]$noise} else if (result[x,]$dist=="log") {
+            exp(result[x,]$noise[[1]])} else {plogis(result[x,]$noise[[1]])}
+          })
         matrices <- map(1:timesteps, function(x) {
           matrix(map_dbl(result$natural.noise, x), byrow = T, ncol = stages)
         })
@@ -201,15 +205,15 @@ matrix_model <- function(data, initialPop, timesteps, covMatrix = NULL,
         }
       }
     } else if(survivalOverflow == "scale") {
-      inputs$noise <- colored_multi_rnorm(timesteps, inputs$mean.trans, inputs$sd.trans,
-                                          inputs$autocorrelation, covMatrix) %>% split(col(.))
-      result <- left_join(dat, inputs, by = c("mean", "sd", "autocorrelation", "dist", "zero", "ref"))
+      unique$noise <- colored_multi_rnorm(timesteps, unique$mean.trans, unique$sd.trans,
+                                          unique$autocorrelation, covMatrix) %>% split(col(.))
+      result <- dat[unique, on = .(mean, sd, autocorrelation, dist, zero, ref), allow.cartesian = TRUE]
       result$noise[dat$zero==T] <- rep(list(rep.int(0, timesteps)), sum(dat$zero==T))
       # checking for >1 probability
-      result <- result %>% rowwise() %>% mutate(
-        natural.noise = ifelse(zero == T, list(noise),
-                               ifelse(dist == "log", list(exp(noise)),
-                                      list(plogis(noise)))))
+      result$natural.noise <- map(c(1:nrow(result)), function(x) {
+        if (result[x,]$zero) {result[x,]$noise} else if (result[x,]$dist=="log") {
+          exp(result[x,]$noise[[1]])} else {plogis(result[x,]$noise[[1]])}
+      })
       matrices <- map(1:timesteps, function(x) {
         matrix(map_dbl(result$natural.noise, x), byrow = T, ncol = stages)
       })
@@ -226,7 +230,9 @@ matrix_model <- function(data, initialPop, timesteps, covMatrix = NULL,
       }
     } else {stop("survivalOverflow must be set to 'redraw' or 'scale'")}
     pop <- projection(initialPop, matrices)
-    pop %>% map(as_tibble, .name_repair = ~ c(paste0("stage", 1:stages))) %>%
-      bind_rows() %>% group_by(timestep = row_number()) %>% nest(data = -timestep) %>%
-      mutate(total = map_dbl(data, sum)) %>% unnest(data)
+    # CONVERT TO DATA.TABLE
+    output <- pop %>% map(as.data.table) %>% rbindlist() %>%
+      .[, `:=`(total = rowSums(.SD), timestep = seq_len(.N))] %>% setcolorder("timestep")
+    names(output) <- c("timestep", map_chr(1:stages, ~paste0("stage", .)), "total")
+    return(output)
 }
diff --git a/README.Rmd b/README.Rmd
index 035f166..2a563f4 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -11,13 +11,15 @@ knitr::opts_chunk$set(
   fig.path = "man/figures/README-"
 )
 library(ggplot2)
+library(Rcpp)
 ```
 
-# colorednoise <img src='man/figures/logo.png' align="right" height="125" />
+# colorednoise <img src='man/figures/hex.png' align="right" height="125" />
 [![Travis_build_status](https://travis-ci.org/japilo/colorednoise.svg?branch=master)](https://travis-ci.org/japilo/colorednoise)
 [![CRAN_version](https://www.r-pkg.org/badges/version/colorednoise)](https://cran.r-project.org/package=colorednoise)
 [![Coverage Status](https://img.shields.io/codecov/c/github/japilo/colorednoise/master.svg)](https://codecov.io/github/japilo/colorednoise?branch=master)
 [![Download_count](https://cranlogs.r-pkg.org/badges/grand-total/colorednoise)](https://CRAN.R-project.org/package=colorednoise)
+[![Paper_doi](https://img.shields.io/badge/doi-10.1111/oik.06438-orange.svg)](https://doi.org/10.1111/oik.06438)
 
 ## Overview
 
diff --git a/README.md b/README.md
index 56b8c56..3d9f363 100644
--- a/README.md
+++ b/README.md
@@ -1,13 +1,14 @@
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
-# colorednoise <img src='man/figures/logo.png' align="right" height="125" />
+# colorednoise <img src='man/figures/hex.png' align="right" height="125" />
 
 [![Travis\_build\_status](https://travis-ci.org/japilo/colorednoise.svg?branch=master)](https://travis-ci.org/japilo/colorednoise)
 [![CRAN\_version](https://www.r-pkg.org/badges/version/colorednoise)](https://cran.r-project.org/package=colorednoise)
 [![Coverage
 Status](https://img.shields.io/codecov/c/github/japilo/colorednoise/master.svg)](https://codecov.io/github/japilo/colorednoise?branch=master)
 [![Download\_count](https://cranlogs.r-pkg.org/badges/grand-total/colorednoise)](https://CRAN.R-project.org/package=colorednoise)
+[![Paper\_doi](https://img.shields.io/badge/doi-10.1111/oik.06438-orange.svg)](https://doi.org/10.1111/oik.06438)
 
 ## Overview
 
diff --git a/cran-comments.md b/cran-comments.md
index dd8d2a8..39bc97e 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,7 +1,7 @@
 ## Test environments
-* local OS X install, R 3.6
-* Ubuntu 16.04 (on travis-ci), R 3.5, 3.6, devel
-* win-builder, R devel
+* local Windows 10 install, R 4.0
+* Ubuntu 16.04 (on travis-ci), R 3.6, 4.0, devel
+* macOS 10.13 (on travis-ci), R 3.6, 4.0, devel
 
 ## R CMD check results
 
diff --git a/man/autocorr_sim.Rd b/man/autocorr_sim.Rd
index 442e1f3..6838eda 100644
--- a/man/autocorr_sim.Rd
+++ b/man/autocorr_sim.Rd
@@ -5,8 +5,17 @@
 \title{Simulate Temporally Autocorrelated Populations for Every Combination of
 Parameters}
 \usage{
-autocorr_sim(timesteps, start, survPhi, fecundPhi, survMean, survSd,
-  fecundMean, fecundSd, replicates)
+autocorr_sim(
+  timesteps,
+  start,
+  survPhi,
+  fecundPhi,
+  survMean,
+  survSd,
+  fecundMean,
+  fecundSd,
+  replicates
+)
 }
 \arguments{
 \item{timesteps}{The number of timesteps you want to simulate. Individuals
diff --git a/man/colored_multi_rnorm.Rd b/man/colored_multi_rnorm.Rd
index d04433b..3666c0f 100644
--- a/man/colored_multi_rnorm.Rd
+++ b/man/colored_multi_rnorm.Rd
@@ -25,9 +25,12 @@ A matrix with as many rows as timesteps and as many columns as mu/sigma/phi valu
 Generates random variables that are correlated to each other and temporally autocorrelated.
 }
 \examples{
-cov <- matrix(c(0.037, 0.044, -0.048, 0.044, 0.247, -0.008, -0.047, -0.008, 0.074), nrow = 3)
+cov <- matrix(c(1, 0.53, 0.73, 0.53, 1, 0.44, 0.73, 0.44, 1), nrow = 3)
 test <- colored_multi_rnorm(100, c(0, 3, 5), c(1, 0.5, 1), c(0.5, -0.3, 0), cov)
 var(test)
-library(dplyr)
-test \%>\% as.data.frame() \%>\% summarize_all(.funs = c("mean", "sd", "autocorrelation"))
+library(data.table)
+as.data.table(test)[, .(V1_mean = mean(V1), V2_mean = mean(V2), V3_mean = mean(V3),
+V1_sd = sd(V1), V2_sd = sd(V2), V3_sd = sd(V3),
+V1_autocorrelation = autocorrelation(V1), V2_autocorrelation = autocorrelation(V2),
+V3_autocorrelation = autocorrelation(V3))]
 }
diff --git a/man/colorednoise.Rd b/man/colorednoise.Rd
index e3e1f7c..e787019 100644
--- a/man/colorednoise.Rd
+++ b/man/colorednoise.Rd
@@ -3,7 +3,6 @@
 \docType{package}
 \name{colorednoise}
 \alias{colorednoise}
-\alias{colorednoise-package}
 \title{\code{colorednoise} package}
 \description{
 Simulate Temporally Autocorrelated Populations
diff --git a/man/figures/README-example-1.png b/man/figures/README-example-1.png
index fe73f2a..b629885 100644
Binary files a/man/figures/README-example-1.png and b/man/figures/README-example-1.png differ
diff --git a/man/figures/README-example-2.png b/man/figures/README-example-2.png
index d3f7159..4d3a3b5 100644
Binary files a/man/figures/README-example-2.png and b/man/figures/README-example-2.png differ
diff --git a/man/figures/hex.png b/man/figures/hex.png
new file mode 100644
index 0000000..d9cc1db
Binary files /dev/null and b/man/figures/hex.png differ
diff --git a/man/figures/logo.png b/man/figures/logo.png
deleted file mode 100644
index 4bec35e..0000000
Binary files a/man/figures/logo.png and /dev/null differ
diff --git a/man/matrix_model.Rd b/man/matrix_model.Rd
index 725ee91..f9fc49e 100644
--- a/man/matrix_model.Rd
+++ b/man/matrix_model.Rd
@@ -4,9 +4,16 @@
 \alias{matrix_model}
 \title{Temporally Autocorrelated Matrix Population Models}
 \usage{
-matrix_model(data, initialPop, timesteps, covMatrix = NULL,
-  colNames = NULL, matrixStructure = NULL, repeatElements = NULL,
-  survivalOverflow = "scale")
+matrix_model(
+  data,
+  initialPop,
+  timesteps,
+  covMatrix = NULL,
+  colNames = NULL,
+  matrixStructure = NULL,
+  repeatElements = NULL,
+  survivalOverflow = "scale"
+)
 }
 \arguments{
 \item{data}{The input data can be one of two formats: a list of three matrices, or a data frame
diff --git a/man/stdev_transform.Rd b/man/stdev_transform.Rd
new file mode 100644
index 0000000..6a7a1fe
--- /dev/null
+++ b/man/stdev_transform.Rd
@@ -0,0 +1,34 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{stdev_transform}
+\alias{stdev_transform}
+\title{Translate Standard Deviation from the Natural Scale to the Log or Logit Scale}
+\usage{
+stdev_transform(mu, sigma, dist)
+}
+\arguments{
+\item{mu}{The mean of the distribution on the natural scale.}
+
+\item{sigma}{The standard devation of the distribution on the natural scale.}
+
+\item{dist}{The distribution to which the standard deviation should be transformed.}
+}
+\value{
+The standard deviation translated to the log or logit scale.
+}
+\description{
+This function changes a given standard deviation so that when a vector of samples is drawn from the given distribution,
+the original standard deviation will be recovered once it is back-transformed from the log or logit scale. In effect,
+the function "translates" a standard deviation from the natural scale to the log or logit scale for the purposes of
+random draws from a probability distribution.
+}
+\examples{
+mean <- 10
+stdev <- 2
+mean_trans <- log(mean)
+stdev_trans <- stdev_transform(mean, stdev, "log")
+draws <- rnorm(50, mean_trans, stdev_trans)
+natural_scale <- exp(draws)
+mean(draws)
+sd(draws)
+}
diff --git a/man/unstructured_pop.Rd b/man/unstructured_pop.Rd
index 448052b..8920246 100644
--- a/man/unstructured_pop.Rd
+++ b/man/unstructured_pop.Rd
@@ -4,8 +4,16 @@
 \alias{unstructured_pop}
 \title{Simulated Time Series of an Unstructured Temporally Autocorrelated Population}
 \usage{
-unstructured_pop(start, timesteps, survPhi, fecundPhi, survMean, survSd,
-  fecundMean, fecundSd)
+unstructured_pop(
+  start,
+  timesteps,
+  survPhi,
+  fecundPhi,
+  survMean,
+  survSd,
+  fecundMean,
+  fecundSd
+)
 }
 \arguments{
 \item{start}{The starting population size.}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index c9787a8..6759c61 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -60,16 +60,16 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// variancefix
-double variancefix(double mu, double sigma, std::string dist);
-RcppExport SEXP _colorednoise_variancefix(SEXP muSEXP, SEXP sigmaSEXP, SEXP distSEXP) {
+// stdev_transform
+double stdev_transform(double mu, double sigma, std::string dist);
+RcppExport SEXP _colorednoise_stdev_transform(SEXP muSEXP, SEXP sigmaSEXP, SEXP distSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< double >::type mu(muSEXP);
     Rcpp::traits::input_parameter< double >::type sigma(sigmaSEXP);
     Rcpp::traits::input_parameter< std::string >::type dist(distSEXP);
-    rcpp_result_gen = Rcpp::wrap(variancefix(mu, sigma, dist));
+    rcpp_result_gen = Rcpp::wrap(stdev_transform(mu, sigma, dist));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -109,7 +109,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_colorednoise_multi_rnorm", (DL_FUNC) &_colorednoise_multi_rnorm, 3},
     {"_colorednoise_cor2cov", (DL_FUNC) &_colorednoise_cor2cov, 2},
     {"_colorednoise_colored_multi_rnorm", (DL_FUNC) &_colorednoise_colored_multi_rnorm, 5},
-    {"_colorednoise_variancefix", (DL_FUNC) &_colorednoise_variancefix, 3},
+    {"_colorednoise_stdev_transform", (DL_FUNC) &_colorednoise_stdev_transform, 3},
     {"_colorednoise_unstructured_pop", (DL_FUNC) &_colorednoise_unstructured_pop, 8},
     {"_colorednoise_projection", (DL_FUNC) &_colorednoise_projection, 2},
     {NULL, NULL, 0}
diff --git a/src/noise.cpp b/src/noise.cpp
index 60837b0..6ab90e2 100644
--- a/src/noise.cpp
+++ b/src/noise.cpp
@@ -85,11 +85,14 @@ arma::mat cor2cov(Rcpp::NumericVector sigma, Rcpp::NumericMatrix corrMatrix) {
 //' @param covMatrix A valid covariance matrix. The number of rows/columns must match the length of the mu, sigma, and phi vectors.
 //' @return A matrix with as many rows as timesteps and as many columns as mu/sigma/phi values.
 //' @examples
-//' cov <- matrix(c(0.037, 0.044, -0.048, 0.044, 0.247, -0.008, -0.047, -0.008, 0.074), nrow = 3)
+//' cov <- matrix(c(1, 0.53, 0.73, 0.53, 1, 0.44, 0.73, 0.44, 1), nrow = 3)
 //' test <- colored_multi_rnorm(100, c(0, 3, 5), c(1, 0.5, 1), c(0.5, -0.3, 0), cov)
 //' var(test)
-//' library(dplyr)
-//' test %>% as.data.frame() %>% summarize_all(.funs = c("mean", "sd", "autocorrelation"))
+//' library(data.table)
+//' as.data.table(test)[, .(V1_mean = mean(V1), V2_mean = mean(V2), V3_mean = mean(V3),
+//' V1_sd = sd(V1), V2_sd = sd(V2), V3_sd = sd(V3),
+//' V1_autocorrelation = autocorrelation(V1), V2_autocorrelation = autocorrelation(V2),
+//' V3_autocorrelation = autocorrelation(V3))]
 //' @export
 // [[Rcpp::export]]
 Rcpp::NumericMatrix colored_multi_rnorm(int timesteps, Rcpp::NumericVector mean, Rcpp::NumericVector sd, Rcpp::NumericVector phi, Rcpp::NumericMatrix covMatrix) {
diff --git a/src/simulation.cpp b/src/simulation.cpp
index dfcdf33..1e5e6e1 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -2,12 +2,28 @@
 #include "noise.h"
 // [[Rcpp::depends(RcppArmadillo)]]
 
-// Variance Fix
-//
-// This function changes the variance so that once it is back-transformed
-// from the logit scale, the original variance is recovered.
+//' Translate Standard Deviation from the Natural Scale to the Log or Logit Scale
+//'
+//' This function changes a given standard deviation so that when a vector of samples is drawn from the given distribution,
+//' the original standard deviation will be recovered once it is back-transformed from the log or logit scale. In effect,
+//' the function "translates" a standard deviation from the natural scale to the log or logit scale for the purposes of
+//' random draws from a probability distribution.
+//' @param mu The mean of the distribution on the natural scale.
+//' @param sigma The standard devation of the distribution on the natural scale.
+//' @param dist The distribution to which the standard deviation should be transformed.
+//' @return The standard deviation translated to the log or logit scale.
+//' @examples
+//' mean <- 10
+//' stdev <- 2
+//' mean_trans <- log(mean)
+//' stdev_trans <- stdev_transform(mean, stdev, "log")
+//' draws <- rnorm(50, mean_trans, stdev_trans)
+//' natural_scale <- exp(draws)
+//' mean(draws)
+//' sd(draws)
+//' @export
 // [[Rcpp::export]]
-double variancefix(double mu, double sigma, std::string dist){
+double stdev_transform(double mu, double sigma, std::string dist){
   if (sigma == 0) {
     return 0;
   }
@@ -68,11 +84,11 @@ double variancefix(double mu, double sigma, std::string dist){
 Rcpp::DataFrame unstructured_pop(int start, int timesteps, double survPhi, double fecundPhi, double survMean, double survSd, double fecundMean, double fecundSd) {
   // These lines generate the temporally autocorrelated random numbers
   double survmu = R::qlogis(survMean, 0, 1, true, false);
-  double survsigma = variancefix(survMean, survSd, "qlogis");
+  double survsigma = stdev_transform(survMean, survSd, "qlogis");
   Rcpp::NumericVector survnoise = colored_noise(timesteps, survmu, survsigma, survPhi);
   Rcpp::NumericVector St = plogis(survnoise);
   double fecundmu = log(fecundMean);
-  double fecundsigma = variancefix(fecundMean, fecundSd, "log");
+  double fecundsigma = stdev_transform(fecundMean, fecundSd, "log");
   Rcpp::NumericVector fecundnoise = colored_noise(timesteps, fecundmu, fecundsigma, fecundPhi);
   Rcpp::NumericVector Ft = exp(fecundnoise);
   // This loop kills some individuals according to St probabilities, creates
diff --git a/tests/testthat/test_noise.R b/tests/testthat/test_noise.R
index 353e632..a6f8da2 100644
--- a/tests/testthat/test_noise.R
+++ b/tests/testthat/test_noise.R
@@ -1,6 +1,6 @@
 library(colorednoise)
 library(purrr)
-library(dplyr)
+library(data.table)
 context("Autocorrelation of colored noise")
 
 test_that("colored noise can produce blue noise", {
@@ -27,7 +27,9 @@ test_that("colored_multi_rnorm can produce red noise", {
   set.seed(989)
   corr <- matrix(c(1, 0.53, 0.73, 0.53, 1, 0.44, 0.73, 0.44, 1), nrow = 3)
   test <- colored_multi_rnorm(100, c(0, 3, 5), c(1, 0.5, 1), c(0.5, 0.5, 0.5), corr) %>%
-    as_tibble() %>% summarise_all(autocorrelation) %>% as.numeric()
+    as.data.table() %>%
+    .[, .(V1_autocorr = autocorrelation(V1), V2_autocorr = autocorrelation(V2), V3_autocorr = autocorrelation(V3))] %>%
+    as.numeric()
   expect_true(all(test > 0.5)==T)
 })
 
@@ -35,6 +37,8 @@ test_that("colored_multi_rnorm can produce blue noise", {
   set.seed(19045)
   corr <- matrix(c(1, 0.53, 0.73, 0.53, 1, 0.44, 0.73, 0.44, 1), nrow = 3)
   test <- colored_multi_rnorm(100, c(0, 3, 5), c(1, 0.5, 1), c(-0.5, -0.5, -0.5), corr) %>%
-    as_tibble() %>% summarise_all(autocorrelation) %>% as.numeric()
+    as.data.table() %>%
+    .[, .(V1_autocorr = autocorrelation(V1), V2_autocorr = autocorrelation(V2), V3_autocorr = autocorrelation(V3))] %>%
+    as.numeric()
   expect_true(all(test < -0.4)==T)
 })
diff --git a/tests/testthat/test_simulations.R b/tests/testthat/test_simulations.R
index acb027f..76b96e4 100644
--- a/tests/testthat/test_simulations.R
+++ b/tests/testthat/test_simulations.R
@@ -1,20 +1,22 @@
 library(purrr)
-library(dplyr)
+library(data.table)
 context("Consistency in population simulations")
 
 test_that("unstructured_pop can produce blue noise populations", {
     test_blue <- rerun(.n = 1000, unstructured_pop(start = 5000, timesteps = 50,
         survPhi = -0.5, survMean = 0.4, survSd = 0.05, fecundPhi = 0,
-        fecundMean = 1.5, fecundSd = 0.2)) %>% map(~mutate(., est_surv = survivors/population,
-        est_fecund = newborns/survivors)) %>% map_dbl(~autocorrelation(.$est_surv))
+        fecundMean = 1.5, fecundSd = 0.2)) %>% map(setDT) %>%
+      map(~.[, `:=`(est_surv = survivors/population, est_fecund = newborns/survivors)]) %>%
+      map_dbl(~autocorrelation(.$est_surv))
     expect_true(mean(test_blue) < -0.2)
 })
 
 test_that("unstructured_pop can produce red noise populations", {
     test_red <- rerun(.n = 1000, unstructured_pop(start = 5000, timesteps = 50,
         survPhi = 0.5, survMean = 0.4, survSd = 0.05, fecundPhi = 0,
-        fecundMean = 1.5, fecundSd = 0.2)) %>% map(~mutate(., est_surv = survivors/population,
-        est_fecund = newborns/survivors)) %>% map_dbl(~autocorrelation(.$est_surv))
+        fecundMean = 1.5, fecundSd = 0.2)) %>% map(setDT) %>%
+      map(~.[, `:=`(est_surv = survivors/population, est_fecund = newborns/survivors)]) %>%
+      map_dbl(~autocorrelation(.$est_surv))
     expect_true(mean(test_red) > 0.2)
 })
 
@@ -33,11 +35,11 @@ test_that("matrix_model can produce cross-correlated autocorrelated populations
   expect_true(sum(test[nrow(test), 2:3]) > sum(test[1, 2:3]))
 })
 
-test_that("output of variancefix transforms back to the natural scale", {
+test_that("output of stdev_transform transforms back to the natural scale", {
   mu <- 0.5
   sigma <- 0.2
   mu.log <- log(mu)
-  sigma.log <- variancefix(mu, sigma, "log")
+  sigma.log <- stdev_transform(mu, sigma, "log")
   rand <- rnorm(1000, mu.log, sigma.log) %>% exp()
   expect_true(sd(rand)>=sigma && sd(rand)<0.25)
 })
diff --git a/vignettes/noise.Rmd b/vignettes/noise.Rmd
index 8c6b860..3182f1c 100644
--- a/vignettes/noise.Rmd
+++ b/vignettes/noise.Rmd
@@ -15,9 +15,10 @@ knitr::opts_chunk$set(
   comment = "#>"
 )
 library("colorednoise")
-library("dplyr")
 library("ggplot2")
 library("purrr")
+library("data.table")
+library("Rcpp")
 ```
 
 ## What is colored noise? 
@@ -65,9 +66,9 @@ The estimates for these seem to be quite accurate. But let's try generating a wh
 
 ```{r generate replicates}
 raw_sims <- cross_df(list(timesteps = 20, phi = c(-0.5, 0, 0.5), mean = 0.3, sd = c(0.5, 0.7, 0.9, 1.1, 1.3, 1.5))) %>% rerun(.n = 30, pmap(., colored_noise))
-labels <- raw_sims %>% map(1) %>% bind_rows()
+labels <- raw_sims %>% map(1) %>% rbindlist()
 noise <- raw_sims %>% map(2) %>% flatten()
-estimates <- data.frame(mu = map_dbl(noise, mean), sigma = map_dbl(noise, sd), autocorrelation = map_dbl(noise, autocorrelation))
+estimates <- data.table(mu = map_dbl(noise, mean), sigma = map_dbl(noise, sd), autocorrelation = map_dbl(noise, autocorrelation))
 sd_range <- cbind(labels, estimates)
 head(sd_range)
 ```
@@ -75,9 +76,12 @@ head(sd_range)
 We can visualize these results by finding the mean and confidence intervals of the autocorrelation estimates for each combination of variance and noise color, and plotting them out.
 
 ```{r plot CIs}
-sd_range %>% group_by(phi, sd) %>% summarize_at(funs(lower.ci = ((function(bar){quantile(bar, probs=c(0.05, 0.95))[[1]]})(.)),
-           upper.ci = ((function(bar){quantile(bar, probs=c(0.05, 0.95))[[2]]})(.)),
-           mean = mean(.)), .vars = vars(autocorrelation)) -> summ
+summ <- sd_range[, .(lower.ci = ((function(bar) {
+    quantile(bar, probs = c(0.05, 0.95))[[1]]
+})(autocorrelation)), upper.ci = ((function(bar) {
+    quantile(bar, probs = c(0.05, 0.95))[[2]]
+})(autocorrelation)), mean = mean(autocorrelation)), keyby = .(phi, 
+    sd)]
 ggplot(summ, aes(x = sd, y = mean)) +
   geom_pointrange(aes(ymin = lower.ci, ymax = upper.ci, color = factor(phi)), size = 0.8) + 
   geom_hline(yintercept = 0, linetype = 2, color = "#C0C0C0") +
@@ -114,15 +118,17 @@ Here we see one example of a population simulated by the function.
 Now let's estimate the autocorrelation in each of these simulated populations. There is a utility function `autocorrelation` in this package that measures temporal autocorrelation at a time lag of 1.
 
 ```{r estimate autocorrelation}
-sims %>% map(~group_by(., survPhi, timesteps)) %>% map(~summarize(., acf.surv = autocorrelation(est_surv, biasCorrection = F))) %>% bind_rows -> estimates
+estimates <- sims %>% map(~.[, .(acf.surv = autocorrelation(est_surv, biasCorrection = F)), keyby = .(survPhi, timesteps)]) %>% rbindlist()
 ```
 
 I find that this format of ggplot displays the estimates by noise color very nicely.
 
 ```{r plotting estimates, fig.height=5, fig.width=8}
-estimates %>% group_by(survPhi, timesteps) %>% summarize_at(funs(lower.ci = ((function(bar){quantile(bar, probs=c(0.05, 0.95), na.rm = T)[[1]]})(.)),
-           upper.ci = ((function(bar){quantile(bar, probs=c(0.05, 0.95), na.rm = T)[[2]]})(.)),
-           mean = mean(., na.rm = T)), .vars = vars(acf.surv)) -> summ2
+summ2 <- estimates[, .(lower.ci = ((function(bar) {
+    quantile(bar, probs = c(0.05, 0.95))[[1]]
+})(acf.surv)), upper.ci = ((function(bar) {
+    quantile(bar, probs = c(0.05, 0.95))[[2]]
+})(acf.surv)), mean = mean(acf.surv)), keyby = .(survPhi, timesteps)]
 # Noise color values for the graph in hexadecimal
 noise = c("#0033FF", "#3366FF", "#6699FF", "#99CCFF", "#FFFFFF", "#FF9999", "#FF6666","#FF0000", "#CC0000")
 ggplot(summ2, aes(x=timesteps, y=mean)) + geom_point(size=8, aes(color=factor(survPhi))) +
@@ -145,9 +151,11 @@ where \phi is the autocorrelation and T is the length of the time series.
 Here is what the autocorrelation estimates look like when I turn on the bias correction in the autocorrelation function.
 
 ```{r bias correction, message = F}
-sims %>% map(~group_by(., survPhi, timesteps)) %>% map(~summarize(., acf.surv = autocorrelation(est_surv, biasCorrection = TRUE))) %>% bind_rows() %>% group_by(survPhi, timesteps) %>% summarize_at(funs(lower.ci = ((function(bar){quantile(bar, probs=c(0.05, 0.95), na.rm = T)[[1]]})(.)),
-           upper.ci = ((function(bar){quantile(bar, probs=c(0.05, 0.95), na.rm = T)[[2]]})(.)),
-           mean = mean(., na.rm = T)), .vars = vars(acf.surv)) -> summ3
+summ3 <- sims %>% map(~.[, .(acf.surv = autocorrelation(est_surv, biasCorrection = TRUE)), keyby = .(survPhi, timesteps)]) %>% rbindlist() %>% .[, .(lower.ci = ((function(bar) {
+    quantile(bar, probs = c(0.05, 0.95), na.rm=T)[[1]]
+})(acf.surv)), upper.ci = ((function(bar) {
+    quantile(bar, probs = c(0.05, 0.95), na.rm=T)[[2]]
+})(acf.surv)), mean = mean(acf.surv)), keyby = .(survPhi, timesteps)]
 ggplot(summ3, aes(x=timesteps, y=mean)) + 
   geom_point(size=8, aes(color=factor(survPhi))) +
   geom_pointrange(aes(ymin=lower.ci, ymax=upper.ci), size=0.8) +