Two-way ANOVA Simulation Example

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:

  1. Iterations
  2. Runtime
  3. Objective value
  4. MSE $m^{-1} \sum_{j=1}^m (\widehat \sigma_j^2 - \sigma_j^2)^2$. Convergence criterion for MM, EM, and FS is set to be the relative change in objective value less than $10^{-8}$.
In [1]:
versioninfo()
Julia Version 0.6.3
Commit d55cadc350 (2018-05-28 20:20 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, skylake)

Prerequisites

To run this document, we need Julia packages

Pkg.add("JLD")
Pkg.add("RCall")

and R packages lme4

install.packages("lme4")

Run simulations

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:

In [2]:
tic()
include("anova_sim.jl")
toc()
WARNING: RCall.jl: Loading required package: Matrix
(c, vcratio) = (5, 0.0)
(c, vcratio) = (5, 0.05)
(c, vcratio) = (5, 0.1)
(c, vcratio) = (5, 1.0)
(c, vcratio) = (5, 10.0)
(c, vcratio) = (5, 20.0)
(c, vcratio) = (10, 0.0)
(c, vcratio) = (10, 0.05)
(c, vcratio) = (10, 0.1)
(c, vcratio) = (10, 1.0)
(c, vcratio) = (10, 10.0)
(c, vcratio) = (10, 20.0)
(c, vcratio) = (20, 0.0)
(c, vcratio) = (20, 0.05)
(c, vcratio) = (20, 0.1)
(c, vcratio) = (20, 1.0)
(c, vcratio) = (20, 10.0)
(c, vcratio) = (20, 20.0)
(c, vcratio) = (50, 0.0)
(c, vcratio) = (50, 0.05)
(c, vcratio) = (50, 0.1)
(c, vcratio) = (50, 1.0)
(c, vcratio) = (50, 10.0)
(c, vcratio) = (50, 20.0)
Julia Version 0.6.3
Commit d55cadc350 (2018-05-28 20:20 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, skylake)
elapsed time: 233.135299544 seconds
Out[2]:
233.135299544

Display simulation results

To show summary statistics on the number of iterations (Table 1 in manuscript):

In [3]:
showwhat = 1
include("anova_sim_print.jl")
\begin{tabular}{rrrrrr}
\toprule
$\sigma^2_{1}/\sigma^2_{e}$ & Method &
\multicolumn{4}{c}{$c=$ \# observations per combination} \\
\cline{3-6}
& & 5 & 10 & 20 & 50  \\
\midrule
0.00 & MM & 143.12(99.76) & 118.26(62.91) & 96.26(50.61) & 81.10(33.42) \\
& EM & 2297.72(797.95) & 1711.70(485.92) & 1170.06(365.48) & 788.10(216.60) \\
& FS & 25.64(11.20) & 21.10(7.00) & 16.46(4.37) & 13.88(2.88) \\
[0.25cm]
0.05 & MM & 121.86(98.52) & 69.38(50.23) & 55.88(37.34) & 29.50(18.80) \\
& EM & 1464.26(954.27) & 538.04(504.42) & 254.90(253.86) & 104.98(157.97) \\
& FS & 16.78(9.13) & 12.62(6.22) &  9.68(3.22) &  8.10(1.34) \\
[0.25cm]
0.10 & MM & 84.74(59.33) & 62.98(50.48) & 40.46(31.43) & 25.86(18.79) \\
& EM & 985.46(830.49) & 360.32(462.62) & 157.70(231.91) & 68.26(107.85) \\
& FS & 15.20(10.10) & 10.58(5.92) &  8.58(3.56) &  7.50(1.72) \\
[0.25cm]
1.00 & MM & 31.04(33.27) & 29.60(27.66) & 25.32(25.39) & 24.90(20.76) \\
& EM & 130.18(299.03) & 161.14(290.23) & 64.20(135.38) & 84.88(137.88) \\
& FS &  6.62(4.72) &  6.32(3.64) &  5.12(1.87) &  5.36(1.50) \\
[0.25cm]
10.00 & MM & 29.80(35.42) & 34.16(38.25) & 28.82(28.44) & 20.90(14.28) \\
& EM & 115.94(274.33) & 177.30(301.71) & 80.12(155.67) & 75.02(127.38) \\
& FS & 12.72(5.14) & 12.86(4.94) & 11.66(3.95) & 11.76(3.66) \\
[0.25cm]
20.00 & MM & 30.10(32.92) & 32.72(39.02) & 23.70(21.20) & 19.62(15.67) \\
& EM & 148.04(318.40) & 85.86(180.28) & 61.74(140.84) & 37.36(83.89) \\
& FS & 18.76(7.51) & 17.40(5.21) & 17.22(5.67) & 16.28(5.03) \\
\bottomrule
\end{tabular}

To show summary statistics on the run times (Table 2 in manuscript):

In [4]:
showwhat = 3
include("anova_sim_print.jl")
\begin{tabular}{rrrrrr}
\toprule
$\sigma^2_{1}/\sigma^2_{e}$ & Method &
\multicolumn{4}{c}{$c=$ \# observations per combination} \\
\cline{3-6}
& & 5 & 10 & 20 & 50  \\
\midrule
0.00 & MM & 11.46(7.77) & 10.06(5.29) & 11.93(6.35) & 10.44(3.99) \\
& EM & 189.32(71.32) & 148.20(48.13) & 147.87(49.97) & 96.28(24.97) \\
& FS & 34.27(33.47) & 24.89(8.55) & 23.70(14.15) & 20.46(4.54) \\
& lme4 & 25.84(12.10) & 22.32(1.25) & 27.34(4.06) & 36.14(5.59) \\
[0.25cm]
0.05 & MM &  9.79(7.72) &  6.19(4.22) &  6.87(4.37) &  4.45(2.20) \\
& EM & 116.03(75.57) & 47.72(45.35) & 30.60(29.88) & 14.23(19.68) \\
& FS & 19.18(10.23) & 15.37(7.48) & 12.78(4.06) & 12.39(2.35) \\
& lme4 & 22.76(1.96) & 24.88(2.60) & 28.72(3.10) & 47.34(16.29) \\
[0.25cm]
0.10 & MM &  7.07(4.78) &  6.29(4.94) &  5.14(3.72) &  3.95(2.23) \\
& EM & 78.96(66.19) & 35.48(45.81) & 19.53(27.71) &  9.67(13.56) \\
& FS & 17.36(11.26) & 14.44(9.00) & 12.08(6.31) & 11.47(2.40) \\
& lme4 & 22.66(1.83) & 28.90(8.70) & 30.16(4.43) & 44.58(4.89) \\
[0.25cm]
1.00 & MM &  2.66(2.61) &  3.22(2.91) &  3.57(3.15) &  3.85(2.50) \\
& EM & 10.71(23.93) & 15.88(27.52) &  8.35(16.26) & 11.34(16.65) \\
& FS &  7.88(5.44) &  9.10(4.95) &  7.12(2.42) &  8.46(2.27) \\
& lme4 & 23.12(1.75) & 30.22(9.37) & 29.96(4.47) & 42.82(8.32) \\
[0.25cm]
10.00 & MM &  2.48(2.72) &  3.24(3.19) &  3.84(3.35) &  3.35(1.71) \\
& EM &  9.66(22.02) & 15.98(26.57) & 10.24(18.78) & 10.27(15.40) \\
& FS & 15.19(6.05) & 16.39(6.11) & 15.81(5.15) & 18.14(5.46) \\
& lme4 & 35.02(3.83) & 47.12(8.10) & 63.24(15.33) & 102.78(34.49) \\
[0.25cm]
20.00 & MM &  2.57(2.49) &  3.13(3.53) &  3.13(2.44) &  3.07(1.81) \\
& EM & 12.28(25.71) &  8.44(16.89) &  8.01(17.12) &  5.47(9.76) \\
& FS & 22.09(8.53) & 22.03(6.14) & 23.08(7.21) & 23.99(7.38) \\
& lme4 & 37.34(12.91) & 50.24(8.59) & 63.62(17.39) & 91.14(28.39) \\
\bottomrule
\end{tabular}

To show summary statistics on the objective values:

In [5]:
showwhat = 5
include("anova_sim_print.jl")
\begin{tabular}{rrrrrr}
\toprule
$\sigma^2_{1}/\sigma^2_{e}$ & Method &
\multicolumn{4}{c}{$c=$ \# observations per combination} \\
\cline{3-6}
& & 5 & 10 & 20 & 50  \\
\midrule
0.00 & MM & -176.67(7.94) & -353.59(10.19) & -713.42(14.90) & -1776.40(25.02) \\
& EM & -176.68(7.94) & -353.60(10.18) & -713.43(14.90) & -1776.41(25.02) \\
& FS & -176.67(7.94) & -353.59(10.19) & -713.42(14.90) & -1776.40(25.02) \\
& lme4 & -176.67(7.94) & -353.59(10.19) & -713.42(14.90) & -1776.40(25.02) \\
[0.25cm]
0.05 & MM & -181.06(7.24) & -365.39(10.92) & -722.16(15.36) & -1794.26(25.96) \\
& EM & -181.06(7.24) & -365.39(10.92) & -722.16(15.36) & -1794.26(25.96) \\
& FS & -181.06(7.24) & -365.39(10.92) & -722.16(15.36) & -1794.26(25.96) \\
& lme4 & -181.06(7.24) & -365.39(10.92) & -722.16(15.36) & -1794.26(25.96) \\
[0.25cm]
0.10 & MM & -185.10(6.77) & -368.83(10.28) & -726.21(13.35) & -1813.88(20.24) \\
& EM & -185.10(6.77) & -368.83(10.28) & -726.21(13.35) & -1813.88(20.24) \\
& FS & -185.10(6.77) & -368.83(10.28) & -726.21(13.35) & -1813.88(20.24) \\
& lme4 & -185.10(6.77) & -368.83(10.28) & -726.21(13.35) & -1813.88(20.24) \\
[0.25cm]
1.00 & MM & -204.05(8.18) & -392.35(10.94) & -754.95(13.93) & -1831.11(22.63) \\
& EM & -204.05(8.18) & -392.35(10.94) & -754.95(13.93) & -1831.12(22.63) \\
& FS & -204.05(8.18) & -392.35(10.94) & -754.95(13.93) & -1831.11(22.63) \\
& lme4 & -204.05(8.18) & -392.35(10.94) & -754.95(13.93) & -1831.11(22.63) \\
[0.25cm]
10.00 & MM & -233.65(7.90) & -416.56(11.79) & -777.85(15.65) & -1862.05(28.29) \\
& EM & -233.65(7.90) & -416.56(11.79) & -777.85(15.65) & -1862.05(28.29) \\
& FS & -233.65(7.90) & -416.56(11.79) & -777.85(15.65) & -1862.05(28.29) \\
& lme4 & -233.65(7.90) & -416.56(11.79) & -777.85(15.65) & -1862.05(28.29) \\
[0.25cm]
20.00 & MM & -242.21(8.20) & -424.68(10.11) & -795.63(15.56) & -1864.07(24.56) \\
& EM & -242.21(8.20) & -424.68(10.11) & -795.63(15.56) & -1864.07(24.56) \\
& FS & -242.21(8.20) & -424.68(10.11) & -795.63(15.56) & -1864.07(24.56) \\
& lme4 & -242.21(8.20) & -424.68(10.11) & -795.63(15.56) & -1864.07(24.56) \\
\bottomrule
\end{tabular}

To show summary statistics on the MSE:

In [6]:
showwhat = 7
include("anova_sim_print.jl")
\begin{tabular}{rrrrrr}
\toprule
$\sigma^2_{1}/\sigma^2_{e}$ & Method &
\multicolumn{4}{c}{$c=$ \# observations per combination} \\
\cline{3-6}
& & 5 & 10 & 20 & 50  \\
\midrule
0.00 & MM &  0.06(0.04) &  0.03(0.02) &  0.02(0.02) &  0.02(0.01) \\
& EM &  0.06(0.04) &  0.03(0.02) &  0.02(0.02) &  0.02(0.01) \\
& FS &  0.06(0.04) &  0.03(0.02) &  0.02(0.02) &  0.02(0.01) \\
& lme4 &  0.03(0.03) &  0.02(0.01) &  0.01(0.01) &  0.01(0.01) \\
[0.25cm]
0.05 & MM &  0.08(0.03) &  0.06(0.02) &  0.04(0.02) &  0.03(0.01) \\
& EM &  0.08(0.03) &  0.06(0.02) &  0.04(0.02) &  0.03(0.01) \\
& FS &  0.08(0.03) &  0.06(0.02) &  0.04(0.02) &  0.03(0.01) \\
& lme4 &  0.05(0.02) &  0.05(0.02) &  0.04(0.02) &  0.03(0.01) \\
[0.25cm]
0.10 & MM &  0.10(0.04) &  0.08(0.03) &  0.06(0.02) &  0.06(0.03) \\
& EM &  0.10(0.04) &  0.08(0.03) &  0.06(0.02) &  0.06(0.03) \\
& FS &  0.10(0.04) &  0.08(0.03) &  0.06(0.02) &  0.06(0.03) \\
& lme4 &  0.08(0.03) &  0.07(0.02) &  0.06(0.02) &  0.06(0.03) \\
[0.25cm]
1.00 & MM &  0.49(0.19) &  0.50(0.22) &  0.45(0.24) &  0.51(0.23) \\
& EM &  0.49(0.19) &  0.50(0.22) &  0.45(0.24) &  0.51(0.23) \\
& FS &  0.49(0.19) &  0.50(0.22) &  0.45(0.24) &  0.51(0.23) \\
& lme4 &  0.49(0.19) &  0.50(0.22) &  0.45(0.24) &  0.51(0.23) \\
[0.25cm]
10.00 & MM &  5.31(2.49) &  5.16(2.05) &  5.31(3.45) &  4.88(2.23) \\
& EM &  5.31(2.49) &  5.16(2.05) &  5.30(3.46) &  4.88(2.22) \\
& FS &  5.31(2.49) &  5.17(2.06) &  5.31(3.46) &  4.89(2.22) \\
& lme4 &  5.31(2.49) &  5.17(2.06) &  5.31(3.46) &  4.89(2.23) \\
[0.25cm]
20.00 & MM &  9.86(4.12) &  9.98(3.83) &  9.50(4.07) &  8.84(3.44) \\
& EM &  9.86(4.12) &  9.97(3.83) &  9.50(4.07) &  8.84(3.44) \\
& FS &  9.87(4.12) &  9.99(3.83) &  9.51(4.07) &  8.85(3.45) \\
& lme4 &  9.87(4.13) &  9.99(3.83) &  9.51(4.07) &  8.85(3.45) \\
\bottomrule
\end{tabular}