Last updated: 2017-01-15

Code version: 94a3347615361216b39b8a9333bfa8354794e0ff

The purpose of this analysis is to check the performance of the non-zero mode option.

First, we load the necessary libraries.

library(ashr)

Simple simulation

I conjecture that the mean of the posterior means should be close to the optimum for the mode… maybe even equal to it. (That is, this would constitute a fixed point for the update. We aren’t explicitly using that directly in the current implementation; the uniform mixture uses optim to do it numerically; the normal mixture uses a true EM I think…)

check_mode = function(betahat, sebetahat, mixcompdist) {
  z.ash = ash(betahat,sebetahat,mixcompdist = mixcompdist,mode = "estimate")
  average.posteriormean = mean(get_pm(z.ash))
  fitted.mode           = comp_mean(get_fitted_g(z.ash))[1]

  # Refit to get g.
  z.ash1 = ash(betahat - fitted.mode,sebetahat,mixcompdist = mixcompdist)
  g      = get_fitted_g(z.ash1)
  loglik = get_loglik(z.ash1)

  loglik.down          = ash(z - fitted.mode - 0.01,1,g = g)$loglik
  loglik.up            = ash(z - fitted.mode + 0.01,1,g = g)$loglik
  loglik.posteriormean = ash(z - average.posteriormean,1,g = g)$loglik

  return(list(fitted.mode = fitted.mode,
              average.posteriormean = average.posteriormean,
              loglik = c(loglik,loglik.down,loglik.up,loglik.posteriormean)))
}

set.seed(100)
z = rnorm(1000) + 3
print(check_mode(z,1,mixcompdist = "uniform"))
$fitted.mode
[1] 3.019

$average.posteriormean
[1] 3.018

$loglik
[1] -1449 -1449 -1449 -1449
print(check_mode(z,1,mixcompdist = "normal"))
$fitted.mode
[1] 3.017

$average.posteriormean
[1] 3.017

$loglik
[1] -1449 -1449 -1449 -1449
print(check_mode(z,1,mixcompdist = "halfuniform"))
$fitted.mode
[1] 3.074

$average.posteriormean
[1] 3.024

$loglik
[1] -1449 -1449 -1449 -1449

An additional experiment:

set.seed(100)
beta     = rexp(1000)
betahat  = beta + rnorm(1000,0,0.1)
z.ash.hu = ash(betahat,0.1,mixcompdist="halfuniform",outputlevel=4,method="shrink")
Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended
z.ash.pu = ash(betahat,0.1,mixcompdist="+uniform",outputlevel=4,method="shrink")

z.ash.hu2    = ash(betahat - 0.2,0.1,mixcompdist = "halfuniform",outputlevel = 4,
                   method = "shrink")
Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended
Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced

Warning in log(m$pi): NaNs produced
z.ash.pu2    = ash(betahat - 0.1,0.1,mixcompdist = "+uniform",outputlevel = 4,
                   method = "shrink")
z.ash.hu.nzm = ash(betahat,0.1,mixcompdist = "halfuniform",mode = "estimate",
                   method = "shrink")
Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended
Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended
Warning in min(gradient(matrix_lik) + prior[1] - 1, na.rm = TRUE): no non-
missing arguments to min; returning Inf
Warning in stats::optimize(test.op, interval = c(min(betahat),
max(betahat))): NA/Inf replaced by maximum positive value
Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended

Warning in check_args(mixcompdist, df, prior, optmethod, gridmult,
sebetahat, : Use of halfuniform without nullbiased prior can lead to
misleading local false sign rates, and so is not recommended
print(z.ash.hu$loglik,digits = 8)
[1] -1051.1457
print(z.ash.hu.nzm$loglik,digits = 8)
[1] -1048.6462

Note to self: check that the normal version works too; check nonzeromodeEMobj as it doesn’t seem needed.

Session information

sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.2

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] knitr_1.15.1 ashr_2.0.5  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.8       magrittr_1.5      REBayes_0.63     
 [4] MASS_7.3-45       doParallel_1.0.10 pscl_1.4.9       
 [7] SQUAREM_2016.8-2  lattice_0.20-34   foreach_1.4.3    
[10] stringr_1.1.0     tools_3.3.2       parallel_3.3.2   
[13] grid_3.3.2        htmltools_0.3.5   iterators_1.0.8  
[16] assertthat_0.1    yaml_2.1.14       rprojroot_1.1    
[19] digest_0.6.10     Matrix_1.2-7.1    codetools_0.2-15 
[22] evaluate_0.10     rmarkdown_1.3     stringi_1.1.2    
[25] Rmosek_7.1.3      backports_1.0.4   truncnorm_1.0-7