dataesr/scientific-paragraphs-categorization huggingface dataset:
4] [55]  In spite of the urgent need for efficient sharing data system providing access to trial data from around the world, it will take some time before user-friendly systems that respect the rights of all stakeholders  54 55  are in place. Access to the IPD from older trials may never be available. Even if IPD becomes widely available, we recommend at least the following, relevant to our findings: ► Mandatory systematic review (involving a comprehensive search for eligible trials) before meta-analysis; ► A comprehensive search for published trial results should include multiple bibliographic databases and reports in all languages; ► Central registration of all trials is essential for an appropriately comprehensive search of all trials (published and unpublished) eligible for a systematic review; ► Systematic reviews should include evaluation of individual trial RoB, since biased trials can lead to a biased meta-analysis and use tools adapted to the type of data (ie, aggregated data or IPD).
The authors thank the referee for valuable suggestions. R. C. was partially supported by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) with ERC Grant Agreement nr. 615722 MOTMELSUM, and thanks the Labex CEMPI (ANR-11-LABX-0007-01). I. H. was partially supported by the SFB 878 of the Deutsche Forschungsgemeinschaft. F. L. was partially supported by ANR-13-BS01-0006 (Valcomo), by ANR-15-CE40-0008 (Défigéo), by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement nr. 246903 NMNAG and by the Institut Universitaire de France. M. R. was partially supported by ANR-15-CE40-0008 (Défigéo).
New standards have been recently defined and implemented enabling a reliable exchange of pharmacometric models across software tools, and facilitating collaborative drug and disease modeling and simulation (M&S) activities. Among these, the Standard Output (SO) has been proposed as the tool-independent exchange and storage format for typical M&S results. The SO integration within the Drug Disease Model Resource (DDMoRe) interoperability framework (IOF) has already shown its potential for enabling effective data flow across modeling tasks and facilitating information retrieval.
DEVELOPING A STANDARD OUTPUT FORMAT TO EFFECTIVELY SUPPORT PHARMACOMETRIC WORKFLOW INTERACTIONS 2] [3] [4] [5]  There is generally a lack of consistency across these tools in model definition, data format, numerical methods, and output of results. This limits the potential reuse of models and methods across tools by the analyst, who would then need to understand varying model and tool-specific requirements to correctly translate and recode models and data. This lack of standard languages and interoperability across tools also significantly limits the integration of tools into seamless pharmacometric workflow as well as the sharing of existing knowledge through shared code. 7] [8] [9]  Among these, the SO has been proposed as the tool-independent format for storing output typically produced from pharmacometric M&S tasks and workflows to guide the model building through assessment of the goodness-of-fit between the model and the dataset and the appropriateness of the underlying model assumptions, as well as to allow the reproducibility of the results.
POPULATING SO The SO has been designed to allow flexibility in populating its elements. Sections or elements can be omitted if these are not available or not appropriate for a given target software. In addition, elements can be populated through the defined XML structure (e. g. , as inline data or by referring to external files). Unique data types were defined for each SO element reusing some of the structures in the PharmML schema (e. g. , matrices and datasets) and distribution specifications from ProbOnto (e. g. , posterior distributions). 8, 9 r data stored in a tabular form, the columnType attribute was defined to support column descriptors (e. g. , individual parameter, population parameter, and random effect) that improve usability and interpretation of stored values. As an example, the attribute of a column relative to a variability parameter would be specified to clarify whether this reports a variance, standard deviation, covariance, or correlation measure.
Numerical experiments For the following numerical experiments, the code and data are available at https://github. com/nbrosse/normalizingconstant. We first experiment our algorithm to compute the logarithm of the normalizing constant of a multivariate Gaussian distribution in dimension d ∈ {10, 25, 50}, of mean 0 and inverse covariance matrix diag(2, 1 ⊗(d-1) ). We set = µ = 0. 1. The number of phases M of the algorithm and the variances σ 2 i M -1 i=0 are chosen accordingly to the formulas  (34)  and  (36). For each phase of the algorithm, the step size γ i is set equal to 10 -2 (m i + L i ) -1 , the burn-in period N i to 10 4 and the number of samples n i to 10 5 where m i , L i are defined in  (13). We carry out 10 independent runs of the algorithm and compute the boxplots in Figure  1.
For each phase, the step size γ i is set equal to 10 -2 (κ i σ 2 i m i )/(dL 2 i ), the burn-in period N i to 10 3 (κ i γ i ) -1 and the number of samples n i to 10 where m i , L i , κ i are defined in (  13 ) and  (14). The experiments are repeated 10 times and the boxplots for each model M are plotted in Figure  2. Note that for this Gaussian model, the log evidence is known and displayed by the red points in Figure  2. With the same parameters for the algorithm, we run 10 independent runs at each phase to measure the variability of each estimator πi (g i ) defined in  (9). The result is plotted in Figure  3  for the model M 1. The last estimator πM-1 (g M -1 ) is much higher, which underlines the specificity of the last phase in the algorithm. The logistic regression is performed on the Pima Indians dataset 1.
4 and 3. 6 in the paper. The first author is also indebted to Denis-Charles Cisinski and Gijs Heuts for discussions which improved his understanding of shuffles. We also like to thank our home universities and NWO for supporting our mutual visits. Finally, we would like to thank the referees for their careful comments on an earlier version of this paper. The first author is partially funded by the ANR grants ANR-13-BS02-0005-02 CATHRE and ANR-16-CE40-0003 ChroK.
Perspectives In future works, the method will be applied to extreme hydrological flow problems (e. g. [13, 30] ). In particular, we plan to apply the approach to the framework of Ensemble Kalman filters (EnKF)  [11]  for large datasets. The EnKF is a recursive Bayesian estimation technique for dynamical models of the form X (k+1) = HX (k) + L (k) , Y (k+1) = BX (k+1) + E (k) , where X (0) , L (k) , and E (k) are independent Gaussian vectors. Note that the equation above is the same as  (1). To estimate the posterior distribution of X (k+1) , the Kalman filter requires the inversion of the covariance matrix C Y (k+1) at each iteration of the discrete dynamical system. However, in the EnKF, C X (k+1) is estimated using a Monte-Carlo estimator with a sample size that can be much lower than the total number of observations n.
As a consequence, the covariance of the forward state X (k+1) is low-rank and we showed in this paper that a low number of projections of the observations are enough to recover the mutual information between the estimated distribution of the state X (k+1) and the observations Y (k+1). Additional challenges arise when the datasets are high dimensional. Considering the mutual information based technique, the problem could be first tackled using high performance computing. Given that we a priori know an upper bound on the number of projections that is already low, we only need an efficient matrix product computation (e. g. [10] ) to implement the algorithm from  [1]. Further developments are required to appropriately use these approaches in a streaming environment. or equivalently, log f X (x | Y = y) = log f Y (y | X = x) + log f X (x) + k 0 = - 1 2 (y -Bx -m E ) T C -1 E (y -Bx -m E ) + (x -m X ) T C -1 X (x -m X ) + k 1 , where k 0 and k 1 are constants.
For γ 1 varying from 0 to 1 by step of 0. 05 (the value 0 is exluded for DG models to save time), the true intensity λ is estimated by Monte-Carlo methods. For each set of parameters m realizations of the model are generated on the square [-2R, 1 + 2R] 2 and then clipped to the unit square. That strategy is detailed and justified by  Baddeley and Nair (2012). Specifically, the number of points in each realization is averaged to obtain the estimated intensity and its standard error. The simulation results for the Strauss models with β = 50 or 100 were obtained by  Baddeley and Nair (2012) , where m = 10000 realizations were generated and the exact simulation algorithm was used, implemented in the R function rStrauss of the spatstat package (see  Baddeley et al. (2015) ). For the Strauss models with β = 200, SHC models, PS and PSHC models, we generate m = 1000 replications and use the rmh function in the spatstat package which implements a Metropolis-Hastings algorithm.
Even if we use 10 6 iterations of the algorithm, the results may be slightly biased. For the DG models, the R package spatstat provides an exact simulation algorithm (function rDiggleGraton) and for such models we generate 10000 replications when β = 200 and R = 0. 025, 0. 05 and when β = 50 and R = 0. 15. We used 1000 replications when β = 200 and R = 0. 075 to save time. • All results can be found in Figures  3, 4  and 5. Plots provide the same information: we depict intensity approximation λ based on different methods in terms of γ 1. The dashed curve represents the Poisson-saddlepoint approximation imsart-ejs ver. 2014/10/16 file: appIntDPP. tex date: December 1, 2017 proposed by  Baddeley and Nair (2012)  and detailed in Section 3. 1. The solid curve is the DPP approximation we propose in this paper and is given by (3.
Acknowledgements The authors are sincerely grateful to  Adrian Baddeley  and  Gopalan Nair  for sharing the Monte-carlo replications produced in Baddeley and Nair (2012) we used to compare the Poisson-saddlepoint approximation and the DPP approximation (Figure 3 (a),(b) and (d)). The research of J-F. Coeurjolly is supported by the  Natural Sciences and Engineering Research Council of Canada.
Some examples Using the FreeFem++ software  [20] , we test some numerical examples. In all the examples below, we work on the square Ω = [0, 1] × [0, 1] in R 2 and use the discretization size 36 × 36 × 9 for the spatial-time variable. We test the examples for costs of the form c(x, y) = inf ξ    1 0 L(ξ(t), ξ(t)dt : ξ ∈ Lip([0, 1]; Ω), ξ(0) = x, ξ(1) = y    , with L(x, v) = k(x)|v| 2 , k ∈ C(Ω), k(x) > 0 for all x ∈ Ω, v ∈ R 2. For this cost, the last projection in the Step 2 (the projection on K x ) is converted to a problem on R and the latter is computed easily by the bisection method. The optimal active submeasures and the optimal displacement are given in Figure  1. Timestep 0 and timestep 9 (ρ 0 and ρ 1 ) are optimal active submeasures of the source and the target, respectively.
(2016)  detected shifts in the optimal value of a trait evolving as an OU. Uyeda and Harmon (2014)  and  Bastide et al. (2016)  detect all shifts for a given number of shifts and use either rjMCMC or penalized likelihood to select the number of shifts. By contrast,  Ingram and Mahler (2013)  uses a stepwise procedure, based on AIC, to detect shifts sequentially, stopping when adding a shift does not improve the criteria anymore. Extensions from univariate to multivariate shifts are more recent. It should be noted that all methods assume that shifts affect all traits simultaneously. Given known shift locations and a multivariate OU process,  Bartoszek et al. (2012)  was the first to develop a likelihood-based method (package mvSLOUCH) to estimate both matrices of multivariate evolutionary rates and selection strengths.
Clavel et al. (2015)  soon followed with mvmorph, a comprehensive package covering a wide range of multivariate processes. Detection of shifts in multivariate traits is more involved and both  Ingram and Mahler (2013)  and  Khabbazian et al. (2016)  make the simplifying assumption that all traits are independent, conditional on their shared shifts. Ingram and Mahler (2013)  then proceed with the same stepwise procedure as in the univariate case whereas  Khabbazian et al. (2016)  uses a lasso-regression to detect the shifts and a phylogenetic BIC (pBIC) criterion to select the number of shifts.
Our method is fast and can handle missing data. It also proved accurate in a large scale simulation study and was able to find back known shift locations in re-analysis of public datasets. Finally, we show that the standard practice of decorrelating traits using phylogenetic principal component analysis (pPCA) before using a method designed for independent traits can be misleading in the presence of shifts. The article is organized as followed. We present the model and inference procedure in Section Model, the theoretical bias of pPCA in the presence of shifts in Section pPCA and Shifts, the simulation study in Section Simulations Studies, the re-analysis of the New World Monkeys and Greater Antillean Anolis lizards datasets in Section Examples and discuss the results and limitations of our method in Section Discussion.
Note that, when the root is taken random, everything happens as if we added a fictive branch above the root with length ρ (α). The length of this branch increases when α goes to zero. We emphasize that only the distribution of the observed traits Y is preserved and not the distribution of the complete dataset X. As a consequence, ancestral traits at internal nodes cannot be directly inferred using this representation. Still, the equivalence recasts inference of R and W(α)∆ in the scOU model into inference of the same parameters in a much simpler BM model, albeit on a tree with rescaled branch lengths i (α). Note that the rescaling depends on α, which needs to be inferred separately. See the discussion (Section Interpretation Issues) for further analysis of this re-scaling.
-The EM algorithm is known to be very sensitive to the initialization. Following  Bastide et al. (2016) , we take advantage of the linear formulation (2) to initialize the shifts position using a lasso penalization  (Tibshirani 1996). This initialization method is similar to the procedure used in 1ou  (Khabbazian et al. 2016). See Appendix EM Inference for more details. Missing Data. -EM was originally designed to handle missing data. As a consequence, the algorithm described above also applies when some traits are unobserved for some taxa. Indeed, the conditional distribution of the missing traits given the observed ones can be derived in the same way as in the E step. However, missing data break down the factorized structure of the dataset and some computational tricks are needed to handle the missing data efficiently (see Appendix EM Inference).
Implementation We implemented the method presented above in the PhylogeneticEM R package (R Core Team 2017), available on the Comprehensive R Archive Network (CRAN). A thorough documentation of its functions, along with a brief tutorial, is available from the GitHub repository of the project (pbastide. github. io/PhylogeneticEM). Thanks to a comprehensive suite of unitary tests, that cover approximately 79% of the code (codecov. io/gh/pbastide/PhylogeneticEM), and that are run automatically on an independent Ubuntu server using the continuous integration tool Travis CI (travis-ci. org), the package was made as robust as possible. The computationally intensive parts of the analysis, such that the upward-downward algorithm of the M step, have been coded in C++ to improve performance (see Section Simulations Studies for a study of the computation times needed to solve problems of typical size).
Simulations Studies Experimental Design General Setting. -We studied the performance or our method using a "star-like" experimental design, as opposed to a full-factorial design. We first considered a base scenario, corresponding to a base parameter set, and then varied each parameter in turn to assess its impact as in  Khabbazian et al. (2016). The base scenario was chosen to be only moderately difficult, so that our method would find shifts most but not all of the time. For the base scenario, we generated one 160-taxon tree according to a pure birth process, using the R package TreeSim Finally, we added three shifts on this phylogeny, with fixed positions (see Figure  2 ). Shift amplitudes were calibrated so that the means at the tips differ by about 1 standard deviation, which constitute a reasonable shift signal  (Khabbazian et al.
2016). Each configuration was replicated 100 times. We then used both our PhylogeneticEM and 1ou package  (Khabbazian et al. 2016)  to study the simulated data. We excluded SURFACE  (Ingram and Mahler 2013)  from the comparison at is (i) quite slow, (ii) assumes the same evolutionary model as 1ou and (iii) was found to achieve worse accuracy than 1ou  (Khabbazian et al. 2016). We used default setting for both methods. For PhylogeneticEM this implies an inference on an automatically chosen grid with 10 α values, on a log scale, Low scaling values correspond to smaller, harder to detect, shifts and high values to larger and easier to detect shifts. -1. 01 1. 03  Selection Strength. -When exploring parameters not related to the shifts, we considered a base number of 3 shifts and a base scaling factor of 1.
To be fair to both methods, we selected parameter ranges leading to similar Kullback divergences, to achieve similar levels of mis-specifications. However, both deviations produce datasets with groups that are also theoretically easier to discriminate compared to model A (see Figure  3 ). Indeed, we can quantify the difficulty of a dataset in terms of group separation by the Mahalanobis distance between the observed data and their expected mean, (phylogenetically) estimated in the absence of shifts:   4 ). Unequal selection strengths (s > 1) and trait correlations (r d > 0) both increase group separation and make it easier to detect shifts. D = Y vec -(1 T Σ d 1) -1 1 T Σ d Y vec 2 Σ d -1 -(np -N NA ) (4 missing data varied from 5% to 50%. We use this measure rather than the classical precision/sensitivity graphs as only the clustering can be recovered unambiguously (see Section Model).
1ou scales less efficiently: it is faster for very small trees (32 taxa) but median running times go up to 20 hours for the large 256-taxon tree. Those long running times were unexpected and higher than the ones reported in  Khabbazian et al. (2016). This discrepancy is partly due to the maximum number of shifts allowed, which strongly impacts the running time of 1ou. (r d = 0. 4) with three "easy" shifts, it becomes harmful and degrades the performances of 1ou when the number or magnitude of the shifts increases (Fig. 9b ). As expected, PhylogeneticEM is unaffected by the pPCA preprocessing, up to numerical issues. We used PhylogeneticEM to re-analyse two publicly available datasets.
New World Monkeys We first considered the evolution of brain shape in New World Monkeys studied by  Aristide et al. (2016). The dataset consists of 49 species on a time-calibrated maximum-likelihood tree. The traits under study are the first two principal components (PC1, PC2) resulting from a PCA on 399 landmarks describing brain shape. We ran PhylogeneticEM on a grid of 30 values for the α parameter. To make this parameter easily interpretable, we report the phylogenetic half-life t 1/2 = ln(2)/α  (Hansen 1997) , expressed in percentage of total tree height. Here, t 1/2 took values between 0. 46 % and 277. 26 %. We allowed for a maximum of 20 shifts. The inference took 17. 56 minutes, parallelized on 5 cores. The model selection criterion suggests an optimal value of K = 4 shifts (Fig. 10 , inset graph).
Lizards We then considered the dataset from  Mahler et al. (2013) , which consists in 100 lizard species on a time-calibrated maximum likelihood tree and 11 morphological traits. We chose this example because of the large number of traits and the high correlation between traits, as all traits are highly correlated (0. 82 < ρ < 0. 97) with snout-to-vent length (SVL). To deal with the correlation between traits,  Mahler et al. (2010 Mahler et al. ( , 2013) )  first performed a phylogenetic regression of all the traits against SVL, retrieved the residuals and then applied a phylogenetic PCA on SVL and the previous residuals, from which they used the first four components (pPC1 to pPC4) for their shift analysis. We first explored how the number of pPCs used can impact the shift detection. Hence we ran PhylogeneticEM 11 times, including 1 to 11 pPCs in the input dataset.
A natural suspect for such instability is the extreme correlation between some traits (0. 996 for tibia and metatarsal lengths), which results in bad conditioning of several matrices that must be inverted. To circumvent this problem, we used the two pseudo-orthogonalization strategies described above, running PhylogeneticEM on the SVL plus residuals dataset, and on the 11 pPCs, with the same parameters as above. Note that all these transformations use a rotation matrix, so that the likelihood and the least squares of the original or of any of the two transformed datasets are the same. Hence, the objective function, as well as the model selection criterion, should remain unchanged. Still, slight differences were found between the maximized likelihood for each pseudo-orthogonalized datasets.
account for the correlation between traits. Still, high correlations between traits may raise strong numerical issues, so PCA can be used as a pseudo-orthogonalization of traits, as well as any other linear distance-preserving transformation that would reduce the correlation between them. This does not dispense of considering the correlation between the transformed traits in the model. The other interest of PCA is to reduce the dimension of the data, which may be desirable when dealing with a large number of traits, such as the original dataset from  Aristide et al. (2016). Since PCA does not correct for the right correlation, we have no clue whether or not the dimension reduction performed by PCA is relevant for shift detection, or if it may remove precisely the direction along which the shifts occur.
It worked quite well, but was much slower. This approach could be extended to the multivariate setting, although with impaired computational burden. Note also that some shifts configuration that are not identifiable in the absence of fossil data become distinguishable with the addition of fossil data. This affects our model selection criterion, which relies on the number of distinct identifiable solutions. Computing this number on a non-ultrametric tree for an OU remains an open problem, and is probaly highly dependent on the topology of the tree.
E step The aim of the E step is to compute the moments of the completed dataset given the observed traits at the tips, namely: E j = E X j Y d , V j = Var X j Y d , C j,pa(j) = Cov X j ; X pa(j) Y d (8) where we dropped the dependency in θ (h) for the sake of legibility, but all these moments are indeed taken with the laws given by the current parameters. We do so thanks to an upward-downward recursion on the tree, as described below. This algorithm can apply to a broad classes of Gaussian processes, provided that the moments of the traits at a child node are of the form: ∀j ∈ 2 , m + n ,      E X j X pa(j) = m j (X pa(j) ) = Q j X pa(j) + r j Var X j X pa(j) = Σ j For a BM, we get Q j = I p , r j = ∆ j and Σ j = j R. A multivariate OU could also be handled, with: Q j = e -A j , r j = (I p -e -A j )β j and Σ j = Γ -e -A j Γe -A T j.