knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
To simulate co-varying multivariate normal data one can use the mvrnorm()
function from the MASS
package. However, this requires input of a covariance matrix and returns a matrix. The holodeck
package provides functions that are "tidy" in the sense that they work with dataframes and the pipe operator (%>%
). It also includes some functions that provide an interface between the Bioconductor package ropls
and the tidyverse.
library(holodeck) library(dplyr) library(purrr) library(mice)
sim_*()
holodeck
provides functions to simulate different kinds of data as columns in a tibble:
sim_cat()
for categorical variablessim_covar()
for multivariate normal numeric datasim_discr()
for multivariate normal data with different means for levels of some grouping variablesim_missing()
for randomly introducing NA
sTo simulate multivariate data you need to start with a dataframe or a tibble. Once you have a dataframe or tibble, the sim_*()
functions add columns onto it.
df <- tibble(Y = rep(c("a", "b"), each = 5)) df %>% sim_covar(n_vars = 5, var = 1, cov = 0.5)
Optionally you can create a tibble with the sim_covar()
or sim_cat()
functions by providing them with the N
argument instead of .data
.
sim_covar(n_obs = 10, n_vars = 5, var = 1, cov = 0.5)
sim_cat()
is a rather simple wrapper that just creates a column of categorical data. Eventually, it will be expanded to allow creation of crossed and nested factors.
sim_cat(n_obs = 10, n_groups = 2)
sim_discr()
simulates covarying data that differs in means between levels of some grouping variable.
df %>% group_by(Y) %>% sim_discr(n_vars = 5, var = 1, cov = 0.1, group_means = c(1, -1))
%>%
One advantage of the holodeck
package is the ability to chain functions together to create complex data covariance structures. You can chain functions together in any order, although the sim_discr()
function requires a grouping variable.
All of the sim_*
functions (besides sim_missing()
) take an optional name argument which names the variables created.
df <- sim_covar(n_obs = 20, n_vars = 5, var = 1, cov = 0.1, name = "low") %>% #5 variables with low covariance sim_covar(n_vars = 5, var = 1, cov = 0.8, name = "high") #5 variables with high covariance
Now we could add a categorical variable, and some variables that discriminate between levels of our categorical variable
df1 <- df %>% sim_cat(n_groups = 2, name = "factor") %>% group_by(factor) %>% sim_discr(n_vars = 5, var = 1, cov = 0.1, group_means = c(-1, 1), name = "discr") %>% ungroup()
Finally, if you want to simulate missing values, you can use sim_missing()
to randomly introduce NAs.
df2 <- df1 %>% sim_missing(prop = 0.1) df2
cov()
creates a covariance matrix with variance on the diagonal. We can visualize it as a heatmap.
library(ggplot2) df1 %>% select(-factor) %>% cov() %>% heatmap(Rowv = NA, Colv = NA, symm = TRUE)
Values are higher for the discriminating variables because the cov
and var
arguments to sim_discr()
only control the covariance and variance within groups.
One reason to simulate multivariate data is to test the effects of different properties of datasets on analysis results. For example, what's the effect of missing data on a statistical analysis? The sim_missing()
function replaces a proportion of values with NA. Let's see how it affects a PLS-DA analysis.
We can chain several sim_*
functions to quickly create a dataframe.
df2 <- sim_cat(n_obs = 40, n_groups = 3, name = "factor") %>% sim_covar(n_vars = 3, var = 1, cov = 0.0, name = "noise") %>% group_by(factor) %>% sim_discr(n_vars = 5, var = 1, cov = 0, group_means = c(-1, 0, 1), name = "signal") %>% sim_discr(n_vars = 5, var = 1, cov = 0, group_means = c(0, 0.5, 1), name = "signal2") %>% ungroup() df2
We can then use map()
from the purrr
package to create many randomly generated datasets using the same specifications, with and without missing values.
set.seed(100) dfs <- map(1:20, ~sim_cat(n_obs = 40, n_groups = 3, name = "factor") %>% sim_covar(n_vars = 3, var = 1, cov = 0.0, name = "noise") %>% group_by(factor) %>% sim_discr(n_vars = 5, var = 1, cov = 0, group_means = c(-1, 0, 1), name = "signal") %>% sim_discr(n_vars = 5, var = 1, cov = 0, group_means = c(0, 0.5, 1), name = "signal2") %>% ungroup())
Alternatively, you could generate one large dataframe (many rows) and take subsets. Either way, you know the "true" properties of the data and can compare to the results of the analyses you test.
We can now map the sim_missing()
function to randomly introduce NAs to the datasets.
set.seed(101) dfs.missing <- map(dfs, ~sim_missing(., prop = 0.05))
And finally, deal with those NAs with multiple imputation with the mice
package.
# this might take a few seconds dfs.imputed <- map(dfs.missing, ~mice(., printFlag = FALSE) %>% complete())
Here, we can compare an example dataset as original, with NAs, and imputed:
head(dfs[[1]]) head(dfs.missing[[1]]) head(dfs.imputed[[1]])
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.