Monte Carlo Methods

Why Do We Do Statistics?

• To study some target quantities in the population
  • Based on a limited sample
• How do we know that a statistics/statistical method gets us to a reasonable answer?
  • Analytic reasoning
  • Simulation

What is MC (in Statistics)?

• Simulate the process of repeated random sampling
  • E.g., repeatedly drawing sample of IQ scores of size 10 from a population
• Approximate sampling distributions
  • Using pseudorandom samples
• Study properties of estimators
  • regression coefficients, fit index
  • compare multiple estimators or modeling approaches

What is it (cont'd)?

• Based on Carsey & Harden (2014):
  • Simulations as experiments
    • Whether there's a "treatment" effect (but not why)
  • Simulations help develop intuition
    • Shouldn't replace analytically and theoretical reasoning
  • Simulations help evaluate substantive theories and empirical results

Examples in the Literature

• Curran, West, & Finch (1996, Psych Methods) studied the performance of the ${\chi }^2$ test for nonnormal data in CFA
• Kim & Millsap (2014, MBR) studied the performance of the Bollen-Stine Bootstrapping method for evaluating SEM fit indices
• MacCallum, Widaman, Zhang, & Hong (1999, Psych Methods) studied sample size requirement for getting stable EFA results
• Maas & Hox (2005, Methodology) studied the sample size requirement for multilevel models

Generating Random Data in R

With MC, one simulates the process of generating the data with an assumed data generating model

• Model: including both functional form and distributional assumptions

rnorm(5, mean = 0, sd = 1)

## [1]  0.8863733  1.6361050 -1.3694538 -1.1621330  1.1365392

rnorm(5, mean = 0, sd = 1)  # number changed

## [1]  0.6607360 -0.7283291  0.5751322  1.4180376 -1.5442535

Setting the Seed

• Most programs use algorithms to generate numbers that look like random
  • pseudorandom
  • Completely determined by the seed
• For replicability, ALWAYS explicitly set the seed in the beginning

set.seed(1)
rnorm(5, mean = 0, sd = 1)

## [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

set.seed(1)
rnorm(5, mean = 0, sd = 1)  # same seed, same number

## [1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

Generating Data From Univariate Distributions

rnorm(n, mean, sd)      # Normal distribution (mean and SD)
runif(n, min, max)      # Uniform distribution (minimum and maximum)
rchisq(n, df)           # Chi-squared distribution (degrees of freedom)
rbinom(n, size, prob)   # Binomial distribution

MC Approximation of $N(0,1)$

library(tidyverse)
set.seed(123)
nsim <- 20  # 20 samples
sam <- rnorm(nsim)  # default is mean = 0 and sd = 1
ggplot(tibble(x = sam), aes(x = x)) + 
  geom_density(bw = "SJ") + 
  stat_function(fun = dnorm, col = "red")  # overlay normal curve in red

Exercise

Parameter vs Estimator

• Estimator/statistic: $T\left(X\right)$, or simply $T$
  • How good does it estimate the population parameter, $\theta$?
• Examples:
  • $T=\overline{X}$ estimates $\theta =\mu$
  • $T=\frac{{\sum }_i(X_i-\overline{X}{)}^2}{N-1}$ estimates $\theta ={\sigma }^2$

Properties of Estimators

• Bias
• Consistency
• Efficiency
• Robustness

What is a Good Estimator?

Sampling Distribution

• What is it?

When to use MC?

• When it's difficult to analytically derive the sampling distribution
  • E.g., indirect effect, fit-indexes; Cohen's $d$, SEs of estimators
• When required assumptions are violated
  • E.g., normality, large sample
  • Model is misspecified
  • Used to check robustness of the estimator

A Simulation Study is an Experiment

Design

Like experimental designs, conditions should be carefully chosen

• What to manipulate? Sample size? Effect size? Why?
  • Based on statistical theory and reasoning
  • E.g., Gauss-Markov theorem: regression coefficients are unbiased with violations of distributional assumptions
• What levels? Why?
  • Needs to be realistic for empirical research
  • Maybe based on previous systematic reviews,
  • Or a small review of your own

Design (cont'd)

Full Factorial designs are most commonly used

Other alternatives include fractional factorial, random levels, etc

• See Skrondal (2000) for why they should be used more often

Generate

• Starts with a statistical data generating model
  • E.g., $Y_i={\beta }_0+{\beta }_1{X}_i+e_i,e_i\stackrel{\text{i.i.d.}}{\sim }N(0,{\sigma }^2)$
    • Systematic (deterministic) component: $X_i$
    • Random (stochastic) component: $e_i$
    • Constants (parameters): ${\beta }_0$, ${\beta }_1$
  • $Y_i$ completely determined by $X_i,e_i,{\beta }_0,{\beta }_1$

Model-Based Simulation

Analyze

Analyze the simulated data using one or more analytic approaches

• Misspecification: study the impact when analytic model omits important aspects of data generating model
  • E.g., ignoring clustering
• Comparison of approaches
  • E.g., Maximum likelihood vs. multiple imputation for missing data handling

Summarize (Evaluation Criteria)

• $\overline{\hat{\theta }}$ = ${\sum }_{i=1}^R{\hat{\theta }}_i/R$
• $\widehat{SD}\left(\hat{\theta }\right)$ = $\sqrt{\frac{{\sum }_{i=1}^R({\theta }_i-\overline{\hat{\theta }}{)}^2}{R}}$

For evaluating estimators:

• Bias
  • Raw: $\overline{\hat{\theta }}-\theta$
  • Relative: Bias / $\theta$
  • Standardized: Bias / $\widehat{SD}\left(\hat{\theta }\right)$
• Relative efficiency (only for unbiased estimators)
  • $RE(\hat{\theta },\tilde{\theta })$ = $\frac{\widehat{Var}\left(\tilde{\theta }\right)}{\widehat{Var}\left(\hat{\theta }\right)}$

Evaluation Criteria (cont'd)

For uncertainty estimators

• SE bias
  • Raw: $\overline{SE\left(\hat{\theta }\right)}-\widehat{SD}\left(\hat{\theta }\right)$
  • Relative: SE bias / $\widehat{SD}\left(\hat{\theta }\right)$

Combining bias and efficiency

• Mean squared error (MSE): $\frac{{\sum }_{i=1}^R({\theta }_i-\theta {)}^2}{R}$
  • Also = ${Bias}^2+\widehat{Var}\left(\hat{\theta }\right)$
  • Root MSE (RMSE) = $\sqrt{MSE}$
• Mainly to compare 2+ estimators

Evaluation Criteria (cont'd)

For statistical inferences:

• Power/Empirical Type I error rates
  • % with $p<\alpha$ (usually $\alpha$ = .05)
• Coverage of $C$% CI (e.g., $C$ = 95%)
  • % where the sample CI contains $\theta$

Results

Just like you're analyzing real data

• Plots, figures
• ANOVA, regression
  • E.g., 3 (sample size) × 4 (parameter values) 2 (models) design: 2 between factors and 1 within factor

Number of Replications

Should be justified rather than relying on rule of thumbs

Why Does MC Work?

• Law of large number
  • ${\sum }_{i=1}^R{T}_i/R\stackrel{p}{\to }\theta$
• When $R$ is large,
  • the empirical distribution $\hat{F}\left(t\right)$ converges to the true sampling distribution $F\left(t\right)$.

Number of Replications (cont'd)

How Good is the Approximation

• Monte Carlo (MC) Error
  • Like standard error (SE) for a point estimate
• For expectations (e.g., bias)
  • MC Error = $\widehat{SD}\left(\hat{\theta }\right)/\sqrt{R}$

E.g., if one wants the MC error to be ≤2.5% of the sampling variability, R needs to be 1 / ${.025}^2$ = 1,600

Number of Replications (cont'd)

For power (also Type I error) and CI coverage,

• MC Error = $\sqrt{\frac{p(1-p)}{R}}$

E.g., with R = 250, and empirical Type I error = 5%,

sqrt((.05 * (1 - .05)) / 250)

## [1] 0.01378405

So R should be increase for more precise estimates

Reporting MC Results

(Boomsma, 2013, Table 1, p. 521)

See Boomsma (2013), Table 2, p. 526 for a checklist

Efficiency tips

• Things that don't change should be outside of a loop
• Initialize place holders when using for-loops
• Vectorize
• Strip out unnecessary computations
• Parallel computing (using the future package)

Other topics not covered

• Error handling
• Assessing convergence
• Debugging
• Interfacing with other software (e.g., Mplus, LISREL, HLM)

Further Readings

Carsey and Harden (2014) for a gentle introduction

Chalmers (2019) and Sigal and Chalmers (2016) for using the R package SimDesign

Harwell, Kohli, and Peralta-Torres (2018) for a review of design and reporting practices

Skrondal (2000), Serlin (2000), and Bandalos and Leite (2013) for additional topics

References

References (cont'd)

References (cont'd)