Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

differentialTest fails with boot = TRUE #128

Closed
amorris28 opened this issue Jan 18, 2022 · 9 comments
Closed

differentialTest fails with boot = TRUE #128

amorris28 opened this issue Jan 18, 2022 · 9 comments

Comments

@amorris28
Copy link

amorris28 commented Jan 18, 2022

I have a dataset that is 24 samples, 12 per treatment. I have aggregated my ASVs at the family level and there are approximately ~300 families. I have tried setting the minimum prevalence to 10% and 20%. I'm testing differential abundance between the two treatments with the variable treat which is a binary variable (p/n).

Here is the command I'm running:

da_analysis <- differentialTest(formula = ~ treat,
                                 phi.formula = ~ treat,
                                 formula_null = ~ 1,
                                 phi.formula_null = ~ treat,
                                 test = "lrt", boot = TRUE,
                                 data = ps_fam_pas5_prev,
                                 fdr_cutoff = 0.05)

When I run this, it appears to print each model to the Console as it runs using print.bbdml. Most of them run fine, except for the occasional model that says this:

Warning in print.bbdml(mod_null) :
  This model is based on a discriminant taxa.
You may see NAs in the model summary because Wald testing is invalid.
Likelihood ratio testing can be used, but valid standard errors cannot be calculated.

After it finishes I get this error message:

Error in differentialTest(formula = ~treat, phi.formula = ~treat, formula_null = ~1,  : 
  All models failed to converge! 

           If you are seeing this, it is likely that your model is overspecified. This occurs when your sample size is not large enough to estimate all the parameters of your model. This is most commonly due to categorical variables that include many categories. 

           Alternatively, double-check your values for the arguments `link`, `phi.link`, and `method` to makes sure that they follow the specified options. 

           To confirm you have fixed the issue, try running a model for a single taxon with bbdml.

I have tried running a number of taxa individually using bbdml as well and each one I've tried has worked:

corncob <- bbdml(formula = ASV2954 ~ treat,
                                 phi.formula = ~ 1,
                                 formula_null = ~ 1,
                                 phi.formula_null = ~ 1,
                                 test = "lrt", boot = TRUE,
                                 data = ps_fam_pas5_prev,
                                 fdr_cutoff = 0.05)

Based on issue #121 and #123 , I tried reinstalling the trust and optimr packages like this:

install.packages('trust', force=TRUE)
install.packages('optimr', force=TRUE)

Still getting the same error message.

I have not tried setting custom init values, but I don't know how to do that.

Any ideas how to get the parametric bootstrap to work?

Thanks!

sessioninfo():

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Manjaro Linux

Matrix products: default
BLAS:   /usr/lib/libblas.so.3.10.0
LAPACK: /usr/lib/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               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    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] microViz_0.7.10.9029 phyloseq_1.36.0      corncob_0.2.0       

loaded via a namespace (and not attached):
 [1] Biobase_2.52.0           tidyr_1.1.4              VGAM_1.1-5               jsonlite_1.7.2          
 [5] splines_4.1.2            foreach_1.5.1            speedyseq_0.5.3.9018     assertthat_0.2.1        
 [9] BiocManager_1.30.16      ROI.plugin.lpsolve_1.0-1 stats4_4.1.2             renv_0.14.0             
[13] GenomeInfoDbData_1.2.6   slam_0.1-50              numDeriv_2016.8-1.1      pillar_1.6.3            
[17] lattice_0.20-45          glue_1.4.2               detectseparation_0.2     XVector_0.32.0          
[21] colorspace_2.0-2         Matrix_1.3-4             plyr_1.8.6               pkgconfig_2.0.3         
[25] microbiome_1.14.0        zlibbioc_1.38.0          purrr_0.3.4              scales_1.1.1            
[29] Rtsne_0.15               tibble_3.1.5             mgcv_1.8-38              generics_0.1.0          
[33] IRanges_2.26.0           ggplot2_3.3.5            ellipsis_0.3.2           BiocGenerics_0.38.0     
[37] survival_3.2-13          magrittr_2.0.1           crayon_1.4.1             fansi_0.5.0             
[41] nlme_3.1-153             MASS_7.3-54              vegan_2.5-7              tools_4.1.2             
[45] registry_0.5-1           data.table_1.14.2        lifecycle_1.0.1          ROI_1.0-0               
[49] stringr_1.4.0            Rhdf5lib_1.14.2          S4Vectors_0.30.2         trust_0.1-8             
[53] munsell_0.5.0            cluster_2.1.2            Biostrings_2.60.2        ade4_1.7-18             
[57] compiler_4.1.2           GenomeInfoDb_1.28.4      rlang_0.4.11             rhdf5_2.36.0            
[61] grid_4.1.2               RCurl_1.98-1.5           iterators_1.0.13         rhdf5filters_1.4.0      
[65] biomformat_1.20.0        rstudioapi_0.13          igraph_1.2.6             bitops_1.0-7            
[69] gtable_0.3.0             codetools_0.2-18         multtest_2.48.0          DBI_1.1.1               
[73] reshape2_1.4.4           R6_2.5.1                 lpSolveAPI_5.5.2.0-17.7  knitr_1.36              
[77] dplyr_1.0.7              seqinr_4.2-8             utf8_1.2.2               permute_0.9-5           
[81] ape_5.5                  stringi_1.7.6            parallel_4.1.2           Rcpp_1.0.7              
[85] vctrs_0.3.8              tidyselect_1.1.1         xfun_0.26            
@bryandmartin
Copy link
Collaborator

Hi @amorris28 ,

Thanks for posting this issue! I am having a bit of trouble re-producing this issue. Can you please try a few things for me and let me know what the output is? That will help me diagnose the problem here.

  • Try running the example in the documentation, and let me know if that works or what the error is
  • Try running your example again, but with a very small value of B, e.g., B = 5 and let me know if the error persists

To set custom init values, you just need to provide a vector of values. For example, you could just do init = rep(1, p) where p is your number of variables.

Keep me updated! I'll work through this I just need a bit more information to figure it out.

@amorris28
Copy link
Author

Okay, so I ran the example for differentialTest. Worked with boot = F and boot = T, B = 5. Was taking a long time on default B.

I tried running my dataset again with B = 5 and it finished, but no models are significant.

I tried again with B = 5, inits = rep(1, 4) and I got this error:

Error in differentialTest(formula = ~treat, phi.formula = ~1, formula_null = ~1,  : 
  All models failed to converge! 

           If you are seeing this, it is likely that your model is overspecified. This occurs when your sample size is not large enough to estimate all the parameters of your model. This is most commonly due to categorical variables that include many categories. 

           Alternatively, double-check your values for the arguments `link`, `phi.link`, and `method` to makes sure that they follow the specified options. 

           To confirm you have fixed the issue, try running a model for a single taxon with bbdml.

@bryandmartin
Copy link
Collaborator

Hi @amorris28 ,

Is your data available to be shared? If so, would you be willing to send it to me? I'm having trouble diagnosing this and it would be helpful to debug looking at the same data you are encountering the problem with.

Bryan

@mmanus
Copy link

mmanus commented Feb 23, 2022

Hi @amorris28 and @bryandmartin - I've been following this conversation (and similar issues posted by others), as I'm also running into similar problems. Are there any updates?

So far I can't even get things to work when setting B=5 and trying to use custom init values. Regardless of the value I use for p in 'init=rep(1,p)' I receive the error:

"Error in differentialTest(formula = ~source4, phi.formula = ~1, formula_null = ~1, :
argument 10 matches multiple formal arguments"

Thanks in advance for any insight,
Melissa

@amorris28
Copy link
Author

Hey, @bryandmartin.

Just wanted to check on this issue. I emailed you my phyloseq object on 2/6, but I don't know if you receive it. Let me know if you didn't. Planning on submitting this manuscript soon with test = 'wald', boot = FALSE, but would happily rerun with bootstrapping, since that seems to be the recommended approach.

A

@adw96
Copy link
Collaborator

adw96 commented Apr 22, 2022

Hi @bryandmartin -- let me know if you would like help with corncob maintenance. Pauline, for example, will be my RA again soon (🥳) and I'm sure she'd be happy to help with some of these questions.

Hope you're well!

Cheers,

Amy

@emi-panzetta
Copy link

Hi
I am having a similar issue. The analysis is at the Family level, 3 samples belong to one group versus 3 samples in the other group.
OPTION 1: boot FALSE-> 9 significant taxa

da_analysis_donors_f_b<- differentialTest(formula= ~Prog,
phi.formula = ~Prog,
formula_null= ~1,
phi.formula_null= ~1,
test='Wald', boot=FALSE,
data=FMT_f_human.ps,
fdr_cutoff=0.05)

summary(da_analysis_donors_f_b)
Length Class Mode
p 56 -none- numeric
p_fdr 56 -none- numeric
significant_taxa 9 -none- character
significant_models 9 -none- list
all_models 56 -none- list
restrictions_DA 1 -none- character
restrictions_DV 1 -none- character
discriminant_taxa_DA 23 -none- character
discriminant_taxa_DV 23 -none- character
data 1 phyloseq S4
full_output 0 -none- NULL

OPTION 2: boot=TRUE, B=5, --> 0 significant taxa

da_analysis_donors_f_b<- differentialTest(formula= ~Prog,
phi.formula = ~Prog,
formula_null= ~1,
phi.formula_null= ~1,
test='Wald', boot=TRUE,B=5,
data=FMT_f_human.ps,
fdr_cutoff=0.05)

summary(da_analysis_donors_f_b)
Length Class Mode
p 56 -none- numeric
p_fdr 56 -none- numeric
significant_taxa 0 -none- character
significant_models 0 -none- list
all_models 56 -none- list
restrictions_DA 1 -none- character
restrictions_DV 1 -none- character
discriminant_taxa_DA 23 -none- character
discriminant_taxa_DV 23 -none- character
data 1 phyloseq S4
full_output 0 -none- NULL

OPTION 3: BOOT=TRUE, FAILS
da_analysis_donors_f_b<- differentialTest(formula= ~Prog,
phi.formula = ~Prog,
formula_null= ~1,
phi.formula_null= ~1,
test='Wald', boot=TRUE,
data=FMT_f_human.ps,
fdr_cutoff=0.05)

Error in if (is.nan(val) || any(phi <= sqrt(.Machine$double.eps)) || any(phi >= :
missing value where TRUE/FALSE needed

--
Even if this would run with B=1000, and there are no significant differences between the 2 groups, would you report what you find with and without bootstrapping or just the bootstrapping results?

devtools::session_info()
Session info ──────────────────────────────────────────────────────
setting value
version R version 4.2.2 (2022-10-31)
os macOS Ventura 13.1
system aarch64, darwin20
ui RStudio
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/New_York
date 2023-03-03
rstudio 2022.12.0+353 Elsbeth Geranium (desktop)
pandoc 2.19.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)

─ Packages ──────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
ade4 1.7-20 2022-11-01 [1] CRAN (R 4.2.0)
ape 5.6-2 2022-03-02 [1] CRAN (R 4.2.0)
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.2.0)
backports 1.4.1 2021-12-13 [1] CRAN (R 4.2.0)
Biobase 2.58.0 2022-11-07 [1] Bioconductor
BiocGenerics 0.44.0 2022-11-07 [1] Bioconductor
biomformat 1.26.0 2022-11-07 [1] Bioconductor
Biostrings 2.66.0 2022-11-07 [1] Bioconductor
bitops 1.0-7 2021-04-24 [1] CRAN (R 4.2.0)
boot 1.3-28.1 2022-11-22 [1] CRAN (R 4.2.0)
breakaway * 4.8.4 2023-02-09 [1] Github (adw96/breakaway@d81b179)
broom 1.0.2 2022-12-15 [1] CRAN (R 4.2.0)
bslib 0.4.2 2022-12-16 [1] CRAN (R 4.2.0)
cachem 1.0.6 2021-08-19 [1] CRAN (R 4.2.0)
callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.0)
car * 3.1-1 2022-10-19 [1] CRAN (R 4.2.0)
carData * 3.0-5 2022-01-06 [1] CRAN (R 4.2.0)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.2.0)
cli 3.6.0 2023-01-09 [1] CRAN (R 4.2.0)
cluster 2.1.4 2022-08-22 [1] CRAN (R 4.2.2)
codetools 0.2-18 2020-11-04 [1] CRAN (R 4.2.2)
colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.2.0)
corncob * 0.3.1 2022-12-05 [1] CRAN (R 4.2.0)
crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.0)
data.table 1.14.6 2022-11-16 [1] CRAN (R 4.2.0)
DBI 1.1.3 2022-06-18 [1] CRAN (R 4.2.0)
dbplyr 2.3.0 2023-01-16 [1] CRAN (R 4.2.0)
detectseparation 0.3 2022-08-26 [1] CRAN (R 4.2.0)
devtools * 2.4.5 2022-10-11 [1] CRAN (R 4.2.0)
digest 0.6.31 2022-12-11 [1] CRAN (R 4.2.0)
DivNet * 0.4.0 2023-02-09 [1] Github (adw96/DivNet@477270c)
doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.2.0)
dplyr * 1.1.0 2023-01-29 [1] CRAN (R 4.2.0)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.0)
evaluate 0.20 2023-01-17 [1] CRAN (R 4.2.0)
fansi 1.0.4 2023-01-22 [1] CRAN (R 4.2.0)
farver 2.1.1 2022-07-06 [1] CRAN (R 4.2.0)
fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.0)
forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.2.0)
foreach * 1.5.2 2022-02-02 [1] CRAN (R 4.2.0)
fs 1.6.1 2023-02-06 [1] CRAN (R 4.2.0)
gargle 1.3.0 2023-01-30 [1] CRAN (R 4.2.0)
generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.0)
GenomeInfoDb 1.34.7 2023-01-18 [1] Bioconductor
GenomeInfoDbData 1.2.9 2023-01-22 [1] Bioconductor
ggplot2 * 3.4.0 2022-11-04 [1] CRAN (R 4.2.0)
glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0)
googledrive 2.0.0 2021-07-08 [1] CRAN (R 4.2.0)
googlesheets4 1.0.1 2022-08-13 [1] CRAN (R 4.2.0)
gridExtra * 2.3 2017-09-09 [1] CRAN (R 4.2.0)
gtable 0.3.1 2022-09-01 [1] CRAN (R 4.2.0)
haven 2.5.1 2022-08-22 [1] CRAN (R 4.2.0)
highr 0.10 2022-12-22 [1] CRAN (R 4.2.0)
hms 1.1.2 2022-08-19 [1] CRAN (R 4.2.0)
htmltools 0.5.4 2022-12-07 [1] CRAN (R 4.2.0)
htmlwidgets 1.6.1 2023-01-07 [1] CRAN (R 4.2.0)
httpuv 1.6.8 2023-01-12 [1] CRAN (R 4.2.0)
httr 1.4.4 2022-08-17 [1] CRAN (R 4.2.0)
igraph 1.3.5 2022-09-22 [1] CRAN (R 4.2.0)
IRanges 2.32.0 2022-11-07 [1] Bioconductor
iterators 1.0.14 2022-02-05 [1] CRAN (R 4.2.0)
jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.0)
jsonlite 1.8.4 2022-12-06 [1] CRAN (R 4.2.0)
knitr 1.41 2022-11-18 [1] CRAN (R 4.2.0)
labeling 0.4.2 2020-10-20 [1] CRAN (R 4.2.0)
later 1.3.0 2021-08-18 [1] CRAN (R 4.2.0)
lattice 0.20-45 2021-09-22 [1] CRAN (R 4.2.2)
lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.0)
lme4 * 1.1-31 2022-11-01 [1] CRAN (R 4.2.0)
lpSolveAPI 5.5.2.0-17.9 2022-10-20 [1] CRAN (R 4.2.0)
lubridate 1.9.1 2023-01-24 [1] CRAN (R 4.2.0)
magrittr * 2.0.3 2022-03-30 [1] CRAN (R 4.2.0)
MASS 7.3-58.1 2022-08-03 [1] CRAN (R 4.2.2)
Matrix * 1.5-3 2022-11-11 [1] CRAN (R 4.2.0)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.0)
mgcv 1.8-41 2022-10-21 [1] CRAN (R 4.2.2)
mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.0)
minqa 1.2.5 2022-10-19 [1] CRAN (R 4.2.0)
modelr 0.1.10 2022-11-11 [1] CRAN (R 4.2.0)
multtest 2.54.0 2022-11-07 [1] Bioconductor
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.0)
mvnfast 0.2.7 2021-05-20 [1] CRAN (R 4.2.0)
nlme 3.1-161 2022-12-15 [1] CRAN (R 4.2.0)
nloptr 2.0.3 2022-05-26 [1] CRAN (R 4.2.0)
numDeriv 2016.8-1.1 2019-06-06 [1] CRAN (R 4.2.0)
permute 0.9-7 2022-01-27 [1] CRAN (R 4.2.0)
phyloseq * 1.42.0 2022-11-01 [1] Bioconductor
pillar 1.8.1 2022-08-19 [1] CRAN (R 4.2.0)
pkgbuild 1.4.0 2022-11-27 [1] CRAN (R 4.2.0)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0)
pkgload 1.3.2 2022-11-16 [1] CRAN (R 4.2.0)
plyr 1.8.8 2022-11-11 [1] CRAN (R 4.2.0)
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.0)
processx 3.8.0 2022-10-26 [1] CRAN (R 4.2.0)
profvis 0.3.7 2020-11-02 [1] CRAN (R 4.2.0)
promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.0)
ps 1.7.2 2022-10-26 [1] CRAN (R 4.2.0)
purrr * 1.0.1 2023-01-10 [1] CRAN (R 4.2.0)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0)
Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.2.0)
RCurl 1.98-1.9 2022-10-03 [1] CRAN (R 4.2.0)
readr * 2.1.3 2022-10-01 [1] CRAN (R 4.2.0)
readxl * 1.4.2 2023-02-09 [1] CRAN (R 4.2.0)
registry 0.5-1 2019-03-05 [1] CRAN (R 4.2.0)
remotes 2.4.2 2021-11-30 [1] CRAN (R 4.2.0)
reprex 2.0.2 2022-08-17 [1] CRAN (R 4.2.0)
reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.2.0)
rhdf5 2.42.0 2022-11-07 [1] Bioconductor
rhdf5filters 1.10.0 2022-11-07 [1] Bioconductor
Rhdf5lib 1.20.0 2022-11-07 [1] Bioconductor
rlang 1.0.6 2022-09-24 [1] CRAN (R 4.2.0)
rmarkdown 2.20 2023-01-19 [1] CRAN (R 4.2.0)
ROI 1.0-0 2020-08-31 [1] CRAN (R 4.2.0)
ROI.plugin.lpsolve 1.0-1 2021-06-15 [1] CRAN (R 4.2.0)
rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.0)
rvest 1.0.3 2022-08-19 [1] CRAN (R 4.2.0)
S4Vectors 0.36.1 2022-12-07 [1] Bioconductor
sass 0.4.4 2022-11-24 [1] CRAN (R 4.2.0)
scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.0)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0)
shiny 1.7.4 2022-12-15 [1] CRAN (R 4.2.0)
slam 0.1-50 2022-01-08 [1] CRAN (R 4.2.0)
speedyseq * 0.5.3.9018 2023-02-09 [1] Github (mikemc/speedyseq@ceb941f)
stringi 1.7.12 2023-01-11 [1] CRAN (R 4.2.0)
stringr * 1.5.0 2022-12-02 [1] CRAN (R 4.2.0)
survival 3.5-0 2023-01-09 [1] CRAN (R 4.2.0)
tibble * 3.1.8 2022-07-22 [1] CRAN (R 4.2.0)
tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.2.0)
tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.0)
tidyverse * 1.3.2 2022-07-18 [1] CRAN (R 4.2.0)
timechange 0.2.0 2023-01-11 [1] CRAN (R 4.2.0)
trust 0.1-8 2020-01-10 [1] CRAN (R 4.2.0)
tzdb 0.3.0 2022-03-28 [1] CRAN (R 4.2.0)
urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.0)
usethis * 2.1.6 2022-05-25 [1] CRAN (R 4.2.0)
utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.0)
vctrs 0.5.2 2023-01-23 [1] CRAN (R 4.2.0)
vegan 2.6-4 2022-10-11 [1] CRAN (R 4.2.0)
VGAM 1.1-7 2022-07-06 [1] CRAN (R 4.2.0)
withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0)
writexl * 1.4.2 2023-01-06 [1] CRAN (R 4.2.0)
xfun 0.36 2022-12-21 [1] CRAN (R 4.2.0)
xml2 1.3.3 2021-11-30 [1] CRAN (R 4.2.0)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.0)
XVector 0.38.0 2022-11-07 [1] Bioconductor
yaml 2.3.6 2022-10-18 [1] CRAN (R 4.2.0)
zlibbioc 1.44.0 2022-11-07 [1] Bioconductor

[1] /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library

─────────────────────────────────────────────────────────────────────

@adw96
Copy link
Collaborator

adw96 commented Oct 25, 2023

Hi @amorris28 -- thanks for your patience! Bryan has moved on to a new position and I'm helping to maintain the package going forward.

Are you able to send me the data that you sent Bryan? We'll do our best to give you an answer soon. You can find my email in the header at http://statisticaldiversitylab.com/.

I'm very sorry for the delay, and I hope your manuscript submission has gone well.

Amy

@adw96
Copy link
Collaborator

adw96 commented Oct 27, 2023

Thanks so much for sharing your data with me via email, @amorris28. The default B = 1000 took a long time so I ran your code with B=4.

I could reproduce an error, but can diagnose it as a user-input error. Here's your original code:

da_analysis_mini_break <-  differentialTest(formula = ~ treat,
                                            phi.formula = ~ treat,
                                            formula_null = ~ 1,
                                            phi.formula_null = ~ treat,
                                            test = "lrt", boot = TRUE,
                                            data = prune_taxa(taxa_sums(ps) > 5000, ps), # to speed up
                                            fdr_cutoff = 0.05,
                                            B=5) # also to speed up

and here is code that runs:

da_analysis_mini_break_LRT <-  differentialTest(formula = ~ treat,
                                            phi.formula = ~ treat,
                                            formula_null = ~ 1,
                                            phi.formula_null = ~ treat,
                                            test = "LRT", boot = TRUE,
                                            data = prune_taxa(taxa_sums(ps) > 5000, ps), # to speed up
                                            fdr_cutoff = 0.05,
                                            B=5) # also to speed up

(Notice the capitalization of "LRT".)

This issue recently came to our attention in #148 and so actually the amazing @svteichman actually wrote a helpful pull request #151 which now gives more information ("Alternatively, double-check your values for the arguments link, phi.link, and method to makes sure that they follow the specified options."). I apologize that this feature didn't occur to us earlier -- it probably would have helped you in Jan 2022!

I have further refined the error message in #155 -- actually I did this independently of your issue as I'm extending the testing options 😸 -- "If no errors were thrown by bbdml, and the issue is instead in differentialTest, no output is shown." I hope future users will find this helpful.

I believe this resolves the original issue so I'm going to close -- please feel free to reopen if your concern is not fully addressed.

To others on this thread also finding this issue with differentialTest() please

  1. update to the latest dev version (with install_github or similar)
  2. double check the spelling of your input argument options
  3. (if problems persist) open a new issue detailing your problem

Thanks all for your patience and feedback -- it helps us grow better 🌽 !

@adw96 adw96 closed this as completed Oct 27, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants