Genetics Simulation Example

This note instructs how to reproduce the genetics (Table 3) simulation example in the manuscript. This simulation compares the performance among MM, EM, FS (Fisher scoring) and lme4 under model: $$ \mathbf{y} \sim \text{Normal}(\mathbf{X} \beta, \sigma_a^2 \Phi + \sigma_e^2 \mathbf{I}). $$ Here the covariates $\mathbf{X}$ and kinship matrix $\Phi$ are from the example in Option 29 of the Mendel software. In simulation $\beta = [33, 6.0]$ and $\sigma_e^2$ is fixed at 2.0. $\sigma_a = r \cdot \sigma_e^2$, where $r$ varies in {0, 0.05, 0.1, 1, 10, 20}.

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)

Prerequisite

To run this document, we need Julia packages

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

and R packages MixedModels and lme4

install.packages("lme4")
install.packages("devtools")
devtools::install_github("jrklasen/MixedModels")

Run simulations

The script mendel29_sim.jl will run the simulations with $r = \sigma_a^2 / \sigma_e^2$ varying in {0, 0.05, 0.1, 1, 10, 20}, and nsim=50 replicates. It saves the results into a file mendel29_sim_result.jld:

In [2]:
tic()
include("mendel29_sim.jl")
toc()
WARNING: RCall.jl: Loading required package: lme4
Loading required package: Matrix
Loading required package: pedigreemm
vcratio = 0.0
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
vcratio = 0.05
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
vcratio = 0.1
vcratio = 1.0
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
vcratio = 10.0
WARNING: RCall.jl: Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
vcratio = 20.0
WARNING: RCall.jl: Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
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: 530.471809008 seconds
WARNING: RCall.jl: Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,  :
  failure to converge in 10000 evaluations
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
WARNING: RCall.jl: Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
   Hessian is numerically singular: parameters are not uniquely determined
Out[2]:
530.471809008

Display simulation results

To obtain Table 3 in manuscript:

In [3]:
include("mendel29_sim_print.jl")
\begin{tabular}{rrrrr}
\toprule
$\sigma^2_{a}/\sigma^2_{e}$ & Method &
Iteration & Runtime (ms) & Objective  \\
\midrule
 0.00 & MM & 198.02(102.23) & 133.61(822.67) & -375.59(9.63)  \\
 & EM & 1196.10(958.51) & 29.71(12.34) & -375.60(9.64)  \\
 & FS &  7.60(3.07) & 19.34(33.77) & -375.59(9.63) \\
& lme4 & -- & 401.02(142.04) & -375.59(9.64)  \\
[0.25cm]
 0.05 & MM & 185.86(99.41) & 17.26(1.76) & -377.39(10.52)  \\
 & EM & 1227.62(1030.07) & 29.82(12.74) & -377.40(10.52)  \\
 & FS &  7.84(2.74) & 14.97(1.55) & -377.39(10.52) \\
& lme4 & -- & 425.04(144.00) & -377.39(10.52)  \\
[0.25cm]
 0.10 & MM & 169.24(99.75) & 16.97(1.59) & -378.40(11.44)  \\
 & EM & 924.80(912.23) & 26.06(11.26) & -378.41(11.45)  \\
 & FS &  7.32(2.75) & 15.06(1.38) & -378.40(11.44) \\
& lme4 & -- & 435.14(128.87) & -378.40(11.44)  \\
[0.25cm]
 1.00 & MM & 58.96(23.69) & 15.53(0.75) & -409.54(10.90)  \\
 & EM & 105.10(79.65) & 15.49(0.96) & -409.54(10.90)  \\
 & FS &  5.80(1.05) & 14.66(0.89) & -409.54(10.90) \\
& lme4 & -- & 493.14(52.80) & -409.54(10.90)  \\
[0.25cm]
 10.00 & MM & 110.00(63.13) & 16.22(1.12) & -532.48(8.77)  \\
 & EM & 642.48(1470.38) & 22.32(18.37) & -532.57(8.75)  \\
 & FS & 14.98(5.21) & 14.78(0.97) & -531.72(8.92) \\
& lme4 & -- & 2897.12(15006.38) & -532.48(8.77)  \\
[0.25cm]
 20.00 & MM & 110.52(34.81) & 16.07(0.91) & -590.87(7.15)  \\
 & EM & 1014.22(1775.40) & 27.03(22.33) & -590.89(7.15)  \\
 & FS & 17.72(3.13) & 14.79(0.93) & -588.46(7.27) \\
& lme4 & -- & 5059.24(20692.67) & -590.79(7.15)  \\
\bottomrule
\end{tabular}