statdivlab / corncob Goto Github PK

Count Regression for Correlated Observations with the Beta-binomial

R 99.48% Dockerfile 0.52%

corncob's Introduction

corncob

Count Regression for Correlated Observations with the Beta-binomial

corncob is an R package for modeling relative abundance and testing hypotheses about the effect of covariates on relative abundance. The corncob methodology was specifically developed for modelling microbial abundances based on high throughput sequencing data, such as 16S or whole-genome sequencing.

Installation

To install the corncob package, use the code below to download the development version from Github.

# install.packages("remotes")
remotes::install_github("statdivlab/corncob")
library(corncob)

Docker

Instead of installing corncob to your local system, you can use corncob via the pre-compiled Docker image: quay.io/fhcrc-microbiome/corncob.

Use

The vignettes demonstrate example usage of all main functions (they go over the same analysis, one with phyloseq objects and one without phyloseq objects). Please file an issue if you have a request for a tutorial that is not currently included. You can see the vignette by using the following code:

library(corncob)
# Use this to view the vignette in the corncob HTML help
help(package = "corncob", help_type = "html")
# Use this to view the vignette as an isolated HTML file
utils::browseVignettes(package = "corncob")

Note that R does not allow variable names to start with numbers. Sometimes, when going directly from QIIME2 to phyloseq objects, taxa names will be a large string starting with numbers. To clean these taxa names for use with corncob, use clean_taxa_names(my_phyloseq_object), see ?clean_taxa_names for details.

Documentation

We additionally have a pkgdown website that contains pre-built versions of our function documentation, vignette, and a version of the vignette without phyloseq.

Citation

If you use corncob for your analysis, please cite our manuscript:

Bryan D. Martin, Daniela Witten, and Amy D. Willis. (2020). Modeling microbial abundances and dysbiosis with beta-binomial regression. Annals of Applied Statistics, Volume 14, Number 1, pages 94-115.

An open-access preprint is available on arXiv here.

Bug Reports / Change Requests

If you encounter a bug or would like make a change request, please file it as an issue here.

corncob's People

Contributors

Stargazers

Watchers

corncob's Issues

any(phi >= : missing value where TRUE/FALSE needed

Hi @bryandmartin,

I have 4 animals that were sampled different ways on 3 different days. I want to know how the sample types differ. My data for one animal looks like this:

     AnimalID       Sample_Type    Day
1  Animal_2549             Feces  Day_7
2  Animal_2549             Feces  Day_9
3  Animal_2549             Feces Day_11
4  Animal_2549       Grab Sample  Day_7
5  Animal_2549       Grab Sample  Day_9
6  Animal_2549       Grab Sample Day_11
7  Animal_2549   Liquid Strained  Day_7
8  Animal_2549   Liquid Strained  Day_9
9  Animal_2549   Liquid Strained Day_11
10 Animal_2549 Liquid Unstrained Day_11
11 Animal_2549 Liquid Unstrained  Day_9
12 Animal_2549             Solid  Day_7
13 Animal_2549             Solid  Day_9
14 Animal_2549             Solid Day_11
15 Animal_2549      Stomach Tube  Day_9
16 Animal_2549      Stomach Tube Day_11
17 Animal_2549      Stomach Tube  Day_7

I am running the following command to look for deferentially abundant taxa in sample types while controlling for the effects of day and animal sampled.

ps51719 <- readRDS("ps_51719.rds")
set.seed(1)
fullAnalysis_060319 <- differentialTest(formula = ~ Sample_Type + AnimalID + Day, phi.formula = ~ 1,formula_null = ~ AnimalID + Day, phi.formula_null = ~ 1, test="Wald",boot=TRUE, data = ps51719, fdr_cutoff = 0.05)

This produces the following error. Although, I still got a list of significant_taxa so it seems as though it didn't stop differentialTest from running.

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

Have you seen this error? Do I have to have something in the phi.formula is that why its giving the error? I believe the error started after I set boot=TRUE, but have't been able to track down the line that is causing the problem. Thanks for your help!

I am using the following versions:

sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 LTS

Matrix products: default
BLAS:   /afs/genomecenter.ucdavis.edu/software/R/3.6.0/lssc0-linux/lib/R/lib/libRblas.so
LAPACK: /afs/genomecenter.ucdavis.edu/software/R/3.6.0/lssc0-linux/lib/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US       LC_NUMERIC=C         LC_TIME=en_US
 [4] LC_COLLATE=en_US     LC_MONETARY=en_US    LC_MESSAGES=en_US
 [7] LC_PAPER=en_US       LC_NAME=C            LC_ADDRESS=C
[10] LC_TELEPHONE=C       LC_MEASUREMENT=en_US LC_IDENTIFICATION=C

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

other attached packages:
[1] corncob_0.1.0   phyloseq_1.24.2

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5      reshape2_1.4.3        purrr_0.3.2
 [4] splines_3.6.0         lattice_0.20-38       rhdf5_2.24.0
 [7] colorspace_1.3-2      stats4_3.6.0          enrichwith_0.2
[10] mgcv_1.8-28           survival_2.44-1.1     rlang_0.3.4
[13] pillar_1.4.0          glue_1.3.1            BiocGenerics_0.26.0
[16] bindrcpp_0.2.2        trust_0.1-7           foreach_1.4.4
[19] plyr_1.8.4            bindr_0.1.1           stringr_1.4.0
[22] zlibbioc_1.26.0       Biostrings_2.48.0     munsell_0.5.0
[25] gtable_0.3.0          VGAM_1.1-1            lpSolveAPI_5.5.2.0-17
[28] codetools_0.2-16      Biobase_2.40.0        permute_0.9-4
[31] IRanges_2.14.12       biomformat_1.8.0      parallel_3.6.0
[34] Rcpp_1.0.1            scales_1.0.0          vegan_2.5-2
[37] S4Vectors_0.18.3      jsonlite_1.6          XVector_0.20.0
[40] ggplot2_3.1.1         stringi_1.4.3         dplyr_0.7.6
[43] numDeriv_2016.8-1     grid_3.6.0            ade4_1.7-13
[46] tools_3.6.0           magrittr_1.5          lazyeval_0.2.2
[49] tibble_2.1.1          cluster_2.0.8         crayon_1.3.4
[52] ape_5.1               pkgconfig_2.0.2       MASS_7.3-51.4
[55] Matrix_1.2-17         data.table_1.11.6     assertthat_0.2.1
[58] brglm2_0.5.1          iterators_1.0.10      Rhdf5lib_1.2.1
[61] R6_2.4.0              multtest_2.36.0       nnet_7.3-12
[64] igraph_1.2.2          nlme_3.1-127          compiler_3.6.0

(Possible) Vignette correction?

Under the heading "Analysis for Multiple Taxa", the second, third, and fourth paragraphs and the following code reads (emphasis mine):

Next, we use the differentialTest command. We specify the covariates of our model using formula and phi.formula as before, except we no longer include the response term because we are testing multiple taxa. We also specify which covariates we want to test for by removing them in the formula_null and phi.formula_null arguments.

The difference between the formulas and the null version of the formulas will be the variables that are tested. In this case, as when we examined the single taxon, we will be testing the coefficients of DayAmdmt for both the expected relative abundance and the overdispersion.

We set fdr_cutoff to be our controlled false discovery rate.

set.seed(1)
da_analysis <- differentialTest(formula = ~ DayAmdmt,
                                 phi.formula = ~ DayAmdmt,
                                 formula_null = ~ 1,
                                 phi.formula_null = ~ DayAmdmt,
                                 test = "Wald", boot = FALSE,
                                 data = soil,
                                 fdr_cutoff = 0.05)

If I'm reading this right, the code is testing for differential abundance across DayAmdmt, controlling for the effect of DayAmdmt on dispersion. If that's the case, I suggest correcting/updating the vignette.

Visualizing corncob output

Hi there! Our lab is really interested in using GraPhlAn to visualize differential abundances of our microbiome dataset. I'm having a hard time parsing the output from corncob into GraPhlAn. Would you be willing to include a function that parses the corncob output for GraPhlAn? I know this might be a lot to ask, but I think it would be a very useful pipeline for investigators doing microbiome research.

For example, Figure 1 of this paper is exactly what I'm hoping to visualize. Using a cladogram, the authors can compare increases and decreases in taxa (and different levels as well) from 2 experimental groups. They use LEfSe to create the input to GraPhIAn. I would really like to use the output from corncob, so your help is much much appreciated!

Add to unit testing

phyloseq object in bbdml
Group argument in plot, RA
Group argument in plot, AA
print summary, with covariates
print summary, without covariates
qbetabinom, above 0
qbetabinom, below 0
waldtest

Perfect separation becomes an issue with logit link

See perfectSep.Rda in singularHessian. Would happen if only one count is observed.

Use older version of .rda so it doesn't require newer version of R

See r-lib/devtools#1912

Tree warning?

Hello!

I'm currently trying corncob out on my data and noticed that if the physeq object contains a tree it has the following warning:

In prune_taxa(taxa, phy_tree(x)) :
prune_taxa attempted to reduce tree to 1 or fewer tips.
tree replaced with NULL.

Is this intentional or?

[question] how to get statistics for individual coefficients or contrasts

Hi,

I was browsing through the docs but could not find any information how I would go about getting the coefficients for several terms in bbdml or differential tests along with the associated test statistics and p value.

For instance if my design would be abundance ~ var1 + var2 how can I get the fitted coefficients for each variable or a particular contrast (for instance coef1 + coef2) similar to DESeq2 or voom?

Thanks!

genInits for bstar

With large covariate values (such as BMI in tutorial), can generate initializations that don't work as well with trust method.

corncob and phyloseq objects: invalid class for otu_table

Hello,

I have run corncob in both a unsubsetted and subsetted phyloseq object
ps_WWTP
ps_subset2

ps_WWTP_RA = transform_sample_counts(ps_WWTP, function(x) x / sum(x))
ps_subset2_RA = transform_sample_counts(ps_subset2, function(x) x / sum(x))

data_phylum <- ps_subset2_RA %>% phyloseq::tax_glom("Phylum")

corncob <- bbdml(formula = OTU1 ~ 1, phi.formula = ~ 1, data = data_phylum)

I received the following warning:
Error in validObject(.Object) : invalid class “otu_table” object: OTU abundance data must have non-zero dimensions.

I thought this was due to the transformations as the class of the otu_table when from integer to double, even though R is supposed to handle them similarly.

This also may be a phyloseq issue, but any info will help. Thanks!

Upgrading otu_to_taxonomy for improved plotting

Hi @bryandmartin,

Thanks for making the new plotting function for differentialTest() it look great! Unfortunately since my ASVs have been assigned down to species and I didn't use taxa_glom the plots get taken up by a big long name. What if the otu_to_taxonomy function was changed to include an argument for what taxonomy level you want returned?

I came up with this:

otu_to_taxonomy <- function(OTU, data, level) {
  if(missing(level)) {
    return(apply(tax_table(data)[OTU,], 1, function(x) {paste(stats::na.omit(x), collapse = '_')}))
  } else {
  return(apply(tax_table(data)[OTU,level], 1, function(x) {paste(stats::na.omit(x), collapse = '_')}))
  }
}

So the default would still have the same output that it currently does, but you could call specific taxonomy like this:

>otu_to_taxonomy(OTU = da_analysis_sub$significant_taxa, data = ps_sub, level= c("Species","Genus"))
                                       ASV_682                                        ASV_429 
"Lachnospiraceae_Lachnospiraceae_AC2044_group"                "Lachnospiraceae_Acetatifactor" 
                                       ASV_460                                        ASV_164 
                             "Lachnospiraceae"               "Lachnospiraceae_Acetitomaculum" 
                                       ASV_472                                       ASV_1021 
                             "Lachnospiraceae"                              "Lachnospiraceae"

To add this to the plotting functions I added level as an argument to lines 27 and 31 of the plot_differentialTest.R file. Like this:

plot.differentialTest <- function(x, level, ...) {
  signif_taxa <- x$significant_taxa
  if ("phyloseq" %in% class(x$data)) {
    if (!is.null(x$data@tax_table)) {
      signif_taxa <- otu_to_taxonomy(signif_taxa, x$data,level)
    }
  }

This really improved the legibility of my graph. Let me know what you think!

Error in installation (* creating vignettes ... ERROR)

Hello,
Probably a local issue, but figured it would be useful for others, in case they see the same.
In trying to install CornCob and encountered an error in creating vignettes for whatever reason.

* preparing ‘corncob’:
* checking DESCRIPTION meta-information ... OK
* installing the package to build vignettes
* creating vignettes ... ERROR
Error: processing vignette 'corncob-intro.Rmd' failed with diagnostics:
Failed to compile corncob-intro.tex.
Execution halted
Installation failed: Command failed (1)

Installing TinyTex solved it.

tinytex::install_tinytex()

application of corncob to longitudinal data

Hi all,

thanks for corncob. I've enjoyed using this package so far.

I have a conceptual question - is this approach applicable to longitudinal data? If, for example, I wanted to look for differential variability in my replicate time series, is it appropriate to perform the test at every individual time point to test for this?

In your preprint you state accommodating longitudinal data could involve adding random-effects terms to the link functions for the beta-binomial model? Is this something you are planning to incorporate into corncob?

thanks again

-shane hogle

Adding a contrasts argument

This is a great tool and I'm looking forward to using it.

Would there be a way to incorporate a contrast argument into differentialTest? This would be useful for testing the effects of levels within a factor, but would also facilitate repeated measures analyses. I was able to apply the "trick" for 2-way repeated measures data (group x condition) where group-specific condition effects are tested with subjects nested within groups, as described in the DESeq2 User Guide. However, contrast tests are necessary to determine whether the condition effect differs between groups.

Perhaps if there was a way to extract the variance-covariance matrix for the regression coefficients in each model, then a Wald-test between different coefficients could be used, as employed by DESeq2.

Issue including co-variate in phi formula but not formula

Hi Bryan,

I've opened up a separate issue as I have run into an error when supplying the same co-variate to formula vs phi.formula. When I run the covariate in formula everything goes as expected, however, when I put it into the phi.formula I get the follow:


height <- corncob::differentialTest(formula= ~ 1,
                         phi.formula = ~ 1,
                         phi.formula_null = ~ PM_STANDING_HEIGHT_AVG,
                         formula_null = ~ 1,
                         test="Wald", data=phylo_test,
                         fdr_cutoff = 0.1,
                         boot=F)


Error: $ operator is invalid for atomic vectors

Any idea what may be happening here?

Thanks,
Jacob Nearing

Hessian is negative of numerical hessian

Add check for -Inf inits to re-initialize

Error in separator(x = x, y = y, linear_program = control$linear_program ; unexpected result from lpSolveAPI for primal test

Hi Bryan! Thanks again for developing and maintaining such a fantastic package.

My issue is "Error in separator(x = x, y = y, linear_program = control$linear_program, :
unexpected result from lpSolveAPI for primal test". Reproducible example below.

sep_issue1 <- structure(list(sample = c("S112", "S118", "S121", "S126", "S127", 
                          "S131", "S137", "S141", "S143", "S144", "S146", "S153", "S155", 
                          "S158", "S186", "S187", "S193", "S194", "S195", "S198", "S200", 
                          "S205", "S224", "S225", "S227", "S229", "S231", "S235", "S239", 
                          "S243", "S244", "S248", "S249", "S258", "S263", "S264", "S268", 
                          "S269", "S278", "S286", "S288", "S290", "S297", "S299", "S304", 
                          "S305", "S312", "S317", "S327", "S330", "S332", "S334", "S342", 
                          "S344", "S347", "S349", "S350", "S353", "S354", "S356", "S365", 
                          "S367", "S370", "S372", "S378", "S380", "S381", "S384", "S389", 
                          "S393", "S395", "S397", "S398", "S420", "S421", "S422", "S424", 
                          "S432", "S434", "S436", "S439", "S441", "S446", "S451", "S456", 
                          "S457", "S458", "S462", "S463", "S468", "S470", "S472", "S475", 
                          "S476", "S478", "S479", "S483", "S484", "S485", "S487", "S488", 
                          "S495", "S497", "S51", "S517", "S524", "S526", "S529", "S53", 
                          "S530", "S536", "S54", "S543", "S547", "S550", "S551", "S552", 
                          "S558", "S561", "S563", "S577", "S584", "S589", "S59", "S591", 
                          "S598", "S60", "S600", "S602", "S604", "S614", "S617", "S618", 
                          "S628", "S638", "S640", "S643", "S646", "S655", "S77", "S80", 
                          "S88", "S92"), W = c(4, 8, 4, 14780, 4199, 29841, 9485, 333, 
                                               6078, 38979, 0, 2, 1827, 1212, 6425, 5779, 20316, 2, 1, 9077, 
                                               26661, 8, 11266, 3, 7237, 6, 35070, 4, 12584, 791, 0, 7, 14, 
                                               4, 71179, 9763, 13, 7, 10011, 746, 3, 8, 233, 3, 2, 370, 5, 16, 
                                               2789, 8307, 25184, 1, 647, 3, 6, 6, 41558, 5, 8, 3, 26490, 5, 
                                               3, 7, 6764, 6618, 18948, 1006, 2, 0, 3, 16215, 0, 4, 46, 5503, 
                                               4, 1, 6, 37954, 0, 5, 3244, 0, 9, 13, 77997, 8733, 3817, 2846, 
                                               417, 19, 16, 21, 4, 5272, 10, 3, 561, 47336, 28303, 5, 19980, 
                                               34, 1480, 21367, 5072, 1, 10203, 4436, 13507, 35445, 697, 2, 
                                               21, 8, 9829, 10141, 654, 0, 6917, 0, 0, 20323, 269, 7525, 5505, 
                                               3, 2, 18, 0, 1530, 0, 4003, 1, 0, 14669, 2, 48, 9372, 7286, 29666, 
                                               8), Participant = c("S112", "S118", "S121", "S126", "S127", "S131", 
                                                                   "S137", "S141", "S143", "S144", "S146", "S153", "S155", "S158", 
                                                                   "S186", "S187", "S193", "S194", "S195", "S198", "S200", "S205", 
                                                                   "S224", "S225", "S227", "S229", "S231", "S235", "S239", "S243", 
                                                                   "S244", "S248", "S249", "S258", "S263", "S264", "S268", "S269", 
                                                                   "S278", "S286", "S288", "S290", "S297", "S299", "S304", "S305", 
                                                                   "S312", "S317", "S327", "S330", "S332", "S334", "S342", "S344", 
                                                                   "S347", "S349", "S350", "S353", "S354", "S356", "S365", "S367", 
                                                                   "S370", "S372", "S378", "S380", "S381", "S384", "S389", "S393", 
                                                                   "S395", "S397", "S398", "S420", "S421", "S422", "S424", "S432", 
                                                                   "S434", "S436", "S439", "S441", "S446", "S451", "S456", "S457", 
                                                                   "S458", "S462", "S463", "S468", "S470", "S472", "S475", "S476", 
                                                                   "S478", "S479", "S483", "S484", "S485", "S487", "S488", "S495", 
                                                                   "S497", "S51", "S517", "S524", "S526", "S529", "S53", "S530", 
                                                                   "S536", "S54", "S543", "S547", "S550", "S551", "S552", "S558", 
                                                                   "S561", "S563", "S577", "S584", "S589", "S59", "S591", "S598", 
                                                                   "S60", "S600", "S602", "S604", "S614", "S617", "S618", "S628", 
                                                                   "S638", "S640", "S643", "S646", "S655", "S77", "S80", "S88", 
                                                                   "S92"), X1 = c(1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 1, 0, 
                                                                                  0, 0, 0, 0, 1, 2, 0, 0, 1, 1, 0, 2, 1, 1, 2, 2, 1, 2, 2, 0, 0, 
                                                                                  0, 1, 2, 0, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 2, 
                                                                                  1, 1, 0, 1, 2, 1, 2, 2, 0, 1, 2, 1, 1, 1, 2, 1, 2, 2, 0, 1, 2, 
                                                                                  2, 2, 1, 2, 0, 0, 1, 2, 2, 0, 2, 2, 1, 2, 2, 2, 0, 2, 2, 1, 1, 
                                                                                  1, 2, 2, 1, 2, 2, 2, 2, 0, 2, 2, 1, 2, 1, 2, 2, 2, 1, 2, 2, 2, 
                                                                                  0, 2, 0, 2, 2, 1, 1, 2, 2, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 
                                                                                  0), X2 = c(69, 70, 69, 69, 70, 69, 70, 69, 69, 69, 69, 70, 70, 
                                                                                             68, 70, 70, 70, 70, 71, 69, 69, 70, 69, 70, 69, 69, 70, 69, 71, 
                                                                                             70, 69, 69, 70, 69, 69, 69, 69, 70, 70, 69, 70, 70, 69, 70, 70, 
                                                                                             70, 70, 70, 70, 70, 71, 70, 70, 70, 71, 70, 70, 69, 70, 70, 71, 
                                                                                             70, 70, 71, 71, 71, 70, 71, 69, 70, 70, 71, 70, 69, 70, 71, 70, 
                                                                                             71, 69, 71, 70, 70, 68, 71, 70, 71, 70, 70, 71, 70, 70, 70, 70, 
                                                                                             71, 70, 71, 71, 71, 69, 71, 70, 69, 70, 69, 71, 70, 71, 69, 70, 
                                                                                             70, 70, 69, 71, 70, 70, 71, 70, 71, 71, 69, 70, 70, 71, 69, 71, 
                                                                                             71, 69, 71, 71, 71, 71, 71, 71, 71, 71, 71, 70, 71, 71, 69, 70, 
                                                                                             69, 69), X3 = c(38.27625652, 20.18474465, 42.16368883, 40.66852255, 
                                                                                                             20.93232779, 15.69924584, 33.04317458, 38.27625652, 33.94027434, 
                                                                                                             29.90332541, 20.4837779, 37.82770664, 35.13640735, 28.55767576, 
                                                                                                             33.1926912, 26.0158931, 34.38882422, 26.61395961, 19.73619477, 
                                                                                                             33.34220783, 19.8857113999999, 23.32459382, 19.13812826, 31.39849168, 
                                                                                                             33.64124108, 23.02556056, 27.21202612, 17.94199524, 18.09151187, 
                                                                                                             35.13640735, 15.84876247, 28.10912588, 29.90332541, 17.04489548, 
                                                                                                             27.660576, 31.69752493, 29.00622565, 37.22964013, 27.36154275, 
                                                                                                             28.55767576, 20.63329453, 21.97894417, 30.20235866, 21.23136104, 
                                                                                                             23.17507719, 26.31492636, 34.98689073, 36.78109025, 34.98689073, 
                                                                                                             32.59462469, 42.61223871, 20.18474465, 31.99655819, 17.79247862, 
                                                                                                             31.54800831, 27.95960926, 32.29559144, 21.53039429, 21.97894417, 
                                                                                                             29.00622565, 35.73447386, 12.11084679, 15.10117933, 28.85670902, 
                                                                                                             26.61395961, 15.84876247, 18.54006175, 14.20407957, 23.17507719, 
                                                                                                             15.54972921, 22.42749406, 16.44682897, 29.90332541, 26.16540973, 
                                                                                                             21.23136104, 23.32459382, 24.52072683, 19.43716152, 23.32459382, 
                                                                                                             43.06078859, 27.660576, 21.38087767, 21.53039429, 26.76347624, 
                                                                                                             37.82770664, 22.27797743, 21.23136104, 18.2410285, 20.63329453, 
                                                                                                             28.85670902, 30.05284204, 25.71685985, 25.4178266, 19.8857113999999, 
                                                                                                             24.96927672, 15.84876247, 31.39849168, 19.43716152, 13.15746318, 
                                                                                                             24.81976009, 24.67024346, 28.70719239, 30.05284204, 27.06250949, 
                                                                                                             23.02556056, 26.46444299, 33.34220783, 19.8857113999999, 24.37121021, 
                                                                                                             24.96927672, 15.84876247, 37.82770664, 28.10912588, 24.22169358, 
                                                                                                             13.75552969, 49.34048692, 32.59462469, 38.27625652, 4.485498811, 
                                                                                                             20.4837779, 24.67024346, 30.9499418, 18.83909501, 49.34048692, 
                                                                                                             30.05284204, 25.71685985, 23.62362707, 23.32459382, 20.03522802, 
                                                                                                             18.98861163, 39.17335628, 22.87604394, 25.56734322, 34.98689073, 
                                                                                                             41.71513894, 41.41610569, 30.80042517, 12.85842993, 26.76347624, 
                                                                                                             26.0158931, 24.07217695, 32.89365795, 28.40815914), X4 = c(95.70894677, 
                                                                                                                                                                        47.85447339, 43.147476, 77.66545681, 90.21744983, 148.2704175, 
                                                                                                                                                                        72.17395986, 52.56147077, 106.6919407, 79.23445594, 91.78644896, 
                                                                                                                                                                        32.94898168, 65.11346379, 117.674934599999, 47.06997382, 92.57094852, 
                                                                                                                                                                        131.7959267, 70.60496073, 103.5539424, 80.0189555, 163.1759093, 
                                                                                                                                                                        117.674934599999, 75.31195812, 54.1304699, 65.89796335, 111.398937999999, 
                                                                                                                                                                        94.13994765, 313.7998255, 152.977414899999, 90.21744983, 285.5578412, 
                                                                                                                                                                        98.84694503, 54.1304699, 157.6844123, 76.88095724, 65.11346379, 
                                                                                                                                                                        72.95845943, 72.95845943, 55.69946902, 73.74295899, 89.43295026, 
                                                                                                                                                                        183.5728979, 66.68246292, 135.7184245, 70.60496073, 106.6919407, 
                                                                                                                                                                        72.95845943, 42.36297644, 80.80345506, 85.51045245, 74.52745855, 
                                                                                                                                                                        105.907441099999, 61.19096597, 120.8129328, 74.52745855, 53.34597033, 
                                                                                                                                                                        80.0189555, 80.0189555, 127.8734289, 131.0114271, 71.3894603, 
                                                                                                                                                                        498.157223, 97.2779459, 121.5974324, 79.23445594, 185.141897, 
                                                                                                                                                                        98.06244547, 168.6674062, 69.03596161, 271.436849, 97.2779459, 
                                                                                                                                                                        78.44995637, 77.66545681, 79.23445594, 73.74295899, 83.15695375, 
                                                                                                                                                                        128.6579285, 203.185387, 123.9509311, 67.46696248, 50.99247164, 
                                                                                                                                                                        99.63144459, 112.183437599999, 38.44047862, 49.42347251, 98.84694503, 
                                                                                                                                                                        225.9358744, 89.43295026, 123.166431499999, 38.44047862, 43.147476, 
                                                                                                                                                                        122.3819319, 61.97546553, 81.58795463, 59.62196684, 279.281844699999, 
                                                                                                                                                                        92.57094852, 138.0719232, 87.07945157, 57.26846815, 58.83746728, 
                                                                                                                                                                        164.7449084, 34.5179808, 174.1589031, 107.4764402, 88.6484507, 
                                                                                                                                                                        62.7599651, 127.0889293, 80.80345506, 55.69946902, 240.8413661, 
                                                                                                                                                                        45.5009747, 80.80345506, 94.13994765, 261.2383547, 76.09645768, 
                                                                                                                                                                        54.91496946, 61.97546553, 136.5029241, 180.4348997, 135.7184245, 
                                                                                                                                                                        96.49344634, 157.6844123, 61.19096597, 85.51045245, 75.31195812, 
                                                                                                                                                                        57.26846815, 94.13994765, 102.7694428, 257.3158569, 71.3894603, 
                                                                                                                                                                        72.17395986, 90.21744983, 53.34597033, 49.42347251, 65.89796335, 
                                                                                                                                                                        55.69946902, 160.037911, 69.03596161, 43.147476, 112.9679372, 
                                                                                                                                                                        61.19096597, 80.80345506), X5 = c(79.2, 72.6, 79.2, 81.4, 79.2, 
                                                                                                                                                                                                          83.6, 79.2, 77, 81.4, 83.6, 81.4, 81.4, 81.4, 116.6, 88, 79.2, 
                                                                                                                                                                                                          88, 83.6, 77, 77, 90.2, 129.8, 81.4, 79.2, 79.2, 77, 68.2, 167.2, 
                                                                                                                                                                                                          72.6, 74.8, 114.4, 105.6, 99, 92.4, 96.8, 81.4, 81.4, 74.8, 79.2, 
                                                                                                                                                                                                          101.2, 77, 92.4, 77, 88, 92.4, 90.2, 79.2, 83.6, 74.8, 74.8, 
                                                                                                                                                                                                          83.6, 72.6, 81.4, 79.2, 72.6, 59.4, 101.2, 103.4, 81.4, 83.6, 
                                                                                                                                                                                                          88, 83.6, 127.6, 92.4, 99, 107.8, 74.8, 94.6, 112.2, 83.6, 77, 
                                                                                                                                                                                                          79.2, 96.8, 81.4, 105.6, 114.4, 90.2, 81.4, 77, 83.6, 83.6, 88, 
                                                                                                                                                                                                          99, 88, 72.6, 92.4, 110, 92.4, 88, 114.4, 88, 57.2, 99, 101.2, 
                                                                                                                                                                                                          96.8, 92.4, 94.6, 112.2, 77, 79.2, 96.8, 88, 90.2, 77, 101.2, 
                                                                                                                                                                                                          88, 79.2, 169.4, 83.6, 123.2, 105.6, 81.4, 96.8, 83.6, 129.8, 
                                                                                                                                                                                                          96.8, 90.2, 74.8, 83.6, 138.6, 94.6, 70.4, 107.8, 77, 81.4, 88, 
                                                                                                                                                                                                          88, 81.4, 92.4, 129.8, 77, 81.4, 72.6, 74.8, 77, 74.8, 79.2, 
                                                                                                                                                                                                          77, 83.6, 94.6, 77, 83.6, 83.6), X6 = c(31.69752493, 57.86293466, 
                                                                                                                                                                                                                                                  31.54800831, 52.33081946, 62.34843348, 41.86465557, 74.60879689, 
                                                                                                                                                                                                                                                  29.45477553, 47.09773752, 55.02211875, 68.47861518, 62.19891685, 
                                                                                                                                                                                                                                                  42.31320545, 40.07045605, 41.86465557, 60.85326721, 51.58323633, 
                                                                                                                                                                                                                                                  29.3052589, 49.7890368, 52.77936935, 42.31320545, 32.14607481, 
                                                                                                                                                                                                                                                  27.21202612, 60.10568407, 52.77936935, 65.7873159, 55.91921851, 
                                                                                                                                                                                                                                                  30.35187529, 68.32909856, 37.22964013, 65.7873159, 43.35982184, 
                                                                                                                                                                                                                                                  36.18302374, 32.14607481, 27.21202612, 43.35982184, 60.40471732, 
                                                                                                                                                                                                                                                  42.31320545, 68.32909856, 36.78109025, 37.52867339, 47.69580403, 
                                                                                                                                                                                                                                                  69.37571495, 71.91749761, 34.38882422, 58.46100117, 37.22964013, 
                                                                                                                                                                                                                                                  45.15402137, 54.57356887, 27.06250949, 47.84532065, 27.95960926, 
                                                                                                                                                                                                                                                  49.7890368, 59.20858431, 39.62190617, 34.98689073, 42.91127196, 
                                                                                                                                                                                                                                                  56.06873514, 53.82598573, 71.16991447, 42.76175533, 43.06078859, 
                                                                                                                                                                                                                                                  50.83565319, 72.21653086, 50.38710331, 33.79075771, 77.74864606, 
                                                                                                                                                                                                                                                  33.79075771, 46.64918764, 47.09773752, 65.18924939, 57.56390141, 
                                                                                                                                                                                                                                                  36.63157362, 71.46894773, 41.41610569, 21.82942755, 54.27453562, 
                                                                                                                                                                                                                                                  61.89988359, 60.2552006999999, 39.02383966, 24.22169358, 61.30181709, 
                                                                                                                                                                                                                                                  36.33254037, 25.56734322, 53.37743585, 42.76175533, 40.51900593, 
                                                                                                                                                                                                                                                  55.17163538, 32.74414132, 48.29387053, 48.29387053, 40.66852255, 
                                                                                                                                                                                                                                                  48.74242042, 28.55767576, 26.46444299, 46.05112113, 25.11879334, 
                                                                                                                                                                                                                                                  56.66680165, 44.85498811, 37.08012351, 50.83565319, 87.46722682, 
                                                                                                                                                                                                                                                  37.22964013, 43.35982184, 48.74242042, 40.21997267, 31.69752493, 
                                                                                                                                                                                                                                                  20.33426128, 48.89193704, 57.41438478, 30.9499418, 18.68957838, 
                                                                                                                                                                                                                                                  33.79075771, 51.4337197, 35.73447386, 34.38882422, 42.46272208, 
                                                                                                                                                                                                                                                  35.88399049, 42.91127196, 22.72652731, 59.80665082, 49.04145367, 
                                                                                                                                                                                                                                                  69.6747482, 26.16540973, 36.482057, 25.56734322, 53.37743585, 
                                                                                                                                                                                                                                                  30.20235866, 44.55595486, 63.84359975, 55.321152, 58.31148455, 
                                                                                                                                                                                                                                                  52.33081946, 51.28420307, 53.37743585, 24.81976009, 42.76175533, 
                                                                                                                                                                                                                                                  39.32287291, 47.99483728, 63.54456649, 47.09773752, 44.85498811, 
                                                                                                                                                                                                                                                  55.76970189), X7 = c(1.83182999999999, 1.53152999999999, 2.7027, 
                                                                                                                                                                                                                                                                       1.5015, 3.42341999999999, 5.34534, 2.16215999999999, 1.32132, 
                                                                                                                                                                                                                                                                       2.37237, 1.92192, 1.68167999999999, 1.47147, 2.1021, 8.1081, 
                                                                                                                                                                                                                                                                       1.65165, 2.1021, 2.49249, 2.82282, 1.86185999999999, 2.61260999999999, 
                                                                                                                                                                                                                                                                       2.82282, 2.73273, 2.1021, 2.43243, 1.74174, 2.94294, 2.25225, 
                                                                                                                                                                                                                                                                       4.05405, 2.52252, 2.34234, 6.006, 2.4024, 2.91290999999999, 6.006, 
                                                                                                                                                                                                                                                                       1.68167999999999, 1.62162, 1.77177, 1.83182999999999, 2.43243, 
                                                                                                                                                                                                                                                                       2.37237, 1.89189, 2.67267, 2.22222, 3.54354, 2.04204, 1.83182999999999, 
                                                                                                                                                                                                                                                                       2.13213, 2.52252, 2.31230999999999, 2.76275999999999, 4.05405, 
                                                                                                                                                                                                                                                                       3.66365999999999, 2.13213, 5.52551999999999, 2.22222, 1.74174, 
                                                                                                                                                                                                                                                                       2.07207, 3.24324, 3.09309, 2.79279, 1.56156, 4.56456, 8.4084, 
                                                                                                                                                                                                                                                                       3.003, 1.2012, 4.2042, 1.77177, 4.8048, 3.21320999999999, 3.54354, 
                                                                                                                                                                                                                                                                       2.67267, 2.16215999999999, 3.75375, 2.01201, 3.36335999999999, 
                                                                                                                                                                                                                                                                       1.89189, 4.41441, 3.42341999999999, 2.31230999999999, 1.68167999999999, 
                                                                                                                                                                                                                                                                       3.63363, 2.94294, 2.91290999999999, 0.930929999999999, 1.92192, 
                                                                                                                                                                                                                                                                       2.13213, 3.21320999999999, 3.21320999999999, 3.18318, 1.83182999999999, 
                                                                                                                                                                                                                                                                       2.1021, 2.97297, 1.38137999999999, 4.08408, 2.04204, 4.11411, 
                                                                                                                                                                                                                                                                       2.97297, 4.14414, 3.48348, 1.62162, 1.62162, 3.9039, 2.85285, 
                                                                                                                                                                                                                                                                       2.1021, 2.52252, 2.01201, 2.1021, 1.95195, 1.5015, 0.87087, 3.33333, 
                                                                                                                                                                                                                                                                       1.5015, 3.6036, 2.64264, 1.02102, 2.85285, 2.49249, 2.13213, 
                                                                                                                                                                                                                                                                       3.96395999999999, 5.16516, 3.33333, 1.89189, 4.14414, 1.23122999999999, 
                                                                                                                                                                                                                                                                       2.76275999999999, 3.42341999999999, 2.79279, 2.82282, 4.74474, 
                                                                                                                                                                                                                                                                       7.8078, 1.92192, 2.19219, 1.98197999999999, 1.59159, 1.65165, 
                                                                                                                                                                                                                                                                       2.25225, 1.53152999999999, 3.75375, 1.65165, 2.43243, 2.43243, 
                                                                                                                                                                                                                                                                       1.74174, 2.16215999999999), X8 = c(202.2428448, 212.2969824, 
                                                                                                                                                                                                                                                                                                          200.3093568, 257.9272992, 235.4988384, 181.747872, 294.6635712, 
                                                                                                                                                                                                                                                                                                          186.7749408, 233.5653504, 237.4323264, 250.1933472, 266.0479488, 
                                                                                                                                                                                                                                                                                                          215.003865599999, 203.4029376, 204.5630304, 245.552975999999, 
                                                                                                                                                                                                                                                                                                          251.7401376, 160.479504, 203.016239999999, 240.9126048, 197.215776, 
                                                                                                                                                                                                                                                                                                          169.3735488, 136.5042528, 248.646556799999, 237.819024, 254.4470208, 
                                                                                                                                                                                                                                                                                                          235.885535999999, 194.1221952, 257.153904, 206.883216, 274.555296, 
                                                                                                                                                                                                                                                                                                          206.4965184, 182.9079648, 162.0262944, 158.9327136, 208.4300064, 
                                                                                                                                                                                                                                                                                                          247.486464, 221.5777248, 259.8607872, 185.614848, 170.146944, 
                                                                                                                                                                                                                                                                                                          220.8043296, 272.621808, 271.0750176, 164.34648, 242.8460928, 
                                                                                                                                                                                                                                                                                                          203.016239999999, 221.5777248, 249.419952, 173.2405248, 250.5800448, 
                                                                                                                                                                                                                                                                                                          147.718483199999, 225.058003199999, 225.8313984, 201.082752, 
                                                                                                                                                                                                                                                                                                          174.4006176, 212.2969824, 218.484144, 224.284608, 287.7030144, 
                                                                                                                                                                                                                                                                                                          218.870841599999, 253.286928, 195.282288, 288.4764096, 216.9373536, 
                                                                                                                                                                                                                                                                                                          169.7602464, 270.68832, 161.252899199999, 196.0556832, 222.35112, 
                                                                                                                                                                                                                                                                                                          248.259859199999, 209.203401599999, 189.481824, 270.3016224, 
                                                                                                                                                                                                                                                                                                          178.2675936, 135.34416, 232.4052576, 255.220416, 243.619488, 
                                                                                                                                                                                                                                                                                                          227.3781888, 145.0116, 235.885535999999, 174.4006176, 143.8515072, 
                                                                                                                                                                                                                                                                                                          247.0997664, 190.255219199999, 209.9767968, 209.5900992, 165.5065728, 
                                                                                                                                                                                                                                                                                                          208.0433088, 212.2969824, 199.149264, 205.3364256, 143.4648096, 
                                                                                                                                                                                                                                                                                                          146.1716928, 221.9644224, 167.0533632, 227.3781888, 169.3735488, 
                                                                                                                                                                                                                                                                                                          172.853827199999, 208.4300064, 336.8136096, 181.747872, 221.191027199999, 
                                                                                                                                                                                                                                                                                                          209.203401599999, 192.188707199999, 182.1345696, 132.2505792, 
                                                                                                                                                                                                                                                                                                          207.6566112, 225.4447008, 174.787315199999, 156.2258304, 177.880896, 
                                                                                                                                                                                                                                                                                                          216.1639584, 186.0015456, 233.178652799999, 206.4965184, 205.3364256, 
                                                                                                                                                                                                                                                                                                          152.745552, 151.5854592, 248.646556799999, 228.5382816, 264.1144608, 
                                                                                                                                                                                                                                                                                                          208.816704, 191.0286144, 149.6519712, 211.9102848, 159.3194112, 
                                                                                                                                                                                                                                                                                                          189.8685216, 271.0750176, 260.6341824, 226.218096, 221.5777248, 
                                                                                                                                                                                                                                                                                                          235.1121408, 257.153904, 186.0015456, 202.6295424, 170.5336416, 
                                                                                                                                                                                                                                                                                                          208.4300064, 241.2993024, 209.203401599999, 214.617168, 235.885535999999
                                                                                                                                                                                                                                                                       ), X9 = c(142.3422, 95.4954, 106.3062, 81.081, 113.5134, 102.7026, 
                                                                                                                                                                                                                                                                                 90.09, 82.8828, 88.2882, 104.5044, 81.081, 100.9008, 90.09, 198.197999999999, 
                                                                                                                                                                                                                                                                                 100.9008, 95.4954, 104.5044, 102.7026, 90.09, 120.7206, 104.5044, 
                                                                                                                                                                                                                                                                                 149.5494, 102.7026, 113.5134, 104.5044, 86.4864, 88.2882, 201.8016, 
                                                                                                                                                                                                                                                                                 95.4954, 111.7116, 218.0178, 131.5314, 88.2882, 147.7476, 93.6936, 
                                                                                                                                                                                                                                                                                 90.09, 106.3062, 106.3062, 95.4954, 192.7926, 120.7206, 109.9098, 
                                                                                                                                                                                                                                                                                 108.107999999999, 82.8828, 102.7026, 102.7026, 122.5224, 108.107999999999, 
                                                                                                                                                                                                                                                                                 91.8918, 106.3062, 124.3242, 120.7206, 115.3152, 84.6846, 97.2972, 
                                                                                                                                                                                                                                                                                 97.2972, 97.2972, 176.5764, 109.9098, 73.8738, 104.5044, 70.2702, 
                                                                                                                                                                                                                                                                                 158.5584, 104.5044, 109.9098, 117.117, 73.8738, 138.7386, 140.5404, 
                                                                                                                                                                                                                                                                                 95.4954, 100.9008, 93.6936, 100.9008, 99.0989999999999, 122.5224, 
                                                                                                                                                                                                                                                                                 156.7566, 104.5044, 97.2972, 136.9368, 99.0989999999999, 100.9008, 
                                                                                                                                                                                                                                                                                 106.3062, 158.5584, 86.4864, 84.6846, 117.117, 135.135, 129.7296, 
                                                                                                                                                                                                                                                                                 120.7206, 133.3332, 117.117, 99.0989999999999, 154.9548, 149.5494, 
                                                                                                                                                                                                                                                                                 135.135, 104.5044, 136.9368, 115.3152, 104.5044, 86.4864, 91.8918, 
                                                                                                                                                                                                                                                                                 124.3242, 142.3422, 118.9188, 99.0989999999999, 95.4954, 113.5134, 
                                                                                                                                                                                                                                                                                 255.8556, 75.6756, 154.9548, 153.153, 91.8918, 129.7296, 86.4864, 
                                                                                                                                                                                                                                                                                 162.162, 129.7296, 108.107999999999, 131.5314, 79.2792, 236.0358, 
                                                                                                                                                                                                                                                                                 117.117, 84.6846, 138.7386, 97.2972, 100.9008, 126.126, 100.9008, 
                                                                                                                                                                                                                                                                                 115.3152, 117.117, 297.297, 95.4954, 106.3062, 93.6936, 95.4954, 
                                                                                                                                                                                                                                                                                 104.5044, 118.9188, 106.3062, 102.7026, 102.7026, 99.0989999999999, 
                                                                                                                                                                                                                                                                                 106.3062, 102.7026, 90.09), X10 = c(22.2, 12.3, 12.2, 15.6, 11.7, 
                                                                                                                                                                                                                                                                                                                     7, 17.4, 21.4, 12.3, 16.4, 21.1, 27.6, 16.5, 10.5, 20.1, 10.1, 
                                                                                                                                                                                                                                                                                                                     11.7, 14.4, 11.4, 17.1, 9.65, 11.6, 11, 11, 15.7, 9.28, 15.8, 
                                                                                                                                                                                                                                                                                                                     8.59, 16.8, 15.2, 10, 6.85, 13.5, 8.5, 16.2, 30.2, 15.7, 13.1, 
                                                                                                                                                                                                                                                                                                                     14.1, 9.34, 20.5, 10.1, 15.1, 10.9, 15.5, 20.7, 23.3, 13.2, 15.6, 
                                                                                                                                                                                                                                                                                                                     11.3, 24.6, 8.54, 20.6, 6.75, 13, 9.36, 12.9, 26.6, 9.67, 8.13, 
                                                                                                                                                                                                                                                                                                                     16.8, 7.81, 7.46, 14.5, 21.3, 9.84, 19.2, 5.64, 9.55, 15.4, 8.53, 
                                                                                                                                                                                                                                                                                                                     7.67, 9.9, 21.6, 10.2, 8.45, 9.78, 10.1, 10.3, 19.6, 8.2, 11.8, 
                                                                                                                                                                                                                                                                                                                     15, 13.1, 26.3, 8.9, 14.2, 11.1, 6.1, 14.2, 14.2, 12.7, 9.11, 
                                                                                                                                                                                                                                                                                                                     8.82, 12.3, 10.7, 9.87, 6.48, 9.61, 10.9, 21.8, 9.52, 15.3, 15.2, 
                                                                                                                                                                                                                                                                                                                     8.91, 10.7, 17.7, 4.39, 12, 8.8, 8.17, 11.9, 10.5, 7.3, 8.29, 
                                                                                                                                                                                                                                                                                                                     18.3, 15.1, 18.8, 16.3, 4.69, 8.67, 11.9, 3.68, 19, 14.5, 12.3, 
                                                                                                                                                                                                                                                                                                                     7.34, 12.1, 6.43, 11.7, 28, 14.8, 11.8, 11.9, 10.5, 14.4, 34.8, 
                                                                                                                                                                                                                                                                                                                     7.49, 11.7, 24.8, 11.3, 19, 17.8), X11 = c(1.15, 2.16, 1.13, 
                                                                                                                                                                                                                                                                                                                                                                0.39, 0.68, 26.51, 0.82, 1.9, 0.56, 1.27, 7.68, 1.5, 0.6, 6.48, 
                                                                                                                                                                                                                                                                                                                                                                0.51, 2.26, 0.71, 0.21, 3.84, 3, 1.52, 1.28, 7.03, 13.56, 1.15, 
                                                                                                                                                                                                                                                                                                                                                                5.67, 1.49, 1.34, 1.36, 4.11, 4.64, 2.88, 1.38, 1.55, 1.44, 0.63, 
                                                                                                                                                                                                                                                                                                                                                                0.94, 1.12, 7.05, 0.82, 1.24, 1.26, 1.54, 0.67, 0.49, 0.13, 3.14, 
                                                                                                                                                                                                                                                                                                                                                                0.83, 0.49, 2.26, 5.68, 1.07, 0.69, 11.37, 0.69, 10.06, 0.42, 
                                                                                                                                                                                                                                                                                                                                                                2.35, 0.85, 4.48, 1.18, 21, 3.26, 13.2, 0.9, 1.84, 1.47, 1.32, 
                                                                                                                                                                                                                                                                                                                                                                0.65, 2.84, 0.41, 10.87, 2.28, 2.37, 7.5, 3.48, 2.53, 3.2, 2.23, 
                                                                                                                                                                                                                                                                                                                                                                11.82, 0.85, 2, 38.95, 0.77, 0.9, 1.21, 3.53, 0.69, 2.48, 2.38, 
                                                                                                                                                                                                                                                                                                                                                                0.34, 2.52, 0.69, 1.31, 7.64, 1.78, 1.9, 1.21, 0.58, 99.22, 1.46, 
                                                                                                                                                                                                                                                                                                                                                                2.23, 0.81, 0.93, 2.62, 1.37, 1.36, 3.85, 1.25, 1.08, 4.91, 0.99, 
                                                                                                                                                                                                                                                                                                                                                                7.93, 0.64, 12.36, 1.62, 5.21, 1.1, 27.15, 1.33, 0.38, 1.86, 
                                                                                                                                                                                                                                                                                                                                                                3.54, 0.34, 1.81, 3.29, 4.55, 1.53, 3.54, 6.25, 0.94, 2.04, 1.09, 
                                                                                                                                                                                                                                                                                                                                                                14.92, 1.26, 0.71, 2.77, 2.4, 3.52, 1.69, 0.83, 2.98, 1.02), 
               X12 = c(12.2, 23, 72.1, 4.75, 25.4, 29.6, 11.3, 8.96, 22.9, 
                       14.5, 42.6, 6.16, 8.49, 33.4, 11.4, 28.5, 11.5, 15.1, 27, 
                       18.8, 23.4, 23.2, 17.1, 14.8, 30.3, 28.4, 21.7, 57.5, 16.1, 
                       35.2, 40.3, 48.2, 14.6, 45.2, 8.79, 14.7, 7.33, 39.4, 44.9, 
                       19.1, 8.53, 20.1, 10.9, 19.9, 21.1, 7.01, 36, 28, 38.5, 37.6, 
                       57.1, 7.02, 6.21, 105, 22.5, 7.06, 23.3, 36.3, 31.7, 28.9, 
                       23.3, 42.1, 37.5, 91.8, 5.47, 36, 16.2, 49.1, 35, 26, 19.6, 
                       18.2, 33.8, 10.9, 49.9, 21.6, 77.7, 37.8, 13.1, 18.5, 32.9, 
                       29.7, 48.9, 2.17, 6.98, 29.5, 29.8, 25.6, 40.4, 13, 10.4, 
                       35.8, 7.11, 21.1, 43.5, 25.7, 28, 29.1, 15.8, 3.59, 7.28, 
                       16.5, 17, 14.1, 25.1, 14.8, 14.6, 32.4, 5.28, 14.1, 12.6, 
                       3.92, 31.8, 20.6, 57.7, 5.49, 12.2, 22.9, 10.3, 36.6, 31.4, 
                       5.6, 31.4, 3.58, 21.2, 53.6, 49.4, 21.4, 34.7, 36.4, 12.5, 
                       15.4, 12.5, 21.3, 8.12, 13.2, 10.1, 23.4, 15.1, 28.5, 28.2, 
                       7.14, 14.7), X13 = c(0, 1.8, 0.1, 0.7, 0.7, 0.2, 0.8, 0.3, 
                                            0.2, 0.7, 0.4, 0.5, 0, 0.4, 0.2, 0.2, 0.3, 0, 0.1, 0.3, 0.1, 
                                            0.6, 0.6, 0.3, 0.9, 0.4, 0.2, 0.8, 56.6, 1.3, 0.6, 1.6, 0.6, 
                                            0.1, 0.4, 0.2, 0.6, 0, 0.6, 0, 0.1, 0.4, 0.1, 1, 3, 0.8, 
                                            0.9, 0.1, 2, 0.1, 0.1, 4.7, 0, 0.5, 0, 0.1, 0.2, 0.6, 0.7, 
                                            1.5, 0.8, 0.9, 0.8, 0.9, 0.5, 0.9, 0.4, 0.9, 1.5, 0.4, 2.3, 
                                            0.4, 0.9, 0.5, 0.4, 0, 2.5, 0.3, 0.7, 1.5, 0.2, 0.1, 0.4, 
                                            0.3, 1.1, 0, 0.4, 0, 0.1, 0.3, 0.3, 0, 0.8, 1.3, 0.4, 0.5, 
                                            1.2, 1, 0.4, 0.3, 0.5, 0, 0, 0.5, 0.3, 0, 0.2, 0.9, 0.7, 
                                            0, 0, 0.4, 0.4, 0.2, 1.1, 0.3, 1.2, 0.2, 0.2, 0.7, 0.2, 1.7, 
                                            0.3, 0.6, 0.5, 0.4, 0.3, 1.2, 0.6, 0.5, 0.4, 0.2, 0.1, 0.9, 
                                            0.5, 0.6, 0.2, 0.3, 0, 0.5, 0.4, 0.2, 0.9), X14 = c(65, 74, 
                                                                                                68, 63, 114, 71, 58, 75, 71, 94, 72, 63, 62, 97, 52, 85, 
                                                                                                56, 80, 71, 72, 71, 97, 74, 65, 68, 64, 64, 110, 65, 69, 
                                                                                                65, 64, 62, 60, 65, 72, 69, 58, 64, 54, 59, 77, 65, 70, 57, 
                                                                                                54, 73, 92, 76, 182, 112, 59, 67, 85, 63, 48, 61, 78, 58, 
                                                                                                81, 72, 60, 95, 132, 54, 69, 74, 72, 70, 85, 66, 67, 72, 
                                                                                                73, 56, 60, 69, 71, 64, 81, 69, 70, 57, 58, 54, 64, 60, 57, 
                                                                                                61, 47, 60, 72, 59, 55, 55, 74, 67, 93, 70, 61, 70, 49, 77, 
                                                                                                72, 75, 62, 59, 81, 76, 69, 55, 72, 59, 64, 89, 72, 50, 58, 
                                                                                                74, 68, 67, 89, 75, 60, 77, 69, 53, 76, 61, 90, 66, 76, 59, 
                                                                                                44, 60, 69, 59, 64, 72, 63, 69, 63, 70), total = c(6188331, 
                                                                                                                                                   11432135, 9071216, 8856043, 7600919, 7310208, 1793383, 7276145, 
                                                                                                                                                   9605078, 13603589, 7893161, 14845141, 13972354, 16163738, 
                                                                                                                                                   11098252, 15797133, 16727169, 14325755, 14169519, 14849646, 
                                                                                                                                                   12321433, 11356091, 6141549, 7367357, 5546644, 7184518, 6675717, 
                                                                                                                                                   5388676, 5559357, 9755207, 9321726, 9023634, 2803159, 9562367, 
                                                                                                                                                   6882326, 5514297, 8460187, 6701281, 8116523, 2877017, 8124274, 
                                                                                                                                                   8091025, 1995441, 7200738, 9781964, 9913532, 6988329, 8092894, 
                                                                                                                                                   9142106, 4765359, 11357280, 6690537, 2293559, 7354137, 7336240, 
                                                                                                                                                   8521415, 8509766, 8027941, 10796926, 10565994, 7530208, 8511576, 
                                                                                                                                                   11255767, 11206524, 7456083, 5966828, 9071099, 8483974, 8161814, 
                                                                                                                                                   5576187, 8496265, 6056579, 7885259, 8506625, 16497934, 9169591, 
                                                                                                                                                   10478059, 13309231, 17638017, 13806592, 16827390, 18019971, 
                                                                                                                                                   6622133, 11867686, 8551369, 5839951, 21206918, 15032242, 
                                                                                                                                                   10882909, 11157812, 2291965, 16656335, 3348443, 17769050, 
                                                                                                                                                   15360414, 10409531, 18071928, 17095556, 2950090, 13968352, 
                                                                                                                                                   13400391, 14778644, 13675474, 7875209, 14793548, 12932461, 
                                                                                                                                                   13340486, 18289402, 9782114, 12962250, 15451028, 8229705, 
                                                                                                                                                   14663503, 14253617, 3289372, 16912942, 8727386, 12024105, 
                                                                                                                                                   16372540, 12686343, 8195228, 9090207, 11006956, 8046474, 
                                                                                                                                                   13301476, 11472474, 20224012, 10786250, 12684928, 12972200, 
                                                                                                                                                   5471791, 3037373, 6146142, 4207054, 3678159, 7310788, 4678621, 
                                                                                                                                                   4992290, 4573959, 8524345, 6191269, 7374788, 7583572)), row.names = c(NA, 
                                                                                                                                                                                                                         -143L), class = c("tbl_df", "tbl", "data.frame"))
sep_issue1 %>% 
  corncob::bbdml(formula = cbind(W, total - W) ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + 
                   X9 + X10 + X11 + X12 + X13 + X14,
                 phi.formula = ~ 1,
                 data = .)

Plots

Reorder by color and sample, add facet, shape

Intuition adding covariates in formula

Hello there,

I am having difficulties to understand in my dataset how I should specify the formula for testing for differential abundance. I have a marine multiyear time series with two fractions differentiated (let’s say free living FL bacteria and particle attached AT).

I want to test for differential abundance between the fractions. As there are monthly seasonal changes as covariates, should include them in the formula? How should I write it down?

Can I put directly only this?

test <- differentialTest(formula = ~fraction, 
                         formula_null = ~1,
                         phi.formula = ~fraction,
                         phi.formula_null = ~fraction,
                         data = diff.phy, 
                         test = 'Wald', 
                         boot = FALSE,
                         fdr_cutoff = 0.05)

Here in my mind I am only testing for abundance differentiation (not variation) for the fraction variable.

Or should I specify in both sites that month is a covariate?

test <- differentialTest(formula = ~fraction + month, 
                         formula_null = ~ month,
                         phi.formula = ~fraction + month,
                         phi.formula_null = ~ fraction + month,
                         data = diff.phy, 
                         test = 'Wald', 
                         boot = FALSE,
                         fdr_cutoff = 0.05)

Additionally, in the parameter interpretation section with the soil dataset from Dr. Thea Whitman you state ‘ we can see that the DayAmdmt21 abundance coefficient is negative and statistically significant`. How should I interpret in an intuitive way this coefficient of change? A value of -0.16 related to a relative abundance perspective what does it mean?

Thank you for your time and to develop this package, it is really useful

Add catch for infeasible initial values

Change manually with printed out warning. Happens when phi is numerically equal to 1

DifferentialTest giving on a subset of samples gives "Error: $ operator is invalid for atomic vectors"

Hi @bryandmartin ,
Really enjoying Corncob so far! I'm running into the classic error that I can't figure out.

Error: $ operator is invalid for atomic vectors

Summary: I've got a phyloseq object with 4 groups across 2 time points (wk0, wk6).

> ps
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 2239 taxa and 63 samples ]
sample_data() Sample Data:       [ 63 samples by 34 sample variables ]
tax_table()   Taxonomy Table:    [ 2239 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 2239 tips and 2238 internal nodes ]
refseq()      DNAStringSet:      [ 2239 reference sequences ]

What I'm interested in is looking at each group's change separately across time.
For example, I can successfully subset my phyloseq object like so:

> subset_samples(ps, na_removed_group=="VWR")
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 2239 taxa and 12 samples ]
sample_data() Sample Data:       [ 12 samples by 34 sample variables ]
tax_table()   Taxonomy Table:    [ 2239 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 2239 tips and 2238 internal nodes ]
refseq()      DNAStringSet:      [ 2239 reference sequences ]

So, with corncob (0.1.0):

cc.VWR <- differentialTest(formula = ~ time,
                                  phi.formula = ~ time,
                                  formula_null = ~ 1,
                                  phi.formula_null = ~ 1,
                                  test = "Wald", boot = F,
                                  data = subset_samples(ps, na_removed_group=="VWR"),
                                  fdr_cutoff = 0.05)
#which gives me the error
Error: $ operator is invalid for atomic vectors

This problem only seems to occur with this particular group. Using the same formula but subsetting with any of the other 3 groups runs successfully without any errors. If I run the whole data set together, it also works just fine. There is nothing special about this VWR group as far as I can see.
Any thoughts?

Add zero counts to filter_discriminant

metabolomics tutorial

Not now (!) but let's find a metabolomics dataset and put together a vignette showing people how to use CORNCOB in this setting.

Vignette missing?

I just installed corncob, and the vignette seems to be missing / not rendered during the install

> library(corncob)
> browseVignettes(package = "corncob")
No vignettes found by browseVignettes(package = "corncob")

phyloseq as suggests not depends

Hi Bryan -- phyloseq can be very difficult to install, and it's not really critical for any of the functionality of corncob. Could it please be downgraded from a required dependence to suggested?

This would greatly speed up my installation 🎉

p-values between differentialTest() and bbdml()

Hi, thank you for developing corncob.
As I was going through the vignette, I tried to look at the direction of the differential abundance of the taxa by putting them through bbdml() individually (as an attempt to bypass #15 ). I noticed that OTU.19 was reported as a significant taxa:

otu_to_taxonomy(OTU = da_analysis$significant_taxa, data = soil)

      OTU.19  
3.323423e-03

However, when I put it through bbdml():

test <- bbdml(formula = OTU.19 ~ DayAmdmt,
      phi.formula = ~ DayAmdmt,
      data = soil)

summary(test)

Call:
bbdml(formula = OTU.19 ~ DayAmdmt, phi.formula = ~DayAmdmt, data = soil)

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
mu.(Intercept)  -1.64987    0.04930 -33.467  < 2e-16 ***
mu.DayAmdmt21    0.02926    0.05321   0.550  0.58672    <-- ?
phi.(Intercept) -5.24564    0.35403 -14.817 8.87e-15 ***
phi.DayAmdmt21  -1.79881    0.50455  -3.565  0.00133 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-likelihood: -281.05

It appears I've got a non-significant result. This raises the question:

Shouldn't differentialTest() be more stringent due to FDR correction? (As seen for OTU.4206 in the vignette).
Could I still use the coefficient 0.02926 to say that OTU.19 is differentially more abundant in DayAmdmt21?

Thanks again for your time

Kevin

Curl error when trying to install package via install_github()

I get this curl error when trying to install the package:

devtools::install_github("bryandmartin/corncob")

Error in curl::curl_fetch_memory(url, handle = h) : 
Failed to connect to api.github.com port 443: Connection refused

genInits not finished

metadata in `soil`

Add the actual names for the metadata in soil. Here's the key:

Amdmt: 0 (no additions), 1 (biochar additions), 2 (fresh biomass additions).
Day: 0 (initial sampling point), 1, ( 12 days after treatment additions), and 2 (82 days after treatment additions).

Handling taxa/OTUs which have singular Hessian matrices?

I'm exploring the output from differentialTest() and I've observed that, for a small number of taxa, the coefficients table isn't fully populated. More specifically, it returns "Estimate" values for each of the terms, but does not populate the "Std. Error", "t value", or "pvalue" fields. For example: this is the output for one of my OTUs:

Applying bbdml() to just the single taxa and then passing the output to waldt() reveals that the issue is due to a singlular Hessian matrix. In this case, waldt() will still create the coefficients table, but that table will only have values for the Estimates.

NOTE: I'm working with a very small dataset (in terms of read counts), so it might just be the particular material I'm working with, with much larger (more realistic) read counts there may be less of a chance of encountering this issue.

I thought I'd bring this to @bryandmartin and @adw96's attention.

I don't know if there's a way to get around the non-invertable matrix problem (possibly?) but I would encourage the developers to include in the documentation and/or error/warning message a discussion/information about what this output means.

Add testthat checks

differentialTest

Want not just that they are significant, but direction / plot of significant taxa.

Change plot error bars to (2.5%-97.5%)

phi.init no longer used in main function

To remove

inits length check

If inits are wrong length should add a check so it fails right away

Interpreting Effect Sizes

Hi Bryan,

Thanks for the work on the tool. I was wondering if you had any suggestions on how to interpret the effect sizes calculated from corncob. Are there any current guidelines for this?

I find that looking at effect sizes and significance tends to be more useful then significance alone, however, I'm not sure how to interpret the magnitude of the effect sizes reported by corncob's lot.differentialTest (I'm guessing its a plot of the coefficients of the model?) For instance if modelling height vs microbe abundance would a value of -0.02 indicate that for every unit increase in height, relative abundance of that microbe goes down by 0.02?

Thanks
Jacob Nearing

Ellipsis argument in plot not being passed properly

Error in trust::trust -- parinit not feasible

Hi @bryandmartin! I've got a dataset that is breaking trust::trust:

Error in trust::trust(objfun, parinit = theta.init, W = W, M = M, X = X.b, :
parinit not feasible

library(tidyverse)
breaks_corn <- structure(list(W = c(22012, 59233, 223268, 122825, 53763, 29208, 
                     20568, 175146, 242484, 3069, 29495, 461774, 277620, 425460, 47622, 
                     268909, 133373, 232651, 201, 81882, 81301, 6125, 132758, 143665, 
                     74543, 34721, 10149, 4575, 2217, 85642, 289783, 35408, 10529, 
                     36338, 83723, 75024, 18172, 12845, 118999, 17063, 170357, 45213, 
                     6455, 72889, 114416, 338174, 18268, 99339, 134994, 74336, 4335, 
                     6503, 34263, 146457, 30725, 54379, 62320, 5825, 573, 418272, 
                     22355, 450960, 266391, 239118, 95671, 17589, 109710, 151699, 
                     1076, 2278, 202113, 157465, 94524, 308358, 60273, 70820, 97379, 
                     235883, 150704, 466330, 660503, 75217, 1873, 558812, 113987, 
                     68933, 205555, 280552, 19630, 1847, 13782, 79365, 20071, 350286, 
                     91333, 515, 801971, 346819, 16109, 40779, 92535, 78280, 113517, 
                     32868, 138647, 14270, 117617, 47739, 107495, 12350, 77869, 37222, 
                     248932, 68421, 49951, 2481, 132092, 38095, 58736, 252037, 169843, 
                     39361, 39701, 33712, 332291, 205788, 84891, 36467, 52626, 22417, 
                     10488, 54524, 89814, 129946, 157795, 41366, 151332, 253242, 46999, 
                     216132, 73530, 40310, 2421), M = c(6188331, 11432135, 9071216, 
                                                        8856043, 7600919, 7310208, 1793383, 7276145, 9605078, 13603589, 
                                                        7893161, 14845141, 13972354, 16163738, 11098252, 15797133, 16727169, 
                                                        14325755, 14169519, 14849646, 12321433, 11356091, 6141549, 7367357, 
                                                        5546644, 7184518, 6675717, 5388676, 5559357, 9755207, 9321726, 
                                                        9023634, 2803159, 9562367, 6882326, 5514297, 8460187, 6701281, 
                                                        8116523, 2877017, 8124274, 8091025, 1995441, 7200738, 9781964, 
                                                        9913532, 6988329, 8092894, 9142106, 4765359, 11357280, 6690537, 
                                                        2293559, 7354137, 7336240, 8521415, 8509766, 8027941, 10796926, 
                                                        10565994, 7530208, 8511576, 11255767, 11206524, 7456083, 5966828, 
                                                        9071099, 8483974, 8161814, 5576187, 8496265, 6056579, 7885259, 
                                                        8506625, 16497934, 9169591, 10478059, 13309231, 17638017, 13806592, 
                                                        16827390, 18019971, 6622133, 11867686, 8551369, 5839951, 21206918, 
                                                        15032242, 10882909, 11157812, 2291965, 16656335, 3348443, 17769050, 
                                                        15360414, 10409531, 18071928, 17095556, 2950090, 13968352, 13400391, 
                                                        14778644, 13675474, 7875209, 14793548, 12932461, 13340486, 18289402, 
                                                        9782114, 12962250, 15451028, 8229705, 14663503, 14253617, 3289372, 
                                                        16912942, 8727386, 12024105, 16372540, 12686343, 8195228, 9090207, 
                                                        11006956, 8046474, 13301476, 11472474, 20224012, 10786250, 12684928, 
                                                        12972200, 5471791, 3037373, 6146142, 4207054, 3678159, 7310788, 
                                                        4678621, 4992290, 4573959, 8524345, 6191269, 7374788, 7583572
                     ), X = c(95.70894677, 47.85447339, 43.147476, 77.66545681, 90.21744983, 
                              148.2704175, 72.17395986, 52.56147077, 106.6919407, 79.23445594, 
                              91.78644896, 32.94898168, 65.11346379, 117.674934599999, 47.06997382, 
                              92.57094852, 131.7959267, 70.60496073, 103.5539424, 80.0189555, 
                              163.1759093, 117.674934599999, 75.31195812, 54.1304699, 65.89796335, 
                              111.398937999999, 94.13994765, 313.7998255, 152.977414899999, 
                              90.21744983, 285.5578412, 98.84694503, 54.1304699, 157.6844123, 
                              76.88095724, 65.11346379, 72.95845943, 72.95845943, 55.69946902, 
                              73.74295899, 89.43295026, 183.5728979, 66.68246292, 135.7184245, 
                              70.60496073, 106.6919407, 72.95845943, 42.36297644, 80.80345506, 
                              85.51045245, 74.52745855, 105.907441099999, 61.19096597, 120.8129328, 
                              74.52745855, 53.34597033, 80.0189555, 80.0189555, 127.8734289, 
                              131.0114271, 71.3894603, 498.157223, 97.2779459, 121.5974324, 
                              79.23445594, 185.141897, 98.06244547, 168.6674062, 69.03596161, 
                              271.436849, 97.2779459, 78.44995637, 77.66545681, 79.23445594, 
                              73.74295899, 83.15695375, 128.6579285, 203.185387, 123.9509311, 
                              67.46696248, 50.99247164, 99.63144459, 112.183437599999, 38.44047862, 
                              49.42347251, 98.84694503, 225.9358744, 89.43295026, 123.166431499999, 
                              38.44047862, 43.147476, 122.3819319, 61.97546553, 81.58795463, 
                              59.62196684, 279.281844699999, 92.57094852, 138.0719232, 87.07945157, 
                              57.26846815, 58.83746728, 164.7449084, 34.5179808, 174.1589031, 
                              107.4764402, 88.6484507, 62.7599651, 127.0889293, 80.80345506, 
                              55.69946902, 240.8413661, 45.5009747, 80.80345506, 94.13994765, 
                              261.2383547, 76.09645768, 54.91496946, 61.97546553, 136.5029241, 
                              180.4348997, 135.7184245, 96.49344634, 157.6844123, 61.19096597, 
                              85.51045245, 75.31195812, 57.26846815, 94.13994765, 102.7694428, 
                              257.3158569, 71.3894603, 72.17395986, 90.21744983, 53.34597033, 
                              49.42347251, 65.89796335, 55.69946902, 160.037911, 69.03596161, 
                              43.147476, 112.9679372, 61.19096597, 80.80345506)), row.names = c(NA, 
                                                                                                -143L), class = c("tbl_df", "tbl", "data.frame"))
breaks_corn %>%
  bbdml(formula = cbind(W, M - W) ~ X,
        phi.formula = ~ X,
        data = .)

Help appreciated!

Relative abundance flag?

Hi Bryan,

Corncob takes absolute values and then converts counts to relative abundance within the corncob or differentialTest function correct? Is there an argument already or might I suggest an enhancement to flag that you are giving the corncob or the differentialTest function relative abundances rather than absolute? I was looking to filter out relatively low abundance organisms to speed up initial testing and double check things before I run all the ASVs. I can reduce the data set in other ways, just thought a flag for relative abundance would improve functionality. Thoughts? Looking forward to the paper!

Cheers,

Remove term absolute abundances from vignette

Rename as "counts" or "Read counts" instead

option to define order of levels in bbmdl?

Hi @bryandmartin,

This might be a redundant question - Is there a way to determine which level for a mapping variable is tested against the rest (i.e., control)? right now the model orders alphanumerically.

Thanks,
Gil

Suppress printed warnings in STAMPS vignette

Excessive warnings for discriminant taxa in differentialTest with bootstrap methods

Happens within call to doBoot if bootstrap simulations are discriminant

Unable to identify Differentially abundant or variable taxa?

So I'm going through the vignette for corncob and for some reason calling
$DV or $DA on the results of differentialTest is not possible?
Has that been removed or?

Best,
N

Check zero overdispersion with logit

seq_depth <- rpois(20, lambda = 1000)
my_counts <- rpois(20, lambda = 50)
my_covariate <- cbind(rep(c(0,1), each = 10))
colnames(my_covariate) <- c("X1")

test_data <- data.frame("W" = my_counts, "M" = seq_depth, my_covariate)

bbdml(formula = cbind(W, M - W) ~ X1,
      phi.formula = ~ X1,
      data = test_data,
      link = "logit",
      phi.link = "logit")

FDR correction of both DA and DV

Thanks for developing corncob, it seems like a great tool!

I noticed today that the FDR correction is for both the DA and DV p-values in differentialTest. I was only interested in testing for DA in the particular analysis I was running and I realized I was being unnecessary stringent by correcting for the DV values as well. I think others will likely be confused by this as well so I thought I'd mention it in case it could be clarified for future releases!

Also, I would have thought since testing for DA and DV were totally different hypothesis tests that the p-values from these tests shouldn't be corrected together anyway (but that's likely something I'm just not understanding).

Thanks,

Gavin

statdivlab / corncob Goto Github PK

corncob's Introduction

corncob

Installation

Docker

Use

Documentation

Citation

Bug Reports / Change Requests

corncob's People

Contributors

Stargazers

Watchers

Forkers

corncob's Issues

Recommend Projects

Recommend Topics

Recommend Org