Monte Carlo studies are a common tool in statistics and related fields. They are used for everything from the evaluation of the finite sample properties of new statistical methods to the generation of probability distributions for risk management.
The MonteCarlo package for the R language provides tools to create simulation studies quickly and easily and it also allows to summarize the results in LaTeX tables. This vignette gives details on the implementation of the MonteCarlo package and presents examples for its application.
To run a simulation study, the user has to nest both - the generation of a sample and the calculation of the desired statistics from this sample - in a single function. This function is passed to MonteCarlo(). No additional programming is required.
The idea behind this approach is to allow the user full control and flexibility with regard to the design of the Monte Carlo experiment. It also makes MonteCarlo() very versatile. Finally, it is very intuitive. The user formulates his experiment as if he/she was only interested in a single draw.
A First Example
The best way to get working with the MonteCarlo package is to look at an example. Here we evaluate the performance of a standard t-test of the hypothesis \(H_0: \mu=0\) vs \(H_1: \mu\neq0\). We are interested to see how the size and power of the test change with the sample size (n), the distance from the null hypothesis (loc - for location) and the standard deviation of the distribution (scale).
The test statistic is given by
\(t=\frac{\bar x}{\hat \sigma}\),
where \(\bar x\) and \(\hat \sigma\) are the arithmetic mean and the estimated standard deviation. The sample is generated from a normal distribution.
To conduct this analysis, we proceed as follows. First, we load the MonteCarlo package.
Then we define the following function.
#########################################
## Example: t-test
# Define function that generates data and applies the method of interest
ttest<-function(n,loc,scale){
# generate sample:
sample<-rnorm(n, loc, scale)
# calculate test statistic:
stat<-sqrt(n)*mean(sample)/sd(sample)
# get test decision:
decision<-abs(stat)>1.96
# return result:
return(list("decision"=decision))
}
As discussed above, ttest() is formulated in a way as if we only want to generate a single test decision. The arguments of the function are the parameters we are interested in.
Our ttest() function carries out 4 steps:
- Draw a sample of n observations from a normal distribution with mean \(\mu=loc\) and standard deviation \(\sigma=scale\).
- Calculate the t-statistic.
- Determine the test decision.
- Return the desired result in form of a list.
We then define the combinations of parameters that we are interested in and collect them in a list. The elements of this list must have the same names as the parameters for which we want to supply grids.
# define parameter grid:
n_grid<-c(50,100,250,500)
loc_grid<-seq(0,1,0.2)
scale_grid<-c(1,2)
# collect parameter grids in list:
param_list=list("n"=n_grid, "loc"=loc_grid, "scale"=scale_grid)
To run the simulation, the function ttest() and the parameter grid (param_list) are passed to MonteCarlo(), together with the desired number of Monte Carlo repetitions (here nrep=1000).
# run simulation:
MC_result<-MonteCarlo(func=ttest, nrep=1000, param_list=param_list)
There is no further coding required. All the mechanics of the Monte Carlo experiment are handled by the MonteCarlo() function.
Calling summary produces a short information on the simulation.
## Simulation of function:
##
## function(n,loc,scale){
##
## # generate sample:
## sample<-rnorm(n, loc, scale)
##
## # calculate test statistic:
## stat<-sqrt(n)*mean(sample)/sd(sample)
##
## # get test decision:
## decision<-abs(stat)>1.96
##
## # return result:
## return(list("decision"=decision))
## }
## <bytecode: 0x000000000bda0e68>
##
## Required time: 20.51 secs for nrep = 1000 repetitions on 1 CPUs
##
## Parameter grid:
##
## n : 50 100 250 500
## loc : 0 0.2 0.4 0.6 0.8 1
## scale : 1 2
##
##
## 1 output arrays of dimensions: 4 6 2 1000
As one can see from the summary, the simulation results are stored in an array of dimension c(4,6,2,1000). The Monte Carlo repetitions are collected in the last dimension of the array.
To summarize the results in a reasonable way and to include them as a table in a paper or report, we have to represent them in a matrix. This is handled by the MakeTable() function, that stacks the submatrices collected in the array in the rows and columns of a new and larger matrix and prints the result in the form of LaTeX code to generate the desired table.
To determine in which order the results are stacked in rows and columns, we supply the function arguments rows and cols that are vectors of the names of the parameters in the order in which we want them to appear in the table, sorted from the inside to the outside.
# generate table:
MakeTable(output=MC_result, rows="n", cols=c("loc","scale"), digits=2, include_meta=FALSE)
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrrrrrrrrrrrr }
## \hline\hline\\\\
## scale && \multicolumn{ 6 }{c}{ 1 } & & \multicolumn{ 6 }{c}{ 2 } \\
## n/loc & & 0 & 0.2 & 0.4 & 0.6 & 0.8 & 1 & & 0 & 0.2 & 0.4 & 0.6 & 0.8 & 1 \\
## & & & & & & & & & & & & & & \\
## 50 & & 0.06 & 0.28 & 0.81 & 0.99 & 1.00 & 1.00 & & 0.05 & 0.12 & 0.28 & 0.56 & 0.80 & 0.94 \\
## 100 & & 0.04 & 0.49 & 0.97 & 1.00 & 1.00 & 1.00 & & 0.06 & 0.19 & 0.53 & 0.83 & 0.98 & 1.00 \\
## 250 & & 0.05 & 0.88 & 1.00 & 1.00 & 1.00 & 1.00 & & 0.05 & 0.36 & 0.89 & 1.00 & 1.00 & 1.00 \\
## 500 & & 0.05 & 0.99 & 1.00 & 1.00 & 1.00 & 1.00 & & 0.05 & 0.59 & 1.00 & 1.00 & 1.00 & 1.00 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ decision }
## \end{table}
To change the ordering, just change the vectors rows and cols. As can be seen below, the results change accordingly.
# generate table:
MakeTable(output=MC_result, rows=c("n","scale"), cols="loc", digits=2, include_meta=FALSE)
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrrrrrr }
## \hline\hline\\\\
## scale & n/loc & & 0 & 0.2 & 0.4 & 0.6 & 0.8 & 1 \\
## & & & & & & & & \\
## \multirow{ 4 }{*}{ 1 } & 50 & & 0.06 & 0.28 & 0.81 & 0.99 & 1.00 & 1.00 \\
## & 100 & & 0.04 & 0.49 & 0.97 & 1.00 & 1.00 & 1.00 \\
## & 250 & & 0.05 & 0.88 & 1.00 & 1.00 & 1.00 & 1.00 \\
## & 500 & & 0.05 & 0.99 & 1.00 & 1.00 & 1.00 & 1.00 \\
## & & & & & & & & \\
## \multirow{ 4 }{*}{ 2 } & 50 & & 0.05 & 0.12 & 0.28 & 0.56 & 0.80 & 0.94 \\
## & 100 & & 0.06 & 0.19 & 0.53 & 0.83 & 0.98 & 1.00 \\
## & 250 & & 0.05 & 0.36 & 0.89 & 1.00 & 1.00 & 1.00 \\
## & 500 & & 0.05 & 0.59 & 1.00 & 1.00 & 1.00 & 1.00 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ decision }
## \end{table}
Now we can simply copy the LaTeX code and add it to our paper, report or presentation. That is all. The user can focus on the design of the experiment that he is interested in and on the interpretation of the results.
We will now discuss the details of MonteCarlo() and MakeTable() and then give further examples.
MakeFrame()
It is often usefull to visualize the results of simulation studies with help of density plots, box plots and so on. To facilitate this, the MonteCarlo package does not provide its own suite of plotting functions. Instead, the function MakeFrame() converts the simulation results contained in the output of the MonteCarlo() function into a data.frame that is suitable to be manipulated and plotted with powerful standard tools such as dplyr, ggplot2 and of course base graphics. The structure of the dataframe is such that each row contains the parameters used and the results of the simulation for one repetition. This is illustrated below.
ttest<-function(n,loc,scale){
sample<-rnorm(n, loc, scale)
stat<-sqrt(n)*mean(sample)/sd(sample)
return(list("stat"=stat))
}
MC_result<-MonteCarlo(func=ttest, nrep=1000, param_list=param_list)
df<-MakeFrame(MC_result)
head(df)
## n loc scale stat
## 1 50 0.0 1 0.06306346
## 2 100 0.0 1 -0.14677467
## 3 250 0.0 1 1.24365601
## 4 500 0.0 1 -0.69706901
## 5 50 0.2 1 2.08853104
## 6 100 0.2 1 0.67448822
If func returns a list with more than one element, then each of the list elements has its own column in the dataframe
One can now use the dataframe to conduct further analyses. For example, we can plot estimated densities for the distribution of the t-test with different sample sizes n.
library(dplyr)
library(ggplot2)
tbl <- tbl_df(df)
ggplot(filter(tbl, loc==0, scale==1)) + geom_density(aes(x=stat, col=factor(n)))
Further Examples: Comparison of Two Estimators
In the previous sections we focused on the t-test as an example for the use of MonteCarlo() and MakeTable(). Here, we consider a second example: the evaluation of estimators. In particular, we compare the mean() and the median() as estimators for the expected value of a Gaussian random variable. This example serves to illustrate the use of the optinal arguments of the MakeTable() function.
Multiple Outputs
First, we specify the function of interest: mean_vs_median() (defined below) generates a sample of n normal distributed random variables with standard deviation scale. It then calulates the mean() and the median() and returns both in a list. Therefore, in contrast to the t-test example used above, this function returns a list with two elements.
mean_vs_median<-function(n,scale){
# generate sample
sample<-rnorm(n, 0, scale)
# calculate estimators
mean_sample<-mean(sample)
median_sample<-median(sample)
# return results
return(list("mean"=mean_sample, "median"=median_sample))
}
n_grid<-c(50, 250, 500)
scale_grid<-c(1, 2, 4)
param_list=list("n"=n_grid, "scale"=scale_grid)
# run simulation:
erg_mean_median<-MonteCarlo(func=mean_vs_median, nrep=1000, param_list=param_list)
Like before, we pass the function and the parameter grid to MonteCarlo() and then feed the results to MakeTable().
MakeTable(output=erg_mean_median, rows="n", cols="scale", digits=2, include_meta=FALSE)
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrr }
## \hline\hline\\\\
## n/scale & & 1 & 2 & 4 \\
## & & & & \\
## 50 & & -0.01 & 0.01 & -0.01 \\
## 250 & & 0.00 & 0.00 & 0.00 \\
## 500 & & 0.00 & 0.00 & 0.00 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ mean }
## \end{table}
##
##
##
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrr }
## \hline\hline\\\\
## n/scale & & 1 & 2 & 4 \\
## & & & & \\
## 50 & & -0.01 & 0.01 & 0.00 \\
## 250 & & 0.00 & 0.00 & 0.00 \\
## 500 & & 0.00 & 0.00 & 0.00 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ median }
## \end{table}
As one can see, the default behavior of MakeTable() is to print separate tables for each value in the list returned by func.
Merge Several Outputs Into One Table
To include both, the estimates for the mean and median in the same table, we simply include the string “list” in either rows or cols. Here, we choose to include it as the last element of cols, to have the results for mean() and median() next to each other.
MakeTable(output=erg_mean_median, rows=c("n"), cols=c("scale","list"), digits=2, include_meta=FALSE)
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrrrrrr }
## \hline\hline\\\\
## list && \multicolumn{ 3 }{c}{ mean } & & \multicolumn{ 3 }{c}{ median } \\
## n/scale & & 1 & 2 & 4 & & 1 & 2 & 4 \\
## & & & & & & & & \\
## 50 & & -0.01 & 0.01 & -0.01 & & -0.01 & 0.01 & 0.00 \\
## 250 & & 0.00 & 0.00 & 0.00 & & 0.00 & 0.00 & 0.00 \\
## 500 & & 0.00 & 0.00 & 0.00 & & 0.00 & 0.00 & 0.00 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ mean, median }
## \end{table}
Generation of Tables for Partial Parameter Grids
For large parameter grids, one might want to focus on a subset of the results. To achive this, just include the optional argument partial_grid. Here, the grids for n and scale both had three values and we only want to use the largest and the smallest values in the grids. We therefore supply a list with the names of the parameters and vectors specifying the positions of the values that we want to keep in the grid.
# use partial_grid
MakeTable(output=erg_mean_median, rows="n", cols=c("scale","list"), digits=2,
partial_grid=list("n"=c(1,3), "scale"=c(1,3)), include_meta=FALSE)
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrrrr }
## \hline\hline\\\\
## list && \multicolumn{ 2 }{c}{ mean } & & \multicolumn{ 2 }{c}{ median } \\
## n/scale & & 1 & 4 & & 1 & 4 \\
## & & & & & & \\
## 50 & & -0.01 & -0.01 & & -0.01 & 0.00 \\
## 500 & & 0.00 & 0.00 & & 0.00 & 0.00 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ mean, median }
## \end{table}
Collapse
Since both estimators - the mean() and the median() are unbiased in this setup, comparing them by their bias might not be the most interesting analysis. Instead, we might be interested in other aspects of their distributions - for example their standard deviations. To modify the function used by MakeTable() to summarize the result from the nrep repetitions of the experiment, a list can be passed to collapse that contains the names of the functions that are supposed to be used to collapse the realisations of the estimators to a meaningful table. Here we use sd(), since we know that mean() is the more efficient estimator in this setup.
MakeTable(output=erg_mean_median, rows="n", cols=c("scale", "list"), digits=2, collapse=list("sd", "sd"), include_meta=FALSE)
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrrrrrr }
## \hline\hline\\\\
## list && \multicolumn{ 3 }{c}{ mean } & & \multicolumn{ 3 }{c}{ median } \\
## n/scale & & 1 & 2 & 4 & & 1 & 2 & 4 \\
## & & & & & & & & \\
## 50 & & 0.14 & 0.29 & 0.56 & & 0.17 & 0.35 & 0.69 \\
## 250 & & 0.06 & 0.13 & 0.26 & & 0.08 & 0.16 & 0.33 \\
## 500 & & 0.04 & 0.09 & 0.18 & & 0.05 & 0.11 & 0.23 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ mean, median }
## \end{table}
Now, it is easy to see, that even though the bias terms are very similar, the mean() has a smaller variance than the median().
Non-Nested Data Generating Processes
One of the key concepts behind the design of MonteCarlo() is that it is generally assumed that the samples are generated from a data generated process that nests all the situations of interest. This is also the case in the examples discussed so far, where the data is generated from a normal distribution and only the standard deviation of the distribution and the number of observations change.
However, this does not mean that MonteCarlo() cannot handle non-nested data generating processes. As an example, we consider the case where we want to compare the standard deviation of the mean() and the median() as we did above, but for samples from a normal distribution and from a uniform distribution. These are clearly non-nested. To deal with this situation, we include DGP as an additional argument in mean_vs_median() and create a case distinction, so that the value of DGP determines from which distribution the sample is generated. We then include the possible values of DGP in an additional parameter grid.
mean_vs_median<-function(n,scale,DGP){
# generate sample
if(DGP=="normal"){sample<-rnorm(n, 0, scale)}
if(DGP=="uniform"){
b<-scale/(2*sqrt(1/12))
sample<-runif(n, -b, b)
}
# calculate estimators
mean_sample<-mean(sample)
median_sample<-median(sample)
# return results
return(list("mean"=mean_sample, "median"=median_sample))
}
n_grid<-c(50, 250, 500)
scale_grid<-c(1, 2, 4)
DGP_grid<-c("normal", "uniform")
param_list=list("n"=n_grid, "scale"=scale_grid, "DGP"=DGP_grid)
Note that the width of the uniform distribution is specified so that both the Gaussian and the uniform distribution have a standard deviation of scale.
# run simulation:
erg_mean_median<-MonteCarlo(func=mean_vs_median, nrep=1000, param_list=param_list)
After running the simulation, we obtain the following results.
MakeTable(output=erg_mean_median, rows=c("n","DGP"), cols=c("scale", "list"), digits=2, collapse=list("sd", "sd"), include_meta=FALSE)
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrrrrrrr }
## \hline\hline\\\\
## & list && \multicolumn{ 3 }{c}{ mean } & & \multicolumn{ 3 }{c}{ median } \\
## DGP & n/scale & & 1 & 2 & 4 & & 1 & 2 & 4 \\
## & & & & & & & & & \\
## \multirow{ 3 }{*}{ normal } & 50 & & 0.14 & 0.29 & 0.56 & & 0.17 & 0.36 & 0.69 \\
## & 250 & & 0.06 & 0.13 & 0.26 & & 0.08 & 0.16 & 0.33 \\
## & 500 & & 0.04 & 0.09 & 0.18 & & 0.06 & 0.11 & 0.23 \\
## & & & & & & & & & \\
## \multirow{ 3 }{*}{ uniform } & 50 & & 0.15 & 0.29 & 0.57 & & 0.24 & 0.48 & 0.98 \\
## & 250 & & 0.06 & 0.12 & 0.24 & & 0.11 & 0.21 & 0.42 \\
## & 500 & & 0.05 & 0.09 & 0.19 & & 0.08 & 0.15 & 0.31 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ mean, median }
## \end{table}
We can see that the standard deviation of the median() is far larger for the uniform distribution whereas that of the mean() remains nearly unchanged.
This example shows that MonteCarlo() is very versatile. Depending on the specification of func, it is possible to deal with non-nested DGPs. Similarly, multivariate outputs could be handled by splitting them in several elements and returning them seperately.