dataesr/scientific-paragraphs-categorization huggingface dataset:
Bootstrap Confidence Intervals In this section, we use the parametric bootstrap method to approximate the distribution of MLEs of the GDUSLN model parameters. Then, we can use the bootstrap distribution to estimate confidence intervals of each parameter for the fitted GDUSLN distribution. Let Θ be a MLE on the set of parameters of interest Θ = (α, µ, σ) using a given dataset {x 1 , x 2 ,. , x n }. The bootstrap is a method to estimate the distribution of statistic Θ by getting a random sample Θ * 1 , Θ * 2 ,. , Θ * B for Θ based on B random samples that are drawn with replacement from {x 1 , x 2 ,. , x n }, see  [14]. The bootstrap sample Θ * 1 , Θ * 2 ,. , Θ * B can be used to construct bootstrap confidence intervals for the parametric set Θ = (α, µ, σ) of the GDUSLN distribution. Thus, using the following formulae, we calculate the 100 × (1δ)% bootstrap confidence intervals for parameters: α ∈ α ∓ z δ/2 ŝe α,boot , µ ∈ μ ∓ z δ/2 ŝe µ,boot , σ ∈ σ ∓ z δ/2 ŝe σ,boot , where z δ denotes the δ th percentile of the bootstrap sample and, for θ ∈ {α, µ, σ}, ŝe θ,boot = 1 B B ∑ b=1 θ * b - 1 B B ∑ b=1 θ * b 2.
Cancer Survival Data First, we utilize the dataset from  [15]  as a biological dataset, which represents an uncensored univariate dataset comprised of the remission times (in months) of a random sample of 128 bladder cancer patients. The descriptive measures of the real dataset, which include sample size (n), minimum (min), first quartile (Q 1 ), median (Md), third quartile (Q 3 ), maximum (max), and inter-quartile range (IQR) are given in Table  9. We also investigate the empirical hazard rate function for the biology dataset using the idea of a total time on test (TTT) plot. The TTT plot is a graph being used to distinguish between several types of aging as displayed in the hazard rate shapes. The common shapes of the hazard rate possess constant, increasing, decreasing, bathtub, and upside-down bathtub shapes, and can be identified by using the TTT plot by the following methods:.
• A plot around the diagonal indicates a constant hazard rate, that is, the failure times can be considered exponentially distributed. For more about the TTT plot, see details in  [16]. The TTT plot is drawn by plotting T i n = i ∑ r=1 x r:n + (n -i)x i:n n ∑ r=1 x r:n against i/n, where i = 1, 2,. , n and x r:n , r = 1, 2,. , n are the order statistics of the sample. Thus, the plot in Figure  4  indicates that this dataset represents an upside-down bathtub shaped hazard rate function. This case is covered by the characteristics of the GDUSLN distribution. To show the potential advantage of the GDUSLN distribution, the following distributions are considered for comparison:.
• The modified Weibull (MoW) distribution (see  [20] ) with pdf f (x) = µσ x α µ-1 exp x α µ + ασ 1 -e (x/α) µ , x > 0, α, µ, σ > 0. • The Weibull distribution with pdf f (x) = α σ x σ α-1 e -(x/σ) α , x > 0, α, σ > 0. We compare the competitive models to the proposed models using the following statistical tools: negative log-likelihood (-log L), Kolmogorov-Smirnov (KS), Cramér-von Misses (W * ), Anderson-Darling (A * ) statistics, Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) values. Tables  10  and 11  display the corresponding MLEs and goodness-of-fit (GOF) statistics of the considered distributions corresponding to the bladder cancer dataset. From these tables, we see that the GOF statistics values of the GDUSLN distribution are smaller than those of the other compared distributions.
It can also be noted that the optimization algorithm possesses successful convergence as indicated in Section 5. 1. The empirical cdf and quantile-quantile (Q-Q) plots for the real dataset are given in Figure  5. This figure shows some nice-shaped curves for those empirical and fitted functions. Thus, we conclude that the GDUSLN distribution is the most suitable distribution for this dataset compared to that of the other distributions. Now, the Hessian matrix corresponding to bladder cancer dataset is obtained as It is observed that the determinant value of the observed information matrix (|J( Θ)|) is non-zero, and hence satisfies the non-singularity condition of the information matrix. Now, Table  12  provides the 95 percent asymptotic confidence intervals for the GDUSLN parameters. Next, we focus on estimating the parameters of the GDUSLN distribution using the Bayesian procedure based on the above discussed univariate bladder cancer survival dataset.
In the context of Bayesian estimation, the analysis was performed using the Metropolis-Hastings algorithm of the MCMC method with 1000 iterations. For comparing Bayes estimates with the MLEs, both the estimates of the GDUSLN parameters for the real dataset are given in Table  13. The numerical computations on Bayesian estimation are done using RStudio software. In this subsection, for the considered dataset, we utilize the computed MLEs to construct the 95 percent bootstrap confidence intervals for the parameters α, µ, and σ. Based on the GDUSLN distribution, we simulate 1001 samples of the same size as the real dataset, with true values of the parameters chosen as MLEs of the respective parameters. We calculate the MLEs α * b , μ * b and σ * b , for b ∈ {1, 2,. , 1001} for each sample obtained.
Likelihood Ratio Test We also utilized the likelihood ratio (LR) test for comparing the GDUSLN distribution, which has an additional parameter α with the LN distribution based on the above discussed bladder cancer survival dataset. The LR statistic for comparing the nested model H 0 : LN against H A : GDUSLN is LR = -2 log likelihood under the null hypothesis likelihood in the whole parameter space which asymptotically follows a chi-square distribution having d degrees of freedom, d being the number of additional parameters in the GDUSLN model. By using this result and standard statistical tables, we can obtain critical values for the LR test statistics for the given bladder cancer dataset. Table  15  includes the LR statistic and the corresponding p-value. Given, the values of test statistic and the associated p-value, we reject the null hypothesis for the above discussed bladder cancer dataset and conclude that the GDUSLN distribution provides a significantly better representation than the LN distribution.
Stanford Heart Transplant Data In this application, we validate the prominence of the GDUSLN regression model by applying it to the real dataset, the renowned Stanford heart transplant data. The dataset is given in  [21] , which can also be found in the R package p3state. msm. The goal of this study is to investigate the survival times (y i ) of patients with covariates x 1 -year of acceptance to the program, x 2 -age of patient (in years), and x 3 -previous surgery status (1 = yes, 0 = no). In this study, the transplant indicator is used as the censoring variable.
Results Using the GDUSLN Regression Model The fitted non-linear regression model is given by x i = exp(τ 0 + τ 1 v 1 + τ 2 v 2 + τ 3 v 3 + σz i ), where the response variable x i is observed follows a random variable following the GDUSLN distribution. In Table  16 , we compare the performance of the GDUSLN regression model with that of the LN regression model, as well as the summaries due to the real dataset, which include estimates of all parameters, negative log-likelihood (-l(ψ)), and the value of AIC.
Further, we adopted the parametric bootstrap technique to obtain confidence intervals for the model parameters. More importantly, we introduced a parametric regression model and a Bayesian regression method based on the new distribution. Simulation studies were conducted to analyze the performance of ML and Bayesian estimates of the GDUSLN parameters and they confirm their consistency. The usefulness of the new model was illustrated by two applications of real datasets, which are related to the field of biology and used goodness-of-fit tests. The novel model consistently outperforms previous models in the literature in terms of fitting. We anticipate that the suggested model would find a wider range of applications in the modeling of positive real-world datasets, that is, not only in the area of biology but also in many other areas such as physics, astronomy, engineering, survival analysis, hydrology, economics, and so on.
Inferring such models raises both statistical and computational issues, many of which were solved in recent contributions using variational techniques and convex optimization tools. The PLN model turns out to be a versatile framework, within which a variety of analyses can be performed, including multivariate sample comparison, clustering of sites or samples, dimension reduction (ordination) for visualization purposes, or inferring interaction networks. This paper presents the general PLN framework and illustrates its use on a series a typical experimental datasets. All the models and methods are implemented in the R package PLNmodels, available from cran. r-project. org.
This includes multivariate regression in its simplest form, but also multivariate sample comparison via linear discriminant analysis (LDA), modelbased clustering using mixture models, dimension reduction via principal component analysis (PCA:  Chiquet et al. , 2018) , and network inference  (Chiquet et al. , 2019). All these analyses are implemented in the R package PLNmodels, available from cran. r-project. org. Because it involves a latent layer, the inference of this model raises a series of computational issues, which can be circumvented via a variational approximation  (Blei et al. , 2017). The rest of this section is devoted to a review of existing methods, the precise definition of the PLN model and how it differs from or is similar to other methods. Section 2 provides a series of examples illustrating how the PLN model can be adapted to tackle some specific questions (sample comparison, clustering, dimension reduction and visualization, or network inference) using various extensions summarized in Figure  1.
The PLN model introduced here is yet another instance of the rich family of GLLVM which takes the middle road. Our goal is to develop a generic and versatile umbrella framework under which one can perform various tasks: ordination  (Chiquet et al. , 2018) , classification, group prediction, network inference  (Chiquet et al. , 2019) , etc. Unlike  Hui et al. (2017) , and in line with  Ovaskainen et al. (2017) , we allow the user to carefully constrain the model parameters: each set of constraints leading to a different task. In contrast to HMSC, we do not specify a grand overarching model but rather several easier to estimate variations around a central model. In addition, we want to analyze large datasets, similarly to  Hui et al. (2017) , and thus rely on the efficient inference machinery of variational approximation rather the slower MCMC.
ADAPTING THE PLN FRAMEWORK TO DIFFERENT TASKS We now introduce a series of extensions of the basic PLN model defined in Equation (1). As explained above, the PLN model deals with both abiotic effects through the mean vector µ i Each method will be accompanied with a specific example. LDA will be used to compare the bacterial communities collected in different body sites of dairy cattle. Model-based clustering will be applied to the microbiota of leaves from several oaks and will prove to be able to recover, in a blind way, the tree of origin of each leaf. Dimension reduction using PCA will be applied to the same dataset for visualization purposes and to exhibit which species contribute most to the community structure. The last example is an attempt to reconstruct the interaction network of fish species from the Barents Sea.
Cow Holobionts To illustrate PLN-LDA, we consider the dataset introduced in Mariadassou et al. (  2020 ) which consists of n = 256 bacterial communities sampled on three body sites (nose, mouth, and vagina) of 45 primiparous Prim'Holstein dairy cattle. The cattle comes from two divergent lineages, each sampled twice (1 month and 3 months after first delivery). Communities were sequenced using the hypervariable V3-V4 region of 16S rRNA as marker-gene. Sequences were cleaned and analyzed with DADA2  (Callahan et al. , 2016)  to create p = 1, 077 Amplicon Sequence Variants (ASVs). The aim is to assess whether the communities living in different body sites are different and how they differ. We analyzed this data set by running PLN-LDA with time and lineage as covariates and body site as class.
Parametrization of the Covariance in PLN-Mixture Models When using parametric mixture models like Gaussian mixture models, it is not recommended to consider general covariance matrices k with no special restriction, especially when dealing with a large number of species. Indeed, the total number of parameters to estimate in the model can become prohibitive: in the general case, a PLN-mixture model with K components like in (4) has K × (p + p(p + 1)/2) model parameters, plus the K ×2(n×p) variational parameters. To reduce the computational burden and avoid over-fitting two different, more constrained parameterizations of the covariance matrices of each component are currently implemented in the PLNmodels package (on top of the general form of k ): diagonal covariances: k = diag(d k ) (2Kp parameters), spherical covariances: k = σ 2 k I (K(p + 1) parameters).
The diagonal structure assumes that, given the group membership of a site, all species abundances are independent. The spherical structure further assumes that all species have the same biological variability. In particular, in both parametrization, all observed covariations are caused only by the group structure. For readers familiar with the mclust R package  (Fraley and Raftery, 1999) , which implements Gaussian mixture models with many variants of covariance matrices of each component, the spherical model corresponds to VII (spherical, unequal volume) and the diagonal model to VVI (diagonal, varying volume and shape). Using constrained forms of the covariance matrices enables PLN-mixture to provide a clustering even when the number of sites n remains of the same order, or smaller, than the number of species p.
Oaks Powdery Mildew To illustrate PLN-mixture, we consider the dataset introduced in Jakuschkin et al. (  2016 ) which consists of microbial communities sampled on the surface of n = 116 oak leaves. Communities were sequenced with both the hypervariable V6 region of 16S rRNA as marker-gene for bacteria and the ITS1 as marker-gene for fungi. Sequences were cleaned, clustered at the 97% identity level to create OTUs and only the most abundant ones were kept (see  Jakuschkin et al. , 2016  for details) resulting in a total of p = 114 OTUs (66 bacterial ones and 48 fungal ones). One aim of this experiment is to understand the association between the abundance of the fungal pathogenic species E. alphitoides, responsible for the oak powdery mildew, and the other species. Furthermore, the leaves were collected on three trees with different resistance levels to the pathogen, which we call the tree susceptibility.
We use this example to assess the ability of model-based clustering to recover, without feeding this information to the model, the existence of groups of leaves with different origins. We analyzed this data set by running PLN-mixture for a number of component varying from 1 to 6. We selected the final number of components with a variant of the Integrated Classification Likelihood (ICL:  Biernacki et al. , 2000) , tuned for our own PLN framework. Since the abundances were measured separately for fungi and bacteria, we define a different offset term o ij for each OTU type to take into account the differences in sampling effort and marker genes. Offsets are still computed as the log-total sums of reads, including those of filtered out OTUs, for each OTU type. The results of our analysis with PLN-mixture are displayed in Figure  4A  shows the evolution of the approximated loglikelihood (which is strictly increasing with the number of components, as expected) and the evolution of the ICL criterion, which suggest a model with 4 or 5 components.
Oaks Powdery Mildew We continue the analysis of the oaks data set begun for PLNmixture. As seen before, there is an obvious structure in the data explained by the tree susceptibility. Knowing this, one may be interested in exhibiting a remaining dependence structure that is not induced by the tree susceptibility. In Figure  5A , we represent the biplot for the first two principle components obtained with PLN-PCA, not accounting for the tree susceptibility, after selecting the best possible rank q thank to the ICL criterion. As expected after our clustering study with PLN-mixture, the factorial map exhibits a strong structure where individuals are spread into three groups corresponding to the level of susceptibility of the trees where the leaves were sampled. We also projected the 10 species with the highest contribution to the first two principal components.
Barents Fish We illustrate the use of the PLN-network model on the Barents fish dataset. We focus on the way the inferred network is modified when introducing covariates in the model. To this aim, we fitted the PLN-network model with (a) no covariates, (b) two environmental covariates (temperature and depth), and (c) all covariates (i. e. , the previous two plus the geographical location) using a common λ-grid (with 20 values spaced equally in log scale between λ min = 0. 03 to λ max = 15. 17) for al models. Figure  6  (top right panel) shows that the number of edges increases as the penalty decreases, as expected. It also shows that, for any penalty, the number of edges decreases as (plain lines) the richness of the number of edges increases (c > b > a) and that most edges recovered in the full (c) model are also recovered in the partial models (a, black dotted curve) and (b, blue dotted curve).
Simulation Settings We simulated count data according to a PLN model with the following parameters: number of samples n ∈ {50, 100, 500, 1000, 10000}, number of species p ∈ {20, 200}, number of covariates d ∈ {2, 5, 10}, sampling effort N ∈ {low, medium, high} (calibrated to roughly correspond to total sums of counts per sample of 10 4 , 10 5 , and 10 6 ), and noise level σ 2 ∈ {0. 2, 0. 5, 1, 2}. These parameters cover values typically observed in real datasets and range from very hard (n = 50, p = 200, d = 10, N = low) to ridiculously easy (n = 10, 000, p = 20, d = 2, N = high). For each of the 360 parameter combinations, hereafter referred to as simulation setup, we generated a variance matrix as σ jk = σ 2 ρ |j-k| , with ρ = 0. 2, a design matrix X with 1s in the first column (intercept) and all other entries sampled from a standard Gaussian distribution (we also centered all columns but the first one to avoid interplay between the intercept and the covariates), a regression coefficient matrix with all entries sampled from a centered Gaussian distribution with variance 1/d.
PLNmodels Package All the variants of the PLN model presented in this paper are available as an R/C++ package PLNmodels, distributed on the CRAN CRAN. R-project. org/package=PLNmodels. The package comes with a set of accompanying functions and methods for visualization and diagnostic. It also relies on the user-friendly GLM-like syntax to define all the models, so that users familiar with (generalized) linear models will feel at home. The development version is available on github. com/pln-team/PLNmodels and all models are fully documented as vignettes available from the package website https://pln-team. github. io/PLNmodels/. The Barents and oaks mildew data sets analyzed in sections 2 and 3 are included in the package. The natural competitors to PLNmodels are (i) package Hmsc of  Tikhonov et al.
(2020) , implementing the hierarchical modeling framework of  Ovaskainen et al. (2017) , and (ii) package gllvm of  Niku et al. (2019b) , implementing the generalized latent variable models of  Niku et al. (2017)  and boral of Hui (  2016 ), implementing the latent variable model of  Hui et al. (2015)  in a Bayesian framework. We believe that all these tools are complementary, each having advantages over the others. The different variants available in PLNmodel (Discriminant Analysis; Dimension Reduction; Sparse Network Reconstruction; Mixture models) is indeed an asset, as well as our fast variational algorithm. gllvm benefits from the variety of available distributions for the counts (Poisson, NB, ZIP) but performs neither clustering nor network inference. Finally the Bayesian framework and careful decomposition of the effect adopted in Hmsc gives access to a fine analysis of the model parameters with posterior distributions but is limited to medium size problems.
Such a fine modeling will requires specific developments in each variant, but a first approximation can also be done at less effort by adding observation weights in all PLN models tuned by the estimated probabilities for a count to be a structural zero for a given species in a given site. On the computational aspects, we plan to use optimization tools from machine learning (stochastic gradient descent algorithms and variants, GPU, distributed computing, autodifferentiation) to propose an implementation of our models that allow to fit large datasets. The different PLN variants that we introduced can address problems involving hundreds of species and sites on a routine basis, and can go up to a few thousand species or sites. Using the tools mentioned above would allow us to gain 1 or even 2 orders of magnitude in terms of speed.
3. 1. 4. A First Example As a first illustration of the use of the PLN model for analyzing abundance data, we consider the dataset introduced by Fossheim et al. (2006) (and re-analyzed by Greenacre, 2013; Greenacre and Primicerio, 2014), which consists in the abundances of p = 30 fish species measured in n = 89 sites of the Barents sea between April and May 1997. The species under study are sensitive to environmental drivers (temperature, water depth,. ) but are also related by trophic interactions. A first aim is to study the relative contribution of each type of interaction (abiotic vs. biotic) to the structure of the community.
The campaign primarily relied on measurements of chemical species aboard the Deutsches Zentrum für Luft-und Raumfahrt (DLR) Falcon aircraft based in Miri, Sarawak, Malaysia, which were complemented by ship-and ground-based observations  (Pfeilsticker et al. , 2013). By sampling convective outflows under the influence of high CHBr 3 coastal emission zones, several of the SHIVA flights were particularly well designed to document the impact of deep convection on the bromoform distribution in the upper troposphere. We selected the case study from the SHIVA campaign on the afternoon of 19 November 2011. Figure  1  shows hourly maps of brightness temperature contours measured by the 11 µm channel IR108 aboard the MTSAT-2 satellite (http://database. rish. kyoto-u. ac. jp/arch/ctop/index_ e. html, last access: 15 November 2021;  Nishi et al.
The other convective cell, Obs_Conv2, initiated at 06:00 UTC (14:00 LT) and was located at about 4. 3 • N and 114. 4 • E. It later produced an anvil of convective outflow that developed and moved north-westward similarly to Obs_Conv1. It lasted several hours and started to decay from 10:00 UTC (18:00 LT). The Falcon flight on the afternoon of 19 November 2011 was aimed at sampling the outflow of the Obs_Conv1 system. Krysztofiak et al. (2018)  identified from humidity data and webcam images the times when the aircraft flew within the convective outflow (i. e. in cloudy conditions) and when it was in cloud-free conditions. This information has been used here to show in Fig. 1  (09:00 UTC panel) where the flight sampled cloud-free air or cloudy air. The Falcon aircraft sampled the Obs_Conv1 system multiple times but also flew within the Obs_Conv2 system at ∼ 12.