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:

Iterations
Runtime
Objective value
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}$.

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:

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

Display simulation results¶

To obtain Table 3 in manuscript:

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}