first-style=long-short ATE short = ATE, long = average treatment effect ATT
short = ATT, long = average treatment effect of the treated CART short = CART,
long = Classification And Regression Tree DAG short = DAG, long = directed
acyclic graph IPTW short = IPTW, long = Inverse Probability of Treatment
Weighting GBSG short = GBSG, long = German Breast Cancer Study Group MFP short
= MFP, long = Multivariable Fractional Polynomial MOB short = MOB, long =
Model-based Recursive Partitioning for Subgroup Analysis MSE short = MSE, long
= Mean Squared Error palmtree short = palmtree, long = partially additive
(generalized) linear model tree RFS short = RFS, long = Relapse-free survival
PI short = PI, long = permutation importance

# Statistical Plasmode Simulations - Potentials, Challenges and
Recommendations

Nicholas Schreck
Biostatistics Division
German Cancer Research Center
D-69120 Heidelberg, GERMANY
<EMAIL_ADDRESS>&Alla Slynko11footnotemark: 1
Department of Statistics and Actuarial Science
University of Waterloo
Ontario, N2L 3G1, CANADA
&Maral Saadati
Freelance Statistician
Saadati Solutions
D-68526 Ladenburg, GERMANY
&Axel Benner
Biostatistics Division
German Cancer Research Center
D-69120 Heidelberg, GERMANY
Both authors contributed equally.

###### Abstract

Statistical data simulation is essential in the development of statistical
models and methods as well as in their performance evaluation. To capture
complex data structures, in particular for high-dimensional data, a variety of
simulation approaches have been introduced including parametric and the so-
called plasmode simulations. While there are concerns about the realism of
parametrically simulated data, it is widely claimed that plasmodes come very
close to reality with some aspects of the “truth” known. However, there are no
explicit guidelines or state-of-the-art on how to perform plasmode data
simulations. In the present paper, we first review existing literature and
introduce the concept of statistical plasmode simulation. We then discuss
advantages and challenges of statistical plasmodes and provide a step-wise
procedure for their generation, including key steps to their implementation
and reporting. Finally, we illustrate the concept of statistical plasmodes as
well as the proposed plasmode generation procedure by means of a public real
RNA dataset on breast carcinoma patients.

##  1 Introduction

Data availability is a crucial issue that arises in the context of statistical
model development and validation, inference derivation, introduction of
statistical concepts and many others 1, 2, 3. In some cases, especially for
high-dimensional data (HDD) where the number of features is substantially
larger than the number of observations, the number of data samples and data
sets is not large enough to properly and reliably perform all required tasks.
To overcome that deficiency, alternative data generation approaches such as
the generation of artificial data are required. The generated data should
match as close as possible the real-life data underlying the research question
of interest, in particular with respect to its probabilistic structure as well
as possible dependencies.

A common approach for data generation is statistical data simulation. In this
paper, we interpret data simulation as a data generation procedure that
follows a data-generating process (DGP), with marginal distributions and a
dependence structure as its basic components. For explanatory and prediction
models being in focus of the present paper, data generation procedures should
necessarily include some steps on the generation of outcome variable(s). To
this end, outcome-generating models (OGMs) are usually utilized to generate
the outcomes using the available covariate information. Parameters of the OGM
are to be estimated from real data, taken from literature or just set based on
investigator’s choice. Examples for outcome data generation can be found in
Reeb and Steibel 4, Franklin et al. 5, Schulz et al. 6, Atiquzzaman et al. 7
and many others. Shmueli 8 depicts the application of OGMs both for explaining
and prediction purposes.

For simulated data, the "truth" is assumed to be known a priori, at least to a
some extent 9. That "truth" can be represented by simulation parameters or
prespecified effect sizes, and is used to reliably evaluate the obtained
estimates or predictions.

The most detailed practical introduction to simulation studies that includes
structural approaches for their planning and reporting as well as the
discussion on the appropriate performance measures is presented in Morris et
al. 2. In particular, the authors introduce a coherent terminology on
simulation studies, data-generating mechanisms, and provide guidance on coding
simulation studies. Most prominently, they introduce the Aims, Data-generating
mechanisms, Estimands, Methods, and Performance (ADEMP) criteria as a guidance
for the planning and performing of simulation studies.

For those researchers who use the results of simulation studies without being
familiar with the entire simulation process, the discussion presented in
Boulesteix et al. 10 can be of great assistance. The paper not only contains
many useful examples and applications, but also describes basic principles of
simulation studies, gives insights into sampling variability and data
generating processes, and demonstrates the role of statistical simulations in
health research.

Two common approaches for data generation are parametric and plasmode
simulation, with the first approach being the most extensively studied and
widely used one. Parametric simulations assume that the parametric stochastic
model used to generate data is realistic and representative, with parameters
of interest estimated from real data, derived from the literature or even set
up by the user, in order to model specific scenarios 1, 2. Plasmode data
generation usually begins by resampling covariate information from the
original real data 11, 12. External (parametric) "truth" such as effect sizes
or model parameter values in explanatory or prediction models can then be
added to the covariable data sets to define the relationship between the
covariables and the outcome. With the OGM being a part of the plasmode data
generation procedure, the resulting plasmode simulation can then be viewed as
a semi-parametric data generation procedure. Of note, parametrically simulated
data may often be considered to be purely artificial, whereas plasmode data is
claimed to reflect reality in the most close way 9, 11, 13.

In the present paper, we provide an extensive literature overview of
parametric and plasmode simulations. In particular, we discuss the difference
between biological and statistical plasmodes. We address advantages and
challenges of parametric and statistical plasmode simulation approaches in
various contexts and provide step-by-step recommendations for the generation
of statistical plasmodes.

At this point, we have no intention to demonstrate the superiority of a
plasmode simulation over a parametric simulation or vice versa. We aim to
analyze advantages and challenges of both data generation methods, and to
illustrate the usefulness of plasmode simulations as a complement to and
possible extension of parametric simulation studies.

The paper is organized as follows: Section 2 discusses parametric and plasmode
simulations, compares their characteristics and provides an extensive
literature review on both data generation methods. Section 3 analyzes
statistical plasmodes in more detail by discussing their challenges. Section 4
provides recommendations for planning, performing and reporting of statistical
plasmode simulations. In Section 5, we present a numerical example to
illustrate the application of such a data generation approach. Section 6
concludes with a discussion.

##  2 From Parametric to Plasmode Simulation Studies

This section provides a comparative introduction to parametric and plasmode
simulations. In particular, we start with a description of the main properties
of parametric simulations as one of the most well-established types of data
generation. Then we move on to plasmode simulations which are often claimed to
be a more close-to-reality approach for data generation.

###  2.1 Parametric Data Simulation

In cases where an underlying DGP is defined in closed form and represented by
a parametric stochastic model, we speak of parametric simulations.

The main asset of parametric simulations is their flexibility in terms of the
chosen DGP. That is, one can easily generate a variety of independent data
sets by varying the assumptions imposed on the DGP and its crucial parameters.
As a result, the simulated data sets can cover many complex but relevant
scenarios as well as extreme situations that do not reflect reality. This
feature becomes particularly important whenever we aim to analyze the
behaviour and performance of different statistical methods. In addition, the
knowledge of the DGP makes the corresponding parametric simulation more
transparent and plausible.

Parametric simulations depreciate the sample size issue as an unlimited amount
of data can be generated by means of a particular DGP. For instance, when
applied to simulation of continuous random variables, an “infinite” number of
distinct data points can be generated without much effort.

The corresponding DGP represents a cornerstone of any parametric simulation.
However, the existence of an appropriate model for the DGP that best fits some
underlying real data cannot be taken for granted. On the other hand, even with
a DGP available, the conclusions based on parametric simulations might be
limited or even biased by the parameters of the chosen DGP model.

Obviously, the quality of a parametric simulation depends on the level of our
comprehension of the underlying processes, distributions and possible
dependencies. For instance, Vaughan et al. 13 state that simulations might not
fully reflect the complexity of the biological data that originates from
nonrandom mating, recombination, hot spots, and other genetic mechanisms.
Boulesteix et al. 10 provides a similar statement and claims that many
simulation studies are too simplified to describe the complexity of the real
life data and thus may lead to inaccurate or even misleading findings.

When generating new data, we strive to preserve not only the marginal
distributions but also the underlying dependence structure. In this context,
the specification of appropriate dependence metrics, such as correlation,
emerge. One of the options for such a specification is estimation of the
dependence structure from the data at hand. Computationally, such an
estimation can be very expensive and time-consuming, especially in case of
large data sets. For parametric simulations, this computational issue is one
of possible reasons for making independence assumptions on certain variables
that then leads to a block diagonal structure for the corresponding
correlation matrix 11. Furthermore, a multivariate normal distribution
provides the simplest model for multivariate covariate data with a pre-
specified mean and correlation structure 1. For generating non-normal
correlated data, diverse copulas or the extended Fleischman power method can
be utilized 14. Obviously, certain concerns about the accuracy and realism of
the underlying modelling assumptions emerge for all such approaches.

One of the assumptions imposed in the context of parametric simulations is
that the "truth" must be known a priori. Such an assumption may not hold in
both low- and high-dimensional situations. However, high-dimensionality may
even exacerbate this issue making parametric simulations inapplicable. For
instance, the “truth” about the set of biological markers truly associated
with a given outcome may be unknown 15.

Within the framework of a parametric simulation, the underlying dependence
structure has to be completely specified in advance. The specification of such
a dependence structure for high-dimensional data may become a challenge.
Possible reasons can be not only computational efficiency issue, but also
spurious correlations 16, sparsity of the data, some nonlinear or even hidden
dependencies and other issues related to large covariance matrices; for more
discussion and examples see Fan and Li 17, Johnstone and Titterington 18 and
the references therein. Pitfalls in the specification of the dependence
structure may then lead to false research discoveries and incorrect
statistical inferences.

Dimensionality reduction and feature extraction play pivotal roles and are
often fundamental in many high-dimensional settings 17. However, it is not
obvious how a parametric simulation may impact the findings of those
procedures considering possible non-representativeness of parametrically
generated covariate data sets.

Altogether, in many cases parametric simulations may turn out to be infeasible
for high-dimensional data generation as it is not obvious how such simulations
would cope with features of high-dimensionality.

A number of papers share our concerns on the applicability of parametric
simulations for the generation of high-dimensional data sets. For instance,
Gadbury et al. 12 question the applicability of standard simulations performed
in a high-dimensional experiment where hundreds of hypotheses are to be
tested. Also Franklin et al. 5 sees the application of ordinary simulation
methods as an issue when comparing high-dimensional variable selection
strategies. In particular, the authors point out that the performance of those
strategies depend "[…] on the information richness and complexity of the
underlying empirical data source" and it is doubtful whether a parametric
simulation is able to capture the richness and complexity of the information.

To overcome, at least partially, the limitations of parametric simulations
outlined above, plasmode simulations have been introduced 11, 13, 5.

###  2.2 Plasmode Data Simulation

The term "plasmode" has been first introduced in Cattell and Jaspers 19, with
a plasmode data set defined as " […] a set of numerical values fitting a
mathematico-theoretical model". In their seminal work, the authors emphasized
that the certainty for the produced plasmode data set to fit the model comes
either because there is a real life experiment producing the data of that kind
or because the simulated data is produced mathematically to fit the functions.
Two different approaches for plasmode generation, performed either in a lab
experiment or by resampling, were also mentioned in Mehta et al. 11. In the
present paper, we refine the discussions provided by these authors, introduce
the concept of statistical plasmodes and analyze their properties. From our
perspective, the classification of plasmodes in biological and statistical
depends on the procedure used for their generation.

Biological plasmodes are those generated " […] by natural biological
processes, under experimental conditions that allow some aspects of the truth
to be known" 13. Such plasmodes may be created, e.g., in a wet lab by
manipulating biological samples as in case of a "spike in" experiment 11, 13.
The latter paper provides a very illustrative introduction to biological
plasmodes. In their detailed definition, the authors state that " […] a
plasmode can be defined as a collection of data that (i) is the result of a
real biological process and not merely the result of a computer simulation;
(ii) has been constructed so that at least some aspect of the "truth" of the
DGP is known”.

A number of research papers such as Mehta et al. 9, Vaughan et al. 13 deal
with "spike in" experiments in microarray expression analysis as an example
for a biological plasmode data set. As part of that experiment, real cases
from one population are randomly assigned to two groups. Then, a known amount
of transcript is added to serve as a positive control. As a result,
distributions and correlations in the generated data are viewed as most
realistic since being taken directly from real data. Besides others, Mehta et
al. 11 discusses application of plasmode data sets in high- dimensional
biology. Vaughan et al. 13 use plasmodes for the estimation of admixture, or
the proportion of an individual’s genome that originated from different
founding populations and thus illustrates the utility of plasmodes in the
evaluation of statistical genetics methodologies. Several authors such as
Sokal et al. 20, Mehta et al. 11 provide helpful insights into generation and
application of biological plasmodes as those to be expected to incorporate
valuable information on biological variation and capture biological reality.
Biological plasmode data sets have also been utilized to evaluate the
performance of statistical methods 21 and their validity 9. Plasmode data sets
were also used to investigate the validity of multiple factor analysis in a
known biological model 20.

Despite their ability to create new and more advanced biological set-ups, e.
g. , by crossing mice 13, sometimes biological plasmodes may become not only
very time-consuming but also require high experimental costs. Researchers
might eventually do not have a lab available to construct biological
plasmodes. In some cases, ethical reasons may also speak against the
construction of biological plasmodes. In all such situations, statistical
plasmodes offer an advantageous alternative.

Statistical plasmodes, being in focus of the present paper, begin with
generation of covariate information performed by applying resampling-based
methods to real data set 22, 4, 5; note that no biologically new samples are
created in the context of such resampling procedure. Further, an appropriate
OGM has to be applied to generate outcomes based on the resampled covariates
5, 23, 24, 25, 26, 27, 28, 7, 29, 30, 31, 32. In cases where the exposure
modeling is also a part of the study 5, 23, 24, 27, 29, some known “truth”
such as e. g. treatment effects can be added manually 5, 28, 26, 7, 29. For
instance, in Franklin et al. 5 the authors create statistical plasmode data
sets by " […] resampling from the observed covariate and exposure data without
modification to preserve the empirical associations among the variables.” In
that paper, the “true” treatment effect and the baseline hazard function are
estimated from the empirical data. Further, the associations between outcome
and the covariates as well as between censoring times and the covariates have
been described by means of two Cox proportional hazard models. Such modeling
approach corresponds to the application of a OGM in our terminology.

According to our interpretation, statistical plasmode simulations utilize
aspects of resampling (when generating the covariate information) as well as
parametric modeling (e.g., application of OGM, modeling of exposure etc.) and
thus can be interpreted as semi-parametric methods.

There are numerous applications of plasmodes generated by certain methods of
data modifications in the literature. Some of those applications are based on
statistical approaches in the sense of our definition. For instance,
Tibshirani 22 utilize plasmodes to assess sample size requirements in
microarray experiments when estimating the false discovery rate and false
negative rate for a list of genes. Gadbury et al. 12 illustrates use of
plasmodes by comparing the performance of 15 statistical methods for
estimating the false discovery rate in data from an high-dimensional
experiment. Elobeid et al. 33 employs plasmode data sets to analyze the
performance of several statistical methods used to handle missing data in
obesity randomized controlled trials. Reeb and Steibel 4 suggest an
interesting application of plasmode data sets to complement the evaluation of
statistical models for RNA-seq data. In their subsequent paper, Reeb et al. 34
then use plasmode data sets to assess dissimilarity measures for sample-based
hierarchical clustering of RNA sequencing data. In Franklin et al. 5,
plasmode-based studies are used for the evaluation of pharmacoepidemiologic
methods in complex healthcare databases. Resampling in combination with
outcome generation by a logistic model to compare the HDD propensity score
method with ridge regression and lasso is used by Franklin et al. 23. Franklin
et al. 24 use plasmode-based studies to compare the performance of propensity
score methods in the context of rare outcomes. In Desai et al. 26, the authors
utilize plasmode data sets to analyze the uncertainty in using bootstrap
methods for propensity score estimation whereas Liu et al. 28 conducts a
plasmode-based study to compare the validity and precision of marginal
structural models estimates using complete case analysis, multiple imputation,
and inverse probability weighting in the presence of missing data on time-
independent and time-varying confounders. The issue of data imputation has
also been addressed in Atiquzzaman et al. 7 where the authors used plasmodes
to compare two imputation techniques when imputing body mass index variable in
osteoarthritis-cardiovascular disease relationship. In Ejima et al. 35, the
authors use statistical plasmodes to assess type I and type II error rates of
analyses commonly used in murine genetic models of obesity. Similarly, Alfaras
et al. 36 resample from the empirical distributions to create plasmode data
sets for murine aging data. Those plasmodes are then utilized to compute type
I error rates and power for commonly used statistical tests without assuming a
normal distribution of residuals. In their most recent study, Hafermann et al.
31 designs a plasmode simulation study to investigate how random forest and
machine learning methods may benefit from external information provided by
prior variable selection studies. Rodriguez et al. 32 evaluate plasmodes as
being useful for preserving the underlying dependencies among hundreds of
variables in real-world data used to evaluate the potential utility of novel
risk prediction models in clinical practice; the authors generate plasmodes
when studying lung transplant referral decisions in cystic fibrosis.

To our understanding, two central steps in the statistical plasmode generation
procedure can be derived, namely:

  * (i)

Generation of the covariate structure by resampling from an original data set

  * (ii)

Outcome generation that includes

    * (ii.1)

Choice of an appropriate outcome generating model (OGM)

    * (ii.2)

Choice of covariate effects either by individual specification or by
estimation based on the original data

    * (ii.3)

Generation of new outcomes by applying the OGM chosen in (ii.1), with the
effects specified in (ii.2) and applied to the covariates generated in (i)

The discussion performed in this section is summarized in Table 1 that
provides a comparative summary of parametric and plasmode simulation studies.

Table 1: Parametric simulations versus statistical plasmodes: Similarities and
differences

Feature |  Parametric Simulations |  Statistical Plasmodes
---|---|---
Data-generating process (DGP) |  DGP is to be specified in advance |  No DGP specification is required
Outcome-generating model (OGM) |  Parameters of a chosen OGM to be estimated from data or derived from literature or set manually |  Parameters of a chosen OGM to be estimated from the original data or derived from literature or set manually
Range of possible scenarios |  Arbitrary scenarios, in particular, extreme and rare scenarios, can be generated |  Only reality bounded to the sample at hand can be generated
Knowledge of "truth" |  "Truth" must be completely known in advance |  At least some "truth" such as effect sizes should be known a priori
Data availability and representativeness |  Irrelevant for simulations based on literature results or previous knowledge |  Crucial, as the simulated data is always limited to the sample at hand
Reality reflection |  Parametric simulations may not be able to capture the complexity of real life data |  Plasmodes are expected to resemble the reality in the most accurate way
High-dimensional data simulations |  Usually time - and cost - consuming. Latent dependencies may also become an issue |  Mostly straightforward, as no estimation of distributions and/or dependencies is required
Small sample sizes |  Essentially uncomplicated, but may become an issue in cases when simulation parameters are to be estimated from the real data at hand |  Difficult due to resampling
Dependence structure |  Becomes a challenge with complex dependencies |  No modeling/estimation of dependence structure is required

That discussion as well as the supporting literature imply that plasmodes
provide an attractive supplement to parametric simulations in data-based
research. In particular, it is expected that plasmode datasets resemble the
reality most closely, especially regarding dependency structures. In the
following, we will analyze plasmodes in detail to examine their strengths and
weaknesses in more detail.

##  3 Challenges of Statistical Plasmode Simulations

Simulations studies are, at least in the scope of the present work, designed
to enable the practical analysis of statistical methods. To this end, data
generation should satisfy several criteria, such as, amongst others, to
provide the basis for subsequent undistorted model comparisons or to enable a
specific covariate dependence structure. Constructing, reporting and
comprehending a parametric simulation study is mostly straightforward and
transparent, as the resulting data is artificially generated in a target-
oriented way. Critical steps in the construction of parametric simulations
include, for instance, the investigator’s choice of the outcome-covariable
association. While this ambiguity is shared by statistical plasmode
simulations, many of the properties of statistical plasmodes are typically
less obvious and verifiable because statistical plasmodes are designed with
the complex task to mimic reality in the closest way while simultaneously
specifying some aspects of the truth. The main advantage of statistical
plasmodes lies in their ability to generate data with specific distributions
and dependence structures without the need for explicit assumptions. The
assumption that statistical plasmodes can faithfully generate data that
closely resemble reality has rarely been questioned. In practice, however,
this assumption can present challenges. For instance, the lack of statistical
analyses or simulations to verify the preservation of dependence structures
can undermine the reliability of the generated data. Consequently, the
advantages attributed to statistical plasmodes can also transform into
challenges. Further potentially critical steps in the construction of plasmode
data include the representativeness of the underlying data and the choice of
the resampling scheme. Below, we theoretically discuss these potential
pitfalls in more detail while also providing corresponding examples from
literature.

###  3.1 Resampling of Covariate Information

In our concept of statistical plasmodes, the simulation is based on the
generation of covariable information by resampling from a real data set. This
has the intention to preserve the characteristics of the original underlying
dataset such as, amongst others, the number and type of covariables and the
corresponding dependence structure, see e. g. Franklin et al. 23, Atiquzzaman
et al. 7, Conover et al. 29. This preservation is primarily achieved through
the use of appropriate resampling techniques. Consequently, the applied
resampling scheme, which consists of specifying the number of generated
datasets ($N$) and the resampling technique, has central importance for the
generated plasmode (covariable) datasets.

Of note, while utilizing resampling, statistical plasmodes are arguably even
more complicated to analyze because of the additional artificial outcome
generation. Consequently, not all established theoretical results concerning
resampling might be transferable to the full plasmode dataset but only to the
plasmode covariable datasets. We use the terminology resampling and bootstrap
interchangeably and indicate the concrete resampling/bootstrapping technique
if necessary. The analysis of resampling methods is almost exclusively
formulated in terms of the asymptotic performance of the bootstrap
distribution $L^{*}$ of an estimator $T$ (e. g. variance, confidence interval)
applied to the empirical distribution of the resampled data 37 e.g.. For
statistical plasmode covariable datasets, the estimator $T$ could be, for
example, some function of the covariance matrix of the covariables
(preservation of correlation structure). When considering the statistical
plasmode procedure as a whole, the estimator of interest $T$ typically
utilizes the artificial outcomes, e. g. $T$ could be the linear predictor in
ridge regression when investigating its performance compared to other models.
The resampling is said to "have worked", if $L^{*}$ convergences weakly to $L$
(the theoretical distribution of $T$) for increasing sample size $n$ of the
underlying data 38 e.g.. Otherwise, one speaks of “bootstrap failure”. In the
following, we discuss the influence of the chosen resampling scheme on the
generated plasmodes in more detail, focusing mainly on the preserveness of the
covariable information.

Number of Plasmode Datasets $N$. The specification of the number of resampled
plasmode datasets $N$ is often performed ad-hoc and potentially leads to
different answers to the same question, in particular if $N$ is specified as
too small 39. In the framework of bootstrap tests, Davidson and MacKinnon 40
propose a pretest procedure for choosing the number of bootstrap samples to
minimize the loss of power due to $N$ being finite. A more general, data-
dependent three-step procedure is proposed by Andrews and Buchinsky 39 who
estimate $N$ to achieve a desired accuracy of the approximation of the
bootstrap to the ideal ($N\rightarrow\infty$) distribution of the estimator of
interest. However, to the best of our knowledge, there is no general guideline
to theoretically specify the number $N$ of datasets to be generated in a data-
independent way (i. e. without already performing the resampling scheme) such
that asymptotic resampling results hold with sufficient accuracy. Moreover,
existing results might not be valid for statistical plasmodes due to the
additional artificial outcome-generation procedure.

In the plasmode literature, $N=500$ 5, 26, 7 e.g. and $N=1000$ 28, 27, 31 e.g.
seem to be popular ad-hoc choices. We have not seen any application where the
choice of $N$ was explicitly justified or the convergence or stability of the
subsequent analyses applied to the plasmode data have been checked for
increasing $N$. In summary, the number of datasets can be critical aspect in
the generation of statistical plasmodes, in particular if convergence of $T$
is not reached. In Section 4 we provide some recommendations for determining
$N$ which we further illustrate in Section 5.

Resampling Technique. Resampling can be performed without replacement as in
the $n$-over-$m$ bootstrap (subsampling with $m<n$) and sample-splitting
(cross-validation) bootstrap. Subsampling draws from the data-generating
process of the original data 41 and has been shown to lead to consistent
estimators under minimal conditions, see Theorem $1$ in Bickel et al. 37, as
long as the subsampling size $m$ and the size of the original dataset $n$ are
appropriately specified. Alternatively, resampling with replacement such as
the $n$-out-of-$n$ bootstrap (also called nonparametric bootstrap) can be
utilized. Resampling schemes based on drawing with replacement draw from
empirical probability distribution derived from the underlying data 41 and
require additional assumptions for consistent estimation, but are more
efficient if the bootstrap "works" 37. However, the nonparametric bootstrap
can fail, for example when the limiting distribution of the estimator has
discontinuities, when estimating extrema and when setting critical values for
some test statistics 37, 42, 43. As a remedy, the $m$-out-of-$n$ bootstrap
(sampling $m\leq n$ with replacement) has been introduced to prevent bootstrap
failure while losing efficiency if the nonparametric bootstrap was consistent.
Sampling fewer than $n$ observations has since been treated as a “cure-all”
method (being asymptotically valid under weak assumptions and not failing)
which has been critically discussed, for instance, in Andrews and Guggenberger
43. A comprehensive overview of resampling techniques is provided, e. g. , in
Bickel et al. 37.

For increasing sample size $n$, estimators based on subsampling and the
$m$-out-of-$n$ bootstrap become more similar as the probability of repeating
observations decreases. Note that, contrary to subsampling, resampling with
replacement allows for $m=n$. The additional requirements for the consistency
of sampling with replacement compared to sampling without replacement mainly
state, informally speaking, that the influence of tied observations on the
bootstrap estimator should be small 38.

In the majority of literature concerning plasmode generation, the $m$-out-
of-$n$ bootstrap has been used 5, 25, 30 e.g., whereas the nonparametric
bootstrap has been used by Rodriguez et al. 32, the sample-split bootstrap by
Gerard 44 and subsampling by Hafermann et al. 31. In some publications we did
not find indications whether resampling was performed with or without
replacement, e. g. in Ju et al. 45, Ripollone et al. 27. All in total, the
type of resampling technique influences asymptotic properties of the
covariables and hence of the plasmode datasets, effects whether the resampling
“has worked” and consequently impacts subsequent analyses on the generated
plasmodes datasets. However, to the best of our knowledge, we have not seen
any application in the literature concerning plasmode data generation in which
the choice of a particular resampling technique has been explicitly justified.

Resampling Size $m$. Using resampling with replacement of size $m$, with
$m\rightarrow\infty$ and $m/n\rightarrow 0$, typically resolves failure of the
$n$-out-of-$n$ bootstrap, but requires the specific choice of $m$ as a key
issue 38. An adaptive rule for the choice of $m$ for subsampling and the
$m$-out-of-$n$ bootstrap in the case of independent observations has been
proposed by Bickel and Sakov 38 and is further illustrated in our example in
Section 5. Informally speaking, if $m$ is in the right range of values, the
bootstrap distributions of the estimator for the similar $m$’s are close to
each other, indicating consistency of the estimator. The rule provides an
adaptive estimator $m^{*}(n)$ and leads to optimal convergence rates of the
estimator irrespective whether the nonparametric bootstrap would work in the
example (then $m^{*}(n)/n\rightarrow 1$ as $n\rightarrow\infty$) or would fail
(then $m^{*}(n)/n\rightarrow 0$).

To the best of our knowledge, only fixed resampling sizes $m$ have been chosen
in the plasmode literature, and we have not observed any explicit
justification of the specific value of $m$. In other words, $m$ appeared to be
chosen arbitrarily. For instance, Hafermann et al. 31 used a selection of
$m^{\prime}s$ ($250,500,1,000,2,000,4,000$) which are small compared to the
number of observations $n=198,895$ while Liu et al. 28 chose $m=500$ for a
dataset with $n=646$ and Atiquzzaman et al. 7 sampled $m=75,000$ out of
$n=84,452$. Interestingly, different authors used different values of $m$
which have been chosen without justification for the same underlying dataset,
as exemplified for the NSAID dataset with $n=49,653$. While $m=30,000$ has
been picked by Franklin et al. 5 and Ripollone et al. 27 used a comparably
large $m=25,000$ as well, Ju et al. 45 chose a much smaller value in $m=1,000$
and Wyss et al. 30 set $m=10,000$. In summary, when applying subsampling or
the $m$-out-of-$n$ bootstrap, the value of $m$ matters for the consistency of
$T$, and should be properly justified and adapted to the underlying data and
estimator(s) of interest.

Covariable Dependence Structure and HDD. Resampling the covariable information
has, amongst others, the aim of preserving the covariable dependence structure
of the underlying dataset; see e. g. Franklin et al. 5, Karim et al. 25,
Conover et al. 29. Under some assumptions such as i. i. d. observations and
finite fourth moments of the covariables, Beran and Srivastava 46 have shown
that the resampled covariance matrix converges to the original covariance
matrix for the nonparametric bootstrap when $n$ increases and the number of
covariables $p$ is fixed (i. e. most HDD situations excluded, see also below).
However, for other resampling schemes similar results have, to the best of our
knowledge, not been shown. For the $m$-out-of-$n$ bootstrap, the optimal $m$
could be estimated with the estimator $T$ specified to represent the
covariable covariance matrix in order to investigate and ensure that the
resampling scheme works (at least for that certain aspect of the data), see
also our example in Section 5. However, other aspects of the covariable
information, such as extreme values, might be more important in some
applications.

A well-discussed issue, in particular in the context of HDD, is the occurrence
of spurious correlations. Amongst others, Fan et al. 16 have shown that
sampling $p$ independent normal $n$-vectors leads to empirical covariance
structures strongly deviating from a diagonal matrix, in particular if $p\gg
n$. However, the risk of spurious correlations is not limited to parametric
simulations. An increasing number of covariates $p$ increases the risk that
the underlying data sample suffers from spurious correlations, which may be
propagated to the generated plasmode datasets by resampling. Further spurious
correlation is likely to distort the empirical covariance structure of the
statistical plasmodes, leading to even stronger deviations to the population
covariance matrix.

For a fixed number of covariables $p$, the bootstrap has been shown to work in
linear models if $p/n$ is small 47. If the number of covariables growths with
the number of observations, Mammen 48 have shown that the bootstrap works for
effect estimates in high-dimensional linear models if $p(n)\rightarrow\infty$
and $p(n)/n\rightarrow 0$ as $n\rightarrow\infty$, and Karoui and Purdom 49
have shown that the confidence intervals of the pairs and residual bootstrap
in linear models are too wide if $p(n)/n\rightarrow c,c\in\mathcal{R}$.

In summary, the goal of preserving the covariable dependence structure can be
used to determine an optimal resampling scheme. In particular in high-
dimensions, the covariance structure could, however, be distorted by spurious
correlations and whether resampling and subsequently statistical plasmodes
work in these scenarios might require additional research.

###  3.2 Representativity of the Underlying Data Sample

One of the main assets attributed to statistical plasmode simulations is that
they are expected to preserve the complex real-world data structure by
resampling the covariable information from a real dataset. Naturally, an
appropriate representative dataset has to be available and constitutes the
basis for the entire plasmode simulation study. Parametric simulations, on the
other hand, can be artificially constructed without requiring representative
data. The data sample is expected to represent the population of interest.
This limits the generalisability of the results of the analyses that the
plasmode simulation study was designed for, which has been acknowledged,
amongst other, by Franklin et al. 5, Liu et al. 28, Atiquzzaman et al. 7.

The data sample should satisfy the assumptions of the applied resampling
technique. As a result, the choice of the resampling technique depends
strongly on the underlying dataset at hand. Standard resampling techniques,
such as discussed above, assume that the observations are independent 37 e.g..
This assumption is violated, for instance, if the observations show clusters,
repeated measures, population structure or longitudinal measurements. In this
context, more sophisticated resampling scheme including block-wise resampling
have to be applied, for which most of the asymptotic results are not
explicitly formulated 38.

Depending on the underlying data and the resampling scheme, the
characteristics of the original dataset to be conserved might not be reflected
by the generated data. This is acknowledged by Karim et al. 25 who state that
“[…] it is possible that important confounders in the empirical study might
not remain important in the plasmode samples”.

In summary, the generated statistical plasmode datasets depend strongly on the
representativeness of the underlying real data and are limited to the
population represented by the data sample. The resampling scheme should be
adaptive to the characteristics of the real data such as population structure,
which have to be identified and reported.

###  3.3 Investigator’s Choice of the “Truth”

The concept of plasmode simulations is mainly based on preserving the complex
but realistic structure of the underlying data while inserting some “truth” by
investigator’s choice. These specification can be manifold in type, and
potentially distort the real-word characteristics of the generated plasmode
data.

Artificial Covariables. Additional to resampling covariable information,
important covariables such as exposure or treatment variables can be
artificially created to model some aspects of the simulation study. For
instance, Franklin et al. 24 and Conover et al. 29 model a binary exposure
variable in relation to confounder variables via logistic regression, while
Rodriguez et al. 32 simulate covariables at a later stage of a longitudinal
study. Naturally, artificially generating covariables does not preserve the
full real-world setting and should be performed with care.

Artificial Outcome Generation. With the covariable information generated (by
resampling or artificially), corresponding artificial outcomes are created, in
our concept of statistical plasmodes, according to some outcome-covariable
association specified by the investigator. A straightforward way to create a
transparent association between resampled covariates and the artificial
outcomes is to utilize regression models specified by the combination of a
link function (type of OGM) and the linear predictor (effect structure).

The OGM determines the type of the artificial outcome (e. g. binary, survival)
and strongly influences subsequent analyses on the generated data. For
instance, if the aim of the study is the performance assessement of several
models, the model closest to the chosen type of OGM has an advantage induced
by the investigator, leading to potentially distorted comparisons. Most
commonly, logistic regression is used for binary outcomes 24, 25, 7 e.g. and
the Cox model for survival endpoints 5, 23, 26 e.g., whereas Rodriguez et al.
32 apply an exponential survival model. Normal linear regression is used by
Liu et al. 28.

Besides the type of OGM, the determination of the effect structure of the
corresponding linear predictor is vital. Parts of the effect vector have been
specified by literature review 29, by sampling from independent standard
normal distributions 45, by estimation on the original dataset 28, 30 or
manually by investigator’s choice 5, 23, 26. Some authors specify the
treatment or exposure effect by hand while estimating the confounder effects
on the original data 25 e.g.. Specifying the value of the effects of
covariables might represent a strong intervention in the generation process of
a realistic dataset. Potential problems include the creation of artificial
outcome-covariate associations and invalidating or even nullifying existing
“real” associations between the covariables and the novel, artificial
outcomes, in particular if the effects are set manually by investigator’s
choice. If the effects are estimated, they depend on the underlying data
sample and estimation uncertainty is ignored. Additionally, the specification
of the type of OGM influences effect estimates and subsequent analyses might
become problematic. For example, an effect vector estimated by a sparse method
will induce advantages of sparse methods in model comparisons on the generated
data.

In summary, a crucial assumption for both parametric and plasmode simulations
is that the chosen outcome generation reflects realistic, natural or
biological associations between outcome and covariables. Whereas the choice of
the OGM and the effect structure is a natural aspect of parametric simulation
studies, plasmodes are often described as closely depicting reality. However,
the outcome data in statistical plasmodes are also artificially created while
inducing some investigator’s choice “truth”. These manipulations of the real-
data are harder to assess and less transparent than in parametric simulations
studies as the structure of the data is typically more complex. In the end,
plasmode generation also leads, at least in part, to artificial data and
constructed associations.

##  4 Statistical Plasmodes: Step-by-step Recommendations

We provide a hands-on overview of our recommendations for the generation and
reporting of our concept of statistical plasmodes in Figure 1. This summary
extends the basic plasmode generation procedure described in Section 2 and
addresses the critical steps discussed in the previous section. We
theoretically discuss our step-by-step procedure below and illustrate its
application in a real data example in Section 5.

Figure 1: Statistical plasmode data generation procedure step-by-step

Step 1: Planning of the Simulation Study. We recommend to clearly formulate
the research problem and to plan the simulation study using the ADEMP criteria
2. Fixing the aims and the data-generating processes aids in the choice
whether statistical plasmodes are needed in the first place or whether a
parametric simulation study might be more appropriate. Additionally, it guides
the choice of the population of interest in Step 2 and potentially the choice
of the resampling technique in Step 4. Importantly, the methods and
performance measures indicate the subsequent choice of the OGM in Step 6. For
instance, if we plan to assess the prediction performance of several models,
we should make sure that the chosen OGM does not bias the subsequent model
comparisons. Also, the choice of the OGM and data generating process determine
the scenarios in which the properties of novel statistical methods can be
empirically assessed in the context of the simulation study.

Step 2: Population of Interest. The population of interest can be of primary
interest and consequently be strongly connected to the aim of the simulation
study and the data-generating processes determined in the previous step, in
particular if methods are developed to deal with populations with certain
characteristics (e. g. many missing values, complex covariance structure,
high-dimensionality). However, we might also be mainly be interested in the
analysis of statistical methods such that the population of interest serves
mainly as an illustration and is not necessarily connected to Step 1. In the
latter case, particular effort should be taken clarify why the chosen
population covers those situations in which the methods under consideration
are claimed to work. In any case, the hypothetical population should be stated
and described as clearly as possible, e. g. the entity of interest,
covariables and population structure. Since the datasets generated by the
statistical plasmodes should be representative of the population of interest,
this step influences several of the following steps, in particular Steps 3, 4
and 6.

Step 3: Representative Sample. It is a central aspect of statistical plasmodes
that the underlying sample is representative of the population of interest
clarified in the previous step, refer also to the discussion in Section 3.2.
Consequently, it is vital to investigate and communicate why the utilized data
sample represents the population of interest and which potential limitations
arise. In particular, it should be stated how the data was sampled, which
covariables are included and what the endpoint of interest is. Additionally,
the sample size should be justified and potential population structures
investigated, as this can influence the choice of the resampling technique,
see Steps 4 and 5 and the discussion in Section 3.1. Note that even if the
sample is representative, the generated plasmode datasets might not be, e. g.
as a result of a poor resampling plan or a outcome generation that distorts
the relationship between the covariables or the outcome-covariable
association.

Steps 4 & 5: Resampling Scheme The resampling scheme consists of the number of
bootstrap samples, the type of resampling technique used and, if applicable,
the justification of the resampling size, see also Section 3.1. It determines,
together with the data sample, the plasmode covariable datasets. Each of the
aspects of the resampling scheme plays a crucial role for the asymptotic
properties of the estimators applied to the data generated by resampling and
to properties such as the preserveness of the covariable correlation
structure, see the discussions in Section 3.1. Unfortunately, the resampling
scheme has to be decided on for each application individually, while keeping
those research aims and properties of the population of interest, that should
be preserved with high priority, in mind. Ensuring the plasmode sets are drawn
from the hypothetical population is only possible by applying subsampling 41,
whereas sampling with replacement draws from the empirical distribution of the
data sample specified in Step 3. However, if the underlying dataset is
representative of the population of interest, drawing with replacement might
become preferable due to its increased efficiency and second-order properties
37. As described in Section 3.1, the nonparametric bootstrap potentially
fails, although this is often impossible to know before the application. To
avoid bootstrap failure, we recommend to utilize the $m$-out-of-$n$ bootstrap
although this might lead to efficiency losses. The optimal resampling size $m$
can be determined by applying the algorithm introduced in Bickel and Sakov 38
while using those properties of priority as estimator in Step 2 of the
optimization algorithm for $m$. It has to be noted that the estimation of $m$
might require high additional computational cost. In many applications it
might be meaningful to opt for some function of the covariable covariance
structure as an estimator, as it is often stressed that the empirical
dependence structure of the original dataset should be preserved.

Step 6: Outcome Generating Model. The choice of the OGM includes the type of
model to determine the association between the resampled covariables and the
novel artificial outcomes, as well as corresponding OGM components such as
effect sizes, for example. The OGM determines the artificial outcomes in type
and value, and is a crucial component for many research questions formulated
in Step 1. Also, it gives the investigator the opportunity to fix some aspects
of the “truth”, see also the discussion in Section 3.3. Special care should be
taken that the OGM does not bias the subsequent analyses that the statistical
plasmode simulations are generated for. To do so, it might be helpful to
investigate the models or methods to be compared in detail and contrast them
with the OGM. For instance, a sparse OGM will most likely support sparse
models in subsequent model comparisons. If the effect structure is chosen in a
sparse way, a sparse model might be more likely to correctly estimate the
effect sizes or perform valid predictions. Additionally, if important
relationships between variables have been detected, the effect structure
should be chosen accordingly to preserve these. For instance, in linear
predictor models, the observed outcome variation depends on the (co-)variances
of the covariables weighted by their corresponding effects, stressing their
influence on the artifical outcomes of the plasmode datasets.

Step 7: Outcome Generation. Each of the $N$ plasmode covariable sampled in
Steps 4 and 5 is combined with the OGM determined in Step $6$ to create $N$
corresponding artificial plasmode outcome vectors.

Step 8: Quality Checks. The quality of the covariables can be assured by
appropriate resampling as described in Steps 4 and 5. It is, however, often
not feasible to compare the original covariable covariance structure with
those of the $N$ statistical plasmode datasets. More research might be
necessary to judge the distance of the original and the generated data. The
original outcome values of the real dataset are, if at all, only explicitly
used to determine the effect structure in Step 6. The quality of the simulated
data could be checked by comparing the generated outcomes of (some of) the
statistical plasmode datasets with the original outcome. The type of
potentially meaningful checks depends on the type of outcome. For continuous
observations, the distributions of the two outcomes could be compared by the
empirical densities or histograms as is done e. g. in Franklin et al. 5.
Additionally, the range of the data should be checked as well as potential
outliers. For categorical (including binary) outcomes, the prevalence of the
classes can be compared, see for example Franklin et al. 5.

Step 9: Reporting. Each step of the statistical plasmode generation should be
justified and reported to enhance reproducibility and transparency of the
proposed data generation procedure. Whenever appropriate, plasmode generation
should follow the scheme presented in Figure 1 and the corresponding
descriptions provided in the present section. Additionally, the research
question determined in Step 1 should be addressed.

##  5 Statistical Plasmodes: A numerical example

The following example has been constructed to illustrate the step-by-step
procedure introduced in the previous section.

Step 1. Assume that we are interested in the aim (A) of investigating the
application of ridge regression 50 and the linear mixed model 51 e.g. in the
context of high-dimensional RNA-expression data with sparse effects on a
normal outcome (data-generating process, D). The estimands (E) are specified
as the parameter vector and the linear predictor in the respective model. We
split the sample once into training and test data ($2:1$), which we deem
sufficient for our illustration purposes. The plasmode datasets are generated
using the training data. Each plasmode dataset of size $m$ and number of
covariates $p$ is analysed (methods, M) using ridge regression of the form

$\displaystyle\ y=\mu
1_{m}+X\beta+\varepsilon,\quad\lVert\beta\rVert_{L_{2}}^{2}\leq\lambda,\quad\varepsilon\sim\mathcal{N}(0,\sigma^{2}I_{m\times
m})$ (1)

via penalized maximum likelihood with cross-validation for $\lambda$ as
implemented in the R- package glmnet Friedman2010, as well as the linear mixed
model in the variance components form

$\displaystyle\ y=\mu
1_{m}+X\beta+\varepsilon,\quad\beta\sim\mathcal{N}(0,\sigma_{\beta}^{2}I_{p\times
p}),\quad\varepsilon\sim\mathcal{N}(0,\sigma_{\varepsilon}^{2}I_{m\times m})$
(2)

with restricted maximum likelihood estimation as implemented in the R-package
sommer 52. Here, $1_{m}$ denotes the $m$-column vector of ones while
$I_{p\times p}$ denotes the identity matrix of dimension $p$. As performance
measures, we utilize the mean absolute bias

$\displaystyle\
\mathrm{MAB}=\frac{1}{p+1}\lVert(\hat{\mu},\hat{\beta})-(\mu,\beta)\rVert_{L_{1}}$
(3)

where $\mu$ and $\beta$ are known as part of the “truth”, and the sample-split
mean squared error of prediction

$\displaystyle\
\mathrm{MSEP}=\frac{1}{m}\lVert\hat{y}-y_{\text{test}}\rVert_{L_{2}}^{2},\quad\hat{y}=\hat{\mu}1_{m}+X_{\text{test}}\hat{\beta}$
(4)

where $y$ corresponds to the artificial outcome in the test split. We estimate
both measures using the mean of the estimates (indexed by superscript $b$) in
the generated $N$ statistical plasmode datasets

$\displaystyle\
\widehat{\mathrm{MAB}}=\frac{1}{N}\sum_{b=1}^{N}\frac{1}{p+1}\lVert(\hat{\mu},\hat{\beta})^{(b)}-(\mu,\beta)\rVert_{L_{1}},\quad\widehat{\mathrm{MSEP}}=\frac{1}{N}\sum_{b=1}^{N}\frac{1}{m}\lVert\hat{y}^{(b)}-y_{\text{test}}\rVert_{L_{2}}^{2}.$
(5)

and visualize the $N$ individual measures via boxplots, see Step 9 and Figures
4 and 5.

Step 2. In the scope of this example, we are interested in the model choice
for high-dimensional RNA-expression data with normal outcomes for female
breast cancer patients which constitutes the population of interest.

Step 3. The data sample underlying the statistical plasmode simulation was
generated by The Cancer Genome Atlas (TCGA) Research Network
(https://www.cancer.gov/tcga). The breast carcinoma (BRCA) cohort which
provides a basis for the following numerical example was last updated on May
31, 2016.

We restrict the publicly available data to $n=1098$ female patients with
breast cancer with cancer tissue, i. e. excluding normal tissue and male
patients. RNAseqV2 gene expression data and clinical data for BRCA were
obtained from the TCGA Data Portal 53 via the R/Bioconductor package
TCGAbiolinks 54, 55, 56. For computational reasons, we choose $p=5000$ out of
the $25828$ available genes at random. The R/Bioconductor package limma 57 has
been utilized to normalize the RNA gene expression data. The expression levels
can be assumed to be measured continuously and they show different shapes and
ranges. This is illustrated in Figure 2 using their empirical distributions at
four randomly chosen genes.

The outcome of interest is age at diagnosis date which can be considered to be
approximately normally distributed, see Figure 3A. While the dataset can be
considered to be representative of a female breast cancer population from the
United States of America, we acknowledge that RNA expression data from other
populations (e. g .different countries) might lead to different results for
our research question.

Step 4. Before the analysis, we set the number of plasmode datasets to be
generated to $N=500$. In the final step 9, we investigate the convergence of
the estimators of the performance measures, see equation (5), in the
statistical plasmode datasets, see also Figure 6.

We choose the $m$-out-of-$n$ bootstrap in order to prevent potential bootstrap
failure but with the potential drawback of losing estimation efficiency.
Performance analysis of resampling method and the estimation of the optimal
$m$ requires the specification of an estimator which is applied to the
generated data. Since resampling in statistical plasmodes is primarily
concerned with the covariate information, already using the performance
measures defined in equation (5) as estimators is not feasible as they require
the subsequent artificial outcome generation. Naturally, there are several
reasonable estimators that could be considered. In this example, we opt for
the covariate dependence structure as the measure of interest because the
majority of publications which applied plasmodes referred to the advantage of
the preserveness of the original covariable dependence structure.

We determine the resampling size $m$ via the algorithm described in Bickel and
Sakov 38. In particular, to adopt that algorithm to our problem formulation,
we specify the sequence of potential $m$’s by setting $q=0.97$, choose the
$L_{2}$-norm of the covariance matrix (of the resampled covariate data) as a
metric, and calculate the resulting empirical distribution functions. We
estimate the covariance matrix using the Ledoit-Wolf linear shrinkage
estimator 58 to obtain a more precise estimate which is necessary because the
covariate data are high-dimensional. The optimal resampling size $m^{*}$ is
the one which minimizes the distance between the distributions of subsequent
$m$’s, where the distance is exemplarily measured by the Wasserstein metric.
The optimal resampling size based on the Wasserstein metric using $100$
iterations resulted in $m^{*}=711$.

We acknowledge that there is variety of optimal resampling sizes $m^{*}$ if
any of the parameters of the algorithms would be changed (such as, amongst
others, estimator, distance metric for empirical distributions and sequence of
potential m’s).

Step 5. We apply resampling with replacement of size $m^{*}=711$ to the matrix
of covariable information to obtain $N=500$ statistical plasmode covariable
datasets. As we have determined the subsampling size with optimality criterion
as the $L_{2}$-norm of the covariance matrix of the covariables, the empirical
covariance structure of the original dataset should be sufficiently preserved.

Step 6. We choose the LASSO 59 as an appropriate OGM to represent the sparse
effect structure associated with the high-dimensional data as required in Step
1. Additionally, the LASSO most likely does not distort the comparison between
ridge regression and the linear mixed model as both of these methods are
shrinkage methods used to model polygenic effects. The “true” effect structure
for the LASSO is chosen as the vector of estimated effect sizes obtained after
a LASSO had been fit to the original data. The proportion of non-zero
estimated effects was $95.4\%$ ($4768$ vs $232$). This implies that $232$
covariables are selected in the investigator’s choice “truth” while $4768$
genes are given a null effect. The effect sizes of the selected covariables
has a median of $-0.01$ (range $[-2.97,2.05]$).

Step 7. We generate one artificial outcome vector of size $m^{*}=711$ for each
of the $N=500$ plasmode covariate datasets by calculating the linear predictor
based on the combination of the resampled covariable information (step 5) and
the “true” effects (step 6). Thus, we obtain statistical plasmode simulations
based on real covariate information with an investigor’s choice “truth”.

Step 8. The artificial outcomes (of some of) the $N=500$ plasmode datasets are
compared with the original outcomes via histograms in Figure 3. The
distribution of the original and artificial outcomes is very similar in shape
and mean. The range of the original outcomes is larger than the range of the
artificial outcome which can be explained by the outcome generation via
resampled covariables and effects determined by LASSO (sparse and shrunken
effects) which most likely will not lead to more extreme outcome values than
contained in the underlying dataset.

We conclude that the artificial outcome data come close to reality but might
not properly reflect extreme values. The range of the artificial outcomes
could be increased, e. g. , by manually altering some elements of the effect
vector estimated by LASSO (as investigator’s choice of the “truth”). By doing
so, however, we would further alter the association between some of the
covariables and the novel outcomes.

Step 9. Finally, we compare the performance of ridge regression and the linear
mixed model in our statistical plasmode simulations. The MAB of the Ridge
regression is estimated as $0.025$ while the estimate of the MAB of the linear
mixed model is $0.022$, see also equation (5). The sample-split MSEP of ridge
regression is estimated as $18.30$ while the sample-split MSEP of the linear
mixed model is $12.83$, see also equation (5). In Figure 4 and 5 we depict the
estimated values for each plasmode dataset via boxplots. This suggests that in
our generated statistical plasmode simulations, which represent high-
dimensional RNA expression data with sparse known effects and artificial
normal outcomes, the linear mixed model performs superior to ridge regression.

Additionally, we illustrate in Figure 6 the convergence of the performance
measures for increasing number of plasmode dataset. The estimators for MAB and
MSEP for both ridge regression and the linear mixed model seem to have
stabilized at about $300$ generated simulations. Thus, we conclude that the
generated number of statistical plasmode is sufficient to obtain stable
estimates of the performance measures defined in step 1.

Figure 2: Empirical distributions illustrated by histograms (15 breaks each)
and smoothed densities for four genes selected at random to illustrate the
differences in location and shape. Figure 3:  Empirical distributions
illustrated by histograms (15 breaks each) and smoothed densities for A: the
original outcome (age at diagnosis) and B-D: vs. artificial outcomes of three
plasmode datasets selected at random.  Figure 4: Boxplots of the mean absolute
bias for both the linear mixed model and ridge regression in $500$ statistical
plasmode datasets. Figure 5: Boxplots of the sample-split mean squared error
of prediction for both the linear mixed model and ridge regression in $500$
statistical plasmode datasets. Figure 6: Convergence of the performance
measures MAB and MSEP for increasing number $N$ of statistical plasmode
datasets.

##  6 Conclusions and Outlook

Many simulation studies impose relatively strong assumptions regarding the
nature of randomness in the data and its dependence structure. Mostly of
theoretical kind, those assumptions primarily rely on the assumptions inherent
in the statistical models applied to generate the data. Since not all
assumptions can be justified in applied settings, the corresponding simulation
studies may not be able to capture biologically meaningful relationships and
thus result in misleading conclusions and research findings.

To avoid (at least some of) those issues, plasmode data sets are considered as
an alternative data generation approach. While parametric simulations are
known to provide only a partial representation of reality 4, plasmodes have
been declared to generate data that resemble reality in the closest way 9
e.g.. Highly appreciated for their ability to generate most realistic data,
plasmode do not impose any specific model assumptions on their data generation
process. Thus, no assumptions need to be justified to address the
applicability of plasmodes. Nevertheless, a number of assumptions such as the
representativeness of the underlying data sample have to be verified in order
to guarantee the reliability of the generated plasmode data.

Plasmodes can accommodate unknown features such as dependence structure,
distributions, and others, in particular, in the case of high-dimensional
data. We recall that in case of parametric simulations most of those
quantities are to be specified in advance. All in total, plasmode data sets
may provide an attractive supplement to parametric simulations and can be
applied in order to increase the reliability of the obtained research results.

In the present paper, we first discuss the concept of statistical plasmodes as
those created by resampling of covariate information from empirical data at
hand and subsequent outcome generation using an appropriate outcome-generating
model. This is what distinguishes them from biological plasmodes which are
usually created by conducting lab experiments. We interpret statistical
plasmodes as an intermediate step between the parametric and nonparametric
simulations, with the parametric component represented by the chosen outcome-
generating model. After the introduction of statistical plasmodes, we discuss
their main advantages and challenges and propose a step-by-step scheme for
their generation and reporting. That scheme is then illustrated by means of a
numerical example. All discussions in the present paper are presented in the
context of prediction and explanatory models.

Plasmodes are bounded to the sample they are based on, and thus cannot produce
the same variety of different scenarios as parametric simulations do. In this
context, questions on the data availability and representativeness arise. In
particular, even if plasmodes offer a flexible data generation procedure which
creates realistic data, the representativeness of the generated data still
substantially depends on the representativeness of the underlying real data
set. To address this limitation, some authors such as Ejima et al. 35 assume
that the empirical data at hand represents the entire population of interest.
Of course, such an assumption cannot be satisfied in each particular
situation.

Spurious correlations are another issue closely related to the question of
representativeness. Although plasmodes do not specify the underlying
dependence structure explicitly, they do reproduce it to a certain extent
while generating new data. Thus, if the sample at hand does not adequately
represent the population of interest, the existing spurious correlations may
be increased or even distorted for the generated plasmode data sets. As a
result, the corresponding generated dependence structure will not represent
the real one.

Statistical plasmodes as introduced in the present paper incorporate features
from both parametric simulations and resampling approaches, and, as a result,
inherit the strengths and weaknesses of each data generation method. On one
hand, statistical plasmodes offer the advantage of creating more realistic
data by generating covariate information through resampling techniques. On the
other hand, they may also introduce certain challenges with respect to the
subsequent model comparisons, as compared to purely parametric simulations 30.
Statistical plasmodes enable to control and manipulate certain aspects of the
“truth" through the use of parametric OGMs, which can be advantageous over
pure resampling methods. Nevertheless, asymptotic results established for
resampling techniques may not be directly applicable to statistical plasmodes.

Our discussion points out several interesting options for future research.
First, basic expectations placed on the plasmodes are related to their ability
to preserve real data distributions, the underlying dependence structure and,
as a result, the existing empirical associations. Those expectations are to be
guaranteed by resampling from the observed covariate data at hand, without any
additional data modification. However, it is not obvious how the choice of a
particular resampling technique and specification of its parameters (such as
the subsampling proportion in case of the subsampling technique) might impact
the robustness of the obtained data generation results, e. g. , in the context
of spurious correlations or sparse data. Additionally, the calculation of the
optimal $m$ might require high computational costs. A closer analysis of these
impacts are possible topics for future research.

Second, a data generation method is considered to be realistic if it reflects
the real data structure and the existing dependencies in the most accurate
way. Thus, appropriate distance measures need to be specified in advance and
also included into the reporting step of the data generation procedure. Such
measures can then be used to measure the closeness of the generated plasmode
data set to the underlying real data set. The choice of an appropriate
distance measure, as well as the robustness of the plamode generation
procedure with respect to that choice, can also be an interesting research
topic.

Finally, the outcome-generating models present the major obstacle for
plasmodes to become a purely non-parametric data generation approach. In the
future we intend to analyze the impact of an OGM on the performance of the
plasmode data generation procedure and to construct examples where the
replacement of a parametric OGM with a non-parametric one improves the
obtained data generation results. It is also of great interest to address
possible “plasmode failure" for data sets generated through statistical
plasmodes.

In total, our paper presents a comprehensive analysis of statistical plasmode
simulations, discusses their potentials and central challenges and provides
step-by-step recommendations for their generation. Our future research aims to
address (at least some of) the pitfalls in the most close way to potentially
provide more understanding and further novel insights into statistical
plasmode generation.

## Acknowledgements

The authors would like to thank Jörg Rahnenführer, Andrea Bommert and Marieke
Stolte for helpful discussions.

## References

  * Burton et al. 2006 A. Burton, D. G. Altman, P. Royston, and R. L. Holder.  The design of simulation studies in medical statistics.  _Statistics in Medicine_ , 25:4279–4292, 2006.
  * Morris et al. 2019 T. P. Morris, I. R. White, and M. J. Crowther.  Using simulation studies to evaluate statistical methods.  _Statistics in Medicine_ , 38:2074–2102, 2019.
  * Bin et al. 2019 R. D. Bin, A.-L. Boulesteix, A. Benner, N. Becker, and W. Sauerbrei.  Combining clinical and molecular data in regression prediction models: insights from a simulation study.  _Briefings in Bioinformatics_ , 00((0)):1–17, 2019.
  * Reeb and Steibel 2013 P. D. Reeb and J. P. Steibel.  Evaluating statistical analysis models for rna sequencing experiments.  _frontiers in Genetics_ , 4:1–9, 2013.
  * Franklin et al. 2014 J. M. Franklin, S. Schneeweiss, J. M. Polinski, and J. A. Rassen.  Plasmode simulation for the evaluation of pharmacoepidemiologic methods in complex healthcare databases.  _Computational Statistics and Data Analysis_ , pages 219–226, 2014\.
  * Schulz et al. 2017 A. Schulz, D. Zöller, S. Nickels, M. E. Beutel, M. Blettner, P. S. Wild, and H. Binder.  Simulation of complex data structures for planning of studies with focus on biomarker comparison.  _BMC Medical Research Methodology_ , 17, 2017.
  * Atiquzzaman et al. 2020 M. Atiquzzaman, M. E. Karim, J. Kopec, H. Wong, M. A. Vera, and A. H. Anis.  Using external data to incorporate unmeasured confounders: A plasmode simulation study comparing alternative approaches to impute body max index in a study of the relationship between osteoarthritis and cardiovascular disease.  _Journal of Statistical Research_ , 54(2):131–145, 2020\.
  * Shmueli 2010 G. Shmueli.  To explain or to predict?  _Statistical Science_ , 25:289–310, 2010.
  * Mehta et al. 2004 T. Mehta, M. Tanik, and D. B. Allison.  Towards sound epistemological foundations of statistical methods for high-dimensional biology.  _Nature genetics_ , 36:943–947, 2004.
  * Boulesteix et al. 2020 A.-L. Boulesteix, R. H. Groenwold, M. Abrahamowicz, H. Binder, M. Briel, R. Hornung, T. P. Morris, J. Rahnenführer, and W. Sauerbrei.  Introduction to statistical simulations in health research.  _BMJ Open_ , 2020.
  * Mehta et al. 2006 T. Mehta, S. O. Zakharkin, G. L. Gadbury, and D. B. Allison.  Epistemological issues in omics and high-dimensional biology: give the people what they want.  _Physiol Genomics_ , 28:24–32, 2006.
  * Gadbury et al. 2008 G. L. Gadbury, Q. Xiang, L. Yang, S. Barnes, G. P. Page, and D. B. Allison.  Evaluating statistical methods using plasmode data sets in the age of massive public databses: an illustration using false discovery rates.  _PLoS Genetics_ , 6, 2008.
  * Vaughan et al. 2009 L. K. Vaughan, J. Divers, M. A. Padilla, D. T. Redden, H. K. Tiwari, D. Pomp, and D. B. Allison.  The use of plasmodes as a supplement to simulations: A simple example evaluating individual admixture estimation methodologies.  _Computational Statistics and Data Analysis_ , 53:1755–1766, 2009.
  * Headrick and Sawilowsky 1999 T. C. Headrick and S. S. Sawilowsky.  Simulating correlated multivariate nonnormal distributions: extending the fleschman power method.  _Psychometrika_ , 64((1)):25–35, 1999.
  * Azuero et al. 2012 A. Azuero, D. T. Redden, H. K. Tiwari, S. G. Asmelash, and C. J. Piyathilake.  A simple distribution-free algorithm for generating simulated high-dimensional correlated data with an autoregressive structure.  _Commun Stat Simul Comput_ , 41(1):89–98, 2012.
  * Fan et al. 2014 J. Fan, F. Han, and H. Liu.  Challenges of big data analysis.  _National Science Review_ , 1((2)):293–314, 2014\.
  * Fan and Li 2006 J. Fan and R. Li.  Statistical challenges with high dimensionality: feature selection in knowledge discovery.  _Proceedings of the International Congress of Mathematicians, Madrid, Spain_ , pages 595–622, 2006.
  * Johnstone and Titterington 2009 I. M. Johnstone and D. M. Titterington.  Statistical challenges of high-dimensional data.  _Phil. Trans. R. Soc. A_ , 367:4237–4253, 2009.
  * Cattell and Jaspers 1967 R. Cattell and J. Jaspers.  A general plasmode (no. 30-10-5-2) for factor analytic exercises and research.  67-3, 1967.
  * Sokal et al. 1980 R. R. Sokal, F. J. Rohlf, E. Zang, and W. Osness.  Reification in factor analysis: a plasmode based on human physiology-of-exercise variables.  _Multivariate Behavioral Research_ , 4, 1980.
  * Irizarry et al. 2003 R. A. Irizarry, B. Hobbs, F. Collin, Y. D. Beazer-Barclay, K. J. Antonellis, U. Scherf, and T. P. Speed.  Exploration, normalization, and summaries of high density oligonucleotide array probe level data.  _Biostatistics_ , 4:249–264, 2003.
  * Tibshirani 2006 R. Tibshirani.  A simple method for assessing sample sizes in microarray experiments.  _BMC Bioinformatics_ , 7, 2006.
  * Franklin et al. 2015 J. M. Franklin, W. Eddings, R. J. Glynn, and S. Schneeweiss.  Regularized regression versus the high-dimensional propensity score for confounding adjustment in secondary database analyses.  _American Journal of Epidemiology_ , 182((7)):651–659, 2015.
  * Franklin et al. 2017 J. M. Franklin, W. Eddings, P. C. Austin, E. A. Stuart, and S. Schneeweiss.  Comparing the performance of propensity score methods in healthcare database studies with rare outcomes.  _Statistics in Medicine_ , 36:1946–1963, 2017.
  * Karim et al. 2018 M. E. Karim, M. Pang, and R. W. Platt.  Can we train machine learning methods to outperform the high-dimensional propensity score algorithm.  _Epidemiology_ , 2018.
  * Desai et al. 2019 R. J. Desai, R. Wyss, Y. Abdia, S. Toh, M. Johnson, H. Lee, S. Karami, J. M. Major, M. Nguyen, S. V. Wang, J. M. Franklin, and J. J. Gagne.  Evaluating the use of bootstrapping in cohort studies conducted with 1:1 propensity score matching - a plasmode simulation study.  _Pharmacoepidemiol Drug Saf_ , 28:879–886, 2019.
  * Ripollone et al. 2019 J. E. Ripollone, K. F. Huybrechts, K. J. Rothman, R. E. Ferguson, and J. M. Franklin.  Evaluating the utility of coarsened exact matching for pharmacoepidemiology using real and simulated claims data.  _Practice of Epidemiology_ , 189((6)):613–622, 2019.
  * Liu et al. 2019 S.-H. Liu, S. A. Chrysanthopoulou, Q. Chang, J. N. Hunnicutt, and K. L. Lapane.  Missing data in marginal structural models: A plasmode simulation study comparing multiple imputation and inverse probability weighting.  _Medical Care_ , 57((3)):237–243, 2019.
  * Conover et al. 2021 M. M. Conover, K. J. Rothman, T. Stürmer, A. R. Ellis, C. Poole, and M. J. Funk.  Propensity score trimming mitigates bias due to covariate measurement error in inverse probability of treatment weighted analyses: A plasmode simulation.  _Statistics in Medicine_ , 40:2101–2112, 2021.
  * Wyss et al. 2021 R. Wyss, S. Schneeweiss, M. van der Laan, S. D. Lendle, C. Ju, and J. M. Franklin.  Using super learner prediction modeling to improve high-dimensional propensity score estimation.  _Epidemiology_ , 2021.
  * Hafermann et al. 2022 L. Hafermann, N. Klein, G. Rauch, M. Kammer, and G. Heinze.  Using background knowledge from preceding studies for building a random forest prediction model: A plasmode simulation study.  _Entropy_ , 24, 2022.
  * Rodriguez et al. 2022 P. J. Rodriguez, D. L. Veenstra, P. J. Heagerty, C. H. Goss, K. J. Ramos, and A. Bansal.  A framework for using real-world data and health outcomes modeling to evaluate machine learning based risk prediction models.  _Value Health_ , 25(3):350–358, 2022.
  * Elobeid et al. 2009 M. A. Elobeid, M. A. Padilla, T. McVie, O. Thomas, D. W. Brock, B. Musser, K. Lu, C. S. Coffey, R. A. Desmond, M.-P. St-Onge, K. M. Gadde, S. B. Heymsfield, and D. B. Allison.  Missing data in randomized clinical trials for weight loss: Scope of the problem, state of the field, and performance of statistical methods.  _PLoS One_ , 4((8)), 2009.
  * Reeb et al. 2015 P. D. Reeb, S. J. Bramardi, and J. P. Steibel.  Assessing dissimilarity measures for sample based hierarchical clustering of rna sequencing data using plasmode datasets.  _PLOS Genetics_ , 2015.
  * Ejima et al. 2020 K. Ejima, A. W. Brown, D. L. Smith, U. Beyaztas, and D. B. Allison.  Murine genetic models of obesity: type i error rates and the power of commonly used analyses as assessed by plasmode-based simulation.  _International Journal of Obesity_ , 44:1440–1449, 2020\.
  * Alfaras et al. 2021 I. Alfaras, K. Ejima, C. V. L. Teixeira, C. D. Germanio, S. J. Mitchell, S. Hamilton, L. Ferrucci, N. L. Price, D. B. Allison, M. Bernier, and R. de Cabo.  Empirical versus theoretical power and type i error (false-positive) rates estimated from real murine aging research data.  _Cell Reports_ , 36(7), 2021.
  * Bickel et al. 1997 B. J. Bickel, F. Götze, and W. R. van Zweet.  Resampling fewer than n observations: gains, losses, and remedies for losses.  _Statistica Sinica_ , 7:1–31, 1997.
  * Bickel and Sakov 2008 B. J. Bickel and A. Sakov.  On the choice of m in the m out of n bootstrap and confidence bounds for extrema.  _Statistica Sinica_ , 18:967–985, 2008.
  * Andrews and Buchinsky 2000 D. W. K. Andrews and M. Buchinsky.  A three-step method for choosing the number of bootstrap repetitions.  _Econometrica_ , 68:23–51, 2000.
  * Davidson and MacKinnon 2000 R. Davidson and J. G. MacKinnon.  Bootstrap tests: How many bootstraps?  _Econometric Reviews_ , 19:55–68, 2000.
  * Politis et al. 1999 D. Politis, J. Romano, and M. Wolf.  _Subsampling_.  Springer, 1999.
  * Bickel and Ren 2001 B. J. Bickel and J.-J. Ren.  The bootstrap in hypothesis testing.  _Lecture Notes-Monograph Series, State of the Art in Probability and Statistics_ , 36:91–112, 2001.
  * Andrews and Guggenberger 2010 D. W. K. Andrews and P. Guggenberger.  Asymptotic size and a problem with subsampling and with the m out of n bootstrap.  _Econometric Theory_ , 26:426–468, 2010.
  * Gerard 2020 D. Gerard.  Data-based rna-seq simulations by binomial thinning.  _BMC Boinformatics_ , 21:206–220, 2020.
  * Ju et al. 2019 C. Ju, R. Wyss, J. M. Franklin, S. Schneeweiss, J. Häggström, and M. J. van der Laan.  Collaborative-controlled LASSO for constructing propensity score-based estimators in high-dimensional data.  _Statistical Methods in Medical Research_ , 28:1044–1063, 2019.
  * Beran and Srivastava 1985 R. Beran and M. S. Srivastava.  Bootstrap tests and confidence regions for functions of a covariance matrix.  _The Annals of Statistics_ , 13:95–115, 1985.
  * 47 P. J. Bickel and D. A. Freedman.  Bootstrapping regression models with many parameters.  In _A Festschrift for Erich L. Lehmann in Honor of his Sixty-fith Birthday_.
  * Mammen 1993 E. Mammen.  Bootstrap and wild bootstrap for high dimensional linear models.  _The Annals of Statistics_ , 21:255–285, 1993.
  * Karoui and Purdom 2018 N. E. Karoui and E. Purdom.  Can we trust the bootstrap in high-dimensions? the case of linear models.  _Journal of Machine Learning Research_ , 19:1–66, 2018\.
  * Hoerl and Kennard 1970 A. E. Hoerl and R. W. Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  _Technometrics_ , 12:55–67, 1970.
  * Searle et al. 1992 S. R. Searle, G. Casella, and C. E. McCulloch.  _Variance Components_.  Wiley Interscience, New Jersey, 1992.
  * Covarrubias-Pazaran 2017 G. Covarrubias-Pazaran.  Genome assisted prediction of quantitative traits using the R package sommer.  _PLoS ONE_ , 11:1–15, 2017.
  * Weinstein et al. 2013 J. Weinstein, E. Collisson, and et. al.  The cancer genome atlas pan-cancer analysis project.  _Nature Genetics_ , 45:1113–1120, 2013.  URL https://doi.org/10.1038/ng.2764.
  * Colaprico et al. 2015 A. Colaprico, T. C. Silva, C. Olsen, L. Garofano, C. Cava, D. Garolini, T. Sabedot, T. M. Malta, S. M. Pagnotta, I. Castiglioni, M. Ceccarelli, G. Bontempi, and H. Noushmehr.  Tcgabiolinks: An r/bioconductor package for integrative analysis of tcga data.  _Nucleic Acids Research_ , 2015.  doi: 10.1093/nar/gkv1507.  URL http://doi.org/10.1093/nar/gkv1507.
  * Silva et al. 2016 Silva, T. C, Colaprico, Antonio, Olsen, Catharina, D’Angelo, Fulvio, Bontempi, Gianluca, Ceccarelli, Michele, Noushmehr, and Houtan.  Tcga workflow: Analyze cancer genomics and epigenomics data using bioconductor packages.  _F1000Research_ , 5, 2016.
  * Mounir et al. 2019 Mounir, Mohamed, Lucchetta, Marta, Silva, T. C, Olsen, Catharina, Bontempi, Gianluca, Chen, Xi, Noushmehr, Houtan, Colaprico, Antonio, Papaleo, and Elena.  New functionalities in the tcgabiolinks package for the study and integration of cancer data from gdc and gtex.  _PLoS computational biology_ , 15(3):e1006701, 2019.
  * Ritchie et al. 2015 M. E. Ritchie, B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, and G. K. Smyth.  limma powers differential expression analyses for RNA-sequencing and microarray studies.  _Nucleic Acids Research_ , 43(7):e47, 2015.  doi: 10.1093/nar/gkv007.
  * Ledoit and Wolf 2004 O. Ledoit and M. Wolf.  A well-conditioned estimator for large-dimensional covariance matrices.  _Journal of Multivariate Analysis_ , 88(2):365–411, 2004.  doi: https://www.sciencedirect.com/science/article/pii/S0047259X03000964.
  * Tibshirani 1996 R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society B_ , 58:267–288, 1996.