Set a global default and user-overrides for the maxdf parameter in huge.mb

I'm estimating networks using small samples sizes (10-20) and low OTU richness (<200). I start by adjusting lambda.min.ratio so that stars/pulsar are able to successfully complete sparsity selection and output a network. Then, I try to increase nlambda to improve stability (getStability(graph) near 0.05). However, each time I increase nlambda beyond a point (varies by network, but not higher than 50), stars/pulsar return an empty matrix. My understanding is that nlambda controls the number of times stars/pulsar sample the lambda path (all values of lambda between lambda.min.ratio and 1). Why would increasing the number of lambda values examined cause stars/pulsar to fail when lambda.min.ratio is the same and stars/pulsar are able to select a network at smaller nlambda? Perhaps I am misunderstanding nlambda or perhaps my data are too noisy? I'm able to improve stability by adjusting lambda.min.ratio further and keeping nlambda low. Below is an example from one of my datasets.

Summary of Phyloseq Object

BM.fragment.list[["BM.REGUA"]] phyloseq-class experiment-level object otu_table() OTU Table: [ 162 taxa and 22 samples ] sample_data() Sample Data: [ 22 samples by 28 sample variables ] tax_table() Taxonomy Table: [ 162 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 162 tips and 160 internal nodes ]

Set Pulsar Parameters

pargs<- list(seed=30010, rep.num=50, ncores=2)

Original Attempt at Network

seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=5e-3, nlambda=50, pulsar.params=pargs) getStability(seBM.REGUA) 0.042795 seBM.REGUA$select$stars$summary [1] 0.00000000 0.01716752 0.02599298 0.03428142 0.04279500 0.05075881 0.05852226 [8] 0.06562417 0.07306078 0.07947427 0.08551123 0.09147114 0.09816167 0.10332705 [15] 0.10812435 0.11379864 0.11891407 0.12431427 0.12845150 0.13389808 0.13965162 [22] 0.14462964 0.15349351 0.16050021 0.16632712 0.17323199 0.17987332 0.19066207 [29] 0.20544106 0.21202509 0.21757865 0.22154215 0.22903716 0.24647926 0.24713160 [36] 0.25737118 0.26442635 0.27172481 0.27824958 0.28186219 0.28858867 0.29670265 [43] 0.30500299 0.31137563 0.31749185 0.32746065 0.33611410 0.33907130 0.34439887 [50] 0.34942020 seBM.REGUA$lambda [1] 1.000000000 0.897511876 0.805527568 0.722970560 0.648874664 0.582372717 [7] 0.522686430 0.469117279 0.421038329 0.377886901 0.339157981 0.304398316 [13] 0.273201104 0.245201236 0.220071021 0.197516355 0.177273274 0.159104869 [19] 0.142798510 0.128163358 0.115028136 0.103239118 0.092658335 0.083161956 [25] 0.074638843 0.066989248 0.060123646 0.053961686 0.048431254 0.043467626 [31] 0.039012711 0.035014371 0.031425814 0.028205041 0.025314359 0.022719938 [37] 0.020391414 0.018301537 0.016425846 0.014742392 0.013231472 0.011875403 [43] 0.010658316 0.009565965 0.008585567 0.007705648 0.006915911 0.006207112 [49] 0.005570957 0.005000000

Increase nlambda by 10 to improve stability, network selection fails

seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=5e-3, nlambda=60, pulsar.params=pargs) Applying data transformations... Selecting model with pulsar using stars... Fitting final estimate with mb... done Warning messages: 1: 1 job had warning: "1 parallel function call did not deliver a result" 2: In pulsar(data = X, fun = match.fun(estFun), fargs = args, seed = 30010, : Optimal lambda may be smaller than the supplied values getStability(seBM.REGUA) 0 seBM.REGUA$select$stars$summary [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [41] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 seBM.REGUA$lambda [1] 1.000000000 0.914112171 0.835601062 0.763833101 0.698229135 0.638259750 [7] 0.583441006 0.533330525 0.487523924 0.445651553 0.407375509 0.372386911 [13] 0.340403408 0.311166898 0.284441449 0.260011391 0.237679577 0.217265794 [19] 0.198605307 0.181547528 0.165954805 0.151701307 0.138672011 0.126761774 [25] 0.115874480 0.105922273 0.096824839 0.088508763 0.080906938 0.073958017 [31] 0.067605923 0.061799397 0.056491581 0.051639642 0.047204425 0.043150140 [37] 0.039444068 0.036056303 0.032959505 0.030128685 0.027540997 0.025175561 [43] 0.023013287 0.021036725 0.019229927 0.017578310 0.016068547 0.014688455 [49] 0.013426895 0.012273688 0.011219528 0.010255907 0.009375049 0.008569847 [55] 0.007833801 0.007160973 0.006545933 0.005983717 0.005469788 0.005000000

Set nlambda back to 50, adjust lambda.min.ratio to improve stability

seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=1e-2, nlambda=50, pulsar.params=pargs) getStability(seBM.REGUA) 0.04987851 seBM.REGUA$select$stars$summary [1] 0.00000000 0.01699146 0.02581667 0.03391274 0.04184624 0.04987851 0.05764245 [8] 0.06473111 0.07170837 0.07853385 0.08456394 0.09047306 0.10101544 0.10155283 [15] 0.10676040 0.11177315 0.11673436 0.12127463 0.12609916 0.13093438 0.13505124 [22] 0.13974082 0.14445481 0.14884576 0.15381630 0.15879948 0.16497238 0.17003874 [29] 0.17601325 0.18260686 0.18987513 0.19746302 0.21383759 0.21907927 0.22542394 [36] 0.23146299 0.23705814 0.24317730 0.24749636 0.25743658 0.26383645 0.27007782 [43] 0.27579492 0.28320540 0.28999666 0.29574677 0.30125444 0.30896539 0.31419345 [50] 0.32060768 seBM.REGUA$lambda [1] 1.00000000 0.91029818 0.82864277 0.75431201 0.68664885 0.62505519 0.56898660 [8] 0.51794747 0.47148664 0.42919343 0.39069399 0.35564803 0.32374575 0.29470517 [15] 0.26826958 0.24420531 0.22229965 0.20235896 0.18420700 0.16768329 0.15264180 [22] 0.13894955 0.12648552 0.11513954 0.10481131 0.09540955 0.08685114 0.07906043 [29] 0.07196857 0.06551286 0.05963623 0.05428675 0.04941713 0.04498433 0.04094915 [36] 0.03727594 0.03393222 0.03088844 0.02811769 0.02559548 0.02329952 0.02120951 [43] 0.01930698 0.01757511 0.01599859 0.01456348 0.01325711 0.01206793 0.01098541 [50] 0.01000000

Thanks for the report - would it be possible to send me your data via email so I can reproduce this bug? It seems that (due to the small sample size?) there's a toxic subsample that gets generated when nlambda gets large enough - and parallel can't collect the results.

Actually, it should be possible to exclude that result using batch mode pulsar, which is set to ignored failed jobs, at the risk of biasing the StARS estimate.

Just an update on this: I'm able to recreate a warning as a result of the same underlying issue #111. I'm going to have to reach out the huge package developers to see if there's a fix.

Thanks for the update! I don't quite understand the underlying issue in 111. Here is my current understanding: Small sample sizes like I have (and like those discussed in 111) produce a toxic subsample, because the exclusion of a single sample in small datasets can have a large impact on the estimated network (my case) and if sample size is small enough, the subsampling routine might not have enough samples to estimate a network (issue 111).

What sort of fix would prevent this issue? Is the fix also related to issue 197?

If I'm able to adjust lambda.min.ratio to reach reasonable stability, but nlambda is < 50, is the resulting network reliable? Especially given my lambda.min.ratios aren't very small (>1e-3), so the lambda path is well-sampled.

I was able to get batch mode pulsar to run as you suggested above. Increasing nlambda to 100 didn't change the results much compared to nlambda=50 that I report above.

bargs<- list(rep.num=50, seed=30010, conffile="/Library/Frameworks/R.framework/Versions/4.1/Resources/library/pulsar/config/batchtools.conf.parallel.R")

seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=1e-3, nlambda=100, sel.criterion='stars', pulsar.select='batch', pulsar.params=bargs) Applying data transformations... Selecting model with batch.pulsar using stars... Sourcing configuration file '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/pulsar/config/batchtools.conf.parallel.R' ... Created registry in '/var/folders/p0/bxmttvw95r3bl3752w_ggcy00000gp/T/RtmpplnGU0/registry1866deadccb' using cluster functions 'Multicore' Adding 51 jobs ... Submitting 51 jobs in 51 chunks using cluster functions 'Multicore' ... Fitting final estimate with mb...
done Warning messages: 1: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 2: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 3: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 4: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 5: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 6: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 7: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 8: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 9: In parallel::mccollect(jobs, wait = FALSE, timeout = timeout) : 1 parallel job did not deliver a result 10: In batch.pulsar(data = X, fun = match.fun(estFun), fargs = args, : Only c(1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 17, 18, 20, 22, 25, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 40, 41, 42, 43, 44, 46, 47, 48, 49, 50, 51) jobs completed... proceeding anyway

getStability(seBM.REGUA) [1] 0.04937195

sum(getRefit(seBM.REGUA))/2 [1] 150

seBM.REGUA$select$stars$summary [1] 0.00000000 0.01486877 0.02268013 0.02944632 0.03629645 0.04274735 [7] 0.04937195 0.05540974 0.06130339 0.07146418 0.07217287 0.07733592 [13] 0.08222345 0.08632565 0.09671707 0.09708894 0.10080284 0.10359567 [19] 0.10647023 0.11174459 0.11762108 0.11852976 0.11894724 0.12119120 [25] 0.12266188 0.12790886 0.12878798 0.13476994 0.13528121 0.13648366 [31] 0.13980746 0.14128068 0.14411365 0.14691523 0.15062476 0.15300958 [37] 0.15686545 0.16098954 0.16507458 0.16996211 0.17359136 0.17988095 [43] 0.18405102 0.18987241 0.19581605 0.20115171 0.20614361 0.21201281 [49] 0.21698573 0.22268815 0.22649949 0.23815358 0.24384578 0.24920298 [55] 0.25417736 0.25816825 0.26275617 0.26874616 0.27380776 0.27882082 [61] 0.28171327 0.28609427 0.29131097 0.29593647 0.30022442 0.30478789 [67] 0.30915758 0.31258611 0.31629053 0.32100435 0.32572072 0.32975176 [73] 0.33475023 0.33694383 0.34054534 0.34393044 0.34715753 0.34966790 [79] 0.35250049 0.35608558 0.35795805 0.36138475 0.36379768 0.36616426 [85] 0.36799147 0.36965082 0.37050585 0.37190317 0.37329356 0.37509815 [91] 0.37669691 0.37875622 0.37920947 0.38107281 0.38214534 0.38330290 [97] 0.38472249 0.38626761 0.38768464 0.38921735

seBM.REGUA$lambda [1] 1.000000000 0.932603347 0.869749003 0.811130831 0.756463328 0.705480231 [7] 0.657933225 0.613590727 0.572236766 0.533669923 0.497702356 0.464158883 [13] 0.432876128 0.403701726 0.376493581 0.351119173 0.327454916 0.305385551 [19] 0.284803587 0.265608778 0.247707636 0.231012970 0.215443469 0.200923300 [25] 0.187381742 0.174752840 0.162975083 0.151991108 0.141747416 0.132194115 [31] 0.123284674 0.114975700 0.107226722 0.100000000 0.093260335 0.086974900 [37] 0.081113083 0.075646333 0.070548023 0.065793322 0.061359073 0.057223677 [43] 0.053366992 0.049770236 0.046415888 0.043287613 0.040370173 0.037649358 [49] 0.035111917 0.032745492 0.030538555 0.028480359 0.026560878 0.024770764 [55] 0.023101297 0.021544347 0.020092330 0.018738174 0.017475284 0.016297508 [61] 0.015199111 0.014174742 0.013219411 0.012328467 0.011497570 0.010722672 [67] 0.010000000 0.009326033 0.008697490 0.008111308 0.007564633 0.007054802 [73] 0.006579332 0.006135907 0.005722368 0.005336699 0.004977024 0.004641589 [79] 0.004328761 0.004037017 0.003764936 0.003511192 0.003274549 0.003053856 [85] 0.002848036 0.002656088 0.002477076 0.002310130 0.002154435 0.002009233 [91] 0.001873817 0.001747528 0.001629751 0.001519911 0.001417474 0.001321941 [97] 0.001232847 0.001149757 0.001072267 0.001000000

Can you try the hugemb-patch branch? It's hacky work-around but at seems to resolve the error your data. https://github.com/zdk123/SpiecEasi/tree/hugemb-patch

Thank you for this work-around! I was able to get nlambda=50 and nlambda=60 to work without triggering the warning! Simply for testing things out, I also tested other nlambda values. I was able to trigger the warning again by increasing nlamba to 100, but nlambda values between 50 and 90 all seem to work fine! Below is what I did:

Install hugemb-patch and Load SpiecEasi

library(devtools) install_github("zdk123/SpiecEasi", ref="hugemb-patch") library(SpiecEasi)

Try SpiecEasi at initial parameters reported above, nlambda =50, lambda.min.ratio=5e-3 NOTE: The stability and number of edges are different than what I initially reported, which I took as confirmation that I had successfully installed the hugemb-patch

pargs<- list(rep.num=50, seed=30010, ncores=2) seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=5e-3, nlambda=50, pulsar.params=pargs) getStability(seBM.REGUA) [1] 0.042795 sum(getRefit(seBM.REGUA))/2 [1] 147 seBM.REGUA$select$stars$summary [1] 0.00000000 0.01716752 0.02599298 0.03428142 0.04279500 0.05075881 0.05852226 [8] 0.06562417 0.07306078 0.07947427 0.08551123 0.09147114 0.09816167 0.10332705 [15] 0.10812435 0.11379864 0.11891407 0.12431427 0.12845150 0.13389808 0.13965162 [22] 0.14462964 0.15349351 0.16050021 0.16632712 0.17323199 0.17987332 0.19066207 [29] 0.20544106 0.21202509 0.21757865 0.22154215 0.22903716 0.24647926 0.24713160 [36] 0.25737118 0.26442635 0.27172481 0.27824958 0.28186219 0.28858867 0.29670265 [43] 0.30500299 0.31137563 0.31749185 0.32746065 0.33611410 0.33907130 0.34439887 [50] 0.34942020 seBM.REGUA$lambda [1] 1.000000000 0.897511876 0.805527568 0.722970560 0.648874664 0.582372717 [7] 0.522686430 0.469117279 0.421038329 0.377886901 0.339157981 0.304398316 [13] 0.273201104 0.245201236 0.220071021 0.197516355 0.177273274 0.159104869 [19] 0.142798510 0.128163358 0.115028136 0.103239118 0.092658335 0.083161956 [25] 0.074638843 0.066989248 0.060123646 0.053961686 0.048431254 0.043467626 [31] 0.039012711 0.035014371 0.031425814 0.028205041 0.025314359 0.022719938 [37] 0.020391414 0.018301537 0.016425846 0.014742392 0.013231472 0.011875403 [43] 0.010658316 0.009565965 0.008585567 0.007705648 0.006915911 0.006207112 [49] 0.005570957 0.005000000

Increase nlambda to 60

seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=5e-3, nlambda=60, pulsar.params=pargs) getStability(seBM.REGUA) [1] 0.04816097 sum(getRefit(seBM.REGUA))/2 [1] 139 seBM.REGUA$select$stars$summary [1] 0.00000000 0.01665210 0.02570711 0.03384109 0.04816097 0.05008990 0.05745964 [8] 0.06444095 0.07150077 0.07824774 0.08440371 0.09014204 0.10504912 0.11004641 [15] 0.11416388 0.11857213 0.11857213 0.12436764 0.12857738 0.13283227 0.13733352 [22] 0.14075534 0.14487036 0.14906452 0.15344541 0.15917835 0.16329116 0.16877920 [29] 0.18539474 0.18938216 0.19340811 0.20446358 0.20561258 0.21057085 0.21465802 [36] 0.21960684 0.22556086 0.23205300 0.24193642 0.25592418 0.26665109 0.26818876 [43] 0.27338983 0.27786483 0.28442788 0.28964135 0.29587817 0.30155074 0.30769002 [50] 0.31432718 0.31937134 0.32505360 0.32960724 0.33477408 0.34014753 0.34507489 [57] 0.34990777 0.35430449 0.36062032 0.36477914 seBM.REGUA$lambda [1] 1.000000000 0.914112171 0.835601062 0.763833101 0.698229135 0.638259750 [7] 0.583441006 0.533330525 0.487523924 0.445651553 0.407375509 0.372386911 [13] 0.340403408 0.311166898 0.284441449 0.260011391 0.237679577 0.217265794 [19] 0.198605307 0.181547528 0.165954805 0.151701307 0.138672011 0.126761774 [25] 0.115874480 0.105922273 0.096824839 0.088508763 0.080906938 0.073958017 [31] 0.067605923 0.061799397 0.056491581 0.051639642 0.047204425 0.043150140 [37] 0.039444068 0.036056303 0.032959505 0.030128685 0.027540997 0.025175561 [43] 0.023013287 0.021036725 0.019229927 0.017578310 0.016068547 0.014688455 [49] 0.013426895 0.012273688 0.011219528 0.010255907 0.009375049 0.008569847 [55] 0.007833801 0.007160973 0.006545933 0.005983717 0.005469788 0.005000000

Increasing nlambda values from 60 to 90 in increments of 10 all work well Increase nlambda to 100

seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=5e-3, nlambda=100, pulsar.params=pargs) Applying data transformations... Selecting model with pulsar using stars... Fitting final estimate with mb... done Warning messages: 1: 1 job had warning: "1 parallel function call did not deliver a result" 2: In pulsar(data = X, fun = match.fun(estFun), fargs = args, rep.num = 50, : Optimal lambda may be smaller than the supplied values

It looks like the issue is allocating a large enough vector in the underlying C++ code. huge attempts to save space by preallocating up to the theoretical maximum, but this undershoots for small sample size and large nlambda.

Can you please update the hugemb-patch branch again and try:

seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=5e-3, nlambda=100, pulsar.params=pargs, maxdf=ntaxa(BM.fragment.list[["BM.REGUA"]]))

thanks!

Hello! Sorry for the delay!

I deleted my previously posted response because it was inaccurate. After downloading the new patch and playing around with lambda.min.ratio, I was able to get a resolved network when nlambda = 100 using maxdf=ntaxa(BM.fragment.list[["BM.REGUA"]]). Thank you!

seBM.REGUA <- spiec.easi(BM.fragment.list[["BM.REGUA"]], method='mb', lambda.min.ratio=1e-2, nlambda=100, pulsar.params=pargs, maxdf=ntaxa(BM.fragment.list[["BM.REGUA"]])) Applying data transformations... Selecting model with pulsar using stars... Fitting final estimate with mb... done getStability(seBM.REGUA) [1] 0.04543382 sum(getRefit(seBM.REGUA))/2 [1] 127

zdk123 / SpiecEasi

Set a global default and user-overrides for the maxdf parameter in huge.mb #196