Truncated ASH

Introduction

Suppose we have \(m + n\) observations of \((\hat\beta_j, \hat s_j)\) in two groups such that, with a pre-specified \(t\).

\[ \text{Group 1: }(\hat\beta_1, \hat s_1), \ldots, (\hat\beta_m, \hat s_m), \text{with } |\hat\beta_j/\hat s_j| \leq t \]

\[ \text{Group 2: }(\hat\beta_{m+1}, \hat s_{m+1}), \ldots, (\hat\beta_{m+n}, \hat s_{m+n}), \text{with } |\hat\beta_j/\hat s_j| > t \]

For Group 1, we’ll only use the information that for each one, \(|\hat\beta_j/\hat s_j| \leq t\); that is, they are moderate observations. For Group 2, we’ll use the full observation \((\hat\beta_j, \hat s_j)\).

Now we have usual ASH assumptions that for both groups,

\[ \begin{array}{c} \hat\beta_j | \hat s_j, \beta_j \sim N(\beta_j, \hat s_j^2)\\ \beta_j \sim \sum_k\pi_k N(0, \sigma_k^2) \end{array} \]

Now the estimate of \(\pi_k\) from the marginal probability of the data in both groups is

\[ \max_{\{\pi_k\}}\sum\limits_{j=1}^m\log(\sum_k\pi_k(2\Phi(\frac{t\hat s_j}{\sqrt{\sigma_k^2+\hat s_j^2}}) - 1)) + \sum\limits_{l=1}^{n}\log(\sum_k\pi_k N(\hat\beta_{m + l}; 0, \sigma_k^2 + \hat s_{m + l}^2)) \]

where \(\Phi(\cdot)\) is the cdf of \(N(0, 1)\), and \(N(\cdot; 0, \sigma_k^2 + \hat s_{m + l}^2)\) is the pdf of \(N(0, \sigma_k^2+\hat s_j^2)\).

Note that

1. The second part of the objective function is the standard objective function for the regular ASH. The first part of the objective function is new, but it would be as if we are using a different likelihood, so it wouldn’t make any difference numerically.

2. The first part of the objective,

\[ \max_{\{\pi_k\}}\sum\limits_{j=1}^m\log(\sum_k\pi_k(2\Phi(\frac{t\hat s_j}{\sqrt{\sigma_k^2+\hat s_j^2}}) - 1)) \]

and its special case, when we set \(\hat s_j \equiv 1\),

\[ m\log(\sum_k\pi_k(2\Phi(\frac{t}{\sqrt{\sigma_k^2+1}}) - 1)) \]

will move more weight to \(\hat\pi_1\). Customarily, we assume \(\sigma_1 \leq \cdots \leq \sigma_K\), so to maximize this first part alone will be essentially to maximize \(\hat\pi_1\). In other words, if we only have Group 1 and no “extreme” observations in Group 2, the estimated prior of \(\beta\) will be a point mass at \(0\). This is true no matter whether we are using \(z\)-scores (thus making \(\hat s_j \equiv 1\) effectively) or not, although using \(z\)-scores will make the case more obviously.

Implementation

Here is the main truncash function, utilizing autoselect.mixsd, mixIP, normalmix from ashr.

truncash

Simulation

source("../code/truncash.R")

set.seed(777)

sebetahat = rep(1, 1000)
t = qnorm(0.975)

pihat0_ash = pihat0_truncash = mse_ash = mse_truncash = mse_mle = lfsr0.05_ash = lfsr0.05_truncash = c()


for (i in 1:100) {
beta1 = rep(0, 900)
betahat1 = rnorm(length(beta1), beta1, 1)
beta2 = rnorm(100, 0, 1)
betahat2 = rnorm(length(beta2), beta2, 1)
beta = c(beta1, beta2)
betahat = c(betahat1, betahat2)

res.ash = ash(betahat, 1, mixcompdist = "normal", method = "fdr")
res.truncash = truncash(betahat, sebetahat, t)

pihat0_ash[i] = get_fitted_g(res.ash)$pi[1]
pihat0_truncash[i] = get_fitted_g(res.truncash)$pi[1]

mse_ash[i] = mean((get_pm(res.ash) - beta)^2)
mse_truncash[i] = mean((get_pm(res.truncash) - beta)^2)
mse_mle[i] = mean((betahat - beta)^2)

lfsr0.05_ash = mean(get_lfsr(res.ash) <= 0.05)
lfsr0.05_truncash = mean(get_lfsr(res.truncash) <= 0.05)
}

plot(pihat0_ash, pihat0_truncash, xlab = "ash", ylab = "truncash", main = expression(hat(pi)[0]),
     xlim = c(min(c(pihat0_ash, pihat0_truncash)), max(c(pihat0_ash, pihat0_truncash))),
     ylim = c(min(c(pihat0_ash, pihat0_truncash)), max(c(pihat0_ash, pihat0_truncash)))
     )
abline(0, 1, lty = 3, col = "blue")
abline(v = 0.9, lty = 3, col = "blue")
abline(h = 0.9, lty = 3, col = "blue")
points(0.9, 0.9, col = "red", cex = 1.5, pch = 19)

plot(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9), xlab = "ash", ylab = "truncash",
     main = expression(paste("|", hat(pi)[0] - 0.9, "|")), 
     xlim = c(0, max(c(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9)))),
     ylim = c(0, max(c(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9))))
     )
abline(0, 1, lty = 3, col = "blue")
points(0.9, 0.9, col = "red", cex = 1.5, pch = 19)

plot(mse_ash, mse_truncash, xlab = "ash", ylab = "truncash", main = "MSE",
     xlim = c(min(c(mse_ash, mse_truncash)), max(c(mse_ash, mse_truncash))),
     ylim = c(min(c(mse_ash, mse_truncash)), max(c(mse_ash, mse_truncash)))
     )
abline(0, 1, lty = 3, col = "blue")

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.3

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] qvalue_2.4.2      edgeR_3.14.0      limma_3.28.5      SQUAREM_2016.10-1
[5] ashr_2.1.2        workflowr_0.3.0   rmarkdown_1.3    

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.9       knitr_1.15.1      magrittr_1.5     
 [4] splines_3.3.2     REBayes_0.62      MASS_7.3-45      
 [7] munsell_0.4.3     doParallel_1.0.10 pscl_1.4.9       
[10] colorspace_1.2-6  lattice_0.20-34   foreach_1.4.3    
[13] plyr_1.8.3        stringr_1.1.0     tools_3.3.2      
[16] parallel_3.3.2    grid_3.3.2        gtable_0.2.0     
[19] git2r_0.18.0      htmltools_0.3.5   iterators_1.0.8  
[22] assertthat_0.1    yaml_2.1.14       rprojroot_1.2    
[25] digest_0.6.11     Matrix_1.2-7.1    reshape2_1.4.1   
[28] ggplot2_2.1.0     codetools_0.2-15  evaluate_0.10    
[31] stringi_1.1.2     scales_0.4.0      Rmosek_7.1.3     
[34] backports_1.0.5   truncnorm_1.0-7