Parallel and other simulations in R made easy: An end-to-end study

It is shown how to set up, conduct, and analyze large simulation studies with the new R package simsalapar = simulations simplified and launched parallel. A simulation study typically starts with determining a collection of input variables and their values on which the study depends, such as sample sizes, dimensions, types and degrees of dependence, estimation methods, etc. Computations are desired for all com- binations of these variables. If conducting these computations sequentially is too time- consuming, parallel computing can be applied over all combinations of select variables. The final result object of a simulation study is typically an array. From this array, sum- mary statistics can be derived and presented in terms of (flat contingency or LATEX) tables or visualized in terms of (matrix-like) figures. The R package simsalapar provides several tools to achieve the above tasks. Warnings and errors are dealt with correctly, various seeding methods are available, and run time is measured. Furthermore, tools for analyzing the results via tables or graphics are pro- vided. In contrast to rather minimal examples typically found in R packages or vignettes, an end-to-end, not-so-minimal simulation problem from the realm of quantitative risk management is given. The concepts presented and solutions provided by simsalapar may be of interest to students, researchers, and practitioners as a how-to for conducting real- istic, large-scale simulation studies in R. Also, the development of the package revealed useful improvements to R itself, which are available in R 3.0.0.


Introduction
Realistic mathematical or statistical models are often complex and not analytically tractable, thus require to be evaluated by simulation. In many areas such as finance, insurance, or statistics, it is therefore necessary to set up, conduct, and analyze simulation studies. Apart from minimal examples which address particular tasks, one often faces more difficult setups with a complex simulation problem at hand. For example, if a comparably small simulation already reveals an interesting result, it is often desired to conduct a larger study, involving more parameters, a larger sample size, or more simulation replications. However, run time for In the language of mathematics, this can be made precise as follows. Let S t,j denote the value of the jth of d stocks at time t ≥ 0. The value of a portfolio with these d stocks at time t is thus where β 1 , . . . , β d denote weights, typically the number of shares of stock j in the portfolio. Considering the logarithmic stock prices as risk factors, the risk-factor changes are given by X t+1,j = log(S t+1,j ) − log(S t,j ) = log(S t+1,j /S t,j ), j ∈ {1, . . . , d}. ( Assume that all quantities at time point t (interpreted as today) are known, and we are interested in the time point t + 1 (one period ahead, for example one year). The loss of the portfolio at t + 1 can therefore be expressed as = − d j=1 w t,j (exp(X t+1,j ) − 1) that is, in terms of the known weights w t,j (at time t, β j and S t,j , j ∈ {1, . . . , d}, are known), and the unknown risk-factor changes. Value-at-Risk (VaR α ) of L t+1 at level α ∈ (0, 1) is given by where F − L t+1 (y) = inf{x ∈ R : F L t+1 (x) ≥ y} denotes the quantile function of the distribution function F L t+1 of L t+1 (equal to the ordinary inverse F −1 L t+1 if F L t+1 is continuous and strictly increasing; see Embrechts and Hofert (2013) for more details about such functions).
For simplicity, we drop the time index t + 1 in what follows. Let X = (X 1 , . . . , X d ) be the d-dimensional vector of (possibly) dependent risk-factor changes. By Sklar (1959), its distribution function H can be expressed as for a copula C and the marginal distribution functions F 1 , . . . , F d of H. A copula is a distribution function with standard uniform univariate margins; for an introduction to copulas,

Translating the scientific problem to R
To summarize, our goal is to simulate, for each sample size n, dimension d, copula family C, and strength of dependence Kendall's tau τ , N sim times n losses L ki , k ∈ {1, . . . , N sim }, i ∈ {1, . . . , n}, and to compute in the kth of the N sim replications VaR α (L) as the empirical α-quantile of L ki , i ∈ {1, . . . , n}, for each α. Since different α-quantiles can (and should!) be estimated based on the same simulated losses, we do not have to generate additional samples for different values of α, VaR α (L) can be estimated simultaneously for all α under consideration. Table 1 provides a summary of all variables involved in our simulation study, their names in R, L A T E X expressions, type, and the corresponding values we choose. Note that this table is produced entirely with simsalapar's toLatex(varList, ....); see page 6. For the moment, The variable N sim gives the number of simulation ("bootstrap") replications in our study. This variable is present in many statistical simulations and allows one to provide an error measure of a statistical quantity such as an estimator. Because of this special meaning, it gets the type "N", and there can be only one variable of this type in a simulation study. If it is not given, it will implicitly be treated as 1.
frozen: The variable w is a list of length equal to the number of dimensions considered, where each entry is a vector (in our case a value which will be sufficiently often recycled by R) of length equal to the corresponding dimension. Variables such as w (or the marginal quantile functions) remain the same throughout the whole simulation study, but one might want to change them if the study is conducted again. Variables of this type are assigned the type "frozen", since they remain fixed throughout the whole study.
grid: Variables of type "grid" are used to build a (physical) grid. In R this grid is implemented as a data frame. Each row in this data frame contains a unique combination of variables of type "grid". The number of rows n G of this grid, is thus the product of the lengths of all variables of type "grid". The simulation will iterate N sim times over all n G rows and conduct the required computations. Conceptually, this corresponds to visiting each of the N sim × n G rows of a virtual grid (seen as N sim copies of the grid pasted together). The computations for one row in this virtual grid are viewed as one sub-job. In many situations, computing all sub-jobs sequentially turns out to be time-consuming (even after profiling of the code and removing time bombs such as deeply nested 'for' loops). In this situation, we can apply parallel computing and distribute the sub-jobs over several cores of a multi-core processor or several machines (nodes) in a cluster.
inner: Finally, variables of type "inner" are all dealt with within a sub-job for reasons of convenience, speed, load balancing etc. As mentioned before, in our example, α plays such a role since VaR α (L) can be estimated simultaneously for all α under consideration based on the same simulated losses.
As result of a simulation, we naturally obtain an array. This array has one dimension for each variable of type "grid" or "inner", and one additional dimension if N sim > 1. Besides the variable names, their type, and their values, we also define R expressions for each variable. These expressions are later used to label tables or plots when the simulation results are analyzed.

Remark 2.1
As an advantage of our approach based on n.sim in terms of load-balancing, each repeated simulation has the same expected run time. Note, however, that thousands of fast sub-jobs might lead to a comparably large overall run time due to both the waiting times for the jobs to start on a cluster and due to the overhead in communication between the master and the slaves. It might therefore be more efficient to send blocks of sub-jobs (say, 10 sub-jobs) to the same core or node. This feature is provided by the argument block.size in the do*() functions (doLapply(), doForeach(), doRmpi(), doMclapply(), doClusterApply()) presented later.
We are now ready to start writing an R script which can be run on a single computer or on a computer cluster. Since cluster types and interfaces are quite different, we only focus on how to write the R script here 2 . The first task is to implement the variable list presented above.
Note that one actually does not need to specify a type for n.sim or variables of type "frozen", the default chosen is "frozen" unless the variable is n.sim in which case it is "N".
The function getEl() can be used to extract elements of a certain type from a variable list (defaults to all values).
1 > str(getEl(varList, "grid")) # extract "value" of variables of type "grid" List of 4 $ n : num [1:2] 64 256 $ d : num [1:4] 5 20 100 500 $ family: chr [1:2] "Clayton" "Gumbel" $ tau : num [1:2] 0.25 0.5 in, one can submit the script simu.R via bsub -N -W 01:00 -n 48 -R "select[model==Opteron8380]" -R " span[ptile=16]" mpirun -n 1 R CMD BATCH simu.R, for example, where the meaning of the various options is as follows: -N sends an email to the user when the batch job has finished; -W 01:00 submits the job to the one-hour queue (jobs with this maximal wall-clock run time) on the cluster; the option -n 48 asks for 48 cores (one is used as master, 47 as slaves); -R "select[model==Opteron8380]" specifies X86_64 nodes with AMD Opteron 8380 CPUs for the sub-jobs to be run (this is important if run-time comparisons are required, since one has to make sure that the same architecture is used when computations are carried out in parallel); the option -R "span[ptile=16]" specifies that (all) 16 cores (on each node) are used on a single node (that means our job fully occupies 48/16 = 3 nodes); mpirun specifies an Open MPI job which runs only one copy (-n 1) of the program; and finally, R CMD BATCH simu.R is the standard call of the R script simu.R in batch mode. To have a look at the grid for our working example (containing all combinations of variables of type "grid"), the function mkGrid() can be used as follows.

The result of a simulation
Our route from here is to conduct the simulations required for each line of the virtual grid (in parallel). As an important point, note that each computational result naturally consists of the following components:

value:
The actual value. This is can be a scalar, numeric vector, or numeric array whose dimensions depend on variables of type "inner". The computed entries also depend on variables of type "frozen", but they do not enter the result array as additional dimensions.

error:
It is important to adequately track errors during simulation studies. If one computation fails, we lose all results computed so far and thus have to do the work again (fix the error, move the files to the cluster, wait for the simulation job to start, wait for it to fail or to finish successfully in this next trial run etc.). To avoid this, we capture the errors to be able to deal with them after the simulation has been conducted. This also allows us to compute statistics about errors, such as percentages of runs producing errors etc.
warning: Similar to errors, warnings are important to catch. They may indicate nonconvergence of an algorithm (or a maximal number of iterations reached etc.) and therefore impact reliability of the results.

time:
Measured run time can also be an indicator of reliability in the sense that if computations are too fast/slow, there might be a programming error (not leading to an error or warning and thus being detected). For example, if one accidentally switches a logical condition, a large computation may return in almost no time because it simply ended up in the wrong case. If the value computed from this case is not suspicious, and if there were no warnings and errors, then run time is the only indicator of a possible bug in the code. Furthermore, measuring run time is also helpful for benchmarking and assessing the usefulness of a result (even if a computation or algorithm only runs sufficiently fast on a large cluster, it might not be suitable for a notebook and therefore might have limited use overall).
.Random.seed: The random seed right before the user-specified computations are carried out. This is useful for reproducing single results for debugging purposes.
In many simulation studies, also on an academic level, focus is put on value only. We therefore particularly stress all of these components, since they become more and more important for obtaining reliable results the larger the conducted simulation study is. Furthermore, error, warning, and .Random.seed are important to consider especially during experimental stage of the simulation, for checking an implementation, and testing it for numerical stability.
The paradigm of simsalapar is that the user only has to take care of how to compute the value (the statistic the user is most interested in). All other components addressed above are automatically dealt with by simsalapar. We will come back to the latter in Section 2.5, after having thought about how to compute the value for our working example in the following section.

Writing the problem-specific function doOne()
Programming in R is about writing functions. Our goal is now to write the workhorse of the simulation study: doOne(). This function has to be designed for the particular simulation problem at hand and is therefore given here (with Roxygen documentation) instead of being part of simsalapar. doOne() computes the value (a numeric vector here) for the given arguments, that is, the component value. For functions doOne() for other simulations, we refer to the demos of simsalapar, see for example demo(TGforecasts) for reproducing the simulation conducted by Gneiting (2011 To conduct the main simulation, we only need one more function which iterates over all sub-jobs and calls doOne(). There are several options: sequential (see Section 2.6) versus various approaches for parallel computing (see Section 3), for which we provide the do*() functions explained below. Since these functions are quite technical and lengthy, we will present the details in Section 5. For the moment, our goal is to understand the functions they call in order to understand how the simulation works. Figure 1 visualizes the main functions involved in conducting the simulation. These functions break down the whole task into smaller pieces (which improves readability of the code and simplifies debugging when procedures fail).
We have already discussed the innermost, user-provided function doOne( Parallel and other simulations in R made easy: An end-to-end study part of this function deals with correctly setting the seed. It also provides a monitor feature; see Section 5.1 for the details.
As mentioned before, there are several choices available for the outermost layer of functions, depending on whether, and if yes, what kind of parallel computing should be used to deal with the rows of the virtual grid. In particular, simsalapar provides the following functions, see Section 5: doLapply(): a wrapper for the non-parallel function lapply(). This is useful for testing the code with a small number of different parameters so that the simulation still runs locally on the computer at hand.

doForeach():
a wrapper for the function foreach() of the R package foreach to conduct computations in parallel on several cores or nodes. A version specific to our working example based on nested foreach() loops is presented in Section 5.

doRmpi():
a wrapper for the function mpi.apply() or its load-balancing version mpi.applyLB() (default) from the R package Rmpi for parallel computing on several cores or nodes.
doMclapply(): a wrapper for the function mclapply() (with (default) or without loadbalancing) of the R package parallel for parallel computing on several cores (not working on Windows).
doClusterApply(): a wrapper for the function clusterApply() or its load-balancing version clusterApplyLB() (default) of the R package parallel for parallel computing on several cores or nodes.

Remark 2.2
The user of simsalapar can call one of the above functions do*() to finally run the whole simulation study; see Sections 2.6 and 3. To this end, these functions iterate over all sub-jobs and finally call the function saveSim(); see Section 5.1. saveSim() tries to convert the resulting list of lists of length four or five to an array of lists of length four or five and saves it in the .rds file specified by the argument sfile. If this non-trivial conversion fails 4 , the raw list of lists of length four or five is saved instead, so that results are not lost. This behavior can also be obtained by directly specifying doAL=FALSE when calling the do*() functions. To further avoid that the conversion fails, the functions do*() conduct a basic check of the correctness of the return value of doOne() by calling the function doCheck(). This can also be called by the user after implementing doOne() to verify the correctness of doOne(); see, for example, demo(VaRsuperadd).

Running the simulation sequentially: doLapply() based on lapply()
In Sections 3 and 5, we will compare different approaches for parallel computing in R. To make this easier to follow, we start with doLapply(), see Section 5.1, which is a wrapper for the sequential (non-parallel) function lapply() to iterate over all rows of the virtual grid. This sequential approach is often the first choice to try (for a smaller number of parameter combinations) in order to check whether the simulation actually does what it should, for debugging etc. If sequential computations based on lapply() turn out to be too slow, one can easily use one of the parallel computing approaches described in Sections 3 and 5, since they share the same interface.
We now demonstrate the use of doLapply() to run the whole simulation. Note that names is an optional argument to our doOne() and the argument monitor, passed to subjob(), allows progress monitoring.

Parallel computing in R
In the same way that doLapply() wraps around lapply(), simsalapar provides convenient wrapper functions to conduct the same computations (but) in parallel. These different approaches are useful for different kinds of setups, such as different available computer architectures or different specifications of the simulation study considered. Before we go into the details, let us mention that one should only use one of the do*() functions. Mixing several different ways of conducting parallel computations in the same R process might lead to weird errors, conflicts of various kinds, or unreliable results at best.
For conducting computations in parallel with R, one just needs to replace doLapply() above (Section 2.6) by one of its "parallelized" do*() versions listed in Section 2.5. We will take doClusterApply() as an example here and refer to Section 5 for a more in-depth analysis and comparison of the results obtained from these different approaches to those from doLapply() to check their correctness, consistency, and efficiency. 1 > res5 ← doClusterApply(varList, sfile="res5_clApply_seq.rds", 2 doOne=doOne, names=TRUE) Indeed, doClusterApply() produces the same result as doLapply() did above: 1 > stopifnot(doRes.equal(res5, res)) # note: doRes.equal() is part of simsalapar

Data Analysis
After having conducted the main simulation, the final task is to analyze the data and present the results. It seems difficult to provide a general solution for this part of the simulation study. Besides the solutions provided by simsalapar however, it might therefore be required to write additional problem-specific functions. In this case, functions from simsalapar may at least serve as good starting points.
The function getArray(), presented in Section 5.2, is a function from simsalapar which, given the result object of the simulation and one of the components "value" (the default), "error", "warning", or "time" creates an array containing the corresponding results. This is typically more convenient than working with an array of lists, which the object as returned by one of the do*() functions naturally is. For the components being "error" or "warning", the array created contains (by default) boolean variables indicating whether there was an error or warning, respectively. This behavior can be changed by providing a suitable argument FUN to getArray(). Additionally, getArray() allows for an argument err.value, defaulting to NA, for replacing values in case there was an error. As mentioned before, each "value", can be a scalar, a numeric vector, or a numeric array, often with dimnames, e.g., resulting from (the outer product of) variables of type "inner". Note that for conducting the simulation, variables sometimes can be declared as "inner" or "frozen" interchangeably. However, this changes the dimension of the result object for the analysis in the sense that variables of type "inner" appear as additional dimensions in the result array and can thus serve as a proper quantity/dimension in a table or plot, whereas variables of type "frozen" do not.
As a first part of the analysis, we are interested in how reliable our results are. We thus consider possible errors and warnings of the computations conducted. Flat contingency tables (obtained by ftable()) allow us to conveniently get an overview as follows. Parallel and other simulations in R made easy: An end-to-end study Since we neither have warnings nor errors in our numerically non-critical example study, let us briefly consider the run times: In what follows, we exclusively focus on the actual computed values, hence the array val. We apply tools from simsalapar that allow us to create flexible L A T E X tables and sophisticated graphs for representing these results.

Creating L A T E X tables
In this section, we create L A T E X tables of the results. Our goal is to make this process modular and flexible. We thus leave tasks such as formatting of table entries as much as possible to the user. Note that there are already R packages available for generating L A T E X tables, for example the well-known xtable or the rather new tables. However, they do not fulfill the above requirements (and come with other unwanted side effects concerning the table headers or formatting of entries we do not want to cope with). We therefore present new tools for constructing tables with simsalapar. For inclusion in L A T E X documents, only the L A T E X package tabularx, and, due to our defaults following the paradigm of booktabs, the L A T E X package booktabs have to be loaded in the .tex document. Much more sophisticated alignment of column entries for L A T E X tables than we show here (even including units) can be achieved in combination with the L A T E X package siunitx; see its corresponding extensive manual. Note that these packages all come with standard L A T E X distributions.
After having computed arrays of (robust) Value-at-Risk estimates and (robust) standard deviations via 1 > non.sim.margins ← setdiff(names(dimnames(val)), "n.sim") 2 > huber. ← function ( we format and merge the arrays. As just mentioned, we specifically leave this task to the user to guarantee flexibility. As an example, we put the (robust) standard deviations in parentheses and colorize 5 all entries corresponding to the largest level α. Next, we create a flat contingency table from the array of formatted results fres. The arguments row.vars and col.vars of ftable() specify the basic layout of Table 2 below.

> ## format values and mads
1 > ft ← ftable(fres, row.vars=c("family","n","d"), col.vars=c("tau","alpha")) To summarize, using functions from simsalapar and packages from L A T E X, one can create flexible L A T E X tables. If the simulation results become sufficiently complicated, creating L A T E X tables (or at least parts of them) from R reduces a lot of work, especially if the simulation study has to be repeated due to bug fixes, improvements, or changes in the implementation. Note that the table header typically constitutes the main complication when constructing tables. It might still require manual modifications in case our carefully chosen defaults do not suffice. simsalapar provides many other functions not presented here, including the (currently nonexported) functions ftable2latex() and fftable() and the (exported) functions tablines()   and wrapLaTable(). These ingredient functions of the method toLatex.ftable can still be useful if one encounters very specific requirements not covered by toLatex.ftable. More details on the latter can be found in Section 5.2. A crucial step in the development of tablines() was the correct formatting of an ftable without introducing empty rows or columns. For this we introduced four different methods of "compactness" of a formatted ftable which are available in format.ftable() from R version 3.0.0 and for earlier versions in simsalapar.

Graphical analysis
Next we show how simsalapar can be applied to visualize the results of our study. In modern statistics, displaying results with graphics (as opposed to tables) is typically good practice, since it is easier to see the story the data would like to tell us. For example, in a table, the human eye can only compare two numbers at a time, in well-designed graphics much more information is visible.
There are various different approaches of how to create graphics in R, for example, with the traditional graphics package, the lattice, or the ggplot2 package. The most flexible approach is based on grid graphics; see Murrell (2006). In what follows, we apply the function mayplot() (based on grid and graphics via gridBase) from simsalapar for creating a plot matrix (also known as conditioning plot) from an array of values. Within each cell of this plot a traditional graphic is drawn to visualize the results.
In our example study, the strength of dependence in terms of Kendall's tau determines the columns of the matrix-like plot and the copula family determines its rows. In each cell, there is an x and a y axis. For making comparisons easier, one typically would like to have the same limits on the y axes across different rows of the plot matrix. Sometimes it makes sense to have separate scales for y axes in different rows (while still having the same scales for all plots within the same row). This behavior can be determined with the argument ylim (being "global" (the default) or "local") of mayplot(). For our working example, the x axis provides the different significance levels α. We thus naturally can depict three different input variables in such a layout (copula families, Kendall's taus, and significance levels α). The y axis may show point estimates or boxplots of the simulated Value-at-Risk values as given in val.
All other variables (sample sizes n, dimensions d) then have to be depicted in the same cell, visually distinguished by different line types or colors, for example (currently one such variable is allowed; we chose d below by fixing n = 256). If more variables are involved, one might even want to put more variables in one cell, rethink the design, or split different values of a variable over separate plots. N sim , if available, enters the scene through a second label on the right side of graphic.
With mayplot() it is easy to create a graphical result (a pdf file for inclusion in a L A T E X document, for example) 6 . Figures 2 and 3 display the results for n = 256. The former shows boxplots of all the N sim simulated Value-at-Risk estimates VaR α (L), whereas the latter depicts corresponding robust Huber "means" and also demonstrates mayplot() for N sim = 1 or, equivalently, no N sim at all. Overall, we see that a graphic such as Figure 2 is easier to grasp and to infer conclusions from than > mayplot(v256, varList, row.vars="family", col.vars="tau", xvar="alpha",   The function doCallWE() The R package simsalapar provides the following auxiliary function doCallWE() for computing the components value, error, warning, and time as addressed in Section 2.3. It is called from subjob() and based on tryCatch.W.E() which is part of R's demo(error.catching) for catching both warnings and errors. A numeric vector, say s, of length n.sim, providing seeds for each of the n.sim simulation replications, i.e., simulation i receives seed set.seed(s[i]), for i from 1 to n.sim. For a fixed replication i, the seed is the same no matter what row in the (physical) grid is considered. This ensures least variance across the computations for the same replication i.
In particular, it also leads to the same results no matter which variables are of type "grid" or "inner"; see demo(robust.mean) where this is tested. This is important to guarantee since one might want to change certain "inner" variables to "grid" variables due to load-balancing while computing the desired statistics based on the same seed (or generated data from this seed). Clearly, since replication i is guaranteed to get seed s[i] (no matter when the corresponding sub-job is computed relative to all other sub-jobs), this seeding method provides reproducible results.
A list of length n.sim which provides seeds for each of the n.sim simulation replications. In contrast to the case of a numeric vector, this case is meant to be for providing more general seeds. At the moment, seeds for l'Ecuyer's random number generator L'Ecuyer-CMRG can be provided; see l'Ecuyer, Simard, Chen, and Kelton (2002) for a reference and Section 5.3 for how to use it. This seeding method also provides reproducible results.
NA: In this case .Random.seed remains untouched. In contrast to NULL, it is not even generated if it does not exist. Also, the fifth component .Random.seed is not concatenated to the result in this case. In all other cases, it is appended if keepSeed=TRUE. As mentioned before, the default keepSeed=FALSE has been chosen to avoid large result objects. Clearly, seeding method NA typically does not provide reproducible results.
a character string, specifying a certain seeding method. Currently, only "seq" is provided, a convenient special case of the second case addressed above, where the vector of seeds is simply 1:n.sim, and thus provides reproducible results.
If keepSeed=TRUE and seed is not NA, subjob() saves .Random.seed as the fifth component of the output vector (besides the four components returned by doCallWE()). This is useful for reproducing the result of the corresponding call of doOne() for debugging purposes, for example.
The default seeding method in the do*() functions is "seq". This is a comparably simple default which guarantees reproducibility. Note, however, that for very large simulations, there is no guarantee that the random-number streams are sufficiently "apart". For this, we recommend l'Ecuyer's random number generator L'Ecuyer-CMRG; see Section 5.3 for an example.

The function doLapply()
As mentioned before, doLapply() is essentially a wrapper for lapply() to iterate (sequentially) over all rows in the virtual grid, that is, over all sub-jobs. As an important ingredient, saveSim(), explained below, is used to deal with the raw result list.  The functions saveSim() and maybeRead() After having conducted the main simulation with one of the do*() functions, we would like to create and store the result array. It can then be loaded and worked on for the analysis of the study which is often done on a different computer. For creating, checking, and saving the array, simsalapar provides the function saveSim().
If possible, saveSim() creates an array of lists (via mkAL()), where each element of the array is a list of length four or five as returned by subjob(). If this fails, saveSim() simply takes its input list. It then stores this array (or list) in the given .rds file (via saveRDS()) and returns it for further usage. In our working example, the array itself is five-dimensional, the dimensions corresponding to n, d, C, τ , and N sim . c("value", "error", "warning", "time") %in% names (x1) For reading a saved object of a simulation study, simsalapar provides the function maybeRead().
If the provided .rds file exists, maybeRead() reads and returns the object. Otherwise, maybeRead() does nothing (hence the name). This is useful for reading and analyzing the result object at a later stage by executing the same R script containing both the simulation and its analysis 7 .

Select functions for the analysis
The function getArray() As promised in Section 4, we now present the implementation of the function getArray(). This function receives the result array of lists, picks out a specific component of the lists, and returns an array containing these components. This is especially useful when analyzing the results of a simulation.

The method toLatex.ftable and related functions
The ftable method toLatex.ftable for creating L A T E X tables calls several auxiliary functions, 7 Note that the first part of this paper is itself such an example. typically, this is the header of the flat contingency table created by fftable(). head contains a "collapsed" version of head.raw but in a much more sophisticated way. \multicolumn statements for centering of column headings and title rules for separating groups of columns are introduced (\cmidrule if booktabs=TRUE; otherwise \cline). The list component align is a string which contains the alignment of the table entries (as accepted by L A T E X's tabular environment). The default implies that all columns containing row names are left-aligned and all other columns are right-aligned. The component rsepcol is a vector of characters which contain the row separators rsep or, additionally, \addlinespace commands for separating blocks of rows belonging to the same row variables or groups of such. The default chooses a larger space between groups of variables which appear in a smaller column number. In other words, the "largest" group is determined by the variables which appear in the first column, the second-largest by those in the second column etc. up to the second-last column containing row variables. For more details we refer to the source code of tablines() in simsalapar.
Finally, the method toLatex.ftable calls wrapLaTable(). This function wraps a L A T E X table and tabular environment around, which can be put in a L A T E X document.

Function mayplot() to visualize a 5D array
We will now present a bit more details about the function mayplot() for creating matrix-like plots of arrays up to dimension five. Due to space limitations, we only describe mayplot() verbally here and refer to the source code of simsalapar for the exact implementation.

Marius Hofert, Martin Mächler 29
mayplot() utilizes the function grid.layout() to determine the matrix-like layout, including spaces for labels; call mayplot() with show.layout=TRUE to see how the layout looks like. pushViewport() is then used to put the focus on a particular cell of the plot matrix (or several cells simultaneously, see the global y axis label, for example). The focus is released via popViewport(). Within a particular cell of the plot matrix a panel function is chosen for plotting. This is achieved by gridBase. The default panel function is either boxplot.matrix() or lines() depending on whether n.sim exists. We also display a background with grid lines similar to the style of ggplot2. Axes (for the y axis in logarithmic scale using eaxis from sfsmisc) are then printed depending on which cell the focus is on; similar for the row and column labels of the cells, again in ggplot2-style. Due to the flexibility of grid, we can also create a legend in the same way as in the plot. Finally, we save initial graphical parameters with opar <-par(no.readonly=TRUE) and restore them on function exit in order to not change graphical parameters for possible subsequent plots.
Overall, mayplot() is quite flexible in visualizing results contained in arrays of dimensions up to five, see the corresponding help file for more customizations.

Using foreach
The wrapper doForeach() is based on the function foreach() of the package foreach. It allows to carry out parallel computations on multiple nodes or cores. In principle, different parallel backends can be used to conduct parallel computations with foreach(). For example, SNOW cluster types could be specified with registerDoSNOW() from the package doSNOW.

Using foreach with nested loops
The approach we present next is similar to doForeach(). However, it uses nested foreach() loops to iterate over the grid variables and replications; see the vignettes of foreach for the technical details. Since this is context specific, doNestForeach() is not part of simsalapar. Unfortunately, it is not possible to execute statements between different foreach() calls. This would be interesting for efficiently computing those quantities only once which remain fixed in subsequent foreach() loops. Note that this is also not possible for the other methods for parallel computing and thus not a limitation of this method alone.

Using Rmpi
The following wrapper function doRmpi() utilizes only tools from the R package Rmpi for parallel computing on multiple nodes or cores in R via MPI. With load.balancing=TRUE (the default), the load-balancing version mpi.applyLB() is utilized (otherwise mpi.apply()) which sends the next sub-job to a slave who just finished one. if(!is.null(r ← maybeRead(sfile))) return ( Similar as before, we now call doRmpi() for our working example, with seed=NULL, and n.sim=1, respectively. We also show here, that seed=NULL is typically non-reproducible.   To see that doMclapply() and doLapply() yield the same result, let us check for equality of res4 with res. We also check equality of res42 with res02 which shows the same for l'Ecuyer's random number generator. if(!is.null(r ← maybeRead(sfile))) return (

Conclusion
The R package simsalapar allows one to easily set up, conduct, and analyze large-scale simulations studies. The user of our package only has to provide the list of input variables on which the simulation study depends (which can be created with the function varlist()) and the function which computes the desired statistic (or result of the study) for one combination of input variables (termed doOne() here). The user can then choose between different functions to conduct the simulation (sequentially via doLapply() or in parallel via one of doForeach(), doRmpi(), doMclapply(), or doClusterApply()), possibly involving replicates (via a variable of type "N" as our n.sim here). Important aspects of a simulation study such as catching of errors and warnings, measuring run time, or dealing with seeds are automatically taken care of and adjusted easily. Furthermore, simsalapar provides various tools to analyze the results. Besides several useful auxiliary functions, the high-level functions toLatex() and mayplot() can be used to create sophisticated L A T E X tables and matrix-like figures of the results, respectively.
In the first part of the paper (up to and including Section 4), we explained and guided the user/reader through a working example end-to-end, which highlights various of the above steps. More advanced information about simsalapar, including explanations of functions under the hood, tests, and further examples were either addressed in the second part of the paper (Section 5) or can be found in the package itself; see, for example, the demos of simsalapar.