Mosaic Monte Carlo: A new method of simulation design to improve the generalizability of findings

1 Introduction

Monte Carlo simulation studies are the bread and butter of testing the performance of statistical methods. They are often implemented by choosing a number of factors to manipulate – such as sample size, magnitude of correlations, and model size – and then randomly generating data under the corresponding data-generating model(s) across a large number of replications. Research on topics such as pharmacology, structural equation modeling, quality control, and many others use Monte Carlo approaches to simulate models, compute statistical precision and accuracy, and evaluate new methods (Curran, Bollen, Paxton, Kirby, & Chen, 2002; Goutelle et al., 2009; Iranmanesh, Parchami, & Sadeghpour Gildeh, 2022).

Various concerns about the design and implementation of Monte Carlo simulation studies have been raised in recent years. Choosing the appropriate sample size range is one potential challenge (Paxton, Curran, Bollen, Kirby, & Chen, 2001; Skrondal, 2000; Spence, 1983). Other concerns include non-converged results and appropriate selection of software (Paxton et al., 2001; Spence, 1983). Concerns and guidance about best practices regarding the design and implementation of Monte Carlo simulation studies more broadly have been discussed by Morris, White, and Crowther (2019) and Siepe et al. (2024). Morris et al. (2019) propose an ADEMP structure (an acronym for Aims, Data-generating mechanism, Estimands, Methods, and Performance measures) to aid in planning simulation studies, as shown in Table 1. The ADEMP structure is implemented in modern statistical simulation studies to ensure reproducibility. Siepe et al. (2024) expand on this framework tailored to an audience of psychology researchers. Burton, Altman, Royston, and Holder (2006), who defined the essential preparatory procedures for simulations in medical statistics, historically contributed to the shift toward these standardized practices.

Another issue that has garnered attention is the number of replications used in Monte Carlo simulation studies (Schaffer & Kim, 2007). Brooks (2002) states that if not enough replications are utilized in simulations, the results could be misleading. A few studies (Preecha, 2004; Supawan, 2004) suggest guidelines for the number of replications for common statistical analyses such as ANOVA and regression. However, a recent development that has gained traction is the use of so-called Monte Carlo Standard Errors(MCSEs; Koehler, Brown, & Haneuse, 2009). Koehler et al. (2009) propose using these to quantify the uncertainty in Monte Carlo designs. In their article, the authors found that there can be unexpectedly substantial uncertainty in results even with 1000 replications. The authors suggest that several thousand replications may be needed at times in order to obtain good precision of estimates. Of course, the number of replications needed for a particular simulation study will be context-dependent, but MCSEs can be helpful for determining how many replications are adequate. Monte Carlo standard errors are now recommended to be part of the standard output reported in simulation studies (Morris et al., 2019; Siepe et al., 2024).

In contrast, Skrondal (2000) argues that one of the main issues plaguing simulation studies is external validity rather than precision. Boulesteix, Stierle, and Hapfelmeier (2015) state that methodological research remains deeply vulnerable to publication bias, which is the systematic suppression of negative data in favor of over-optimistic conclusions. In addition, standard computational literature frequently suffers from optimization bias, in which authors evaluate their own procedures under non-neutral, overly favorable conditions (Boulesteix, Lauer, & Eugster, 2013). A main contributing factor of this issue is the choice of the data-generating mechanism, which serves as the primary component affecting variation across simulation conditions (Boulesteix et al., 2020).

However, it is common practice in Monte Carlo simulations to choose few sets of population parameters, especially in cases when they are not the topic of interest and not thought to influence the performance of the statistical method being studied. The value of parameters are often chosen arbitrarily or without justification, as in Collins, Schafer, and Kam (2001); Gomer, Jiang, and Yuan (2019); Hu and Bentler (1999), and many more. Alternatively, they may be based on values from a real dataset, as in Ke and Wang (2015). However, neither of these approaches can guarantee that findings will generalize to other data-generating models. Not only is there potential for misspecification of population parameters according to Paxton et al. (2001), but Kulinskaya, Hoaglin, and Bakbergenuly (2021) show that even modifying subtle assumptions in data generation (i.e., such as shifting between fixed and random intercept models) can also affect the apparent performance and ranking of statistical methods.

Morris et al. (2019) explicitly ask whether a chosen data-generating mechanism may favor some methods over others, and how the choice of data-generating mechanism can be checked and justified. However, even when researchers do justify such choices, some areas of study may be particularly sensitive to the choice of parameter values, such as missing data analysis and statistical approaches that use Bayesian methodology. For example, consider the simulation studies in Gomer and Yuan (2021) involving a multiple regression model. In Study 1, regression coefficients were set to either \((0, 0, 0, 0)'\) or \((0, 1, 1, 1)'\) and predictors were set to be correlated or uncorrelated. When comparing two versions of the missing not at random (MNAR) mechanism (focused MNAR and diffuse MNAR 1a), it would appear that one systematically yields more biased results than the other across these conditions. However, in a small simulation study in the supplemental material, the regression coefficients are arbitrarily set to either \((0, 0, 0, 0)'\) or \((0, 2, 2, 2)'\) and the pattern reverses. The difference in coefficient values between Study 1 and the supplemental material appears trivial. It would be hard to anticipate that such a change would be influential, especially since the relationships between variables are all the same. If this paper had been smaller in scope and focused only on Study 1, it would have been easy to draw an entirely wrong conclusion from the findings.

Critically, this is an example in the literature where the choice of parameter values can play a significant role in the overall conclusion of an article, even though they may not seem related to the topic being studied. Methodologists design simulation studies with the belief that simulation factors could have systematic effects on the outcome of interest. However, there may be other scenarios, particularly in missing data analysis and Bayesian contexts, in which small adjustments to parameter values influence simulation results in ways that are nonsystematic. We will show an example of this in a later section. If these factors are not varied in a simulation design then generalizability will be threatened and recommendations become questionable. However, interpreting findings that vary nonsystematically is another difficulty.

Franklin, Schneeweiss, Polinski, and Rassen (2014) and Schreck, Slynko, and Saadati (2024) show an alternative methodological approach to addressing the external validity of simulation studies using plasmode simulations. In a plasmode simulation, real empirical datasets are used to anchor the data-generating mechanism by resampling observed covariate data (Franklin et al., 2014). This preserves the complex, real-world dependence structures that are often difficult to synthesize parametrically (Schreck et al., 2024). An investigator-specified outcome-generating model is then applied to the resampled data to introduce a known “truth,” such as a specific treatment effect or effect size. While plasmode simulation achieves realism by grounding the simulation in a complex empirical reality, its generalizability remains bounded by the specific characteristics and representativeness of the chosen source sample (Schreck et al., 2024).

Nevertheless, it is difficult to test multiple data-generating models and parameters with Monte Carlo simulations due to computational burden and the need to summarize results. For example, let us consider a basic regression Monte Carlo simulation with 4 sample sizes, 4 probability distributions, 4 sets of regression coefficients, and 4 sets of predictor correlations. With 1,000 replications, this simulation requires calculations on 256,000 datasets. If the statistical method being studied requires iteration (i.e., multiple imputation, Bayesian estimation) then this design is already computationally intensive. What’s more, this design still considers a very limited set of data-generating models and still faces issues of generalizability. There are also 256 distinct simulation factors to sift through and summarize in a meaningful way. In many cases, it is simply not feasible to consider a wider range of data-generating models.

Skrondal (2000) suggests expanding experimental factors and using fewer replications in order to address the issue of computational complexity. He also illustrates how to use fractional factorial designs in order to save computational time. However, issues related to summarizing large sets of results and nonsystematic relationships remain. Leigh and Bryant (2015) focus on a so-called “Bayesian style” approach to selecting parameter values in simulation studies instead of grid-based approaches to address issues of computational demand. In their approach, a posterior distribution of parameters is reported conditional on a particular simulation outcome. This also sidesteps issues related to summarizing results. In their approach, parameters are sampled using an algorithm based on importance sampling which selects parameters adaptively rather than over an entire grid. The suggestions in both articles are aimed at solving the generalizability issue of Monte Carlo simulation studies without placing significant burdens on computational complexity.

In this article, we test a new type of simulation design with the goal of improving the generalizability of Monte Carlo simulation studies. We call this method “Mosaic Monte Carlo.” In art, a mosaic is a pattern or an image that is created from many small pieces. In Mosaic Monte Carlo, the idea is to break up the replications of a given simulation condition into smaller subsets, each using a different data-generating model. For example, instead of using the same data-generating model across 1,000 replications, we use 20 data-generating models with 50 replications each. The data-generating models act like the small mosaic pieces in the art metaphor. In the context of a multiple regression model, we might choose to vary the values of the regression coefficients and the predictor correlations. The performance measures calculated from the data-generating models – for example, bias, variance, and type I error rates – can then be summarized and presented for the simulation condition of interest, which is similar to how small mosaic pieces are used to form an overall image or pattern in art.

The primary purpose of Mosaic Monte Carlo is to account for factors that have a nonsystematic impact on the outcome(s) of interest, thus improving generalizability of simulation studies. While researchers often can generate a list of experimental factors that have systematic impacts on their performance measures before conducting their simulation study, it may be more difficult to predict what aspects of data-generating mechanisms lead to nonsystematic impacts. Thus, we propose using Mosaic Monte Carlo as a two-stage procedure involving a screening phase and a simulation phase. The screening phase helps researchers to identify which aspects of data-generating mechanisms have an impact on their performance measures and whether these effects are systematic or nonsystematic. The simulation phase incorporates this information, conducts the simulation study of interest, and summarizes the results.

In a similar vein as Skrondal (2000), we propose using fewer replications; however, we do so across factors that are not experimental conditions in order to retain precision across those dimensions. An advantage of Mosaic Monte Carlo over that of Leigh and Bryant (2015) is that our method does not require knowledge of Bayesian statistics. Mosaic Monte Carlo also is designed to account for nonsystematic relationships between parameter values and the outcome(s) of interest, whereas the approach suggested by Leigh and Bryant (2015) assume the presence of a systematic relationship. An advantage of Mosaic Monte Carlo over plasmode simulations (Franklin et al., 2014; Schreck et al., 2024) is that while both plasmode simulation and Mosaic Monte Carlo share the fundamental goal of improving the generalizability of methodological findings beyond narrow, arbitrary parameter specifications, Mosaic Monte Carlo utilizes a broad, synthetic parameter space, rather than anchoring to one empirical dataset. This strategy mitigates the nonsystematic impacts of arbitrary parameter choices without relying on real-world data availability.

The remainder of this article is organized as follows. First, we describe how to implement Mosaic Monte Carlo in practice. Then, we perform a small demonstration that highlights how Mosaic Monte Carlo can be used in a context with missing data. Next, we present the methods and results for a simulation study that tests the baseline performance of Mosaic Monte Carlo. We also provide an empirical example that uses Mosaic Monte Carlo to replicate findings from a simulation study in the context of structural equation modeling. We conclude with a discussion of our main findings and recommendations.

2 How to implement Mosaic Monte Carlo

In this section, we describe how Mosaic Monte Carlo can be implemented in practice as a two-stage procedure. Some example code is also provided in a github repository which can serve as a reference. As a note, the code is written for clarity rather than computational efficiency. This procedure is applicable across modeling contexts but we have written the guidance under the assumption that at least some of the estimands of interest are parameter estimates and thus that MCSEs can be calculated. However, Mosaic Monte Carlo can still be used for other targets of simulation studies, as we illustrate later in our empirical example.

2.1 Terminology

Throughout the remainder of this paper, we will use terminology consistent with that used in the ADEMP structure developed by Morris et al. (2019). Aspects of data generation that are varied in Mosaic Monte Carlo are called mosaic pieces. In the screening phase, these are referred to as candidate mosaic pieces until the researcher verifies that these aspects of data generation impact results nonsystematically.

2.2 Stage 1: Screening phase

The purpose of the screening phase in Mosaic Monte Carlo is to help researchers distinguish factors that have systematic vs. nonsystematic effects on their simulation results and choose how to allocate replications across data-generating models. The flowchart in Figure 1 breaks this process up into 3 steps. We will discuss each of these steps in the subsections below.

Step 1: Setting up the screening phase

In the first step in Figure 1(a), the researcher makes a series of decisions about the design of his simulation study. Researchers following the ADEMP framework for planning simulations (Morris et al., 2019) will complete all steps before conducting the screening phase of Mosaic Monte Carlo. Then, researchers revisit the data-generating mechanisms step and consider whether these factors are expected to influence any targets of the simulation study. Unless a particular factor is expected to systematically impact results, it should be treated as a candidate mosaic piece and investigated further rather than being held as fixed. When choosing population parameters, the following categories are often appropriate candidate mosaic pieces in Mosaic Monte Carlo:

• Population parameters that are not essential to the research aims. An example of this could be the measurement error of variables in a simulation study examining the impact of model misspecification in the context of confirmatory factor analysis. A second example could be the difficulty parameter of items in a simulation study seeking to identify the breaking point of estimation for long tests given to small sample sizes.
• Population parameters that help address the research aims but whose values can be chosen arbitrarily. An example of this could be the value of regression coefficients in a simulation study examining bias of their estimates in the context of missing data. Another example could be the predictor correlations between variables in a simulation study examining the effects of multicollinearity. In both cases, the specific values of the population parameters are related to the topic of study but the researcher hopes to generalize findings beyond the specific parameter values chosen. In such cases, selecting multiple parameter values helps strengthen the findings.

Depending on the topic of interest, other aspects of data generation that are not population parameters could be candidate mosaic pieces. For example, a researcher examining the performance of a missing-data method with data that is Missing At Random (MAR) would need to choose how to create missing values. She may choose an interval selection method that removes values on a variable \(Y\) if the corresponding value on a second variable \(X\) is above a certain threshold. The value on \(X\) usually is chosen to obtain a specific percentage of missing values, but the researcher could vary which variable serves as the \(X\) variable if there are multiple candidates.

The next decision to be made in step 1 of the screening phase is to choose how to generate the mosaic pieces. In the case of population parameters, these should be generated such that they cover a broad but realistic range of values. A familiar method of doing this would be traditional random sampling. Uniform distributions are a great choice for parameter generation because the resulting parameters will be evenly spread across the range of possible values. In contrast, generating parameters from a normal distribution will result in values that are more tightly clustered towards the center, which is less desirable from a generalizability standpoint. However, other sampling methods such as Latin hypercube sampling and low discrepancy sequences can also be used to generate mosaic pieces which ensure the desired sampling space is adequately represented in its entirety (McKay, Beckman, & Conover, 2000; Tuffin, 1996). Another alternative is to expand plasmode simulation to include many empirical datasets, although such a task may be time-intensive.

If random sampling is used, the corresponding parameters for the distribution itself (analogous to hyperparameters in Bayesian statistics) must be chosen after the generating distribution has been selected. In the case of the uniform distribution, this means specifying meaningful boundaries. For example, correlations should occur on the interval [-1, 1] in order to be valid. If a researcher is interested in the effects of multicollinearity, the interval could be narrowed to [-1, -0.7] and [0.7, 1]. This allows the researcher to specify parameter values according to the topic of interest while still making use of Mosaic Monte Carlo. If a different probability distribution is preferred over the uniform distribution, our general recommendation would be to specify the distribution in such a way that there is as large a variance as possible. This will improve the generalizability of the study. However, care must be taken that this range is still representative of the realistic application space and does not draw in implausible regions. What constitutes a realistic application space will vary by context but could include things such as ensuring that the smallest factor loadings in a structural equation model meet guidelines for interpretability according to the criteria in Stevens (1992).

As mentioned previously, the idea of Mosaic Monte Carlo is not limited to parameter generation, although that is our main focus. When other aspects of simulation studies are used to form the mosaic, choices about how to generate the mosaic pieces may be more systematic. Revisiting our example of the researcher investigating the performance of a missing-data method with MAR data, she may only have 4 variables that could play the role of \(X\). In that case, she may alternate between all 4 possible candidates.

As a note, multiple types of population parameters and/or aspects of the simulation study can be varied simultaneously to create the mosaic. In other words, each data-generating model is likely to be constructed from several mosaic pieces, such as regression coefficients and correlation values. This will help to reduce the computational costs involved in the screening phase.

The last decision the researcher needs to make in this step is the initial number of data-generating models (\(J\)) and replications (\(K\)) to use in the screening simulation. At this point, the researcher should choose a large enough \(K\) to obtain reliable estimates of his performance measures and a large enough \(J\) to reliably calculate a standard deviation and to have sufficient sample size for analyses in step 3. We suggest using a minimum of 10 for both \(J\) and \(K\). However, it is important to note that the starting value of \(K\) will be updated in the next step of the screening phase and that the value of \(J\) used in the screening phase can be increased in the simulation phase.

Step 2: Implementing the screening simulation

In this step shown in Figure 1(b), the researcher conducts a small-scale simulation study that follows the same design of the simulation study he intends to run. This step determines which candidate mosaic pieces actually impact results and which can be treated as fixed later on in the simulation phase. Simulation factors and data-generating mechanisms already known not to interact with candidate mosaic pieces (i.e., based on previous studies in the literature) can be held constant in the screening simulation. In contrast, if these could have possible interaction effects or the effects are unknown, researchers can select several or all values of each factor/data-generating mechanism and run the screening simulation in a fully-crossed design. We acknowledge that this may not be computationally feasible in some applications. In such cases, we recommend researchers select values of factors/data-generating mechanisms that are most likely to elicit effects in the candidate mosaic pieces. For example, this could be the smallest sample size or the largest percentage of missingness.

When programming the simulation study, some care must be taken to ensure that the mosaic is implemented correctly. We will give some brief guidance on how this might be achieved using two alternative strategies. These strategies work in both stages of Mosaic Monte Carlo. The code provided in our supplemental materials provide examples under both strategies.

The first strategy for implementing Mosaic Monte Carlo is to create the mosaic pieces and data-generating models first, and then run the simulation study in a second and separate step. However, the specific parameter values used for candidate mosaic pieces should be saved so that they can be referenced later. For example, these values would be used for calculations of bias.

The second strategy is to create the mosaic pieces and data-generating models while running the simulation study simultaneously. For example, if iterating over looped code, candidate mosaic pieces would be randomly generated inside the loop, incorporated into the data-generating model, and then used to create the datasets. In this strategy, parameter values for candidate mosaic pieces do not need to be saved if all the desired output is also calculated inside the loop.

When conducting the screening simulation, researchers should not only compute the values of performance measures but also calculate their MCSEs for each distinct data-generating model. These will be used to determine if the number of replications is adequate. The MCSEs can be calculated following the guidelines in Koehler et al. (2009) in order to obtain \(J\) separate MCSE values for each performance measure. Additionally, the coefficient of variation (CV) for each performance measure across the \(J\) distinct data-generating models should be calculated. This will be used to determine whether the candidate mosaic pieces actually impact results or if they can be treated as fixed. The formula for the CV in this context is provided below:

where PM stands for performance measure, \(\overline {PM}_j = 1/J\sum _{j=1}^J PM_j\), and the values of the performance measures are calculated in the usual way for each distinct data-generating model.

After the screening simulation has been run, the researcher examines the MCSEs for each performance measure and determines if these values fall below the desired screening threshold. We cannot recommend a specific value for this screening threshold as MCSEs will vary by performance measure and context. However, MCSEs should be small enough to have confidence in using the performance measures to calculate CVs across the \(J\) data-generating models and effect sizes in step 3 of the screening phase. For data-generating models that yield MCSEs exceeding the screening threshold, the number of replications should be increased until acceptable MCSEs are obtained for each performance measure. This means that some data-generating models will have greater numbers of replications than others. This adaptive approach helps direct computational resources more efficiently.

Once acceptable MCSEs have been obtained in all data-generating models for all performance measures, the researcher next determines which candidate mosaic pieces actually impact results using CVs. A CV value less than 0.1 suggests that the performance measure has low variability across data-generating models and thus is relatively unaffected by the candidate mosaic pieces. If this is true for all performance measures, the candidate mosaic pieces can be treated as fixed in the simulation phase of Mosaic Monte Carlo and the screening phase ends without proceeding to step 3. CV values greater than 0.1 suggest that the candidate mosaic pieces do have an impact on results and these should be examined further in step 3 of the screening phase. In cases when performance measures have a target value of 0 (i.e., relative bias), CV values will be unduly inflated because the means are so small. In such cases, the SD of the performance measure can be used instead of the CV with the same 0.1 cutoff.

Step 3: Distinguishing systematic from nonsystematic effects and preparing for the main simulation study

The last step in the screening phase of Mosaic Monte Carlo makes use of effect sizes in order to determine which candidate mosaic pieces have systematic vs. nonsystematic effects on results. As shown in Figure 1(c), this step also helps the researcher decide the numbers of replications and data-generating models to use in the simulation phase of Mosaic Monte Carlo. In this step, the researcher first calculates the effect sizes according to the guidance in Table 2 for each performance measure using the output of the screening simulation. The sample size for such analyses is \(J\), the number of data-generating models. Referring to Table 2, performance measures such as bias and variance can be used as outcome variables in a multiple regression model with candidate mosaic pieces as predictors and \(R^2_{adj}\) serving as a multivariate effect size. If the value of \(R^2_{adj}\) exceeds the suggested cut-off, this suggests that one or more candidate mosaic pieces have a systematic impact on the performance measure and the values of individual regression coefficients should be examined. Candidate mosaic pieces corresponding to large regression coefficients likely have systematic impacts on the results and should be treated as full simulation factors in the simulation phase. In contrast, candidate mosaic pieces corresponding to small regression coefficients likely have nonsystematic impacts on results and thus should be treated as mosaic pieces in the simulation phase.

In this example, what constitutes a “large” vs. “small” regression coefficient depends upon the scale of the performance measure and the candidate mosaic pieces. It is thus up to the discretion of the researcher. However, we note that statistical significance should not be considered when making these decisions, as it is directly influenced by the number of data-generating models and thus is not a reliable criterion. As a note, candidate mosaic pieces that have nonsystematic effects on one performance measure but systematic effects on another should be treated as a full simulation factor in the simulation phase.

After determining which candidate mosaic pieces will remain mosaic pieces in the simulation phase and which will become full simulation factors, the researcher re-evaluates the MCSE values from the output of the screening simulation. Now, the question is whether or not these MCSEs are less than the desired simulation threshold. The simulation threshold should be small enough to obtain adequate precision in the final reporting of simulation results, which is likely to be more strict than the number used for the screening threshold. For data-generating models whose MCSEs already meet this stricter criterion, the same number of replications can be used in the simulation phase. The number of replications should be increased for data-generating models whose MCSEs do not yet meet the simulation threshold until acceptable values are obtained. Once acceptable MCSEs have been obtained for all data-generating models and performance measures, the researcher should decide how many data-generating models can be computationally afforded. In cases when the researcher can afford greater \(J\), MCSEs can be calculated for the new set of data-generating models and compared to the simulation threshold in order to determine appropriate \(K\). Alternatively, the researcher could use the largest \(K\) obtained from the previous data-generating models on the new set.

2.3 Stage 2: Simulation phase

Conducting the simulation study in the second stage of Mosaic Monte Carlo likely involves making minor expansions to the screening simulation in the first stage. Actually, the output from step 3 of the screening phase with finalized \(K\) values can be re-used. The simulation study will need to be run anew for the set of new data-generating models and for any simulation conditions that were omitted from the screening phase.

The final task in the simulation phase of Mosaic Monte Carlo is to summarize the results across data-generating models. There are many possible approaches that could be used and the best approach will depend on the number of simulation conditions and the nature of the mosaic pieces. Some possible approaches include: 1) distributional summaries, 2) aggregation, and 3) frequencies.

Distributional summaries, such as boxplots, can be especially useful for gauging several aspects of performance simultaneously and provide nuanced information. For example, side-by-side boxplots comparing the relative bias of parameter estimates across several estimation methods can illustrate the accuracy of methods by the center of the boxplots. Reliability of each method can be assessed by the range, interquartile range, and presence of outliers. The overlap in boxplots between methods can indicate how advantageous one method is compared to the others. Aggregation can be useful to conserve space when less nuance is required or when the scope is relatively broad, such as constructing overall summaries or making general recommendations. However, care must be taken to aggregate only over mosaic pieces that impact results nonsystematically according to the screening stage. Frequencies can be a useful summary for binary performance measures such as nonconvergence rates or counting the number of times a method was identified as the best-performing method among several candidates. Since many methodological studies seek to answer multiple interrelated questions, it is likely that several summarizing strategies could be used together in one paper.

When aggregation is used as a summarizing strategy, calculations of bias and variance must be adjusted if the corresponding parameters form part of the mosaic. The modified formula for calculating aggregated empirical bias of an estimator \(\widehat {\theta }\) for \(\theta \) is given by

where \(j\) refers to a single data-generating model with \(J\) distinct data-generating models in total (i.e., the number of distinct data-generating models created from the mosaic) and \(k\) denotes a single replication with \(K\) replications in total. The formula for calculating aggregated empirical relative bias is given by

Supposing that \(\overline {\theta }_j = 1/k\sum _{k=1}^K \theta _{jk}\), Equations 1 and 2 will be unbiased estimators of the aggregate (relative) bias as long as all \(\overline {\theta }_j - \theta _j\) are unbiased estimators of bias.

If standard errors or variance are part of the desired output for the simulation study, the corresponding population values should be fixed to a constant to ensure the data-generating models are comparable. To calculate the empirical variance of an estimator \(\widehat {\theta }\) for a \(\theta \) that is varied in the mosaic, the following modified formula can be used:

where \(\widehat {\theta }_j\) is the mean of the estimates across all \(K\) replications for a given data-generating model \(j\), \(\widehat {\theta }_j = \frac {1}{K}\sum _{k=1}^K \widehat {\theta }_{jk}\). In other words, the user should calculate the variance of the estimates across all replications for a given data-generating model. Then, these variances should be averaged across all data-generating models \(J\) to obtain a single value. Since \(1/(K-1)\sum _{k=1}^K(\widehat {\theta }_{jk} - \widehat {\theta }_j)^2\) is the formula for the sample variance of \(\widehat {\theta }_{jk}\) values, Equation 3 inherits unbiasedness and consistency from the sample variance (Casella & Berger, 2001). However, despite these desirable finite-sample properties, each individual sample variance can be considerably noisy for small values of \(K\). Equation 3 is also a mean so it is sensitive to outliers that may arise under such conditions.

Please note that the user should not simply take the variance of the \(\widehat {\theta }_{jk}\) values across all replications and across all data-generating models simultaneously. When \(\theta \) is part of the mosaic, this practice falsely treats the \(\widehat {\theta }_{jk}\) values as though they are comparable across different data-generating models. More formally, and assuming that all data-generating models have the same \(K\) for simplicity,

where \(\overline {\theta }\) is the mean of all \(\widehat {\theta }_{jk}\) values, \(\overline {\theta } = 1/JK \sum _{j=1}^J\sum _{k=1}^K \widehat {\theta }_{jk}\). After some mathematical manipulation expanding the sum of squares, the nonequivalency between the two formulas is even more apparent:

2.4 Considering computational burden

Implementing Mosaic Monte Carlo designs will add computational burden to simulation studies. For researchers who already adapt the number of replications to obtain acceptable MCSEs, the additional burden comes from the screening phase of Mosaic Monte Carlo. We acknowledge that in some applications, it is simply not feasible to follow the procedure outlined in Figure 1 in its entirety. Compromises may need to be made on the number of factors varied in the mosaic, the number of data-generating models, and the number of replications. Compromises can also be made by reducing the levels of simulation factors and/or holding some simulation factors constant. It is important to note that attempts to implement a screening phase will improve the generalizability of findings even if screening is done on a small scale. We recommend to make the screening phase as exhaustive as can be afforded, while acknowledging that this may deviate from the guidance in the flowchart.

3 Demonstration of Mosaic Monte Carlo

This section contains a brief demonstration of Mosaic Monte Carlo and how it can be used when aspects of data generation have nonsystematic impacts on results. We hope this provides readers with additional clarity on the usefulness of Mosaic Monte Carlo and how it can be implemented in practice. The R code for this demonstration is available in a github repository.

Continuing discussion of Gomer and Yuan (2021) as an example, suppose we want to see if one MNAR subtype consistently leads to more biased regression coefficients than the other in a similar setting. Specifically, our research aims are to determine 1) whether diffuse MNAR leads to more biased regression coefficient estimates and 2) how often this occurs. We investigate the relative bias of the regression coefficient associated with missing values (\(\widehat {\beta }_3\) in the original article). We calculate the relative bias of \(\widehat {\beta }_3\) under both focused and diffuse MNAR (corresponding to diffuse MNAR 1a in the original article) when maximum likelihood is used as a missing-data method.

3.1 Stage 1: Screening Phase

We use the same simulation design as Study 1 in Gomer and Yuan (2021) but revisit aspects of data generation for candidate mosaic pieces. Our main concern is that the values of regression coefficients and predictor coefficients may have nonsystematic impacts on the bias of \(\widehat {\beta }_3\) so following the flowchart in Figure 1(a), we treat these as candidate mosaic pieces. We choose to generate regression coefficients from a uniform distribution \(Unif(-5, 5)\) and predictor correlations from \(Unif(-0.6, 0.6)\). To begin with, we will use \(K=10\) replications for each of \(J = 60\) data-generating models. This value of \(J\) was chosen in order to have sufficient sample size for conducting the multiple regression analysis in step 3 of the screening phase assuming that there are 6 predictors (i.e., 3 correlations and 3 regression coefficient population values and following the rule of thumb that there should be 10 observations per predictor).

In step 2 of the screening phase, we run the screening simulation and calculate the MCSEs for the relative bias of \(\widehat {\beta }_3\) for each data-generating model. We also calculate the CV of the relative bias across data-generating models under both focused and diffuse MNAR. Although relative bias has a target of 0, the setting of MNAR subtypes creates levels of relative bias that are noticeably different from 0 and thus will not unduly impact CV values in this case. For the purpose of screening, we will aim to have all MCSEs less than a 0.2 screening threshold.

After adaptively increasing \(K\) for data-generating models whose MCSEs exceeded the 0.2 screening threshold, we obtained acceptable MCSEs for all data-generating models with 310 replications. The CV of the relative bias of \(\widehat {\beta }_3\) was 0.934 under focused MNAR and 4.2 under diffuse MNAR. According to the flowchart in Figure 1(b), this suggests that the candidate mosaic pieces do impact the results under both MNAR subtypes since the CV values are greater than 0.1. We thus proceed to step 3 of the screening phase.

In step 3, we conduct multiple regression analyses using the relative bias of \(\widehat {\beta }_3\) under focused and diffuse MNAR as outcome variables and the parameters \(\beta _1, \beta _2, \beta _3, \rho _{12}, \rho _{13}\), and \(\rho _{23}\) as predictors. The output for this analysis is shown in Table 3. In this table, we can see that the \(R^2_{adj}\) values are quite small and is even negative under diffuse MNAR. With poor model fit, the predictors do little to explain the variance in the bias of \(\widehat {\beta }_3\). This suggests that the candidate mosaic pieces impact the bias of \(\widehat {\beta }_3\) nonsystematically.

The final task is to refine the number of \(J\) and \(K\) for the simulation phase of Mosaic Monte Carlo. Suppose we want to use a simulation threshold of 0.1 for our MCSEs. In this case, there are several data-generating models that need additional replications to meet this stricter cutoff. After adaptively increasing \(K\) for these models, we obtain acceptable MCSEs for all data-generating models with 1010 replications. Although this number may seem large, there was only one data-generating model that required this value.

3.2 Stage 2: Simulation Phase

The output of step 3 can be used for the simulation phase of Mosaic Monte Carlo. Recall that our research aims are to determine 1) whether diffuse MNAR leads to more biased regression coefficient estimates and 2) how often this occurs. To address the first aim, we used the values of relative bias of \(\widehat {\beta }_3\) obtained under focused and diffuse MNAR across 60 data-generating models to construct boxplots and calculate Cohen’s \(d\). The boxplots are distributional summaries that provide nuanced information while Cohen’s \(d\) provides a concise aggregate summary.

The side-by-side boxplots in Figure 2 show that the range of relative bias is larger under diffuse MNAR than under focused MNAR. In fact, the relative bias of \(\widehat {\beta }_3\) rarely exceeds 10% under focused MNAR whereas this occurs more commonly under diffuse MNAR. However, we obtained \(d = 0.2699\). This is a small effect size, suggesting that the bias of \(\widehat {\beta }_3\) is not substantially worse under diffuse MNAR when we examine the trends overall. Taken together, these findings show that the bias of \(\widehat {\beta }_3\) has potential to be much larger under diffuse MNAR than focused MNAR but that this does not occur regularly enough to make a blanket statement that diffuse MNAR leads to worse bias.

To address our second aim, we counted the proportion of times that the bias of \(\widehat {\beta }_3\) was worse under diffuse MNAR compared to focused MNAR. Worse bias was obtained under diffuse MNAR 55% of the time, which translates to a Cohen’s \(h\) of 0.2. This is also a small effect size, suggesting that the bias of \(\widehat {\beta }_3\) is not worse more frequently under diffuse MNAR.

3.3 Summary

In this demonstration, Mosaic Monte Carlo was helpful for illustrating that bias is not always worse under a specific subtype of MNAR. Based on the analyses we performed, we were also able to identify that correlation and regression coefficient values had some nonsystematic impact on whether bias would be worse under diffuse MNAR compared to focused MNAR. Although diffuse MNAR has potential for leading to large levels of bias compared to focused MNAR, this did not occur often enough to warrant a bigger claim that diffuse MNAR is worse than focused MNAR.

4 Methods

In addition to the small demonstration in the previous section, a larger Monte Carlo simulation study was run in R (R Core Team, 2020) to assess possible detrimental impacts of using Mosaic Monte Carlo compared to traditional simulation design. In particular, we sought to determine whether aggregation can safely summarize information related to bias and variance estimates, especially when there is a large number of data-generating models relative to the number of replications. The code for this simulation study is provided in a github repository. We focus on our method’s impact on bias and variance because these are two popular metrics used to evaluate the performance of statistical methods in Monte Carlo simulation studies. Specifically, we sought to address the following research questions:

• RQ1: To what extent does the aggregation strategy of Mosaic Monte Carlo impact the overall bias of statistical estimates in a simulation study?
• RQ2: To what extent does the aggregation strategy of Mosaic Monte Carlo impact the overall variance of statistical estimates in a simulation study?
• RQ3: To what extent does sample size and sampling variability affect the answers to RQ1 and RQ2?

We also evaluate Mosaic Monte Carlo under relatively ideal conditions in which estimates are unbiased. This will establish the baseline performance of Mosaic Monte Carlo.

4.1 Simulation design

The data-generating models used in this study were multiple regression models with 3 predictors of the form:

where the \(\beta _j\) coefficients were treated as candidate mosaic pieces and randomly generated from \(Unif(0.25, 10)\). Here, \(i\) denotes individual data values with \(i= 1, \ldots , n\) and \(j\) denotes a specific data-generating model with \(j = 1, \ldots J\). The residuals were given by \(\varepsilon _{ij} \sim N(0, \sigma ^2_{\varepsilon })\) where \(\sigma ^2_{\varepsilon }\) is set to 9, 25, and 100. The error variance was treated as a full simulation factor in order to make sure that the randomly generated data-generating models are comparable when analyzing the impact of our method on estimates of variance. The sample size \(n\) was also treated as a full simulation factor and set to 50, 75, 100, 200, and 400. The \(X\) variables were treated as fixed across replications and initially generated from a standard multivariate normal distribution using the R package mvtnorm (Genz et al., 2020), with correlation coefficients treated as candidate mosaic pieces and randomly generated from \(Unif(0.1, 0.6)\). Note that this means there is no conceptual difference between the regression coefficients and that they are essentially exchangeable.

The total number of unique data-generating models \(J\) was set to be 1, 10, 25, 50, 100, 1000 and the number of replications for each data-generating model was set to be 1000, 100, 40, 20, 10, and 1 respectively as shown in Table 4. That is, the simulation condition that uses 1 distinct data-generating model was conducted with 1000 replications while the simulation condition that uses 100 distinct data-generating models was conducted with 10 replications for each model. Using this scheme, the total number of replications for each set of simulation conditions was always 1000 but consisted of a mixture of data-generating models. Note that the condition with 1000 replications and 1 distinct data-generating model serves as the control condition, as this corresponds to traditional simulation design. In contrast, the condition with 1 replication and 1000 data-generating models serves as a fully random condition. In Table 4, the ratio of data-generating models to replications increases from left to right, with the most extreme condition tested being 1000 data-generating models with 1 replication each.

The total number of simulation conditions tested was \(n \times J \times \sigma ^2_{\varepsilon }= 5 \times 6 \times 3 = 90\) and each set of conditions (consisting of 1000 replications each) was run 500 times (we will refer to these as “test runs” to distinguish them from replications), yielding \(500\times 1000 = 500,000\) replications in total for each set of simulation conditions.

4.2 Evaluation criteria

Bias

As discussed previously, calculations of bias in Mosaic Monte Carlo under the aggregation summarizing strategy require modified formulas which were provided in Equations 2 and 3. As a reminder, the formula for relative bias is given by

where \(k\) corresponds to a given replication, \(j\) denotes a given data-generating model, and \(\theta \) is the population parameter corresponding to the estimator \(\widehat {\theta }\). Because the setting of the simulation study will yield theoretically unbiased results, we evaluated the performance of Mosaic Monte Carlo by examining the average and maximum absolute relative bias across the 500 test runs. This differs from how Mosaic Monte Carlo would be implemented in practice but allows us to better investigate deviations in performance.

Variance

To examine the impact of Mosaic Monte Carlo via the aggregation summarizing strategy on estimates of variance, we compared the empirical variances of the regression coefficients \(\boldsymbol {\widehat {\beta }}\) to their true variances, where the true variance-covariance matrix is given by

where \(X\) is the \(n\times p+1\) design matrix (consisting of the dataset and an additional column of 1’s) and \(\sigma ^2\) is the true error variance (which is fixed in the data generation process). The empirical variance was computed by calculating the variance of the regression coefficient estimates across all replications and then averaging across the data-generating models, as shown previously in Equation 4.

Monte Carlo standard errors

Monte Carlo standard errors (MCSEs) were calculated for the bias and variance estimates of regression coefficients to assess loss in precision. For bias, the MCSE was calculated for each data-generating model separately by calculating the standard deviation of the bias estimates and dividing by the square root of the number of replications according to the guidelines in Koehler et al. (2009). Bootstrap estimates of MCSEs were obtained for the variance of regression coefficients by drawing 200 bootstrap samples of the regression coefficient estimates for a given data-generating model, calculating their variances, and then taking the standard deviations of the 200 bootstrap values (Koehler et al., 2009).

5 Screening phase

As previously discussed, aggregation is only a viable summarizing strategy when it is done across mosaic pieces which impact results nonsystematically. Thus, we first sought to verify that our candidate mosaic pieces met this criteria by conducting the screening phase of Mosaic Monte Carlo before running our main simulation study. Our main simulation study uses a design with \(J\) and \(K\) as full simulation factors so we will omit parts of step 3 in the screening phase that refine these values for a later simulation phase.

In the screening phase, our candidate mosaic pieces were the values of regression coefficients and predictor correlations. We used an initial \(J = 50\) to ensure adequate sample size for conducting multiple regression analyses and an initial K of 200. We adaptively increased \(K\) for data-generating models whose MCSEs exceeded 20% of the corresponding true variance values. Thus, the MCSE screening thresholds were customized for each data-generating model to reflect true differences in underlying variability. Because regression coefficient estimates are unbiased in this setting, we used the SD instead of CV for calculations involving relative bias as a performance measure. To reduce computational burden, we conducted the screening phase with a fixed \(n=50\) and error variance of 100.

After adaptively increasing K to 500, all data-generating models obtained MCSEs within the desired screening thresholds. The SD for the bias of estimated regression coefficients are all less than the 0.1 threshold for CV values according to the first column of results in Table 5. To examine the impact of the mosaic pieces on variance estimates, we used the ratio of empirical variance to true variance as a performance measure. This quantity is more comparable across data-generating models with differing true variances. However, the CVs for the ratio of empirical to true variances of regression coefficients are just over the 0.1 threshold in Table 5. This suggests that population values of coefficients and correlations do not impact bias of regression coefficient estimates but perhaps could impact their variance estimates. Therefore, we conducted follow-up screening to determine whether these are systematic vs. nonsystematic effects. We used multiple regression with population coefficient and correlation values as predictors and the bias or variance estimates of each regression coefficient. The values of \(R^2_{adj}\) in Table 5 were all less than 0.1, suggesting that the population values of regression coefficients and correlations have nonsystematic impacts on estimates of bias and variance. This implies that aggregating the results of the simulation study is defensible across regression coefficient and correlation mosaic pieces.

6 Results

The results of our simulation study are presented in this section which is organized as follows. We begin with results that show the relative bias of regression coefficients with their MCSEs. Then, we present results regarding the variance of the regression coefficients and their MCSEs. Recall that the regression coefficients \(\widehat {\beta }_1\), \(\widehat {\beta }_2\), and \(\widehat {\beta }_3\) are exchangeable. In most of our tables and figures, the results for \(\widehat {\beta }_2\) and \(\widehat {\beta }_3\) are omitted for brevity. These results can be found in our supplemental materials.

6.1 Bias

The aggregated absolute relative bias was calculated according to Equation 3 for each regression coefficient across the 500 test runs from different random seeds. The results for \(\widehat {\beta }_0\) and \(\widehat {\beta }_1\) are shown in Table 6. The results for \(\widehat {\beta }_2\) and \(\widehat {\beta }_3\) are comparable to \(\widehat {\beta }_1\) and available in our supplemental material.

When there is only one randomized data-generating model with 1000 replications (i.e., under the control condition), the average absolute relative bias for the regression coefficients is minimal. This is true regardless of sample size, with the largest bias reaching 1.5% in a few conditions where the residual variance was 100. For the conditions with larger numbers of data-generating models and fewer replications, the average relative bias also tends to be small. However, in extreme cases when the sample size is 50 or 75 and the residual variance is 100, there are a few conditions where the absolute relative bias can be as large as 2 or 3%. As the number of data-generating models increases and the number of replications decreases, the absolute relative bias also tends to increase in small amounts. For example, when the sample size is 50 and the error variance is 25, the absolute relative bias of \(\widehat {\beta _0}\) for 10 data-generating models with 100 replications each is 0.0096 while the absolute relative bias for 1000 data-generating models with 1 replication each has increased to 0.0109. Compared to the traditional Monte Carlo condition (with an absolute relative bias of 0.0067), the absolute relative bias of \(\widehat {\beta }_0\) increases by approximately 0.003 by including 10 data-generating models and increases by approximately 0.004 by including 1000 data-generating models. The amount of this increase is also influenced by sampling variability (i.e., residual variance) and sample size, with smaller sample sizes and larger sampling variability leading to larger increases in absolute relative bias as the number of data-generating models increase and the number of replications decrease.

In a worst-case-scenario analysis, we also calculated the maximum absolute relative bias in each condition across all 500 test runs. The results in Table 7 show an interesting pattern – in most conditions, traditional Monte Carlo (1 data-generating model with 1000 replications) has the largest maximum biases compared to Mosaic Monte Carlo designs with larger numbers of data-generating models. In conditions where traditional Monte Carlo does not have the largest maximum bias, the largest values tend to occur with smaller numbers of data-generating models (such as with \(n=50\) and \(\sigma _{\varepsilon }^2 = 25\) when the largest maximum bias occurs with 20 data-generating models for \(\widehat {\beta }_1\)). The maximum biases tend to become smaller with greater numbers of data-generating models, although it is not a uniform pattern. This is a result of sampling variability. Occasionally, randomly generated datasets corresponding to a particular data-generating model may over or under represent the extremes of a distribution and lead to inaccurate parameter estimates. Mosaic Monte Carlo designs dilute the impact of these extreme occurrences, whereas traditional Monte Carlo designs offer no protection against this scenario.

To illustrate these differences in sampling variability across the competing Monte Carlo designs, the distribution of the aggregated bias of \(\widehat {\beta }_1\) across designs is shown in the boxplots in Figure 3 when \(n=50\) or \(100\). Each boxplot is constructed from 500 points which correspond to the aggregated bias of \(\widehat {\beta }_1\) calculated from a single test run.

These boxplots demonstrate that across sample sizes and conditions of residual variance, Mosaic Monte Carlo designs lead to smaller ranges of values for the relative bias of \(\widehat {\beta }_1\) across the 500 test runs. Larger sample sizes and smaller values of residual variance lead to greater precision. Traditional Monte Carlo designs with 1 data-generating model and 1000 replications have the smallest interquartile range compared to the Mosaic Monte Carlo designs, but also contain the most and largest outliers and therefore the largest overall ranges. As a reminder, implementing either approach in practice would be comparable to a single test run corresponding to a single point in the boxplots of Figure 3.

Figure 4 contains boxplots of the distribution of relative bias for \(\widehat {\beta }_1\) at larger sample sizes. The trends in these boxplots are similar to the ones in Figure 3: Mosaic Monte Carlo has smaller ranges of relative bias than traditional Monte Carlo designs but has larger interquartile ranges. Greater precision can be obtained with larger sample sizes and smaller values of residual variance. It is worth noting that the number of outliers in the boxplots are a stable feature that depends on the number of data-generating models and are not affected by conditions of sample size or residual variance. Thus, a researcher could not build in protection against possible outliers by adjusting these features. Figures 3 and 4 together illustrate that the probability of inaccuracy is higher for Mosaic Monte Carlo designs but that the magnitude of these inaccuracies can be much less than under traditional Monte Carlo.

Since MCSEs are specific to data-generating models, we present the MCSEs for the relative bias of the regression coefficients in two ways. First, we present the aggregated average MCSEs across all 500 test runs in Table 8. This quantifies the Monte Carlo error that is generally present across the data-generating models in a given condition. Then, we present the maximum of the maximum MCSEs across data-generating models and all 500 test runs in Table 9. That is, the maximum MCSE for a specific data-generating model is recorded and then the maximum of those values across the 500 test runs are shown in the table. This quantifies the maximum Monte Carlo error that is present in a given condition. As a note, the fully random condition with 1000 data-generating models and 1 replication each are not shown because the MCSEs cannot be calculated in conditions with only 1 replication.

Table 8 shows that as the number of data-generating models increase (and the number of replications decrease), there is an increase in the average MCSE values. This is expected since the calculating of MCSEs involves dividing by the square root of the number of replications. In fact, there is about a ten-fold increase in the MCSE from conditions with a single data-generating model to conditions with 100 data-generating models. This means that the differences are mostly in the denominator and that standard deviations of relative bias of the regression coefficients are generally stable across different Mosaic Monte Carlo designs and traditional Monte Carlo. Also as expected, MCSEs tend decrease with larger sample sizes and increase with larger values of residual variance.

Table 9 shows that the maximum MCSEs across Mosaic Monte Carlo designs can be quite different from each other. There is larger than a ten-fold increase between conditions with 1 data-generating model and conditions with 100 data-generating models. This happens because nothing in the table has been aggregated and so extreme occurrences will not be diluted by the Mosaic Monte Carlo procedure. For example, there are 50,000 MCSEs in conditions with 100 data-generating models but only 500 MCSEs in conditions with 1 data-generating model. The extremes of distributions will be better captured when there are 50,000 observations compared to 500. Like Table 8, Table 9 also shows that the maximum MCSEs of the relative bias of regression coefficients decreases with increasing sample size and increases with increasing residual variance.

6.2 Variance

Table 10 contains the aggregated ratios of empirical variance estimates to their true variances for \(\widehat {\beta }_0\) and \(\widehat {\beta }_1\). Values greater than 1 indicate that the empirical variance is larger than the true variance whereas values less than 1 indicate that the empirical variance is smaller than the true variance. Note that the condition with 1000 data-generating models with 1 replication each is omitted from the table; this is because the empirical variance cannot be calculated with only 1 replication.

There are a few trends worth noting in this table. First, the differences between the average empirical and true variances are generally to the third decimal place or less. Most of the exceptions occur regarding \(\widehat {\beta }_1\) when the number of data-generating models is 100 and the sample size is 100 or less, when the differences can be up to the second decimal place. The biggest discrepancy occurs for 100 data-generating models with a sample size of 100 and \(\sigma _{\varepsilon }^2 = 100\), when the ratio of empirical to true variance is 0.9798. Second, all conditions varying the number of data-generating models (i.e., both traditional and Mosaic Monte Carlo designs) tend to underestimate the variance of \(\widehat {\beta }_0\) while overestimating the variance of \(\widehat {\beta }_1\). Another interesting finding from this table is that there does not appear to be a systematic relationship between the number of data-generating models and the ratios of empirical to true variances. For example, the ratio of empirical to true variance of \(\widehat {\beta }_1\) increases from 0.9978 under traditional Monte Carlo (i.e., 1 data-generating model with 1000 replications) to 1.0062 when the number of data-generating models is increased to 25 for \(n = 75\) and \(\sigma _{\varepsilon }^2 = 25\). However, the ratio of empirical to true variance then decreases to 0.9963 when the number of data-generating models is increased to 50 under the same set of conditions. The values in Table 10 show that Mosaic Monte Carlo does not lead to systematic or substantial increases or decreases in variance estimates compared to traditional Monte Carlo designs.

Table 11 shows the MCSEs of the variance estimates of \(\widehat {\beta }_0\) and \(\widehat {\beta }_1\). Unlike the MCSEs of their biases, there is very little change in MCSE values across conditions of sample size or residual variance and so results are shown only for the condition with \(n = 50\) and \(\sigma _{\varepsilon }^2=100\) for the sake of brevity. The aggregated average MCSE of variance estimates in Table 11 show variability across the number of data-generating models that is not explained by the number of replications; there is less than a ten-fold increase in MCSE values as the number of data-generating models is increased from 1 to 100. Interestingly, this suggests that there is actually more variation in variance estimates due to Monte Carlo error in traditional Monte Carlo designs with 1 data-generating model compared to Mosaic designs involving several data-generating models. This pattern is present for both regression coefficients and for both the aggregated average MCSEs and the maximum of the maximum MCSEs.

Since researchers would only run Mosaic Monte Carlo designs a single time in implementation, we also illustrate the differences between empirical and true variance in a single test run to show what results might look like in practice. We will show results for \(\widehat {\beta }_0\) and \(\widehat {\beta }_1\) in conditions with \(n = 100\) and \(\sigma _{\varepsilon }^2 = 100\), as the patterns are the same for the other conditions of sample size and residual variance.

Figure 5 shows the empirical variances of \(\hat {\beta }_0\) for the first test run. The top plot shows the condition for 10 randomized data-generating models with 100 replications each. As we might expect, the differences between the true and empirical variance is very small, since the empirical variance is calculated with 100 replications in the denominator. The differences are less than 1% and the directions of the arrows do not suggest a systematic tendency to either underestimate or overestimate the population variance. The middle plot shows larger differences between the true variance and the empirical variance of the 50 data-generating models, although there are several models whose empirical variance estimates are quite close to the true variance. This pattern is intuitive since the empirical variance is calculated using only 20 replications. We can see an even greater degree of differences between the true variance and empirical variance in the bottom plot, which is the condition for 100 data-generating models with 10 replications each. From Figure 5, it is clear that the empirical variance will have greater differences from the true variance as the number of data-generating models increase. However, if the results are analyzed all together and averaged to form the mosaic, we will see small differences between the average empirical variance and the true variance as in Table 10.

Figure 6 shows the difference between empirical variances and the true variance for the regression coefficient \(\widehat {\beta }_1\). Recall that there is no meaningful difference between the results for \(\widehat {\beta }_1\) and the other two regression coefficients. The true variance fluctuates between data-generating models more than what was seen for the true variance of \(\widehat {\beta }_0\). As a reminder, this fluctuation occurs because we held \(X\) fixed across replications but allowed it to vary across data-generating models. However, the general trends from this plot are similar to what was seen in Figure 5. That is, the differences between the empirical variances and the true variances fluctuate more widely as the number of data-generating models increase.

7 Empirical Example

To demonstrate the performance of Mosaic Monte Carlo in a “real world context”, we used our method to redo part of the simulation study from Gomer et al. (2019). This article developed and compared several effect size measures in the context of structural equation modeling. Our goal is to examine the performance of 6 of these effect size measures under the same conditions of the original simulation study but using Mosaic Monte Carlo with factor loadings and factor correlations as candidate mosaic pieces. We will replicate Tables 4 and 5 from the original article containing results from the hierarchical linear regression and distance of the effect size measures from \(F_0\) (i.e., the population value of the discrepancy function). Since MCSEs are designed for estimators while the original work studied the population values of effect size measures, we will not compute MCSEs in this example. Due to this and and heavy computational burden in the simulation design, we set the number of data-generating models \(J\) to be 100 and the number of replications \(K\) to be 50. These values are the maximum \(J\) and \(K\) that could be afforded that would give a large enough sample size in the regression analyses in step 3 of the screening phase. The code for this empirical example is provided in a github repository.

7.1 Data generation

Data were generated according to a 3 factor measurement model with 9, 15, or 30 manifest variables. The factor correlation values were randomly generated from \(Unif(0.1, 0.7)\) while the factor loadings were randomly generated from \(Unif(0.4, 0.9)\), keeping the first loadings on each factor fixed to 1. These boundaries were chosen in order to obtain interpretable factor loadings according to the criteria in Stevens (1992). The population values for the factor variances and error variances were set to 1. The data distributions were normal, t, elliptical, non-normal, uniform, and mixture distributions according to the same specifications used in Gomer et al. (2019). The sample size was varied from 75, 100, 120, 150, 200, 300, 500, 800, to 1000. In conditions where the model was misspecified, a single cross-loading was omitted whose value ranged from 0.1 to 1 in increments of 0.1. The model was estimated via standard maximum likelihood (ML) or using robust M-estimation. The formulas for the effect sizes we considered are shown in Table 12.

7.2 Evaluation Criteria

Hierarchical linear regression To evaluate the performance of the effect size measures, the original article conducted a hierarchical linear regression to determine the impact of sample size, model size, distribution, and model misspecification on the value of the candidate effect size measures. The reduced model had a single predictor of model misspecification while the full model included predictors for sample size, model size, and distribution. Effect sizes with good performance were considered to be those with large values of \(R^2_{adj}\) and when the change in \(R^2_{adj}\) between the full and reduced models was small; in other words, the values of the effect size measure were mostly determined by the magnitude of model misspecification and not by the extraneous factors of sample size, model size, and distribution.

Distance to \(F_0\) Gomer et al. (2019) also analyzed the performance of the candidate effect size measures by calculating the distance of each effect size measure from \(F_0\), the population value of the discrepancy function. This is a measure of model misspecification. Effect size measures that were closest to \(F_0\) were considered to have good performance.

7.3 Screening phase

Due to computational burden, we performed a limited version of the simulation study in the screening phase. We conducted the screening simulation for models with 9 manifest variables and only considered model misspecification sizes of 0, 0.5, and 1 with a sample size of \(n=200\). However, we retained all 8 data-generating distributions. When the size of model misspecification was 0, the effect size measures had target values close to 0 so the SD was used as criteria in place of the CV for this condition.

In step 2 of the screening phase, we obtained SD and CV values ranging from approximately 0 to 0.406 across conditions of model misspecification size, distribution, and effect size measure. This suggests that the candidate mosaic pieces (i.e., factor loadings and correlations) could impact the values of the effect size measures in some conditions. Thus, we conducted follow-up investigation in step 3 of the screening phase.

We performed multiple regression analyses with the values of effect size measures as outcome variables and the factor loadings and correlations as predictors. The values of \(R^2_{adj}\) from these regressions for each effect size measure are shown in Table 13. As can be seen from the table, the values of \(R^2_{adj}\) are low for all effect size measures. This suggests that factor loadings and correlations have nonsystematic impacts on the values of effect size measures and also that aggregating over these mosaic pieces is defensible.

7.4 Results

The results of the hierarchical regression analysis when implementing the simulation study using Mosaic Monte Carlo are shown in Table 14. Rows with the smallest change in \(R^2_{adj.}\) are bolded. Encouragingly, Mosaic Monte Carlo suggests that the three best effect size measures according to this criteria are \({\mathscr{E}} _1\), \({\mathscr{E}} _3\), and \({\mathscr{E}} _6\). This replicates the results of the original article. The specific values in the table also tend to have small differences from Table 4 in the original article. However, there is one exception: the \(R^2_{adj.}\) values are quite small for \({\mathscr{E}} _1\). When looking at our results further, it appears that there were computational difficulties with this effect size measure because the formula involves the square root of a difference. When the difference is negative, the effect size cannot be calculated. This has a higher chance of occuring under Mosaic Monte Carlo methodology because of random variation in the data-generating process. However, it is impossible to know the characteristics of real datasets in advance so Mosaic Monte Carlo was able to capture a real difficulty of using \({\mathscr{E}} _1\) in practical applications.

To match the results presented in the original article, we re-created Table 5 showing the distances of \({\mathscr{E}} _1\), \({\mathscr{E}} _3\), and \({\mathscr{E}} _6\) from \(F_0\). The results based on our Mosaic Monte Carlo simulation are shown in Table 15. \({\mathscr{E}} _3\) is closest to \(F_0\) in almost every simulation condition, regardless of whether it is estimated with maximum likelihood or M-estimation. This mirrors the findings in the original article.

7.5 Discussion

The goal of replicating the simulation study from Gomer et al. (2019) using Mosaic Monte Carlo was to demonstrate that Mosaic Monte Carlo can be a viable simulation strategy and that it can improve the generalizability of findings. Not only was Mosaic Monte Carlo able to replicate the conclusion of the article and recommend \({\mathscr{E}} _3\) as the best candidate effect size measure, but it was also able to detect problematic features of \({\mathscr{E}} _1\) that failed to be captured originally.

8 Discussion and conclusion

The goal of this paper was to test an approach to Monte Carlo simulation design that breaks up large numbers of replications across factors that nonsystematically affect outcome(s) of interest. In the screening phase, researchers determine which factors have nonsystematic effects and appropriate numbers of data-generating models and replications. The simulation phase involves carrying out the main study of interest. We call this two-stage approach Mosaic Monte Carlo. This method can be used to improve the generalizability of findings from simulation studies compared to traditional designs that hold such factors as fixed. We showed how Mosaic Monte Carlo can be implemented in practice using an example from the literature on missing data. To test the baseline performance of Mosaic Monte Carlo, we examined how estimates of bias (RQ1) and variance (RQ2) were impacted by this approach when using the aggregation summarizing strategy and when the number of data-generating models is large compared to the number of replications. We conducted our study in the context of regression and also investigated how sample size and sampling variability might play a role in estimates of bias and variance (RQ3).

With respect to bias (RQ1), we found that Mosaic Monte Carlo generally leads to small increases in absolute relative bias for the regression coefficients which was also affected by the sample size and residual variance – small sample sizes with larger residual variances lead to increased absolute relative bias. However, traditional Monte Carlo designs had more outliers in performance and had the largest maximum values of absolute relative bias.

Regarding concerns about variance (RQ2), we found that Mosaic Monte Carlo did not systematically impact estimates of variance. Variance was affected mostly by sample size and by residual variance, and large differences between the true and empirical variance estimates occurred even in the control condition. Thus, we believe our method works well according to this criterion.

Our simulation study also showed that traditional Monte Carlo has more outliers in performance than Mosaic Monte Carlo. This means an entire set of results from an individual simulation study could be close or far from the true values with no way to distinguish between the two. It is almost unheard of to consider that an entire set of simulation results could be inaccurate reflections of the underlying data-generating model and yet this occurred multiple times in our main simulation study. In other words, the findings from a traditional Monte Carlo study could be substantially wrong with no way to detect it. Mosaic Monte Carlo helps to insure against this scenario.

The results from our study were encouraging and suggest that it is feasible to conduct Monte Carlo simulation studies by using this Mosaic Monte Carlo approach. This may be particularly useful for topics which may be sensitive to model parameter values such as missing data analysis and studies involving Bayesian methods. The results from the empirical example confirm that Mosaic Monte Carlo may be useful in settings where methods are not thought to be sensitive to the data-generating process, since Mosaic Monte Carlo was able to detect an issue with the formula of one of the candidate effect size measures.

8.1 Limitations

We evaluated the performance of Mosaic Monte Carlo by examining how estimates of bias and variance were impacted, although there may be other criteria of interest in simulation studies. Furthermore, we examined the performance of our method in a case where sample statistics are unbiased. We also only conducted our simulation study in the context of regression with a small number of variables which may be idealistic. We did not evaluate Mosaic Monte Carlo in contexts involving missing data or Bayesian statistics, although these areas were a prime motivation for developing this procedure.

A limitation of Mosaic Monte Carlo itself is that the user must be careful not to vary population parameters that would make models incomparable. In the context of our study, this applied to the error variance. The particular population parameters that should be fixed depends on the model used, the topic of interest, and the summarizing strategy. Importantly, Mosaic Monte Carlo improves the generalizability of simulation results but only within the boundaries studied. There are still risks that findings based on simulation results will fail to generalize to real data with different characteristics.

Furthermore, the guidance we developed for implementing the screening phase of Mosaic Monte Carlo may not fit all applications. As illustrated by our empirical example, compromises may need to made to reduce computational burden and/or adjustments made with respect to the criteria itself depending on the nature of the study.

8.2 Recommendations

We recommend that researchers who use Monte Carlo simulations in their work consider implementing their studies using Mosaic Monte Carlo, especially when studying methods that are sensitive to parameter values such as missing data analysis and Bayesian methodology. There may be unexpected cases in which aspects of data-generation can have nonsystematic influences on results.

References

Appendix: Supplemental Materials

\(n = 50\), \(\sigma _{\varepsilon }^2 = 9\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 9\)
\(n = 50\), \(\sigma _{\varepsilon }^2 = 25\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 25\)
\(n = 50\), \(\sigma _{\varepsilon }^2 = 100\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 100\)

Note. Solid line drawn at 0, dashed and dotted lines denote relative bias of 1% and 2.5% respectively.

\(n = 200\), \(\sigma _{\varepsilon }^2 = 9\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 9\)
\(n = 200\), \(\sigma _{\varepsilon }^2 = 25\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 25\)
\(n = 200\), \(\sigma _{\varepsilon }^2 = 100\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 100\)

Note. Solid line drawn at 0, dashed and dotted lines denote relative bias of 1% and 2.5% respectively.

\(n = 50\), \(\sigma _{\varepsilon }^2 = 9\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 9\)
\(n = 50\), \(\sigma _{\varepsilon }^2 = 25\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 25\)
\(n = 50\), \(\sigma _{\varepsilon }^2 = 100\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 100\)

Note. Solid line drawn at 0, dashed and dotted lines denote relative bias of 1% and 2.5% respectively.

\(n = 200\), \(\sigma _{\varepsilon }^2 = 9\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 9\)
\(n = 200\), \(\sigma _{\varepsilon }^2 = 25\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 25\)
\(n = 200\), \(\sigma _{\varepsilon }^2 = 100\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 100\)

Note. Solid line drawn at 0, dashed and dotted lines denote relative bias of 1% and 2.5% respectively.

\(n = 50\), \(\sigma _{\varepsilon }^2 = 9\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 9\)
\(n = 50\), \(\sigma _{\varepsilon }^2 = 25\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 25\)
\(n = 50\), \(\sigma _{\varepsilon }^2 = 100\)	\(n = 100\), \(\sigma _{\varepsilon }^2 = 100\)

Note. Solid line drawn at 0, dashed and dotted lines denote relative bias of 1% and 2.5% respectively.

\(n = 200\), \(\sigma _{\varepsilon }^2 = 9\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 9\)
\(n = 200\), \(\sigma _{\varepsilon }^2 = 25\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 25\)
\(n = 200\), \(\sigma _{\varepsilon }^2 = 100\)	\(n = 400\), \(\sigma _{\varepsilon }^2 = 100\)

Note. Solid line drawn at 0, dashed and dotted lines denote relative bias of 1% and 2.5% respectively.