Mr.ash2 parameter estimation

Results for 9 simulations runs, using different initiate and update strategy

show_param <- function(tag2){
  f <- lapply(tags, get_files, tag2 = tag2)
  parf <- lapply(f, '[[', "par")
  param <- do.call(rbind, lapply(parf, function(x) t(read.table(x))[2:1,]))
  knitr::kable(param)
}

show_param(tag2s[1])

show_param(tag2s[2])

show_param(tag2s[3])

show_param(tag2s[4])

Regional mr.ash2s PIP overview

Take simulation 1 (NULL; expr-snp; expr-snp) as examples. We use region size 500kb and PIP cut off at 0.5 for SUSIE.

f <- get_files(tag= tags[1], tag2 = tag2s[1])
a <- read.table(f[["rpip"]], header = T)
plot(a$p0, a$rPIP, pch =19, col ='salmon', xlab = "position", ylab= "Sum of PIP")
grid()

PIP scatter plot

mr.ash2s PIP vs. susie PIP.

scatter_plot_PIP<- function(tag2){
  f <- lapply(tags, get_files, tag2 = tag2)
  mrashf <- lapply(f, '[[', "gmrash")
  names(mrashf) <- tags
  
  susief <- lapply(f, '[[', "gsusie")
  names(susief) <- tags

  .tagname <- function(x, flist){
    a <- read.table(flist[[x]], header =T)
    a[, "name"] <- paste0(x, ":", a[, "name"])
    a
  }
  mrashres <- do.call(rbind, lapply(tags, .tagname, flist = mrashf))
  susieres <- do.call(rbind, lapply(tags, .tagname, flist = susief))
 
  res <- merge(mrashres, susieres, by = "name", all = T)
  
  res <- res[complete.cases(res),]
  res <- rename(res, c("PIP" = "mr.ash_PIP", "pip" = "SUSIE_PIP", "pip.null" = "SUSIE_PIP_null") )
  res$ifcausal <- mapvalues(res$ifcausal, 
          from=c(0,1), 
          to=c("Non causal", "Causal"))
  
  fig1 <- plot_ly(data = res, x = ~ mr.ash_PIP, y = ~ SUSIE_PIP, color = ~ ifcausal, 
                 colors = c( "salmon", "darkgreen"))
  
  fig2 <- plot_ly(data = res, x = ~ mr.ash_PIP, y = ~ SUSIE_PIP_null, color = ~ ifcausal, 
                 colors = c( "salmon", "darkgreen"))
  
  fig <- subplot(fig1, fig2, titleX = TRUE, titleY = T, margin = 0.1)
  fig
}

scatter_plot_PIP(tag2s[1])

Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

scatter_plot_PIP(tag2s[2])

scatter_plot_PIP(tag2s[3])

scatter_plot_PIP(tag2s[4])

ROC curve

ROC_plot<- function(tag2){
  f <- lapply(tags, get_files, tag2 = tag2)
  mrashf <- lapply(f, '[[', "gmrash")
  names(mrashf) <- tags
  
  susief <- lapply(f, '[[', "gsusie")
  names(susief) <- tags
  
  gwasf <- lapply(f, '[[', "ggwas")
  names(gwasf) <- tags

  .tagname <- function(x, flist, colnames = NULL){
    a <- read.table(flist[[x]], header =T)
    if (!is.null(colnames)){
      colnames(a) <- colnames
    }
    a[, "name"] <- paste0(x, ":", a[, "name"])
    a
  }
  mrashres <- do.call(rbind, lapply(tags, .tagname, flist = mrashf))
  susieres <- do.call(rbind, lapply(tags, .tagname, flist = susief))
  gwasres <- do.call(rbind, lapply(tags, .tagname, flist = gwasf, 
                                   colnames =  c("chr", "p0",   "p1", "name", "Estimate", "Std.Error", "t-value", "PVALUE")))

  res <- merge(mrashres, susieres, by = "name", all = T)
  res <- merge(res, gwasres, by = "name", all = T)
  
  res <- res[complete.cases(res),]
  res <- rename(res, c("PIP" = "mr.ash", "pip" = "SUSIE", "PVALUE" = "TWAS") )
  res[,"TWAS"] <- -log10(res[, "TWAS"])
  
  roccolors <-  c("red", "green", "blue")
  methods <- c("mr.ash", "SUSIE", "TWAS")
  plot(0, xlim=c(0,1), ylim=c(0,1), col="white", xlab = "FPR", ylab = "TPR")
  for (i in 1:3){
    method <- methods[i]
    bordered <- res[order(res[,method]),] 
    actuals <- bordered$ifcausal == 1
    sens <- (sum(actuals) - cumsum(actuals))/sum(actuals)
    spec <- cumsum(!actuals)/sum(!actuals)
    lines(1 - spec, sens, type = "l", col = roccolors[i])
    abline(c(0,0),c(1,1))
    auc <- sum(spec*diff(c(0, 1 - sens)))
    cat("AUC for ", method, ": ", auc)
  }
  legend(0.6,0.3, legend= methods, col=roccolors, lty=1, cex=0.8)
  grid()
}

ROC_plot(tag2s[1])

AUC for  mr.ash :  0.8532693AUC for  SUSIE :  0.8479238AUC for  TWAS :  0.8474589

ROC_plot(tag2s[2])

AUC for  mr.ash :  0.8532693AUC for  SUSIE :  0.8479238AUC for  TWAS :  0.8474589

ROC_plot(tag2s[3])

AUC for  mr.ash :  0.7802647AUC for  SUSIE :  0.8690825AUC for  TWAS :  0.8674838

ROC_plot(tag2s[4])

AUC for  mr.ash :  0.7769784AUC for  SUSIE :  0.8691713AUC for  TWAS :  0.8674838

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] plyr_1.8.6          tidyr_0.8.3         plotly_4.9.2.9000  
[4] ggplot2_3.3.1       data.table_1.12.7   mr.ash.alpha_0.1-34

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6      highr_0.7         compiler_3.5.1   
 [4] pillar_1.4.4      later_0.7.5       git2r_0.26.1     
 [7] workflowr_1.6.0   tools_3.5.1       digest_0.6.25    
[10] viridisLite_0.3.0 jsonlite_1.6.1    evaluate_0.12    
[13] tibble_3.0.1      lifecycle_0.2.0   gtable_0.2.0     
[16] lattice_0.20-38   pkgconfig_2.0.2   rlang_0.4.6      
[19] Matrix_1.2-15     shiny_1.2.0       crosstalk_1.0.0  
[22] yaml_2.2.0        httr_1.4.1        withr_2.1.2      
[25] stringr_1.4.0     dplyr_1.0.0       knitr_1.20       
[28] htmlwidgets_1.3   generics_0.0.2    fs_1.3.1         
[31] vctrs_0.3.1       tidyselect_1.1.0  rprojroot_1.3-2  
[34] grid_3.5.1        glue_1.4.1        R6_2.3.0         
[37] rmarkdown_1.10    purrr_0.3.4       magrittr_1.5     
[40] backports_1.1.2   scales_1.0.0      promises_1.0.1   
[43] htmltools_0.3.6   ellipsis_0.3.1    xtable_1.8-3     
[46] mime_0.6          colorspace_1.3-2  httpuv_1.4.5     
[49] stringi_1.3.1     lazyeval_0.2.1    munsell_0.5.0    
[52] crayon_1.3.4