This note instructs how to reproduce the two-way ANOVA random effects model (Tables 1 and 2) simulation example in the manuscript. This simulation compares the performance among MM, EM, FS (Fisher scoring) and lme4 under model: $$ y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ijk}, \quad 1 \le i \le a, 1 \le j \le b, 1 \le k \le c, $$ where $\alpha_i \sim \text{Normal}(0, \sigma_1^2)$, $\beta_j \sim \text{Normal}(0, \sigma_2^2)$, $(\alpha \beta)_{ij} \sim \text{Normal}(0, \sigma_3^2)$, and $\epsilon_{ijk} \sim \text{Normal}(0, \sigma_4^2)$ independently.
Comparison metrics include:
versioninfo()
To run this document, we need Julia packages
Pkg.add("JLD")
Pkg.add("RCall")
and R packages lme4
install.packages("lme4")
The script anova_sim.jl
will run the simulations with a=b=5
, c
varying in {5, 10, 20, 50}, vcratio
($\sigma_i / \sigma_4$, $i = 1, 2, 3$) varying in {0, 0.05, 0.1, 1, 10, 20}, and nsim=50
replicates. It saves the results into a file anova_sim_result.jld
:
tic()
include("anova_sim.jl")
toc()
To show summary statistics on the number of iterations (Table 1 in manuscript):
showwhat = 1
include("anova_sim_print.jl")
To show summary statistics on the run times (Table 2 in manuscript):
showwhat = 3
include("anova_sim_print.jl")
To show summary statistics on the objective values:
showwhat = 5
include("anova_sim_print.jl")
To show summary statistics on the MSE:
showwhat = 7
include("anova_sim_print.jl")