David Danan

Virginie Le Rolle
email: virginie.lerolle@univ-rennes1.fr

Arnaud Hubert

Elena Galli

Anne Bernard

Erwan Donal

Alfredo I Hernández
email: alfredo.hernandez@inserm.fr

Validation of a Finite Element Method framework for cardiac mechanics applications

Keywords: Cardiac mechanics, Finite element method, hyperelasticity, Getfem, Benchmark validation

Modeling cardiac mechanics is a particularly challenging task, mainly because of the poor understanding of the underlying physiology, the lack of observability and the complexity of the mechanical properties of myocardial tissues. The choice of cardiac mechanic solvers, especially, implies several difficulties, notably due to the potential instability arising from the nonlinearities inherent to the large deformation framework. Furthermore, the verification of the obtained simulations is a difficult task because there is no analytic solutions for these kinds of problems. Hence, the objective of this work is to provide a quantitative verification of a cardiac mechanics implementation based on two published benchmark problems. The first problem consists in deforming a bar whereas the second problem concerns the inflation of a truncated ellipsoid-shaped ventricle, both in the steady state case. Simulations were obtained by using the finite element software GETFEM++. Results were compared to the consensus solution published by 11 groups and the proposed solutions were indistinguishable. The validation of the proposed mechanical model implementation is an important step toward the proposition of a global model of cardiac electro-mechanical activity.

INTRODUCTION

Modeling of cardiac activity is an active field of research and the state-of-art is wide and very active. Recent projects have provided interesting preliminary results towards the use of individualized, computer-based, human heart models. [START_REF] Talbot | Towards an interactive electromechanical model of the heart[END_REF][START_REF] De Vecchi | Patient specific fluid-structure ventricular modelling for integrated cardiac care[END_REF][START_REF] Mccormick | Simulating left ventricular fluid solid mechanics through the cardiac cycle under lvad support[END_REF] Most of approaches proposed descriptions, at many different levels of detail, of the cardiac electrical activity, [START_REF] Luo | A dynamic model of the cardiac ventricular action potential. ii. afterdepolarizations, triggered activity, and potentiation[END_REF][START_REF] Werner | Electrical excitation propagation in the human heart[END_REF] the excitation-contraction coupling, [START_REF] Kerckhoffs | Homogeneity of cardiac contraction despite physiological asynchrony of depolarization: a model study[END_REF][START_REF] Bovendeerd | Dependence of intramyocardial pressure and coronary flow on ventricular loading and contractility: a model study[END_REF] the mechanical activity [START_REF] May-Newman | Homogenization modeling for the mechanics of perfused myocardium[END_REF][START_REF] Vetter | Three-dimensional stress and strain in passive rabbit left ventricle: a model study[END_REF][START_REF] Kerckhoffs | Timing of depolarization and contraction in the paced canine left ventricle: model and experiment[END_REF][START_REF] Sermesant | An electromechanical model of the heart for image analysis and simulation[END_REF] and the mechano-hydraulic coupling. [START_REF] Verdonck | Fluid-structure interaction modelling of left ventricular filling[END_REF] However, translation of these models into clinical practice remains difficult and new methodological approaches are still needed in order to bring these model-based approaches to the clinical field. The global objective of our project is to propose a model-based method for the analysis of 3D regional cardiac motion and deformation, which will be applicable to Cardiac Resynchronization Therapy (CRT) in the pre-and post-operatory phases.

In order to simulate accurately myocardial deformations, the proposed model should integrate an appropriate description of cardiac mechanics. The ventricular mechanical activity is usually described as a function of its active and passive properties. Active properties are the consequence of the shortening and lengthening of sarcomeres, which are the elementary mechanical contractile elements of myocytes. This mechanical activity is under the influence of an electrical activity, since the variation of calcium concentration during the action potential allows the development of force. Passive properties are mainly related to fiber structure and orientation, collagen properties and metabolic conditions (such as hypoxia or ischemia). This paper mainly focus on the simulation of passive mechanical properties of cardiac tissues.

Passive myocardial properties could be described through specific mechanical constitutive laws. Most of these mechanical laws are hyperelastic, incompressible and anisotropic. [START_REF] Humphrey | Determination of a constitutive relation for passive myocardium: Ii. parameter estimation[END_REF][START_REF] Novak | Regional mechanical properties of passive myocardium[END_REF][START_REF] Hunter | Myocardial constitutive laws for continuum mechanics models of the heart[END_REF] The majority of them have been determined using uniaxial (Mirsky 1976) or biaxial tension tests. [START_REF] Novak | Regional mechanical properties of passive myocardium[END_REF] An empiric law based on the description of sarcomere dynamics has been proposed. [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF] The simulation of these models is often based on finite element methods (FEMs). [START_REF] May-Newman | Homogenization modeling for the mechanics of perfused myocardium[END_REF][START_REF] Vetter | Three-dimensional stress and strain in passive rabbit left ventricle: a model study[END_REF][START_REF] Kerckhoffs | Timing of depolarization and contraction in the paced canine left ventricle: model and experiment[END_REF][START_REF] Sermesant | An electromechanical model of the heart for image analysis and simulation[END_REF] Although this kind of formulation requires higher computational resources than other approaches, it allows a rather detailed description of the myocardium dynamics.

The implementation of FEM-based models implies several difficulties such as the definition of the constitutive law and boundary conditions, the construction of the mesh, the choice of simulation parameters that ensure stability and the validation of simulated deformations. One particular issue is the lack of analytic solutions to the problems. In this context, the use of a published benchmark 17 is a consistent way to validate and verify the accuracy of simulations. Recently, a consortium of several teams 17 provided a benchmark of three problems that require the ability to solve accurately actual features that one may come across in cardiac mechanics: pressure distribution that depends on the deformed configuration , anisotropic material properties or even active contraction of a ventricle. Such a benchmark was solved by 11 different groups and they were able to reach a consensus. The aim of this paper is to provide a validation of the first two benchmark problems by using the finite element software GETFEM++.

The rest of the paper is organized as follows. In Section 2 we introduce some preliminary material to be used in the rest of the paper. In Section 3, we formulate the mathematical model of the hyperelastic problem and provide a description of the benchmark problems we solved. Finally, in Section 4, we present numerical results and confront them with the strain given by the benchmark.

FINITE STRAIN ELASTICITY THEORY AND NOTATIONS

In order to introduce the concepts referred to thereafter, we recall briefly the main aspects of the theory of finite strain elasticity.

Let Ω be the reference configuration and let ϕ be a deformation of Ω that preserves the orientation, that is to say:

ϕ : Ω → Ω ϕ ⊂ R 3 , with det∇ϕ > 0 on Ω. (1)
Now, denote by x ϕ = ϕ(x) the position of a particle x submitted to the deformation ϕ; the displacement field is then defined by

u(x) = ϕ(x) -x, (2)
and we deduce the expression of the gradient of deformation

F F = ∇ϕ = I + ∇u. (3)
To estimate the gap between a deformed configuration and a so-called rigid body motion, the Green-Lagrange strain tensor E is used

E = F T F -I 2 = 1 2 ∇u T ∇u + 1 2 (∇u + ∇u T ). ( 4
)
Note that in the case of linear elasticity, supposedly too restrictive in the context of cardiac mechanics, the second order term 1 2 ∇u T ∇u is neglected.

MECHANICAL PROBLEMS

3.1 Mathematical model

Hyperelastic constitutive law

In order to take into account fiber orientation induced by the anatomical organization of cardiomyocytes, a transversally isotropic constitutive law was chosen. Its nonlinear strain energy function, the so-called Guccione's law, is given by

W (E) = C 2 (e Q -1) (5)
with

Q = b f E 2 11 + b t (E 2 22 + E 2 33 + E 2 23 + E 2 32 ) + b f s (E 2 12 + E 2 21 + E 2 13 + E 2 31 ). ( 6
)
where E ij are components of the Green-Lagrange strain tensor E, defined in (4). Note that, in this configuration, the fibers are in the e 1 -direction.

Strong formulation

We consider a hyperelastic body that occupies the bounded domain Ω, with Γ its boundary, and we denote by ν the unit outward normal on Γ. Since the body is clamped on Γ 1 , the displacement field vanishes there. A volume force of density f 0 acts in Ω, surface tractions of density f 2 act on Γ 2 .

The classical formulation of an incompressible solid in steady-state equilibrium can be stated as Problem P. Find a stress field

Π : Ω × R + → S d and a displacement field u : Ω × R + → R d such that Π = ∂ F Ŵ (F )
in Ω, ( 7)

det(F) = 1 in Ω, (8)
Div Π + f 0 = 0 in Ω, (9)
u = 0 on Γ 1 , (10)
Πν = f 2 on Γ 2 . ( 11
)
Now we shortly describe the physical meaning of relations ( 7)- (11). Equation ( 7) represents the generic form of a hyperelastic constitutive law where Ŵ (F ) = W (E) denote the strain energy function defined at the beginning of this section and Π, the first Piola-Kirchhoff stress tensor. Note that by using advantageously the properties of the second Piola-Kirchhoff stress tensor S = F -1 Π, one can prove that

Π = F ∂ E W (E). (12)
In fact, it is a much easier to work with formulation as the strain energy function used in (5) depends explicitly on E. The expression (8) describes the incompressibility condition, which means the deformation is isochoric throught the static process. Relation ( 9) is the equilibrium equation, where Div represents the divergence operator for tensor-valued function and f 0 , the density of applied volume forces. Condition (10) is the displacement boundary condition. Next, the traction boundary condition (11) states that the stress vector Πν is given on part Γ 2 of the boundary, and is equal to the boundary force of density f 2 .

In fact, one has to take into account this very condition with particular care as the direction of the pressure boundary condition changes with the deformed surface orientation, and its magnitude scales with the deformed area. Therefore, if we denote by p the blood pressure, this so-called follower load expression is given by

f 2 = -p det(F )F -T ν. ( 13
)
From now on, as prescribed in the benchmark, we assume that no volume force acts on the body.

Problems description

In what follows, we present the first two problems, as they are described in the benchmark. 17

Problem 1

The first problem is the deformation of a bar in static.

• Ω = [0, 10 mm] × [0, 1 mm] × [0, 1 mm].

• Constitutive parameters:

Transversely isotropic, C = 2 kPa, b f = 8, b t = 2, b f s = 4.
• Fiber direction: Constant along the long axis, i.e. (1, 0, 0).

• Dirichlet boundary conditions: The left face (x = 0) is fixed in all directions.

• Pressure boundary conditions: A pressure of 0.004 kPa is applied to the bottom face (z = 0).

• The undeformed geometry is defined by the volume between:

-The endocardial surface r s = 7 mm, r l = 17 mm, u ∈ [-π, -arccos 5 17 ], v ∈ [-π, π], -The epicardial surface r s = 10 mm, r l = 20 mm, u ∈ [-π, -arccos 5 20 ], v ∈ [-π, π] -The base plane z = 5 mm which is implicitly defined by the ranges for u.

• Constitutive parameters:

Isotropic, C = 10 kPa, b f = b t = b f s = 1.
• Dirichlet boundary conditions: The base plane (z = 5 mm) is fixed in all directions.

• Pressure boundary conditions: A pressure of 10 kPa is applied to the endocardium.

Numerical methods and stability issues

As stated explicitly in the benchmark, there is no requirement regarding the computational method used to solve the problem. Here, a generalized Newton's method is used to deal with the nonlinearities arising from the hyperelastic constitutive law and the Lagrangian method was considered to enforce the incompressibility condition.

Such an approach leads to the following strain energy function

W Lag = W (E) -p(det(F ) -1) with det(F ) = 1, (15)
with p the hydrostatic pressure. Nevertheless, to anticipate stability issues reported by several groups, [START_REF] Land | Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour[END_REF] we still followed the suggestions of another article [START_REF] Land | Improving the stability of cardiac mechanical simulations[END_REF] and implemented the stabilized version of a strain energy function independent of changes in volume. In order to do so, we introduce the isochoric component of the deformation gradient F = det(F ) -1/3 F ,

which by definition is such that det( F ) = 1 and we deduce that

Ē = F T F -I 2 = det(F ) -2/3 E + (det(F ) -2/3 -1) I 2 . (17)
Then, the expression of the new strain energy is

W stab iso = W ( Ē) -p(det(F ) -1) + κ 2 (det(F ) -1) 2 with det(F ) = 1. (18)
The last term of ( 18) plays the role of a higher order incompressibility penalty term and is supposed to enforce more accurately the incompressibility constraint and avoid to some extent the pathological cases [START_REF] Peyrault | Orientation preservation and newton-raphson convergence in the case of an hyperelastic sphere subjected to an hydrostatic pressure[END_REF] (non-physical and unstable) where det(F ) < 0.

Evaluation of the simulation

In order to evaluate the accuracy of our results, we recall here the definition of the norms used thereafter. Let u be a vector of size N , we have

u 2 = N i=1 u(i) 2 (19)
u ∞ = max i |u(i)|. (20)
For the second problem, we also provide the positions of the endocardial and epicardial apex, before the deformation, and the position interval obtained by the different teams, after the inflation of the ventricle. It will be used later on as one of the criteria to assess the validity of our solution.

Table 1. Position of the apex before and after the deformation.

Endocardial apex Epicardial apex

Reference position (0, 0, -17 mm) (0, 0, -20 mm) Deformed configuration (-z) 25 mm (one group) 26 -27 mm 27.75 -28.75 mm

RESULTS

The aim of this Section is to provide numerical simulations in order to validate the results. Note that the different numerical methods have been implemented in a code which is based on Finite Element Library in C++ under the GNU Public license: GEneric Tools for Finite Elements Methods (GETFEM++) developed by Julien Pommier and Yves Renard. For more details, we refer to http://download.gna.org/getfem/html/homepage/.

Problem 1

For the computation, we used the Augmented Lagrangian method, which require the strain energy provided in (15). We computed the solution with hexahedral elements Q2 for the displacement and Q1 for the multipliers which correspond to a problem of 20197 degrees of freedom (with a single load step) in 1012 CPU time (in seconds) on a computer with Intel Quad core processors (2.8 GHz).

In Figure 1 we plot the deformed bar; just like in the benchmark the reference solution and the computed solution are indistinguishable. To evaluate and analyse our result, we resort to the difference between our solution and the reference solution obtained for the exact same mesh (Table 2). It shows the solution is accurate enough.

Problem 2

For the computation, we used the isochoric stabilized version, which require the strain energy provided in (18). We computed the solution with tetrahedral elements P 2 for the displacement and P 1 for the multipliers which correspond to a problem of 41984 degrees of freedom (with 25 load steps) in 1720 CPU time (in seconds) on a computer with Intel Quad core processors (2.8 GHz).

In Figure 4 we plot the inflated ventricle; just like in the benchmark the reference solution and the computed solution are still hardly distinguishable which is why we provide the difference in norm (Table 3). While we can not deny the results are not as good as in the first problem, the discrepancy among the teams also seems to be enhanced in this very case. Finally, we give in Table 4 the z-coordinate of the endocardial apex and epicardial apex after the inflation. 1 that the obtained value are in an acceptable range in comparison with the other groups who took part in this benchmark. [START_REF] Land | Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour[END_REF]

CONCLUSION

This study presented a validation of cardiac mechanics implementation, based on two benchmark problems. The results, obtained by using the finite element software GETFEM++, were compared to consensus solutions, provided by a consortium of 11 groups. In terms of three-dimensional deformation as visualized for the 2 problems, the proposed solutions are indistinguishable and the computed errors show only slight differences. The validation of the 2 problems implementation is an important step in the process of proposing a global model of cardiac mechanics. Future works will focus on the proposal of a complete model integrating electro-mechanical activations.

Figure 1 . 1 3. 2 . 2 2

 11222 Figure 1. Reference configuration for the problem 1

Figure 2 .

 2 Figure 2. Reference configuration for the problem 2

Figure 3 .

 3 Figure 3. Superposition of the deformed configuration computed for the problem 1 and the reference solution given by the benchmark

u ref -u getf em 2 u

 2 refu getf em ∞ 7 10 -3 mm 4.8 10 -2 mm

Figure 4 .

 4 Figure 4. Superposition of the deformed configuration computed for the problem 2 and the reference solution given by the benchmark

u ref -u getf em 2 u

 2 refu getf em ∞ 4.29 mm 0.43 mm

Table 4 .

 4 Z-Position of the endocardial apex and epicardial given by Getfem in absolute value.

	Endocardial apex Epicardial apex
	26.59 mm	28.26 mm
	One can easily confirm in Table

ACKNOWLEDGMENTS

This work is supported by the French National Research Agency (ANR MAESTRo project).


Alice Abbondanza

Milica Drapsin

Michaela Rusková

Sylvie Dumas

Véronique Bernard

Alena Sumova

Helena Janickova

Nicotinic activation of NPY/AgRP neurons of the arcuate nucleus and its role in stress and feeding

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

In the brain, Neuropeptide Y (NPY) modulates anxiety via amygdalar circuits or regulates energy homeostasis and stress response via hypothalamic nuclei, including the arcuate nucleus (ARC). Within the ARC, the orexogenic NPY/Agouti-related peptide (NPY/AgRP) neurons are responsive to variation of the metabolic state, hormones, GABAergic transmission, and nicotine. Interestingly, nicotine consumption is associated with a reduction of food intake and an increase of the hypothalamo-pituitary-adrenal stress response, whereas NPY release to the hypothalamic paraventricular nucleus (PVN) stimulates feeding behavior. We hypothesize that cholinergic activation of NPY-expressing neurons of the ARC is implicated in the modulation of NPY release. To test this hypothesis, we crossed beta2-flox/flox mice with NPY-IRES-Cre to delete beta2-containing nicotinic acetylcholine receptors (nAChRs) in all the NPY-expressing cells. We used in situ hybridization for NPY and beta2-nAChRs to confirm the presence of nAChRs on NPY/AgRP neurons of the ARC in control animals, and the beta2 deletion in the floxed mutants. Then, we induced stress by restraining and we measured its effect in selected behavioral tasks. Given a dysmorphic association between stress exposure and alterations in feeding, we decided to test male and female mice, both during light and dark regimen. In females, the deletion of beta2-nAChRs increases the susceptibility to daily stressors and impairs the preference for social interaction, while in males it influences metabolism and body weight. The phenotype observed suggests a more complex interplay between the cholinergic system and NPY release, with possible recruitment of other types of nAChRs and compensatory mechanisms that need further investigation.

This work was supported by the project "Grant Schemes at CU"(reg. No. CZ.02.2.69/0.0/0.0/19_073/0016935)


Marjorie Fonnesu

Nicola Kuczewski
email: nicola.kuczewski@ubiv-lyon1.fr

Doubts on the efficacy of outliers correction methods

While the utilisation of different methods of outliers correction has been shown to counteract the inferential error produced by the presence of contaminating data not belonging to the studied population; the effects produced by their utilisation when samples do not contain contaminating outliers are less clear. Here a simulation approach shows that the most popular methods of outliers correction (2 Sigma, 3 Sigma, MAD, IQR, Grubbs and winsorizing) worsen the inferential evaluation of the studied population in this condition, in particular producing an inflation of Type I error and increasing the error committed in estimating the population mean and STD. We show that those methods that have the highest efficacy in counteract the inflation of Type I and Type II errors in the presence of contaminating outliers also produce the stronger increase of false positive results in their absence, suggesting that the systematic utilisation of methods for outliers correction risk to produce more harmful than beneficial effect on statistical inference. We finally propose that the safest way to deal with the presence of outliers for statistical comparisons is the utilisation of non-parametric tests.

Introduction

Outliers have been defined as observations that "deviate so much from other observations as to arouse suspicion that it was generated by a different mechanism" [START_REF] Hawkins | Identification of Outliers[END_REF] or observations that "appears to be inconsistent with the remainder of that set of data" [START_REF] Barnett | Outliers in Statistical Data[END_REF]". The first point highlighted by these definitions is that outliers are generally considered as contaminating data that do not belong to the population of the parameter that is under investigation. As a consequence, the presence of outliers in the analysed sample is seen as an external element worsening the inferential evaluation of the studied population. The second point highlighted by the aforementioned definitions is that the outliers can be identified only relatively to the available data sample; outlier detection is therefore submitted to the inferential errors linked to random sampling. Several methods of outliers detection have been developed among which the most popular are the following: a) the mean plus or minus two-three standard deviation (Sigma 2 and Sigma 3) methods, b) the interquartile range (IQR) method, c) the absolute deviation around the median (MAD) method, d) the Grubbs method and e) the trimming method. [START_REF] Barnett | Outliers in Statistical Data[END_REF][START_REF] Seo | A Review and Comparison of Methods for Detecting Outliers in Univariate Data Sets[END_REF][START_REF] Leys | How to Classify, Detect, and Manage Univariate and Multivariate Outliers, With Emphasis on Pre-Registration[END_REF] for more details on the different methods. Once detected outliers are in general removed for the sample or transformed to be reduced in the limit of the acceptable data, a procedure known as winsorizing. The efficacy of different outliers correction methods has been evaluated and compared mainly by approaches using simulation of data sets [START_REF] Zimmerman | A Note on the Influence of Outliers on Parametric and Nonparametric Tests[END_REF][START_REF] Zimmerman | Increasing the Power of Nonparametric Tests by Detecting and Downweighting Outliers[END_REF][START_REF] Seo | A Review and Comparison of Methods for Detecting Outliers in Univariate Data Sets[END_REF][START_REF] Liao | Outlier Impact and Accommodation Methods: Multiple Comparisons of Type I Error Rates[END_REF][START_REF] Liao | Outlier Impact and Accommodation on Power[END_REF] or by randomly extracting real data from a known population arbitrary divided in 'Normal' and 'Outliers' subjects [START_REF] Osborne | The power of outliers (and why researchers should ALWAYS check for them[END_REF]). These approaches have demonstrated that the presence of contaminating outliers in the analysed sample worsen the inferential evaluation of population, producing an inflation of Type I and Type II error as well as an increase of the error committed in estimating the population effect. They have also highlighted the benefit of the utilisation of the different correction methods in contaminated sample, that in several cases counteract the negative effect produced by the presence of the outliers. While these approaches have focalised their attention to the effect produced by outliers correction in the case where outliers are present in the analysed samples, they do not evaluate the impact produced by the application of these methods on samples that do not contain external contaminating data. Indeed, it is conceivable that some of the data belonging to the studied population would be considered as outliers when screened by the different methods; this as a consequence of the rule of detection (ex. Sigma 2 methods would consider outliers ~5% of the data coming from normally distributed population) as well as a consequence of the inferential error discussed above. From now on such outliers produced by inherent variability will be designated as random sampling outliers (RSO). Given that the presence of RSO in the data set cannot be discriminate from contaminating outliers, the utilisation of outlier correction methods could in several circumstance target mainly or exclusively the former. By using the mean plus or minus sigma method on simulated data coming from normal distribution, Bakker and Wicherts showed that in the absence of contaminating outliers, the removal of RSO produce an increase of the occurrence of Type I when the p-value was evaluated with a t-test statistic [START_REF] Bakker | Outlier removal, sum scores, and the inflation of the Type I error rate in independent samples t tests: the power of alternatives and recommendations[END_REF]. While this report highlight the possible risk of outlier removal the conclusions are limited to a specific experimental condition. The main objective of the present work was therefore to extend the investigation on the impact produced by the correction of RSO on inferential statistic. To this goal we use an approach based on simulated data to address the following questions: 1-What is the proportion of RSO detectable as a function of correction method, samples size and type of data distribution? 2-Does the presence of RSO in the analysed samples affect the inferential evaluation of the studied population? 3-What is the impact of the correction of RSO on inferential statistic? In particular, on the evaluation of population parameters (mean and std), on the occurrence of Type I and Type II error and on the estimation of the real effect when samples belong to different populations. 4-Does an eventual impact of RSO correction on Type I and Type II error depends on the type of statistical test used to compare the samples? 5-Do the different methods of outliers correction have similar impact on inferential statistic depending whether samples present RSO or contaminating outliers?

An important difference compared to Bakker and Wicherts 2014, is the fact that in the present work the average impact produced by the correction of RSO is evaluated on the selected samples population containing RSO and not on the whole samples population. We believe that this approach is more informative for a researcher that, facing to an experimental sample for which some of the data can be labelled as outliers, has to evaluate the risk-benefit ratio of their correction.

Methods

Data simulations were performed with Python 3.7. np.random.normal and np.random.lognormal function were used to generate 100000 random samples, or pairs of samples, for each condition. This number was empirically determined by the procedure depicted in supplementary figure 1. A similar approach was used to determine the number of permutations (600) performed for the permutation test (not shown). These simulations were performed from samples coming from normal and log-normal distributions having µ=0 and σ=1. It should be pointed out that for log-normal distribution µ and σ are not the values for the distribution itself but that of the underlying normal distribution it is derived from (random.lognormal function from numpy), what give a real µ ≃ 1.649 and real σ≃2.16 . Outlier correction methods: Grubbs was imported from Smirnov_grubbs package; Winsorizing was performed with mstats.winsorize with 5% winsorisation; Median absolute deviation used for MAD method was performed with mad_std function (astropy.stats), values x where considered to be outliers when │(x-Median(sample))│/MAD>2.24 . In some of the simulations the more restrictive threshold of 3 for MAD detection was used [START_REF] Leys | Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median[END_REF]. In Winsorizing sigma 2 method data above 2 STD from sample mean were reduced to 2 STD value. All the scripts are available at the Open Science Framework (https://osf.io/gzra8/).

Results

Part 1: Random Sampling Outlier (RSO)

Here RSO are defined as the data present in a sample that, despite belonging to the studied population, are considered to be outliers by the different detection methods. In these simulations therefore all the data are sampled from the target distributions, i.e. no external contaminating data are injected in the samples. It should be stressed the fact that some of the data can be labelled as outliers with respect to a specific sample without necessary be extreme values of the studied population. In the same way extreme values of the studied population could be considered as normal value with respect of a specific sample (some examples of these cases are presented in supplementary figure 4).

Simulation 1: Impact of RSO correction on the estimation of population parameters

Simulation 1 was performed by the procedure descripted in Figure 1 (the script is available at https://osf.io/753yc/). This allows to extract, for the different detection methods, 10000 samples of sizes n having at least one RSO in them. The error in the estimation of population mean and standard deviation was expressed in standard deviations of the population distribution.

The objectives of this simulation were:

1-Determine the proportion of RSO detected by the different methods for different sample sizes 2-Determine whether and how the correction of RSO from the samples modify the estimation of population average and standard deviation

Correction of RSO increase the error in the estimation of population mean and std (samples from normal distribution)

The first objective was to determine the probability to observe at least one RSO in a random sample from normally distributed population (the RSO probability) as a function of sample size and outliers detection method (Figure 2 red lines). With the exception of "Grubbs" methods the RSO probability increase with the sample size eventually becoming 100% for large sample.

As a consequence the difference in RSO probability between different methods is principally observed for the smallest sample (n<100-500), where it is higher when IQR and MAD method are used. "Grubbs" method is particularly inefficient in detecting RSO, especially when sample size increase. It is interesting to note that RSO detectable by "Sigma 3" method are never present when sample size is ≤11 (Fig 2C), such results is in agreement with theoretical prediction [START_REF] Shiffler | Maximum Z Scores and Outliers[END_REF]. The second objective was to determine whether RSO correction impact the estimation of real population mean (µ) and STD(σ). As show in figure 2A and 3A, when samples are not selected based on the presence of outliers the estimation error for the two parameters decrease with the sample size. Similar decrease is observed when samples are selected for the presence of RSO (Fig 2 B-H and 3B-H blue lines). Outliers correction never reduce the average error in estimating µ (Fig 2) and σ (Fig 3). On the contrary most of the correction methods worsen the estimation of this two parameters, especially for the smallest sample size (n<50). The only exception is the "Accommodation sigma 2" method for which outliers correction do not modify the error in estimating µ and produce only a slight increase of the error in σ estimation. Supplementary figure 2 allow to compare the impact of "Sigma 2" and "Accommodation sigma 2" methods on the µ and σ estimation error distributions when sample size =6. It also shows an example of sample correction that brings to an increase of µ and σ evaluation error when "Sigma 2" correction method is used but no relevant change on the evaluation of these two parameters when "Accommodation sigma 2" method is used. Readers interested to run these simulations with other parameters are invited to use the following scripts: https://osf.io/dfeas/ for the error distributions, https://osf.io/t9mvc/ for sample correction example.

Correction of RSO increase the error in the estimation of population mean and std (samples from log-normal distribution)

Compared to normal distribution the probability that a sample data would be considered as an outlier is greatly increased when samples came from log-normal distribution. The RSO probability increase with the sample size for all tested detection It should be noted that when "Accommodation sigma 2" method is used a small decreased of µ estimating error is observed for sample size between 6 and 30 (Fig 5H). Supplementary figure 3 allow to compare the impact of "Sigma 2" and "Accommodation sigma 2" methods on the µ and σ estimation error distributions when sample size =100. It also shows un example of sample correction that brings to an increase of µ and σ evaluation error when "Sigma 2" correction method is used but only a relevant change on the evaluation of σ when "Accommodation sigma 2" method is used.

Conclusions from simulation 1

Simulation 1 shows that:

1-When using current detection methods some of the data belonging to studied population are considered to be outliers (i.e. the tested detection methods produce RSO)

2-The probability of RSO detection increase with the sample size, depend on the detection method and is larger when samples come from log normal than normal distribution.

3-The correction of RSO (removal or winsorisation) worsens the estimation of population mean and standard deviation that can be drawn from the samples.

Simulation 2. Impact of RSO correction on type I error

Simulation 2 was performed by the procedure descripted in Figure 6 (the script is available at https://osf.io/jbuvc/). This allow The objectives of this simulation were:

1-Determine whether the presence of RSO affect the occurrence of Type I error rate.

2-Determine whether and how the correction of RSO from the samples modify the occurrence of type I error.

3-Determine whether the effect of outlier correction depend on the statistical test used for samples comparison.

Correction of RSO increase the occurrence of Type I error (samples from normal distribution)

As shown in figure 7 (uncoloured lines), in samples selected for the presence of RSO the Type I error is, in general, maintained at the σ risk (5%). The only exception is observed when RSO are detectable by "Sigma 3" or "Grubbs" methods. In these cases, an inflation of Type I error is observed, but only for sample size between 6 and 100 and when Mann Whitney test is used.

Correction of RSO produce an inflation of Type I error for all correction methods except for "Accommodation sigma 2" (coloured lines). Type I error inflation after correction is higher for small sample size and is abolished for n≥ 500 when "Sigma 3" and "Grubs" methods are used. Moreover, inflation produced by outlier correction is slightly more important when Permutation and Ttest are used compared to Mann-Whitney test. Supplementary figure 4 show some examples where outliers correction generates false positive results (readers interested to run this simulation with other parameters are invited to use the following scripts: https://osf.io/zhyk3/). We also tested winsorizing correction at 10 and 20 %. In these conditions a strong increase of type I error was produced (Supplementary figure 5).

Conclusions from simulation 2

The results from simulation 2 show that:

1-The presence of RSO can produce a small inflation of Type I error only in the case in which RSO are detectable with "Sigma 3 "and "Grubbs" methods, samples belong to normal distribution, have a size between 6 and 50 and 2-The correction of RSO (removal or winsorisation) produces an inflation of Type I error for all tested method except when the "Accommodation sigma 2" method is used in association with Mann Whitney test.

3-The inflation of type I error produced by RSO correction is more important when samples are collected from lognormal distribution.

4-The inflation of type I error produced by RSO correction is more important when statistical comparison of the sample is made with Permutation test and Ttest compared to Mann Whitney test.

Simulation 3. Impact of RSO correction on type II error

Simulation 3 was performed by the procedure descripted in Figure 9 (the script is available at https://osf.io/jkg75/ ). This allows to extract, for the different detection methods, 10000 pairs of samples of sizes n, belonging to the separate populations and having at least one RSO in them. The occurrence of type II error was calculated by dividing the number of comparison between pairs giving a false negative (FN) result (i.e. when p-value>0.05) by the number of pair (10000). Statistical comparisons were made using Student ttest, Mann-Whitney test and permutation test. In these simulations a statistical power of either 50% or 95% was obtained by the modification of the distributions mean (µ2) from which sample 2 was extracted as a function of the required sample size and distribution types. When samples came from normal distribution the value of µ2 was chosen in order to observe the target power when t-test is used. When samples came from log-normal distribution the value of µ2 was chosen in order to observe the target power when Mann-Whitney test is used (Table 1). The mean of the distribution from which sample 1 was extracted was µ=0. Distributions standard deviation σ=1.

The objectives of this simulation were:

1-Determine whether the presence of RSO affect the occurrence of Type II error rate.

2-Determine whether and how the correction of RSO from the samples modify the occurrence of type II error.

3-Determine whether the effect of outlier correction on type II error depend on the statistical test used for samples comparison and statistical power. were respectively set equal to 0 and 1. For normal distribution, µ2 were calculated using GPower3 (the target power is reached when t-test is used). For log-normal distributions µ2 were empirically determined in order to attend the target power when Mann-Whitney test is used.

µ2

Correction of RSO mainly decrease the occurrence of Type II error when samples come from normal population.

As shown in figure 10

Conclusions from simulation 3

All together the results from simulation 3 show that:

1-While using non-parametric test with normal distribution do not considerably impact the occurrence of Type II error, parametric and permutation test associated to log normal distribution produce a strong inflation of false negative results. 2-The presence of RSO, detectable with the tested methods, do not modify the probability to observe false negatives results (that remain at the target β risk). 3-Outliers correction methods mainly lead to a reduction of Type II error especially from samples belonging to log normal distribution 4-The reduction of Type II error produced by RSO correction is observed for conditions that increase the risk of Type I error in the absence of real effect (compare figures 7 and 8 with figures 10,11,12 and 13).

Simulation 4. Impact of RSO on the estimation of real effect

Simulation 4 was performed by the procedure descripted in Figure 14 (the script is available at https://osf.io/nykbm/). This allows to extract, for the different detection methods, 10000 pairs of samples of sizes n, belonging to the separate populations and having at least one RSO in them. The impact of outlier correction on the estimation of population effect was quantified by comparing the difference of the populations means (µ1-µ2) with the difference of samples means (X̅ 1-X̅ 2). More precisely, the error on effect estimation (EEE) is calculate as the absolute deviation of sample effect from population effect (EEE=abs((µ1-µ2)-(X̅ 1-X̅ 2)). In these simulations a statistical power of either 50% or 95% was obtained by the modification of the distributions mean (µ2) from which sample 2 was extracted as a function of the required sample size and distribution types (table 1).

The objective of simulation 4 is to determine whether outlier correction modify the estimation of population effect that can be drown from the samples.

Correction of RSO increase the error committed to evaluate the population effect from the samples comparison (samples from normal population).

Comparison of figures 15A and 16A show that the error committed in the estimation of real effect from the samples decrease with statistical power. They also show that EEE, expressed as % of modification from real effect, increase with the sample size.

It is important to outline that the absolute EEE do indeed decrease with sample size (see data at https://osf.io/v4d9f/), but to lesser extent that the decrease of real effect needed to keep the statistical power at 50% or 95% (see Table 1). This has the consequence to produce the depicted increase of relative modification whit the increase of sample size. Outlier correction produce an increase of the EEE for the majority of the tested methods (

Correction of RSO mainly decrease the error committed to evaluate the population effect from the samples comparison (samples from normal population).

When sample belongs to log-normal population the error committed in the estimation of real effect is slightly higher than for sample belonging to normal population. The presence of RSO slightly increase the EEE only when RSO are detectable with Sigma 3 method and sample size is small (compare Fig 17A and18 A with 17D and 18 D blue lines, sample size=12). Impact of outlier correction change as function of statistical power and correction method. For low power the estimation error on effect quantification is reduced by outlier correction, whatever the method used (Fig 17). However, for comparisons in high statistical power condition the EEE is slightly reduced when "Sigma 3", "Sigma 2", "Accommodation sigma 2" and "Winsorizing" methods are used and increase with the other tested methods (Fig 18).

Conclusions from simulation 4

All together the results from simulation 4 show that:

1-When data come from normally distributed population, outliers correction methods applied to RSO worsens the estimation real effect that can be drawn from the samples, with the exception of winsorizing and accommodation sigma 2 methods.

2-When data come from population log-normally distributed, outliers correction methods applied to RSO improve the estimation real effect that can be drawn from the samples, if sampling is carried in low powered condition but worsen the EEE in higher powered condition when MAD, IQR and Grubbs methods are used.

Simulation 5. Use of outlier correction methods for p-hacking

As shown before (simulation 2) outlier correction in the absence of real outliers increase the probability of type I error for most of the methods tested. Combination of outliers correction methods could therefore potentially be used as strategy for data dredging. To test the impact of this strategy on the occurrence of Type I error we used the procedure descripted in Figure 19 (the script is available at https://osf.io/7s3gv/). As shown in figure 20 Type I error is strongly inflate when outliers correction is used with a the p-haking strategy. The inflation is greater when samples come from log-normal distribution becoming close to 50% for large sample size.

Conclusions from part 1

The main conclusions are that the presence of extreme values belonging to the studied population, i.e. the RSO, do not or only slightly impact the inferential evaluation of real population from the sample. In particular, for what concern the evaluation of population parameters (mean and STD), Type I and Type II error and real effect estimation. On the other hand, the correction of RSO by most of the methods tested in this study worsen the inferential evaluation of the studied population by increasing the error in the estimation of population parameters (simulation 1) and Type I error (simulation 2). In several condition corrections of RSO can drop the occurrence of Type II error below the expected β risk (simulation 3), however such decrease is often associate with a stronger error in the evaluation of the real effect, in particular for data coming from normal distribution (simulation 4). The impact of the correction change with the method used, with accommodation methods (accommodation sigma 2 and winsorizing) having the lowest negative impact on parameter estimation, Type I error and real effect evaluation when associated to Man Whitney test but not producing reduction of Type II error rate. On the other hand, method having the larger beneficial effect in Type II error (MAD, Sigma 2, IQR) also produce the stronger inflation of Type I error. In general, the negative impact of outlier correction is reduced when statistical comparisons are performed with Mann Whitney test. Finally, our results show that the improper use of outliers correction methods for data dredging can bring to a strong inflation of false positive results (simulation 5).

Part 2: Contaminating outliers

The evaluation of the effect produced by sample contamination with data which do not belong to the studied population (contaminating outliers), was performed through simulations on which extreme values were injected in the analysed sample (similar to the approach used by [START_REF] Bakker | Outlier removal, sum scores, and the inflation of the Type I error rate in independent samples t tests: the power of alternatives and recommendations[END_REF][START_REF] Liao | Outlier Impact and Accommodation Methods: Multiple Comparisons of Type I Error Rates[END_REF][START_REF] Liao | Outlier Impact and Accommodation on Power[END_REF]. In practice "n" data of one of the two samples were substituted by "n" outliers randomly selected from uniform distribution situated between 4 and 8 (simulation of Type I error) or 4 and 5 (simulation of Type II error) STD of the original distribution. Where n= 0,1,2,3,4,5,6,7. When Type II error was evaluated, outliers were injected in the sample coming from the distribution having µ=0. In this way increasing false negative rate.

Conclusions from part 2

The main conclusions are that the presence of extreme values not belonging to the studied population, i.e. the real outliers, negatively impact the inferential evaluation from the sample. All the correction methods tested counteract the increase of the error committed in the evaluation of population mean and STD that is produced by outliers injection; with MAD, IQR and Sigma 2 methods having the highest efficacy. This irrespectively of the fact that sample came from normal or log-normal distribution The methods that show the larger benefit in the presence of contaminating outliers (MAD, IQR and Sigma 2) are also those that produced the higher inflation of Type I error in the absence of contaminating outliers (Fig 9 and10).

Discussion

As shown in the present report, the inefficacy of common methods of outliers correction stems from the fact that none of these is capable to distinguish outliers data belonging to the studied population, i.e. the RSO, from contaminating data. Since RSO occurs with high probability and since their removal worsen the inferential evaluation of the studied population, the systematic utilisation of outliers correction risk to produce more harmful than beneficial effect on statistical inference. This is particularly true for what concern Type I error that is considerably increased by the removal of RSO. This result is in agreement end extend what previously report by Bakker and Wicherts 2014 that found an increase of Type I error when methods based on mean plus or minus sigma were applied to samples coming from normal distribution. In their report however the inflation of false positive results was slightly lower than what reported here and increased whit the sample size instead of decreasing.

This discrepancy is likely due to the fact that in their simulation were also included samples not contains RSO, therefore the inflation Type I error depend both on the fallacy of the correction method as well as on the its probability to detect RSO. We believe that the selective selection of samples presenting RSO, allow a better comparisons of the risk-benefit ratio associated to the different correction methods. Whit this approach we showed that those methods that, in the absence of contaminating outliers, present the lower risk of type I error inflation and the lower risk of inferential evaluation worsening (ex. the Accommodation Sigma 2) also have the lowest efficacy in counteracting Type I and Type II error inflation produced by contaminating outliers. On the contrary those method having the highest efficacy against contaminating outliers (ex. MAD) produce stronger inflation of Type I error in their absence. As suggested by simulation 5, the inflation of false positive results could be further accentuated by the inexistence of specific guidelines on the use of outlier correction methods, that leaves the researcher free to choose or not the method that best fits its own data, or its expectancy on data interpretation [START_REF] Simmons | False-positive psychology: undisclosed flexibility in data collection and analysis allows presenting anything as significant[END_REF][START_REF] Holman | Evidence of Experimental Bias in the Life Sciences: Why We Need Blind Data Recording[END_REF]. Since their diffuse utilization [START_REF] Leys | Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median[END_REF][START_REF] Bakker | Outlier removal, sum scores, and the inflation of the Type I error rate in independent samples t tests: the power of alternatives and recommendations[END_REF] the systematic utilisation of outlier correction procedures investigated here are likely to contribute to the reproducibility crises associated to methodological limits and experimental bias [START_REF] Chambers | The Seven Deadly Sins of Psychology[END_REF][START_REF] Munafò | A manifesto for reproducible science[END_REF],the problem of Type I inflation linked to outlier correction should therefore be explicitly addressed by journal's editorial boards. How to deal with contaminating outlier? As proposed by Leys and coo workers the decision of removing a putative outlier from a data set should not be taken exclusively based on the outcome of a mathematical test (at least not one of those investigated in the present study), "Researchers are recommended to follow a two-step procedure to deal with outliers. First, they should aim to detect the possible candidates by using appropriate quantitative (mathematical) tools. …. Second, they should manage outliers and decide whether to keep, remove, or recode these values, based on qualitative (non-mathematical) information" [START_REF] Leys | How to Classify, Detect, and Manage Univariate and Multivariate Outliers, With Emphasis on Pre-Registration[END_REF].

Although in some specific condition mathematical methods appears to have a high reliability to mainly detect contaminating outliers (ex Sigma 3 method for sample size ≥500, see figures 7D, 8D and 24). The simulation approach used in the present reports could be a useful tool to look for conditions not investigated here or for seeking new mathematical methods more resistant to the inflation Type I error produced by removal of RSO. For the moment the most reliable, non-subjective, way to reduce the impact of contaminating outliers appears to be the utilisation of non-parametric test, at least when data come from

Fig 1 :

 1 Fig 1: Procedure used to simulate the effect of the correction of "Random Sampling Outliers" on the estimation of real population mean and standard deviation. This procedure was separately applied for the different methods of outliers correction, sample size and population distribution (Normal or Log-normal). The error is expressed in population STD. "Number of sampling "and "Number of outliers "were set to zero at the beginning of the simulation (script is available at https://osf.io/753yc/).

Fig 2 :Fig 3 :

 23 Fig 2: Outlier correction (RSO) increase the average error on population mean estimation (data from normal distribution). In A samples are not selected based on the presence of outliers. Note that in C and G sample size start from 12 and 20 respectively. Vertical bars represent 95% CI. All data can be found at https://osf.io/pjcak/

  methods and, when n=100, become 100% whatever methods is used(Fig 4, red line). Even in this case no data further than 3 STD from sample mean are observed when sample size is ≤11 (Fig4 C). The error in estimating µ and σ from sample mean and sample STD is larger compared to sample coming from normal distribution but still decrease with sample size(Fig 4A and Fig 5A). Again outliers correction worsen the estimation of µ and σ but in this case the worsening increase with the sample size.

Fig 4 :Fig 5 :

 45 Fig 4: Outlier correction (RSO) increase the average error on population mean estimation (data from log normal distribution). In A samples are not selected based on the presence of outliers. Note that in C and G sample size start from 12 and 20 respectively. Vertical bars represent 95% CI. All data can be found at https://osf.io/gq63d/

  to extract, for the different detection methods, 10000 pairs of samples of sizes n, belonging to the same population and having at least one RSO in them. The occurrence of type I error was calculated by dividing the number of comparison, within pairs, giving a false positive (FP) result (i.e. when p-value<0.05) by the number of pair (10000). Statistical comparisons were made using Student Ttest, Mann-Whitney test and permutation test.

Fig 6 :

 6 Fig 6: Procedure used to simulate the effect of outlier correction on the occurrence of Type I error for "Random Sampling Outliers". Outliers where separately detected for sample 1 and sample 2. Detection means at least 1 outlier in either sample 1 or sample 2. This procedure was separately applied for the different methods of outliers detection and different type of population distribution (Normal, Log-normal and Exponential). F.P.= false positives. Number of outliers, number of sampling and F.P. were set equal to 0 at the beginning of the simulation. M.W. = Mann Witney test ; Permut.= Permutation test .(Script is available at https://osf.io/jbuvc/).

Fig 7 :

 7 Fig 7: Impact of correction of "Random sampling outlier" on Type I error (data from normal distribution). In A samples are not selected based on the presence of outliers. All data can be found at: https://osf.io/cq4w3/Correction of RSO increase the occurrence of Type I error (Samples from log-normal population)As shown in figure8(uncoloured lines), the presence of RSO do not impact the occurrence of the Type I error, that is maintained at the σ risk (5%). However, outlier correction produces a strong inflation of Type I error, especially when Permutation test and Ttest are used (green and red lines). In these cases, the inflation of Type I error produced by the correction increase with the sample size. A notable exception is the utilization of "Accommodation sigma 2" that, used in association to Mann-Whitney test, keep the Type I error at the σ risk.

Fig 8 :

 8 Fig 8: Impact of correction of "Random sampling outlier" on Type I error (data from log-normal distribution). In A samples are not selected based on the presence of outliers. All data can be found at: https://osf.io/hezmw/

  comparison is made Mann-Witney test (Fig 9 D and F, black line).

Fig 9 :

 9 Fig 9: Procedure used to simulate the effect of outlier correction on the occurrence of Type II error for "Random Sampling Outliers". Outliers where separately detected for sample 1 and sample 2. Detection means at least 1 outlier in either sample 1 or sample 2. This procedure was separately applied for the different methods of outliers detection and different types of population distribution (Normal, Lognormal). F.N.= false negatives. Number of outliers, number of sampling and F.N. were set equal to 0 at the beginning of the simulation. M.W. = Mann Witney test ; Permut.= Permutation test. Script is available at https://osf.io/jkg75/.

  and 11 (uncoloured lines), in samples selected for the presence of RSO the Type II error is maintained at the targets β risk (50% and 95%) when Ttest and Permutation test are used. On the other hand, RSO detectable by "Sigma 3" method reduce occurrence of false negatives when samples obtained in low powered conditions are compared with the Mann-Whitney test(Fig 10D). Outlier correction has a different effects depending on the statistical power. In low powered conditions outlier correction reduce the probability to observe false negative results, in particular when Ttest and Permutation test are used (Fig10). In higher powered conditions outliers correction has more heterogeneous effect, producing either a reduction or an increase of type II error that depend on method and the test used(Fig 11). These modifications are extremely mild (1-2%) except when RSO are detectable by Grubss test and sample size =6, in this case a consistent increase of false positives is observed after correction (Fig11 F).

Fig 10 :

 10 Fig 10: Impact of correction of "Random sampling outlier" on Type II error when statistical power = 50% (data from normal distributions). In A samples are not selected based on the presence of outliers. All data can be found at https://osf.io/p4t9q/

Figure

  Figure 12A and 13A show the inadequacy of Ttest and Permutation test to detect statistically significant differences between samples arising from populations having a log-normal distribution. A strong increase of Type II error is observed in this condition that is slightly lower when permutation test is associated to smaller sample size. In samples selected for the presence of RSO there are not relevant modification of Type II error both when using parametric and non-parametric tests, in the latter case false negative rate remain at the targets β risk (50% and 95%) (uncoloured line of Fig 12 B-H, Fig 13 B-H). Outlier correction has a beneficial impact on the occurrence of false negative results. First it strongly reduces the inflation of type II error produced

Fig 14 :

 14 Fig 14: Procedure used to simulate the effect of outlier correction on the estimation of population difference Outliers where separately detected for sample 1 and sample 2. Detection means at least 1 outlier in either sample 1 or sample 2. This procedure was separately applied for the different methods of outliers detection and different types of population distribution (Normal, Log-normal). Script is available at https://osf.io/nykbm/.

  Fig 15: Impact of correction of "Random sampling outlier" on the estimation of population difference when statistical is power = 50% (data from normal distributions). EEE= error in the estimation of real effect. In A samples are not selected based on the presence of outliers. Vertical bars represent 95% CI. All data can be found at https://osf.io/v4d9f/.

Fig 16 :

 16 Fig 16: Impact of correction of "Random sampling outlier" on the estimation of population difference when statistical is power = 95% (data from normal distributions). EEE= error in the estimation of real effect. In A samples are not selected based on the presence of outliers. Vertical bars represent 95% CI. All data can be found at https://osf.io/xnb7r/

Fig 17 :

 17 Fig 17: Impact of correction of "Random sampling outlier" on the estimation of population difference when statistical is power = 50% (data from log-normal distributions). EEE= error in the estimation of real effect. In A samples are not selected based on the presence of outliers. Vertical bars represent 95% CI. All data can be found at https://osf.io/yrgwk/

Fig 19 :

 19 Fig 19: Procedure used to simulate the impact on type I error produced by the utilisation of outliers correction methods in a data dredging strategy. Seven correction methods are used ( MAD, Sigma 2, Sigma 3, IQR, Grubbs, Winzorizing and Accomodation sigma 2). F.P.= False positives results. "Number of sampling ", "Number of outliers "and "F.P." were set to zero at the beginning of the simulation (script is available at https://osf.io/7s3gv/)

Fig 20 :

 20 Fig 20: Inflation of Type I error produced by the utilisation of outlier correction methods through a data dredging strategy. Samples from normal and log normal distributions. Dashed line represent the σ risk (5%). All data can be found at https://osf.io/xm3nu/

Fig 25 :

 25 Fig 25: Comparison of outlier correction on type I error when MAD method is applied with different threshold. A and B) when sample do not contain contaminating outliers. C and D) in the presence of contaminating outliers. MW= Mann Whitney test

(

  simulation 6). Both before and after outlier correction the inflation of Type I end Type II error produced by outliers injection is lower when samples are compared with non-parametric (Mann Whitney) than parametric (Ttest) test (except the particular case of accommodation methods for Type II error in log-normal distributions, Fig23D). When sample size=20 only MAD, IQR and Sigma 2 methods are able to counteract the inflation of Type I and Type II errors produced by outlier injection, however their efficacy is observed only when sample came from normally distributed population, being often responsible of a further increasing of false positive or false negatives results when data are sampled from log-normal distribution (simulation 7 and 8).

  normal and log-normal distribution. Indeed, as shown in figure 7A, the use of Mann Whitney in the absence of contaminating data keep Type I error at σ risk when is used with normally and log normal distributed data, while Type II error is only slightly increase when samples came from normal distribution (Fig 10 A and 11 A). On the other hand, the inflation of Type I and Type II error produced by contaminating outliers is strongly reduced by the use of Mann Whitney test compared to Ttest for samples coming both from normal or log-normal distribution (Fig 22 and 23).

for normal distribution µ2 for log-normal distribution Real µ2 of the log- normal distribution Sample size

		Power	Power	Power	Power	Power	Power
		50 %	95 %	50 %	95 %	50 %	95 %
	6	1.252	2.3	1.352	2.49	6.37	19.89
	10	0.926	1.7	0.94	1.76	4.22	9.58
	12	0.837	1.54	0.85	1.575	3.86	7.96
	20	0.636	1.17	0.655	1.195	3.17	5.45
	30	0.515	0.95	0.528	0.97	2.79	4.35
	50	0.396	0.73	0.405	0.74	2.47	3.45
	100	0.279	0.51	0.284	0.525	2.19	2.78
	500	0.124	0.23	0.128	0.235	1.87	2.08
	1000	0.0877	0.16	0.09	0.164	1.80	1.94

Table 1 : µ2 values of the simulated distributions allowing to have the required statistical power for

 1

a given sample size. µ1 and Sigma

Acknowledgments:

This work was supported by the CNRS, Inserm, and Lyon 1 University. We wish to thank Christophe Leys for its pertinent comment on the manuscript as well as Nicolas Fourcaud-Trocméand and Samuel Garcia for helping produce the scripts for data simulation.

Simulation 6. Impact of injected outliers and their correction on the estimation of population parameters

This simulation was performed by the procedure depicted in figure 1 with the exception that sample was contaminated with the injected outliers and 200000 sample were analysed regardless the fact that outliers were detected. Only results obtained for sample size equal to 20 are presented. Readers interested to run the simulation with other parameters are invited to use the script available here (https://osf.io/quzm4/).

As expected outlier injection increase the error committed in estimating population mean and STD inferred from sample mean MAD method indubitably showing the higher efficacy. It is important to note that the beneficial effect of outlier correction is present only when the number of outliers in each sample is larger than two (more than 20% of the data are contaminating outliers). On the other hand, when samples come from Log-normal distribution outlier correction further increase the inflation of Type I due to outlier injection. This counterintuitive result, illustrated by the examples in Figure 22 E and F, appears to be due to the fact that contaminating outliers tend to be considered normal data with respect to the sample while at the same time RSO are detected in the sample not having injected outliers.

It's important to note that the inflation of Type I error produced by outlier injection is strongly reduced when samples are compared with non-parametric test with respect to the use of parametric test (compare black triangles in A Vs B and C Vs D).

Simulation 8. Impact of injected outliers and their correction on Type II error

This simulation was performed by the procedure depicted in figure 9 with the exception that sample 1 was contaminated with the injected outliers and 200000 paired of sample were analysed regardless the fact that outliers were detected. The results presented are obtained in a condition in which sample size equal to 20 and statistical power in the absence of outliers is 95%.

Readers interested to run the simulation with other parameters are invited to use the script available here https://osf.io/8ntwy/).

As expected outlier injection in one of the two samples increase the occurrence of Type II error for samples coming both from normal and log-normal distribution (Fig 23 A-D). Again it should be noticed that this inflation is lower when samples are compared with Mann Whitney test (black triangles). When samples came from normal distribution outlier correction counteract the inflation of Type II error associated to outliers injection. Particularly when MAD, IQR and Sigma 2 methods are used. The relative effect of outlier correction is stronger when samples are compared with the parametric test. But this is due to the fact that the Type II inflation in the absence of correction is bigger when Ttest is used. Indeed, the absolute value of Type II error after outliers correction is lower for Mann Whitney than Ttest condition, (use horizontal lines in The results of simulation 7 and 8 show that, when sample size =20, the correction methods that produce the larger benefit in the presence of contaminating outliers are also those that produced the higher inflation of Type I error in their absence (MAD, IQR and Sigma 2 ; Fig 7 and8). In order to look for the existence of experimental conditions for which Type I error is maintained at ~ σ risk in absence of outliers and inferential statistic still benefit for correction of contaminating data, we tested Sigma 3 outlier method for sample size ≥500 (a condition for which no inflation of Type I error is observed in the absence of outliers, fig 7D and8D). As show in figure 24 when data came from normal distribution this method efficiently counteract the inflation of Type I and Type II error produced by outliers injection while a small increase of the two errors is produced when samples came from log-normal distributions. Up to now the simulations using MAD correction method have been performed using a threshold of 2.24. However recent reports suggest that the more restrictive threshold of 3 should be used instead [START_REF] Leys | Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median[END_REF][START_REF] Leys | How to Classify, Detect, and Manage Univariate and Multivariate Outliers, With Emphasis on Pre-Registration[END_REF]. We therefore compared the impact of outlier correction on Type I error for the two thresholds. As shown in figure 25A and 25B the inflation of Type I error in the absence of contaminating outliers is reduced by the utilization of the more restrictive threshold. However, it remains largely above the σ risk, especially for samples coming from log-normal distribution and lower sample sizes when sample came from normal distribution. On the other hand, increasing the threshold reduced the efficacy of MAD method to counteract the inflation of Type I error produced by contaminating outliers ( Supplementary figure 2: A-B Distribution of µ and σ estimation errors, before and after outliers correction. A When "Sigma 2" method is used. B When "Accommodation sigma 2" method is used. C. Example of sample correction that brings to an increase of µ and σ estimation error when "Sigma 2" correction method is used but no relevant change on the estimation of these two parameters when "Accommodation sigma 2" method is used. Sample from normal distribution. Sample size=6. Vertical dot lines represent the means of error distributions. Scripts for simulation A-B and C are respectively available here https://osf.io/dfeas/ and here https://osf.io/t9mvc/. Supplementary figure 3: A-B Distribution of µ and σ estimation errors, before and after outliers correction. A When "Sigma 2" method is used. B When "Accommodation sigma 2" method is used. C. Example of sample correction that brings to an increase of µ and σ estimation error when "Sigma 2" correction method is used but no relevant change on the estimation of these parameters when "Accommodation sigma 2" method is used. Sample from log-normal distribution. Sample size=100. Scripts for simulation A-B and C are respectively available here https://osf.io/dfeas/ and here https://osf.io/t9mvc/.