{'fraction_non_alphanumeric': 0.027141012188231495, 'fraction_numerical': 0.009120994430543225, 'mean_word_length': 4.663009941720946, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parameter space exploration of ecological models

In recent years, we are seeing the formulation and use of elaborate and complex models in ecological studies. The questions related to the efficient, systematic and error-proof exploration of parameter spaces are of great importance to better understand, estimate confidences and make use of the output from these models. In this work, we investigate some of the relevant questions related to parameter space exploration, in particular using the technique known as Latin Hypercube Sampling and focusing in quantitative output analysis. We present the analysis of a structured population growth model and contrast our findings with results from previously used techniques, known as sensitivity and elasticity analyses. We also assess how are the questions related to parameter space analysis being currently addressed in the ecological literature.


Introduction
There is a growing trend in the use of mathematical modelling tools in the study of many areas of the biological sciences. The use of models is essencial as they present an opportunity to address questions that are impossible or impractical to answer in either in purely theoretical analyses or in field or laboratory experiments, and to identify the most important proccesses which should then be investigated by experiments. One compelling example is made by the Individual Based Models (IBM), which represent individuals that move and interact in space, according to some decision-making rules. These models permit a great level of detail and realism to be included, as well as linking multiple levels of complexity in a system.
On the other hand, more realistic models employ a vast selection of input parameters, from temperature and rainfall to metabolic and encounter rates, which may be difficult to accurately measure. Moreover, one may be interested in estimating how much predictions for a model fitted in one place or to one species may be extrapolated to different places or species. While variations in some of those parameters will have negligible impact on the model output, other parameters may profoundly impact the validity of a model's predictions, and it may be impossible to determine a priori which are the most important parameters. In theory, this could be done by evaluating the model at all possible combinations of parameters, however this would require a prohibitive number of model runs, specially considering that a single run of those models may take days to complete. Our challenge then consists in providing the best estimates for the importance of the several parameters, requiring the least number of model runs.
The disciplines of uncertainty and sensivity analysis have been developed in the context of the physical sciences and engineering, and have been greatly developed in the 1980 and 1990 decades [24,33,10,15,14,1,16,47,19,38,17,18,27,34]. More recently, these analyses have been successfully applied to biological models, in order to explore the possible outcomes from the model output, estimate their probability distribution and the dependency of the output on different combinations of parameters, and to assess which parameters require more experimental effort in order to be more confidently estimated. This kind of parameter space exploration is considered a fundamental step prior to using the model in management decisions [2].
One approach to the parameter space exploration, which will be described here, is to generate samples from the parameter space, run the model with these samples, and analyse the qualitative or quantitative differences in the model output. Section 2 will present the sampling techniques, with empashis on the Latin Hypercube method, while section 3 will deal with quantitative analyses. We should emphasize that the analyses tools are not coupled with the sampling techniques used: one can, in principle, use the sampling techniques described and other analyses tools, or apply the analyses described here to a more general class of sampling methods. We also present an example of the sampling and analysis used in section 4. Then, we briefly review some relevant research papers which have used such techniques in the exploration of ecological models in section 5.

Parameter spaces
In order to better pose our questions, we need first to discuss some properties of the parameter space (or PS for short) of our models.
The parameters (or inputs) are quantities x 1 , x 2 , · · · , x m which will be used to run the model. In our discussion, we will assume that all the x i are real valued. These quantities are unknown, and one first challenge is to determine which set of values better fit a model to the available data, which is the subject of linear and nonlinear estimation.
However, the same model may be parametrized in different ways, as discussed by Ross [31]. For example, in population ecology, the logistic growth equation may be represented with two parameters r and K as: However, the same equation may be written as Here, the two parameters are α = r and β = −r/K. While the first equation is far more commonly used in the biological context, both are equivalent, and using one or the other model is simply a matter of choice. There are many other ways of writing this equation, and one of special interest when trying to fit real data is in terms of orthogonal polynomials, such as: Where θ 0 can be calculated from θ 1 and θ 2 . This more complicated equation has several numerical advantages over the previous, as the parameters θ 1 and θ 2 can be estimated with much more accuracy, and will not be correlated (as is the case with α and β, as well as r and K). However, these parameters are hard to interpret in biological terms. These different equations illustrate the existence of interpretable (Eq. 1), defining, or algebraic (Eq. 2) and computing (Eq. 3) parameters. Most of the times, it would be preferable to estimate the values that best fit some data by using computing parameters, and then to transform them to interpretable parameters in order to present the results.
Also, it should be mentioned that the parameter space may be constrained. This will have an impact on some of the available sampling and analysis techniques. The simplest constraint is requiring some parameter to be positive or negative. Also, there may be combinations of values that are meaningless. For example, if we model a community with N individuals and S species, the number of individuals and species, considered on their own, may be any positive number. However, it is clear that the number of species may not be bigger than the number of individuals, which imposes the condition S ≤ N . This condition is called a constraint, and limits the values that the parameter vector may assume.
If we consider the m-dimensional space consisting on all possible combination of values for the parameters, our parameter space will be the subset of this space that respects all our constraints. For example, consider that the parameters we are interested are two angles of a triangle. In this case, the sum of the angles must be less than 180 degrees, a 1 + a 2 < 180•. Clearly, this parameter space is not square, in the sense that, if we define the ranges of the variables a 1 and a 2 independently as (0, 180), not all combinations of parameters will be meaningful. What can be done in this case is to create a new parameterâ 1 , defined aŝ This new parameter varies between 0 and 1, and all combinations ofâ 1 , a 2 are points from our parameter space. Now, care must be exercised after applying such transformations in order to preserve the marginal distributions from the original variables, as exemplified on figure 1. Another related concept that should not be confused with the constraints is the correlation between variables. For example, acidic soils are likely to have a lower cation exchange capacity (CEC), and more alkaline soils are likely to have a larger CEC [39]. Thus, those variables are correlated. Correlations have a profound impact on some analysis, however, they are difficult to measure, and data on correlations are not generally available on the literature. We will return to questions related to correlations in parameter spaces in section 2.2.

Applications of parameter space exploration
Next, we turn our attention to the kind of problems we might want to address with the exploration of the parameter space. First, the simplest case is asking "is there a region of my parameter space where condition X holds?" This condition might be, for example, the extinction or coexistance of species, some pattern of distribution or abundance of species. We also might be interested in mapping where are these regions. In complex models, where several different regions might exist where the qualitative results of the models are very different, we may ask how many of these regions are there, as well as map the frontiers between them.
Another class of problems arises when the model produces some quantitative response, and we are interessed in determining the dependency of this response to the input parameters. For example, when modeling the dynamics of a population, we might want to know how the final population varies with each of the input parameters. In this context of quantitative analysis, the questions are divided in two classes: first, how much the variation of the input parameters is translated into the total variation of the results, which is the topic of uncertainty analysis, and second, how much of the variation in the results can be ascribed to the variation of each individual parameter, which is the topic of sensitivity analysis [15,14]. We will present the techniques and results from both uncertainty and sensitivity analysis in section 3. Also, the model which we wish to analyse can be any function of the input parameters. In particular, there are three classes of models that can be used. First, the model may be a complex mathematical function (for example, Y is defined by a differential equation). Second, the model may be a simulation model, like an IBM. Third, the model may be the result of fitting a statistical model.
All these problems may be formulated in a general way, defining some response from the model Y as a function of the input parameter vector x: In the equation 5, all the quantities are vectors, indicated by the boldface. Here, x = [x 1 , x 2 , . . . , x m ] represent the parameters to the model f , and Y = [y 1 , y 2 , . . . , y n ] represent the some quantitative responses from the model. In some sections, we will discuss the response as a single value y, without loss of generality.
Each of the input parameters x i is associated with a probability distribution D i (x), which represent our degree of knowledge about the values that x i may assume (see figure 2 for examples).
Taken together, all the distributions D i form the joint probability distribution of the parameters, D(x). This function takes into account not only the individual distribution of each parameter, but also all the correlation terms between them.
In very simple models, it may be possible to analitically deduce the behaviour of the model response taken at each point of the joint distribution of parameters. In the general case, however, this is impossible, and a way of investigating the model is to choose some points from the joint distribution and analysing the model at each point. Section 2 will present some strategies for choosing these points. shows the exponential distribution of life expectancy, which can be deduced from theory [5]. Pannel "b" shows a gamma distribution, which can be used to model the number of parasites per host, but has no theoretical derivation [4]. Pannel "c" shows data from an empirical study on the CO 2 uptake from plants [29], and pannel "d" shows an example where no prior information can be used.

Sampling Techniques
There are several strategies that can be used to choose the samples from the parameter space that will be used as input to our model of interest. Here, we will present some of them, along with their limitations, to justify our choice for the Latin Hypercube Sampling, which we will describe in section 2.1. One way of exploring the parameter space is by discretizing every distribution and running the model for every possible combination of values for all parameters. This is called full parameter space exploration, as done in [44], and altough it possesses many advantages, it may become very costly in terms of computer time. In addition, the number of possible combinations increases exponentially with the number of parameter dimensions considered.
To circumvent the exponential increase in the number of samples, it is usual to explore the parameter space in the following fashion: holding all but one parameter constant, we analyse how the output of a model is affected by one parameter dimension at a time (as done in [46]). This analysis is refered to as individual parameter perturbance. This kind of analysis is, however, limited by the fact that the combinations of changed parameters may give rise to complex and unexpected behaviours. Often, algorithms based on individual parameter perturbance are used as a first step in order to discriminate between parameters that may have a substantial impact on the output, and parameters that are less relevant.
Another viable option would be to chose N random samples from the entire space, in order to analyse both the effect of each parameter and the combined effect of changing any combined number of parameters. This sampling scheme is called random sampling, or Monte Carlo sampling, and has been applied to many biological models [20]. One important feature of the Monte Carlo sampling is that its accuracy does not depend on the number of dimensions of the problem, as shown in [22].
Stratified sampling strategies, which are a special case of Monte Carlo sampling, consist in strategies for chosing these random samples while, at the same, making sure that every subdivision (or strata) of the distribution is well represented. As shown on [24], the estimatives of the statistical properties (such as the mean or the variance) of the model output are better represented by stratified random sampling than by simple random sampling (see figure 3 for examples). As we shall see in the next session, the Latin Hypercube sampling is a practical and easy to understand stratified sampling strategy.
Another class of Monte Carlo methods that should be mentioned here is the Markov Chain Monte Carlo (MCMC), which is also used on similar analyses [22]. This methods consists in generating a sequence of points {x (t) } from the parameter space whose distribution converges to the joint probability distribution D(x), and in which each sample x (t) is chosen based on the previous x (t−1) . MCMC methods perform better than LHS methods for estimating the distribution of the model responses, however, they require a number of model runs which is orders of magnitude higher than LHS requirements.

Latin Hypercube: Definition and use
In this section, we describe the Latin Hypercube Sampling, and show how it can be used to efficiently solve the questions posed in section 1.
We also discuss what are the available methods for obtaining the LHS. Firstly, let us define, in the context of statistical sampling, what is a Latin Square: Definition If we divide each side in a square in N intervals, and then take samples from the square, the resulting square will be called Latin if and only if there is exactly one sample in each row and each column.
A Latin Hypercube is simply the generalization of the Latin Square to an arbitrary number of dimensions m. From these definitions, it is clear that the number of samples is fixed as N We will show how to estimate the optimal N in section 2.5, but for now it is relevant to note that N does not depend on the number of parameters considered.
We will now construct the Latin Hypercube. Let's fix our attention in one parameter dimension i of the parameter space. The first step we should take is to divide the range of x i in N equally probable intervals. In order to do so, we will turn our attention to the probability distribution of x i , defined on section 1.1 as D i . Recall that this probability distribution must be chosen in a way that represents our current understanding of the biology of the given system. This function might be estimated by an expert in the field, it might represent a data set from field or laboratory work, or in some cases it may be simply the broadest possible set of parameters, in some cases where the actual values are unknown or experiments are unfeasible (see fig. 2).
In possession of the distribution function D i , we must sample one point from each equally probable interval. There are two approaches used here: it is possible to choose a random value from within the interval (as proposed by [24]), or instead, we can use the midpoint from each interval [16]. As the statistical properties of the generated samples are very similar, we will use the second approach here.
The integral of the distribution function is called the cummulative distribution function F i (x). This function relates the values x that the parameter may assume with the probability p that the parameter is less than or equal to x. We will refer to the inverse of the cummulative distribution function, F −1 i , as the quantile function of the parameter x i , as it associates every probability value p in the range (0, 1) to the value x such that P (x i ≤ x) = p. We divide the range (0, 1) in N intervals of size 1/N , and use this quantile function to determine the x values as the midpoints of each interval. Summarizing, we take the N points, represented as as: The samples from each dimension are subsequently shuffled, to randomize the order in which each value will be used (see example on figure 4). As the samples come from the distributions D i , and are only reordered, their (marginal) distribution will remain that of D i . However, the joint distribution of the parameters is still not well defined. In particular, this simple shuffling may result in some of the parameters to be positively or negatively correlated with each others, which might be undesirable. Some techniques have been developed to eliminate these correlation terms or to impose different correlations between the variables, and will be presented on section 2.2.
It should be noted that, in the mathematical literature, it is usual to refer to a somewhat different object as a Latin Square: this would be a square whose sides are divided in N intervals, and is filled with N different symbols, such that for each row and column there is exactly one occurrence of each symbol, as represented in figure 5.

Algorithms and extensions
As described above, the LH sampling generates an uniform distribution of samples in each parametric dimension. However, there is no guarantee that the correlation between two or more parameters will be zero, and the classical algorithm from McKay [24], described in the previous section, usually produces correlations as high as 0.3 between pairs of factors, which can difficult or even compromise further analyses. In this section, we will present one algorithm designed to take into account the correlation between the parameter variables [16], using a single-switch-optimized sample reordering scheme. We will present the general case of prescribing a correlation matrix, and will also present results for the trivial case of zero correlation terms. Other methods have been proposed to address this problem [47,40,10], including methods that deal with higher-order correlation terms [41] using orthogonal designs. These methods, however, either impose severe restrictions on the number of samples that must be chosen or are too computationally intensive.
In order to obtain the samples with prescribed correlation terms, we define as C i,j the desired m × m correlation matrix between the variables x i and x j , and denote by C * i,j the current correlation between x i and x j . The next step is done iteratively for each parameter dimension, starting with the second one. Suppose that the method has already been applied to i = 1, 2, ...., l − 1, and we will apply it to i = l. The square sum of the errors in the correlations between x l and the anterior parameters is given by Afterwards, we calculate, for each pair of values sampled from the parameter dimension l, what would be the error in the correlation if they were switched. The pair that corresponds to the greater error reduction is then switched, and the procedure is repeated iteratively until the error is acceptably small.

Stochastic models
When dealing with stochastic models, like several relevant individual based models (IBM), the questions presented become complicated by the fact that running the same model with exactly the same parameters might wield largely different results, both quantitative and qualitatively. In this scenario, we must be able to diferentiate the variation in responses due to the variation of the parameters with the variation in response due to stochastic effects.
We will proceed first by defining the variation due to the input parameters as the epistemic uncertainty. This uncertainty arises from the fact that we do now know what are the correct values for a given parameter in a given natural system, and is related to the probability distributions D i , presented in section 1.2. The variation in the behaviour of the model which is caused by stochastic effects is called stochastic uncertainty, and is inherent to the model.
It is important to note that the two uncertainty components are impossible to disentangle in general stochastic models. This has prevented the general analysis of such models until recently. In recent years, studies have shown that the important parameters and their effects can be correctly identified by running such models repeatedly for the same input variables and then averaging the output [35], given that the following conditions are respected: • Sample sizes should be large, relative to the aleatory uncertainty.
• The output values should be unimodal, that is, the output values for a given parameter choice should be clustered around a central value.
• The correct analysis tools should be used (as will be discussed on session 3). We refer to this increase in the number of sample points as a sampling refinement, and as it is done on demand after examining the already generated outputs, we refer to it as an adaptative sampling refinement. This technique can be applied while estimating the Symmetrized Blest Measure of Association (SBMA) which will be used to determine the optimum N size (section 2.5).

Adaptative refinement
Our construction will start by defining a new sequence of values γ i for every parameter dimension, with For every k ∈ [1, N ], and so γ i has 2N elements. Following this, Huntington & Lyrintzis' algorithm [16] is applied on the set of γ i in order to minimize the correlation between them. The extended Latin Hypercube is formed by concatenating the original x i values to the γ i values, and obtaining a cube with 3N samples defined as z i , of which just 2N must be calculated.
The extended Latin Hypercube presents the following desirable properties: • The expected sample mean is equal to the distribution mean for each parameter: z i = D i • The expected sample variance is equal to the distribution variance for each parameter: σ(z i ) = σ(D i ) • The correlation between any two variables is bounded by the following expression:

Measuring the concordance with increasing sample size
We will now turn our attention to the problem of determining the optimal number of model runs we should apply in order to provide a good estimate of which are the relevant parameters for a given model. One way of proceeding is by sistematically increasing the number N of model runs and applying any of the sensitivity analysis techniques, which will be discussed on the following section. If the analyses indicate similar results for consecutive runs, we can presume that increasing the sample size will not yield major changes to the results. All of the sensitivity analyses present us with a list of the parameters that have most influence in the model output. By comparing the resulting lists from two experiments, we can decide to stop increasing N when the lists are sufficiently similar. Our problem then is to determine how similar are two vectors of ranks. In principle, we could apply any distance function to those vectors. However, consider that 3 analyses indicated that the order of the most influential parameters is: H1 = 1 2 3 4 5 6 H2 = 1 2 3 6 4 5 H3 = 2 3 1 4 5 6 By using standard distances (like Spearman's rho or Kendall's tau), we will see the same difference between H1 and H2 and between H1 and H3. On the other hand, in the context of determining the most influential parameters, we would be inclined to see H1 and H2 as more similar than any of them to H3, as the first two preserve the ordering of the three first parameters.
Iman and Conover proposed a correlation coefficient for this problem called Top-Down Correlation Coefficient [18], which is based on Savage Scores. This coefficient, know as TDCC, was extensively used for sensitivity analyses [32]. Another measure of concordance proposed more recently is the Symmetrized Blest Measure of Association (SBMA) [12]. Recent research suggests that estimates for SBMA produce a smaller standard error than TDCC without the assumption that there is no correlation between the variables [23]. Thus, we propose using SBMA as a measure of concordance between analyses from different sample sizes. Defining the ranks from the first sample as R i and the ranks for the second sample as S i , the estimator for the SBMA is: For the techniques that may output negative values, for instance negative correlations, the SBMA must be applied on the absolute values. Otherwise, the parameters which present strong negative effects will be ranked very low, and will not be taken into account by the SBMA.
We will apply SBMA for the PRCC and eFAST techniques (both of which are discussed on section 3) on the example on section 4.

Box 1: Parameter space sampling
In order to adequately explore the parameter space of a model, the following steps should be taken: • The parameters of interest must be identified, in a way that all combinations of parameter values are meaningful.
• The range and distribution of each parameter should be estimated.
• The model to be analysed must be formulated as a function of all the parameters.
• A Latin Hypercube should be generated based on the parameter distributions.
• The correlations between parameters, when present, must be taken into account.
• The model must be evaluated at each combination of parameter values. If the model is stochastic, it should be evaluated repeatedly at each combination.
3 Quantitative output analysis

Uncertainty analysis
The first question we would like to answer, in the context of quantitative analysis, is what is the probability distribution of the response variable y given that we know the join probabilities of the input parameters x (see definitions in section 1.2), which is the subject of uncertainty analysis [15]. This can be done by fitting a density curve to the output y or an empiric cumulative distribution function (ecdf). If there is any theoretical reason to believe that the distribution of y should follow one given distribution, it is possible to fit this function to the actual output data and estimate the distribution parameters. If the joint distribution of the input parameters correspond to the actual probability of some natural system to exhibit some given set of parameter values (as opposed to the case where we have no biologically relevant estimates for some parameters), the estimate represented by the density and ecdf functions approaches the actual distribution that the variable y should present in nature. This functions may be used, for example, to provide confidence intervals on the model responses.
However, this is only the case when the input variables are uncorrelated or when enough correlation terms have been taken into account. Smith ([38]) provides an example where ignoring the correlation terms leads to inaccuracies on the estimation of confidence intervals.
The next reasonable step is to construct and interpret scatterplots relating the result to each input parameter. These scatterplots may aid in the visual identification of patterns, and altough they cannot be used to prove any relationship between the model response and input, they may direct the research effort to the correct analyses. There are extensive reviews of the use of scatterplots to identify the important factors and emerging patterns in sensitivity analyses [19].
We will present here some quantitative analyses tools, aimed at identifying increasingly complex patterns in the model responses. It should be stressed that no single tool will capture all the relations between the input and output. Instead, several tools should be applied to any particular model.

Sensitivity analysis
The question of "what is the effect of some combination of parameters to the model output" may be answered by testing the relation between the parameters and outputs. There are extensive reviews about detecting these relations after generating samples with Latin Hypercubes [32,19], so we will give just a brief overview. We will first note that the methods used must take into account the variation of all the parameters. For example, instead of calculating the correlation between the result and some parameter, partial correlation coefficients should be used, which discount the effect of all other parameters.
These relations between the result and the parameters may be classified, in order of increasing complexity, as: • Linear relation, which can be tested with the Pearson partial correlation coefficient, • Monotonic relation, which can be tested with the Spearman partial correlation coefficient, • Trends in central location, for which the Kruskal-Wallis test may be applied, • Trends in variability, for which the FAST method and Sobol' indices may be used.
Subsections 3.2.1 to 3.2.4 will provide some mathematical background for each method, and section 4 will present examples of use of those tests. We should stress here that the application of one method is not enough to draw conclusions about the relations between the input and output variables, as these techniques test different hypotheses, and have different statistical powers. Instead, every model should be analysed by a combination of techniques, preferably one for each category outlined here.

Linear relation
Under the hypothesis of independence between the central location and dispersion of the model responses, the most straightforward relationship between y and x i is the linear, represented by y ∼ x i . This is case if, every time x i is increased, y increases by aproximately the same ammount. The Pearson correlation coefficient is the commonly used measure to test for a linear correlation: Where σ a is the variance of a and σ ab is the covariance between a and b. The correlation coefficient is a measure of the predicted change in y when x i is changed one unit, relative to its standard deviations, and, as such, approaches ±1 when there is a strong linear relation between the variables. The square of ρ, usually written as R 2 , measures the fraction of the variance in the output that can be accounted for by a linear effect of x i . Is is usual to test the significance of this linear relation by a t-test [11].
Other than examining the individual relationships between the parameters and the output, we can investigate the joint effect of several x i , as y ∼ x 1 +x 2 +· · ·+x m . In this case, the multiple R 2 represent the fraction of the variance on the output due to linear effects of all the x i considered.
However, a measure of ρ close to zero does not mean that no relationship exists between y and x i -for instance, x 2 + y = 1, x ∈ [−1, 1] presents ρ = 0, so clearly other methods might be needed.
The Partial Correlation Coefficient (PCC) between x i and y is the measure of the linear effect of x i on y after the linear effects of the remaining parameters have been discounted. In order to calculate the PCC, first we fit a linear model of x i as a function of the remaining parameters: A corresponding model is done with y: The PCC is calculated as the correlation between the residuals of these two models:

Monotonic relation
Let us refer to each value of y as y k and each value of x i as x ik . The rank transformation of y, represented by r(y k ) can be found by sorting the values y k , and assigning rank 1 to the smallest, 2 to the second smallest, etc, and N to the largest. The rank of x ik , r(x ik ), can be found in a similar way. If there exists a strictly monotonic relation between y and x i , that is, if every time x i increases, y either always increase or always decreases by any positive ammount, it should be clear that the ranks of y and x i present a linear relationship: r(y) ∼ r(x i ).
The correlation between r(y) and r(x i ) is called the Spearman correlation coefficient η yx i . The same analyses presented on section 3.2.1 can also be applied for the rank transformed data, including significance testing and multiple regression.
If the procedure described to calculate the PCC is followed on rank transformed data, that is, if y and x i are rank transformed and fitted as linear models of the remaining parameters, the correlation between the residuals is called PRCC, or Partial Rank Correlation Coefficient. This measure is a robust indicator of monotonic interactions between y and x i , and is subject to significance testing as described in [32]. This measure will perform better with increasing N .

Trends in central location
Even if the relation between y and x i is non monotonic, it may be important and well-defined. The case in which y ∼ x 2 i , x i ∈ (−1, 1) is a common example. This relation may be difficult to visualize, and sometimes may not be expressed analytically. In these cases, the Kruskal-Wallis rank sum test may be used to indicate the presence of such relations [19].
In order to perform the test, the distribution of x i must be divided into a number N test of disjoint intervals. The model response y is then grouped with respect to these intervals, and the Kruskal-Wallis test is used to investigate if the y values have aproximately the same distribution in each of those intervals. A low p-value for this test indicates that the mean and median of y is likely to be different for each interval considered, and thus that the parameter x i have a (possibly non monotonic) relationship with y.
The number of intervals N test is not fixed as any "magical number", and may have a large impact on the test results. It is then recommended that this test should be repeated with different values to obtain a more comprehensive picture of the interactions between x i and y ( fig. 6).

Trends in variability
Other than the central tendency of the results, their dispersal may be dependent on the input parameters. A classical approach which may be used to test whether the dispersal of the output is related to any input parameter is to divide the distribution of x i into a number N test of disjoint intervals and group the model response y with respect to these intervals, as done on the Kruskal-Wallis test. In this case, the ANOVA F statistic can be used to test for equality between the y conditional to each class [19].
A similar approach, which we will use, is to employ the eFAST indices [34,33], which is a variance decomposition method based on the FAST and Sobol' indices. While the Sobol' indices were described in 1969, in Russian, FAST was developed by Cukier et al. in 1973, and both are identical in all but one computation [1].
These methods estimate what fraction of the output variance can be explained by variation in each parameter x i , which is called the first-order sensitivity of x i or main effect of x i . The method estimates as well the fraction which is explained by the higher-order interactions between x i and all other parameters. The sum of all terms related to x i is called the total-order sensitivity of x i .
The eFAST method estimates the main effect of each parameter by chosing a periodic function f i (x i ) for each parameter, where the frequency φ i of each function is distinct, and should, in theory, be incommensurable. Each of this functions is sampled N s times, and a Fourier analysis is applied to the model output. The Fourier coefficients at each frequency φ i is related to the main effect of the variable x i . The total order sensitivity of x i is then calculated as the fraction of the variation which is not explained by the complimentary of x i (that is, all parameters but this one).
There are two things that should be noted here about eFAST indices. The first one is that the eFAST calculation does not involve the LHS sampling scheme, and may require more model evaluations. Also, this method produces small positive total-order sensitivity estimates even for parameters which do not play any role on the model output, as many numeric approximations are involved.

Box 2: Uncertainty and sensitivity
The following questions should be addressed in uncertainty and sensitivy analyses (between parentheses are the section numbers which detail each step): • What is the range and distribution of model responses (3.1)?
• The ecdf plot of model responses suggests any theoretical distribution (3.1)?
• The scatterplots between each variable and the model response suggest any relationship (3.1)?
• Is there any linear or monotonic relation between the model parameters (3.2.1 and 3.2.2)?
• Is there any trend in central location or dispersion of the response (3.2.3 and 3.2.4)?
• The output from the analyses is coherent between different sample sizes (2.5)?
4 Case study: a density-dependent population matrix model

Model description
In this section, we demonstrate the uses and advantages of the methods outlined in the previous sections by performing sensitivity analyses on a density-dependent model of the tropical palm Euterpe edulis (commonly known as palmito juçara). All the data used here was extracted from Silva Matos et al. paper [37], which compared a density independent matrix model of population growth with a density dependent model in which the recruitment of seedlings was affected by the number of seedings and adult trees. Silva Matos provided results, sensitivities and elasticities for the density independent model that can be compared to our findings, and results for mean and maximum values of the density dependent model -but unfortunately, their methods did not allow for a full sensitivity analysis of the density dependent model. We have used the R language to perform the sampling and analysis, with most of the code based on the "sensitivity" package, which implements PRCC analysis (see section 3.2.2) and eFAST analysis (see section 3.2.4), among others. We have also used our implementation of Huntington & Lyrintzis' algorithm to generate zero correlation between LHS samples (see section 2.1). All code used is freely available on the web.
The models analysed are based on a Lefkovitch matrix with seven size classes. The matrix used on the density-independent model is Here, P i is the probability of a tree surviving and remaining in the same class (statis), G i is the probability of a tree surviving and growing to the next class, and F i is the number of offspring produced per reproductive palm.
The dominant eigenvalue of this matrix is related to the predicted population growth rate. We considered the dominant eigenvalue for this matrix as the model output, as usually done on this modelling approach. The density dependent model used the same matrix, but now the growth term of the first size class represented a decreasing function of the population density: Here, N 1 and N 7 represent the number of seedlings and adults per patch. The parameters G m and a represent the maximum transition rate at low densities and the strength of reduction in G 1 with increasing seedling densities. The remaining parameters z and ρ represent the crown area of an adult tree and the plot size (which is fixed as 25m 2 ), and their ratio is related to the reduction of recruitment due to the presence of adults, due to the fact that few seedlings are able to grow underneath the canopy of an adult.
As this model does not produce any static matrix, it is not meaningful to calculate any eigenvalue. Instead, the total population corresponding to the stable population distribution was used as model output.
A naïve approach to estimating the parameter sensitivities of this model would use the stasis, growth and fecundities given. However, this would yield erroneous results, as the probabilities of stasis and growth for a given class are not independent, as P i + G i ≤ 1 for all classes. As discussed on section 1.1, we need to use an alternative parametrization for this model.
We will represent by s i the probability of survival for each class, calculated as s i = P i + G i , and by lowercase g i the probability of growth, calculated as g i = (s i − P i )/s i . Using the notation for complementary probabilities g i = 1 − g i , we can now write the Lefkovitch matrix as: The models have, respectively, 14 and 16 parameters. All analyses have been done with mean and standard deviation calculated from Silva Matos paper, assuming a normal distribution of parameters truncated at the [0, 1] interval for probabilites, and on [0, +∞) for the other parameters. In the case of the density dependence parameters Gm, z and a, only the mean estimate was given on the paper, so conservative values were used for the standard deviations.

Results
First, we have generated Latin Hypercubes consisting of all relevant variables for each model. Then, the models were run for each combination of parameters. By using the SBMA measure of concordance (section 2.5), we have determined that the sample size required for the density independent model is approximately 100, and between 400 and 500 for the density dependent (Table 1).
If we presume that the data collected is representative of our knowledge about each of these parameters, the probability distribution of the model responses can be seen as the probability that the real population of palms exhibit each value of the model output. Figure 7 shows these distributions, which suggest that the population is viable for the vast majority of parameters values in the parameter space considered. Also, the λ calculated from the density-independent model (mean 1.22, standard deviation 0.06), is very close to the value found by Silva Matos (mean 1.24 ± 0.06 se). Considering the density-dependent model, the median stable population predicted (4902 trees in each 25m 2 plot), is comparable to, altough higher than, the population actually measured by the study (1960 ± 560 trees per plot, mean and sd calculated over three years).
We have generated scatterplots between the result from the models and each independent parameter, in order to visually identify the relations between the inputs and outputs (figs. 8 to 11). It is clear from these scatterplots that the fecundity plays a major role on the population dynamics, and may be involved in non-linear interactions. Also, growth probabilities (g i ) have a greater impact on the model output than survival (s i ) on the density-independent model. This is to be contrasted with Silva Matos results, which show all of the elasticities to be aproximately equal for all parameters. In the density-dependent model, the patterns are much more complex. Survival parameters seem to be more influent than growth, and the parameters reducing the recruitment (a and z) show a clear negative effect on the population size. However, there is evidence now for non-linear effects of the parameters, in particular s 1 and F 7 .
These scatterplots show very high dispersion of values, mostly due to the fact that all parameters are being varied between runs. In order to investigate the effect of each parameter on the outputs discounting the effects of the others, we analyse the Partial Rank Correlation Coefficient (PRCC, fig. 12). The PRCC analysis for the density independent model is in agreement with our previous expectations, with F 7 being the most influential parameter, followed by growth probabilities. Survival probabilities follow with low correlations. The density dependent model presents us with some surprises, as now the survival parameter for the smallest and largest size classes jumped to occupy the second and third largest positive PRCCs. The remaining parameters follow the new parameters a and z, which are strongly negatively correlated with the output.
The last analysis we present here is the Extended Fourier Amplitude Sensitivity Test (eFAST, fig. 13). This analysis provides an estimation of the fraction of variation of model output that can be explained by the individual variation of each parameter (which we call first-order sensitivity, or main effect), along with the total     Figure 13: eFAST analysis for the density independent (a) and density dependent (b) models. The bars represent the first and total order estimates for the sensitivity of each parameter in the model output.
variation caused by interaction between that parameter and others (total-order sensitivity). The interaction term for each parameter is the difference between its first and total order sensitivities. The eFAST analysis is substantially more intensive in terms of computer time than the previously mentioned. The analyses presented here required 7392 and 8448 runs of the simulations, respectively. Table 2 shows the SBMA measure of concordance between different sizes. Note that N s reported should be multiplied by the number of parameters in each model to obtain the total number of simulations executed. Notice that, while the main effect D i converge for both model, the total order D t indices are still variable with large sample sizes.  This analysis reveals that the output of the density independent model is mostly explained by first-order relations (which explain 62 % of the output variation), with F 7 and the growth terms being the most important. The importance of the linear terms shouldn't come as a surprise, as the matrix growth model is a linear model. The density-dependent model, which has 66 % of the output variation predicted by linear terms, exhibit more complex interactions between the terms, in particular interactions involving s 7 , g m and a, the survival for adult plants, and the terms associated with competition between seedlings, respectively.

Conclusions
The results from the uncertainty and sensitivity analyses presented show some of the advantages from the methodology described in this work that are unavailable to the usual framework used in ecological studies. First, we have been able to quantify the uncertainty in the assymptotic growth rate (related to λ) and the stable population size related to the uncertainty in the model inputs. Also, we provide a common framework to investigate side-by-side the linear and non-linear matrix models, from which we were able to point out the similarities and discrepancies between the models. Analyses based on matrix elasticities are also unable to investigate the role played by parameters not directly present on the matrix, as the size of the adult trees canopy z in our case. Finally, our approach permits the identification and quantifiation of relative importance of non-linearities and interactions between the input parameters in determining the model's outcome, and allows us to incorporate our previous knowledge about the system in specifying the range and distribution of each input parameter.

Previous use of Latin Hypercube in ecology
The importance of sensitivity and uncertainty analyses in the development and use of ecological models is widely recognized. Searching for the terms "Sensitivity Analysis, Uncertainty Analysis or Parameter Space Exploration" in the Web of Knowledge reports 1199 papers in eight major journals (see fig. 14 and legend for details) since 1971, with 31650 total citations. However, most of these papers rely on full and individual parameter space exploration, which, as discussed in section 2, are not optimal. When restricting those results with the keywords "Latin Hypercube, MCMC, Markov, Monte Carlo", just 120 works show up in the results. Of those, only 13 (about 1%) use Latin Hypercube Sampling [3,6,25,43,45,26,36,7,30,32,28,13,21]. There are also relevant examples of LHS use in other journals [8,9,42].
Also, many of these papers did not explicitly take into account the correlations between parameters. Those who did used mostly Iman and Conover's method [17].
These works have used Latin Hypercubes typically from 10 to 30 dimensions, but ranging from 6 to 143 [3], and the number of simulations ranged from 19 [28] to 2000 [43]. Also, these works are from varied areas within ecology: applied plant ecology [6,43], species richness [13], epidemiology [36] and food chain analysis [7], stressing that the method is useful on varied problems.