binomial thinning ultimate

workflowr

Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  analysis/contrainedclustering.Rmd
    Untracked:  analysis/gsea.Rmd
    Untracked:  analysis/ideas.Rmd
    Untracked:  analysis/methylation.Rmd
    Untracked:  analysis/susie.Rmd
    Untracked:  code/sccytokines.R
    Untracked:  code/scdeCalibration.R
    Untracked:  data/bart/
    Untracked:  data/cytokine/DE_controls_output_filter10.RData
    Untracked:  data/cytokine/DE_controls_output_filter10_addlimma.RData
    Untracked:  data/cytokine/README
    Untracked:  data/cytokine/test.RData
    Untracked:  data/cytokine_normalized.RData
    Untracked:  data/deconv/
    Untracked:  data/scde/

Unstaged changes:
    Modified:   analysis/deconvolution.Rmd
    Modified:   analysis/index.Rmd
    Deleted:    data/mout_high.RData
    Deleted:    data/scCDT.RData
    Deleted:    data/sva_sva_high.RData

Let’s start with the simplest case: we have a count \(Y\) matrix, whose all entries are the same! Then we apply binomial thinning to \(Y\). Half of samples are from group 1 and the rest are from group 2. \(90\%\) of genes have no effects and the rest have effects generated from standard normal distribution.

In the following simulation, we set the number of samples to be \(n=500\), the number of genes \(p=5000\) and entries of \(Y\) vary from 100, 50, 30, 20, 10, 5, 1.

We fit a simple linear regression to each gene and compare the log-likelihood of ash, setting \(\alpha=0\) and \(\alpha=1\), and plot \(\beta\) vs \(\hat\beta\) as well as \(\beta\) vs \(s.e.(\hat\beta)\) for each case.

#'@param Z count matrix, sample by features
#'@param x 1 for group 1, 0 for group 2
#'@param beta effect of fearures,  0 for null.
#'@return W, thinned matrix
bi_thin = function(Z,x,beta){
  
  n=nrow(Z)
  p=ncol(Z)
  
  # group index
  g1 = which(x==1)
  g2 = which(x==0)
  
  
  p2 = 1/(1+exp(beta))
  p1 = 1-p2
  P = matrix(nrow = n,ncol = p)
  P[g1,] = t(replicate(length(g1),p1))
  P[g2,] = t(replicate(length(g2),p2))
  
  W = matrix(rbinom(n*p,Z,P),nrow=n)
  
  W
}

library(limma)

n=500
p=5000
set.seed(12345)

loglik0 = c()
auc0 = c()
loglik1 = c()
auc1 = c()

par(mfrow = c(1,2))
for(ni in c(100,50,30,20,10,5,1)){
  Y = matrix(rep(ni,n*p),nrow=p,ncol=n)

  group_idx = c(rep(1,n/2),rep(0,n/2))
  X = model.matrix(~group_idx)

  b = rnorm(p)
  b[sample(1:p,0.9*p)] = 0
  which_null = 1*(b==0)

  W = bi_thin(t(Y),group_idx,b)

  lmout <- limma::lmFit(object = t(log(W+0.5)), design = X)

  ash0 = ashr::ash(lmout$coefficients[, 2],lmout$stdev.unscaled[, 2]*lmout$sigma,alpha=0)
  loglik0 = c(loglik0,ash0$loglik)
  auc0 = c(auc0,pROC::roc(which_null,ash0$result$lfsr)$auc)

  ash1 = ashr::ash(lmout$coefficients[, 2],lmout$stdev.unscaled[, 2]*lmout$sigma,alpha=1)
  loglik1 = c(loglik1,ash1$loglik)
  auc1 = c(auc1,pROC::roc(which_null,ash1$result$lfsr)$auc)
  
  
  plot((b),lmout$coefficients[, 2],xlab='beta',ylab='estimated',col='grey60',main=paste('Y:',ni))
  abline(a=0,b=1)
  plot(abs(b),lmout$stdev.unscaled[, 2]*lmout$sigma,xlab='absolute value of beta',ylab='sd of beta_hat',main=paste('Y:',ni))
}

knitr::kable(cbind(c(100,50,30,20,10,5,1),loglik0,loglik1),col.names = c('Y','alpha0','alpha1'),
             caption = 'log likelihood')

knitr::kable(cbind(c(100,50,30,20,10,5,1),auc0,auc1),col.names = c('Y','alpha0','alpha1'),
             caption = 'AUC')

So when counts in \(Y\) are large, setting \(\alpha=1\) gives higher likelihood. Estimated \(\beta\)s are close to true \(\beta\)s and apparently, scale of \(\beta\) and \(var(\hat\beta)\) are positively correlated. While when counts in \(Y\) are small, things are opposite.

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] limma_3.38.2

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.2        highr_0.7         pillar_1.3.1     
 [4] plyr_1.8.4        compiler_3.5.1    later_0.7.5      
 [7] git2r_0.26.1      workflowr_1.6.0   iterators_1.0.10 
[10] tools_3.5.1       digest_0.6.18     tibble_2.1.1     
[13] evaluate_0.12     gtable_0.2.0      lattice_0.20-38  
[16] pkgconfig_2.0.2   rlang_0.4.0       Matrix_1.2-15    
[19] foreach_1.4.4     yaml_2.2.0        parallel_3.5.1   
[22] dplyr_0.8.0.1     stringr_1.3.1     knitr_1.20       
[25] pROC_1.13.0       fs_1.3.1          tidyselect_0.2.5 
[28] rprojroot_1.3-2   grid_3.5.1        glue_1.3.0       
[31] R6_2.3.0          rmarkdown_1.10    mixsqp_0.2-2     
[34] purrr_0.3.2       ggplot2_3.1.1     ashr_2.2-39      
[37] magrittr_1.5      whisker_0.3-2     backports_1.1.2  
[40] scales_1.0.0      promises_1.0.1    codetools_0.2-15 
[43] htmltools_0.3.6   MASS_7.3-51.1     assertthat_0.2.0 
[46] colorspace_1.3-2  httpuv_1.4.5      stringi_1.2.4    
[49] lazyeval_0.2.1    doParallel_1.0.14 pscl_1.5.2       
[52] munsell_0.5.0     truncnorm_1.0-8   SQUAREM_2017.10-1
[55] crayon_1.3.4