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.

Packages

library(holodeck)
library(dplyr)
library(purrr)
library(mice)

Simulating data with `sim_*()`

holodeck provides functions to simulate different kinds of data as columns in a tibble:

sim_cat() for categorical variables
sim_covar() for multivariate normal numeric data
sim_discr() for multivariate normal data with different means for levels of some grouping variable
sim_missing() for randomly introducing NAs

To 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)
#> # A tibble: 10 × 6
#>    Y          V1      V2     V3      V4      V5
#>    <chr>   <dbl>   <dbl>  <dbl>   <dbl>   <dbl>
#>  1 a      0.0292  1.02   -0.187  0.147  -1.01  
#>  2 a     -0.913  -0.729   1.25  -0.578  -0.253 
#>  3 a      1.11    0.709   0.266  1.46    0.145 
#>  4 a     -0.475  -0.154   0.487  0.516  -1.08  
#>  5 a     -1.20   -0.361  -0.130 -0.684  -0.950 
#>  6 b      1.27   -0.290   0.850 -0.0649  0.366 
#>  7 b     -0.632  -0.680  -2.61  -0.300   0.0663
#>  8 b      2.00   -0.711   1.41   0.878   0.226 
#>  9 b     -1.09   -0.0561 -0.228 -0.872  -0.847 
#> 10 b     -0.892   0.873   1.32  -0.251   1.20

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)
#> # A tibble: 10 × 5
#>        V1      V2     V3      V4     V5
#>     <dbl>   <dbl>  <dbl>   <dbl>  <dbl>
#>  1  0.373  0.776  -0.809 -0.0727 -0.375
#>  2  1.37   0.649   1.81   0.716   1.94 
#>  3 -1.47  -0.107  -0.634 -1.08    0.193
#>  4  0.376 -0.387   0.356  0.789   0.151
#>  5  0.845  0.438   0.236  0.966   1.08 
#>  6  1.07  -0.541  -1.38  -0.294  -1.84 
#>  7  0.626  1.70   -0.654  0.759   0.949
#>  8 -0.967 -0.0475 -1.14  -0.659  -0.363
#>  9  1.04   2.01    0.606  1.09    1.76 
#> 10  0.176 -0.462   0.914  0.646  -0.629

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)
#> # A tibble: 10 × 1
#>    group
#>    <chr>
#>  1 a    
#>  2 a    
#>  3 a    
#>  4 a    
#>  5 a    
#>  6 b    
#>  7 b    
#>  8 b    
#>  9 b    
#> 10 b

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))
#> # A tibble: 10 × 6
#> # Groups:   Y [2]
#>    Y         V1     V2      V3      V4     V5
#>    <chr>  <dbl>  <dbl>   <dbl>   <dbl>  <dbl>
#>  1 a     -0.532 -0.182  0.812   1.86    1.63 
#>  2 a      1.24   2.04   1.76    1.28    2.37 
#>  3 a     -1.02   0.443  0.389   1.80   -0.964
#>  4 a      2.24   0.934  0.301   1.22    0.363
#>  5 a      0.901  2.80   0.784   0.0205  1.04 
#>  6 b     -0.960 -2.22  -0.964  -0.633  -2.30 
#>  7 b      0.207 -0.809 -0.0428 -0.675  -1.34 
#>  8 b      0.883 -0.875 -1.45   -2.15   -1.61 
#>  9 b      1.10  -1.32  -0.830   0.236  -0.644
#> 10 b      1.42  -0.448 -1.80   -0.772  -1.99

Chaining functions with `%>%`

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
#> # A tibble: 20 × 16
#>    factor    low_1   low_2   low_3  low_4   low_5  high_1   high_2  high_3
#>    <chr>     <dbl>   <dbl>   <dbl>  <dbl>   <dbl>   <dbl>    <dbl>   <dbl>
#>  1 a       2.02    -1.25    0.246   0.360  0.207  -0.270  -0.394   -0.127 
#>  2 a      -0.414   -0.568   0.819  NA     -0.339   0.754   0.622    0.499 
#>  3 a       0.222   NA       1.12    0.118 -1.41    0.750  -0.0761   0.250 
#>  4 a       0.00774 NA      -0.551   2.74  -0.772   0.292   0.0129   0.615 
#>  5 a       0.651   -0.238  -1.18    0.494 -1.44    0.0888  0.907    0.561 
#>  6 a       0.433    0.927  -0.371  -0.623 -0.865  NA      -0.461   -0.872 
#>  7 a      -0.602   -1.12   -1.74   -2.43  -0.135  -1.39    0.399   -0.496 
#>  8 a       0.211   NA       0.407   0.670 -1.42    2.39    2.62     1.69  
#>  9 a      -0.509   -2.48    0.383  -0.252 -0.168  -0.467  -0.0968  -1.39  
#> 10 a      -0.283   -0.360   1.74    0.189 -0.0845  0.206   0.00510 -0.128 
#> 11 b       0.832   -0.311  -0.289  -0.747 -0.752  -1.12   -1.97    -0.560 
#> 12 b       0.233   -0.649  -0.134  NA     -0.608  -0.389  -0.00557 -0.253 
#> 13 b       0.107   -1.86   -0.488  -0.707  1.10   -0.858  -1.10    -1.20  
#> 14 b      NA       -0.499   0.264  -0.608  1.46   -0.0112 -0.649   -0.116 
#> 15 b      -0.434    0.710  -0.788   1.17  -1.23    1.26    0.283   -0.0719
#> 16 b       0.170   -0.550  -0.551  NA      1.25    1.31    1.87     1.99  
#> 17 b       1.08     0.0630 -0.317   0.237 -0.238   0.308  -0.0503   0.475 
#> 18 b      NA        0.244  -0.0940 -0.435 -0.0298  1.18   -0.00967  1.62  
#> 19 b      -0.925   -0.144   1.79    0.167  0.400  -0.154   0.791   -0.731 
#> 20 b      -1.90    -0.838  NA      -1.21   0.193  -0.254  -0.358   NA     
#> # ℹ 7 more variables: high_4 <dbl>, high_5 <dbl>, discr_1 <dbl>, discr_2 <dbl>,
#> #   discr_3 <dbl>, discr_4 <dbl>, discr_5 <dbl>

Visualizing covariance

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.

Example: Effect of missing data

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.

Create datasets

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
#> # A tibble: 40 × 14
#>    factor noise_1 noise_2  noise_3 signal_1 signal_2 signal_3 signal_4 signal_5
#>    <chr>    <dbl>   <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
#>  1 a       0.264  -0.198   1.04      -1.21  -3.01     -3.06    -4.45    -0.511 
#>  2 a       1.35   -1.40   -0.949     -1.31  -2.72     -0.0499  -2.14    -1.30  
#>  3 a       1.22    0.788   0.647     -0.456 -0.394    -0.155   -1.38    -1.64  
#>  4 a       0.480   0.0463 -0.531     -2.46  -0.811    -1.52    -3.24    -1.50  
#>  5 a      -1.56    1.35    0.604     -2.58  -0.262    -0.882    0.108   -1.31  
#>  6 a      -0.405   0.309  -1.06      -0.656 -0.790    -0.114    0.0135  -0.0347
#>  7 a      -0.0862 -1.11    1.33      -0.659 -1.99      0.165   -2.08    -3.54  
#>  8 a      -0.557   1.06    0.111      0.187 -0.0588   -1.07    -3.08    -0.115 
#>  9 a       0.394  -2.23    0.940     -0.541 -0.182     1.18     0.294   -2.86  
#> 10 a      -0.215  -0.926  -0.00246    0.630 -0.00603  -2.00    -2.32    -0.135 
#> # ℹ 30 more rows
#> # ℹ 5 more variables: signal2_1 <dbl>, signal2_2 <dbl>, signal2_3 <dbl>,
#> #   signal2_4 <dbl>, signal2_5 <dbl>

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.

Simulate missing data

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]])
#> # A tibble: 6 × 14
#>   factor noise_1  noise_2 noise_3 signal_1 signal_2 signal_3 signal_4 signal_5
#>   <chr>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
#> 1 a       1.07   -0.00916 -0.493    -0.953    -1.29    0.952   -0.651   -1.23 
#> 2 a      -0.0134 -1.85    -0.292    -0.849    -1.87    0.311    0.856   -0.785
#> 3 a      -0.401  -1.28    -1.48     -0.432    -1.66   -0.458    0.873   -0.748
#> 4 a      -1.82    2.01    -1.12      0.298    -1.18   -1.05    -0.331   -0.514
#> 5 a       1.71    0.597   -0.0753   -1.00      1.04   -2.65    -0.442   -0.836
#> 6 a       0.848   1.13     0.0428   -1.60     -1.37   -3.24    -0.888   -1.41 
#> # ℹ 5 more variables: signal2_1 <dbl>, signal2_2 <dbl>, signal2_3 <dbl>,
#> #   signal2_4 <dbl>, signal2_5 <dbl>
head(dfs.missing[[1]])
#> # A tibble: 6 × 14
#>   factor noise_1  noise_2 noise_3 signal_1 signal_2 signal_3 signal_4 signal_5
#>   <chr>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
#> 1 a       1.07   -0.00916 NA        -0.953    -1.29    0.952   -0.651   -1.23 
#> 2 a      -0.0134 -1.85    -0.292    NA        -1.87    0.311    0.856   -0.785
#> 3 a      NA      -1.28    -1.48     -0.432    -1.66   -0.458    0.873   -0.748
#> 4 a      -1.82    2.01    -1.12      0.298    NA      -1.05    -0.331   -0.514
#> 5 a       1.71    0.597   -0.0753   -1.00      1.04   -2.65    -0.442   -0.836
#> 6 a       0.848   1.13     0.0428   -1.60     -1.37   -3.24    -0.888   -1.41 
#> # ℹ 5 more variables: signal2_1 <dbl>, signal2_2 <dbl>, signal2_3 <dbl>,
#> #   signal2_4 <dbl>, signal2_5 <dbl>
head(dfs.imputed[[1]])
#>   factor    noise_1      noise_2     noise_3   signal_1  signal_2   signal_3
#> 1      a  1.0738532 -0.009162203 -0.07529851 -0.9525532 -1.289870  0.9523219
#> 2      a -0.0134100 -1.851241745 -0.29246857  0.4403128 -1.873263  0.3114562
#> 3      a -0.7215366 -1.284047854 -1.48162821 -0.4322049 -1.656319 -0.4583891
#> 4      a -1.8176185  2.014018985 -1.11967562  0.2977722 -1.394962 -1.0478207
#> 5      a  1.7095029  0.596707289 -0.07529851 -0.9996540  1.039059 -2.6513107
#> 6      a  0.8475474  1.130402793  0.04280330 -1.6014609 -1.370313 -3.2403584
#>     signal_4   signal_5  signal2_1  signal2_2  signal2_3   signal2_4  signal2_5
#> 1 -0.6506898 -1.2266449 -0.9398164  0.7444634 -1.4554263 -0.06208422 -1.0852542
#> 2  0.8561365 -0.7850617  0.2655353 -0.1112986  0.9284334  1.70032445 -0.9680104
#> 3  0.8725502 -0.7483738  2.4863183  0.8013671  0.0386806 -1.23665052  0.9490844
#> 4 -0.3310259 -0.5142761 -0.2309898  0.7444634  0.3386661  1.58683283 -2.8872108
#> 5 -0.4418016 -0.8359875  1.0516597 -0.2346468  0.9491426 -0.14501237 -2.8872108
#> 6 -0.8882406 -1.4122971  1.5272074  1.7493532  2.0554579 -0.13655839 -1.8517622

Simulating Multivariate Data with holodeck