Strange error ashr
Hello,
I was using ashr on simingly non problematic data set (attached).
I got the following error
df <- as.data.frame(read_csv("strange_ash_error.csv"))
ash(df$x, df$sd)
[strange_ash_error.csv](https://github.com/user-attachments/files/18967160/strange_ash_error.csv)
Error in if (!all(nonzero_cols)) { :
missing value where TRUE/FALSE needed
In addition: Warning message:
In log(1 - exp(lpb - lpa)) : NaNs produced
I am a bit surprised that this cause problem given that
hist( df$x/df$sd, nclass=100) seems not to pathologic
here the session Info
sessionInfo()
R version 4.3.3 (2024-02-29 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 26120)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8
[2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8
time zone: Etc/GMT-1
tzcode source: internal
attached base packages:
[1] grid stats4 stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] readr_2.1.5
[2] fsusieR_0.2.94
[3] testthat_3.2.3
[4] TxDb.Hsapiens.UCSC.hg38.knownGene_3.18.0
[5] GenomicFeatures_1.54.4
[6] Gviz_1.46.1
[7] GenomicRanges_1.54.1
[8] GenomeInfoDb_1.38.8
[9] org.Hs.eg.db_3.18.0
[10] AnnotationDbi_1.64.1
[11] IRanges_2.36.0
[12] S4Vectors_0.40.2
[13] Biobase_2.62.0
[14] AnnotationHub_3.10.1
[15] BiocFileCache_2.10.2
[16] dbplyr_2.5.0
[17] BiocGenerics_0.48.1
loaded via a namespace (and not attached):
[1] later_1.4.1 BiocIO_1.12.0
[3] bitops_1.0-9 filelock_1.0.3
[5] tibble_3.2.1 XML_3.99-0.18
[7] rpart_4.1.23 lifecycle_1.0.4
[9] mixsqp_0.3-54 rprojroot_2.0.4
[11] vroom_1.6.5 lattice_0.22-5
[13] ensembldb_2.26.1 MASS_7.3-60.0.1
[15] backports_1.5.0 magrittr_2.0.3
[17] Hmisc_5.2-2 rmarkdown_2.29
[19] yaml_2.3.10 remotes_2.5.0
[21] httpuv_1.6.15 sessioninfo_1.2.3
[23] pkgbuild_1.4.6 cowplot_1.1.3
[25] DBI_1.2.3 RColorBrewer_1.1-3
[27] abind_1.4-8 pkgload_1.4.0
[29] zlibbioc_1.48.2 purrr_1.0.2
[31] AnnotationFilter_1.26.0 biovizBase_1.50.0
[33] RCurl_1.98-1.16 nnet_7.3-20
[35] VariantAnnotation_1.48.1 rappdirs_0.3.3
[37] GenomeInfoDbData_1.2.11 irlba_2.3.5.1
[39] wavethresh_4.7.3 codetools_0.2-19
[41] DelayedArray_0.28.0 xml2_1.3.6
[43] tidyselect_1.2.1 matrixStats_1.5.0
[45] base64enc_0.1-3 GenomicAlignments_1.38.2
[47] ellipsis_0.3.2 Formula_1.2-5
[49] tools_4.3.3 progress_1.2.3
[51] Rcpp_1.0.14 glue_1.8.0
[53] gridExtra_2.3 SparseArray_1.2.4
[55] xfun_0.50 MatrixGenerics_1.14.0
[57] usethis_3.1.0 dplyr_1.1.4
[59] withr_3.0.2 BiocManager_1.30.25
[61] fastmap_1.2.0 latticeExtra_0.6-30
[63] caTools_1.18.3 digest_0.6.37
[65] truncnorm_1.0-9 R6_2.6.1
[67] mime_0.12 colorspace_2.1-1
[69] jpeg_0.1-10 dichromat_2.0-0.1
[71] biomaRt_2.58.2 RSQLite_2.3.9
[73] generics_0.1.3 data.table_1.16.4
[75] rtracklayer_1.62.0 prettyunits_1.2.0
[77] httr_1.4.7 htmlwidgets_1.6.4
[79] S4Arrays_1.2.1 pkgconfig_2.0.3
[81] gtable_0.3.6 blob_1.2.4
[83] XVector_0.42.0 brio_1.1.5
[85] htmltools_0.5.8.1 profvis_0.4.0
[87] smashr_1.2-9 ProtGenerics_1.34.0
[89] scales_1.3.0 png_0.1-8
[91] ashr_2.2-63 knitr_1.49
[93] rstudioapi_0.17.1 tzdb_0.4.0
[95] rjson_0.2.23 checkmate_2.3.2
[97] curl_6.2.0 cachem_1.1.0
[99] stringr_1.5.1 BiocVersion_3.18.1
[101] parallel_4.3.3 miniUI_0.1.1.1
[103] foreign_0.8-86 RcppZiggurat_0.1.6
[105] restfulr_0.0.15 desc_1.4.3
[107] pillar_1.10.1 vctrs_0.6.5
[109] urlchecker_1.0.1 promises_1.3.2
[111] xtable_1.8-4 cluster_2.1.6
[113] htmlTable_2.4.3 evaluate_1.0.3
[115] invgamma_1.1 cli_3.6.3
[117] compiler_4.3.3 Rsamtools_2.18.0
[119] rlang_1.1.5 crayon_1.5.3
[121] SQUAREM_2021.1 LaplacesDemon_16.1.6
[123] interp_1.1-6 fs_1.6.5
[125] stringi_1.8.4 deldir_2.0-4
[127] BiocParallel_1.36.0 munsell_0.5.1
[129] Biostrings_2.70.3 lazyeval_0.2.2
[131] devtools_2.4.5 Matrix_1.6-5
[133] BSgenome_1.70.2 hms_1.1.3
[135] bit64_4.6.0-1 ggplot2_3.5.1
[137] KEGGREST_1.42.0 shiny_1.10.0
[139] SummarizedExperiment_1.32.0 interactiveDisplayBase_1.40.0
[141] Rfast_2.1.4 memoise_2.0.1
[143] RcppParallel_5.1.10 bit_4.5.0.1
@william-denault I'm not quite sure, but I'll note that
fit <- ash(df$x,df$sd,mixcompdist = "normal")
worked for me.
@william-denault @pcarbo There is a floating point issue when evaluating log_comp_dens_conv.unimix on this data set:
Browse[1]> mean(abs(lpa - lpb) < .Machine$double.eps)
[1] 0.5059023
As a result, the log likelihood matrix is filled with NaN and Inf.
@william-denault The data are weird: you have many data points that are very close to zero and also have very small s.e.'s. There could be something wrong with these data.