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 matrix \(Y\), whose all entries are the same! Then we apply binomial thinning to \(Y\).

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.\) Half of samples are from group 1 and the rest are from group 2. About \(90\%\) of genes have no effects and the rest have effects \(\beta_j\) generated from standard normal distribution.

We fit a simple linear regression to each gene with log transformed counts and compare the log-likelihood from 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('count:',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('count:',ni))
}

Past versions of unnamed-chunk-1-1.png

Past versions of unnamed-chunk-1-2.png

Past versions of unnamed-chunk-1-3.png

Past versions of unnamed-chunk-1-4.png

Past versions of unnamed-chunk-1-5.png

Past versions of unnamed-chunk-1-6.png

Past versions of unnamed-chunk-1-7.png

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