Multi-species population indices for sets of species including rare, disappearing or newly occurring species

(MSI) are widely used as ecological indicators and as instruments to inform environmental policies. Many of these indices combine species-specific estimates of relative population sizes using the geometric mean. Because the geometric mean is not defined when values of zero occur, usually only commoner species are included in MSIs and zero values are replaced by a small non-zero value. The latter can exhibit an arbitrary influence on the geometric mean MSI. Here, we show how the compound Poisson and the negative binomial model can be used in such cases to obtain an MSI that has similar features to the geometric mean, including weighting halving and doubling of a species ’ population equally. In contrast to the geometric mean, these two statistical models can handle zero values in population sizes and thus accommodate newly occurring and temporarily or permanently disappearing species in the MSI. We compare the MSIs obtained by the two statistical models with the geometric mean MSI and measure sensitivity to changes in evenness and to population trends in rare and abundant species. Additionally, we outline sources of uncertainty and discuss how to measure them. We found that, in contrast to the geometric mean and the negative binomial MSI, the compound Poisson MSI is less sensitive to changes in evenness when total abundance is constant. Further, we found that the compound Poisson model can be influenced more than the other two methods by trends of species showing a low interannual variance. The negative binomial MSI is less sensitive to trends in rare species compared with the other two methods, and similarly sensitive to trends in abundant species as the geometric mean. While the two new MSIs have the advantage that they are not arbitrarily influenced by rare, newly appearing and disappearing species, both do not weight all species equally. We recommend replacing the geometric mean MSI with either compound Poisson or negative binomial when there are species with a population size of zero in some years having a strong influence on the geometric mean MSI. Further, we recommend providing additional information alongside the MSIs. For example, it is particularly important to give an evenness index in addition to the compound Poisson MSI and to indicate the number of disappearing and newly occurring species alongside the negative binomial MSI.


Introduction
Biodiversity indices are used to synthesise the complex and multidimensional nature of biodiversity (e.g.Mace et al., 2012) into a single figure that measures selected aspects of biodiversity and that is possible to communicate (Duelli and Obrist, 2003;Buckland et al., 2012).The quantified aspects differ between different biodiversity indices and therefore, the choice for a biodiversity index depends on the objectives of a study (McRae et al., 2017;Watermeyer et al., 2021).
Measuring biodiversity and quantifying its change are crucial for taking conservation actions to reduce biodiversity loss (Purvis and Hector, 2000;Loh et al., 2005;Collen et al., 2009;Gregory and van Strien, 2010;Bal et al., 2018).In this way, biodiversity indices support national and global policy decisions (Nicholson et al., 2012).One of the most widely chosen objectives for an index is monitoring of populations because the knowledge of population changes is essential to understand the drivers of population dynamics and for taking any remedial conservation actions (Nichols and Williams, 2006;Roy et al., 2019).Therefore, many biodiversity indices are based on species population sizes such as the Living Planet Index (Loh et al., 2005), the European Butterfly Index (van Swaay et al., 2015), the UK Priority Species Index (Eaton et al., 2015), many Wild Bird Indices (Gregory and van Strien, 2010;Buckland et al., 2012;Gregory et al., 2019), and many national wide indices based on populations or distribution of species.These multi-species indices (MSI) combine relative population sizes (i.e.species indices) across species to measure an average population trend.
Many MSIs like the Living Planet Index, the Wild Bird Index and many country-specific Wild Bird indices are geometric means of relative population sizes of constituent species (Buckland et al., 2005;Gregory et al., 2005;Gregory and van Strien, 2010;Buckland et al., 2012).An MSI based on the geometric mean has three main properties.First, relative changes are weighted equally among species, i.e. if the index of one species declines to 50% ("halves") and the index of another species increases to 200% ("doubles") and everything else stays constant, the MSI stays constant (Buckland et al., 2012;van Strien et al., 2012).Second, when the population index of a single species is zero, the geometric mean, and consequently the MSI, is not defined.Third, random noise in relative population size is typically larger in rare species than in abundant species.Therefore, random variation in rare species can strongly influence the geometric mean MSI (Buckland et al., 2005;Lamb et al., 2009).
Issues arising from the latter properties usually are resolved by selecting commoner species as for Wild Bird Indices in some countries (e. g.Szép et al., 2012) and in the supra-national Wild Bird Indices (Gregory et al., 2005;Gregory et al., 2019).Alternatively, a small value is added to the whole time series (e.g.Loh et al., 2005), zeros are replaced by a small positive value (e.g.Zbinden et al., 2005;O'Brien et al., 2010), or truncated at a small positive value (e.g.Soldaat et al., 2017).In these situations, the value chosen for addition, replacement or truncation has an arbitrary and potentially strong influence on the resulting MSI.
Selecting species according to their commonness or to data availability likely results in a non-representative sample of species included in the MSI (Gregory and van Strien, 2010;Renwick et al., 2012;Buckland and Johnston, 2017;McRae et al., 2017).Representative samples of species according to ecological criteria (Butler et al., 2012;Wade et al., 2014), or including all species using a similar habitat or area (Zbinden et al., 2005) will inevitably contain rare species.Moreover, oncecommon species can become rare and vice versa as time series get longer (Ricklefs and Bermingham, 1999).Therefore, if we include only species into an MSI that were present or even abundant during the reference (first) year, the MSI will be biased towards a negative trend because some species will naturally disappear, and newly occurring species are not included.As a result, the sample of species becomes nonrepresentative.
Alternatively, zero values could be treated as missing values that are either ignored or imputed.When ignoring these values the index represents a cumulative trend over the years averaged across those species present per year, e.g. as done for the Canadian Species Index (Marconi et al., 2021).In consequence, such an index does not allow direct comparison of the relative population size today with the one, e.g. in the reference year, because species composition changes across the time series.Therefore, if long-term comparison of average relative population sizes is an objective of the MSI, zero values should not be ignored.Recently developed workarounds include replacing zeros by truncation according to Soldaat et al. (2017), by predicted values from a population model (e.g.Freeman et al., 2021) or by smoothed species indices (Rosenberg et al., 2019).However, the longer the sequence of zero values, the more difficult it gets to replace them with any of these methods.To summarize, when the goal of the MSI is a comparison of the average relative population size between a specific year (e.g. the most recent year) with a reference year (Fraixedas et al., 2020), zero values cannot be ignored and imputation by a model may be difficult.However, we do not know of any MSI calculation method that allows to have zero values within the set of values from which the average is calculated.
Our aim was to find a method for calculating an MSI that can deal with population sizes of zero, that serves for comparison of average relative population sizes across years, and that can be used as a policy tool.First, we explain why we chose the compound Poisson model and the negative binomial model.We then introduce these two methods and compare them with the standard geometric mean MSI with respect to their sensitivity to changes in evenness (i.e., the distribution pattern of population sizes across species), their weighting of species and sensitivity to trends of single species using virtual and real data.In the discussion section, we explain how uncertainty of species indices can be propagated to the MSI and how different sources of variance may be assessed.Finally, we recommend when and how the new methods may be useful.

Overview and why we chose the compound Poisson and negative binomial model
There are many different paths to go from species-specific population sizes or another species-specific abundance metrics to an MSI representing average relative population sizes (Fig. 1).Most commonly, the geometric mean is applied to relative population sizes (path 1a in Fig. 1) resulting in an average relative population size that weights species equally independent of their abundance.Note that an MSI based on the log-normal model is equivalent to the geometric mean (see 2.2.).Whilst equal weighting of species is generally seen as a desirable property of MSIs, unequal weighting is sometimes applied to address unrepresentative taxonomic or geographic coverage (McRae et al., 2017).Equal weighting across species is also achieved by other methods based on relative population sizes (species indices) like the log-normal growth model (1b), the compound Poisson model (1c), or indices based on smoothed species-specific relative population indices (2).In contrast, most methods based on absolute species-specific population sizes (e.g.counts or estimates; paths 3 and 4) result in MSIs that weight individuals equally or give individuals much more weight compared to species.Such MSIs are proportional or approximately proportional to the arithmetic mean of the population sizes, and therefore mostly represent the trends in the abundant species.However, species identity and the persistence of species in a community is a critical aspect of biodiversity change.Therefore, many biodiversity indices aim at measuring average relative population sizes across species giving relative changes for each species similar weights as is done by the geometric mean (Gregory et al., 2005;Buckland et al., 2011;van Strien et al., 2012).Of the methods applicable to absolute population sizes and considered in this study, only the negative binomial model including additive effects of year and species (path 3b1) resulted in fitted values that weight species equally.
For this study we did not explore more complex model structures that may include linear and non-linear trends, ecological predictors, speciesspecific among-year variances, or interaction between species, because it was not our aim to study mechanisms of biodiversity trends but our aim was simply reporting of changes in biodiversity, as needed for policy.
We finally identified the compound Poisson model fitted to relative population sizes (path 1c) and the negative binomial model fitted to absolute population sizes (path 3b1) as suited for further exploration and comparison with the geometric mean (path 1a).These two approaches were chosen, because they allow the inclusion of zero values, they result in an average relative population size and they do not allow for tuning by the data analyst, such as choosing a degree of smoothing, choosing a threshold, or a value for replacing zero values.

Descriptions of the geometric mean MSI and the two MSIs allowing for population sizes of zero
MSIs are based on predicted, estimated or counted population sizes for species (i) and years (t) for a specific area such as a country.We first describe the three methods assuming (unrealistically) the true population sizes N it are known and we explore in the discussion error propagation and effects of using predicted, estimated or counted population sizes.For the geometric mean and compound Poisson indices, relative population sizes, y it = Nit Ri , relative to a species-specific reference R i were used.The reference R i can be the population in a specific year, e.g.R i = N i2000 , a long-term average or a target population size of species i.To keep comparability across years, we do not change species composition along the time series.Absent species (e.g. in years when colonising species have not arrived yet, or in years after a species has disappeared) were retained at a population size of zero, except for the geometric mean MSI, where zero values and values close to zero were replaced by truncating at 1% of the reference population size.As reference population we used the population size in year 2000 for most species, except for two species with strongly increasing and decreasing population sizes for which we used the year with their maximum population size as recommended by Soldaat et al. (2017).
1) Geometric mean MSI: The geometric mean index is defined as ∑ S i=1 ln(yit ) , where S is the number of species.In words, the geometric mean index for year t is the exponential of the arithmetic mean of the natural logarithms of the relative species-specific population indices in year t.The geometric mean can be derived by assuming a log-normal distribution for y it , i.e. ln ( . Taking the exponential function of the maximum likelihood estimator of the mean equals the geometric mean of y it .Therefore, the geometric mean index I geom t can be calculated by fitting a normal linear model to ln(y it ) with independent means β t for each year: where μ it is a linear predictor.In the simplest case μ it = β t .When including species as a factor, μ it = β t + α i=2,⋯,S .The exponentials of the year-specific intercepts β t referenced to a specific year equals the geometric mean index: , where e.g.R = e β 2000 .Note that the same index I geom t results from the two versions of the linear predictor (including or excluding species as a factor).
2) Compound Poisson MSI: The family of the Tweedie distributions is part of the linear exponential family with a dispersion parameter where the variance is a power of the mean.Tweedie distributions have mass at zero and support on the non-negative reals (Dunn and Smyth, 2005), which means data points can take on real values greater or equal zero.If the power index of the density function is 1, the distribution is a Poisson distribution, if it is 2 the distribution is a gamma distribution.For power indices between 1 and 2, the distributions represent the sum of N gamma distributed variables where N itself follows a Poisson distribution.These distributions have been called compound Poisson (e.g.Smyth and Jørgensen, 2002), compound gamma (e.g.Johnson and Kotz, 1970) or Poisson-gamma distributions (Dunn and Smyth, 2005).We call the distribution compound Poisson distribution following Zhang (2013).Compound Poisson distributions are typical for means of counts (Swallow et al., 2016).Average (e.g., over many sites), expected and relative population sizes are continuous variables that can take on the exact value of zero.When the compound Poisson model is fitted to relative population sizes, y it ∼ compoundPoisson(μ it , p, φ), a geometric meanlike MSI can be obtained.The model parameter p is the power index of the variance function and φ is the dispersion parameter.The natural logarithm is used as a link function and the compound Poisson MSI is obtained as the exponential of the year-specific coefficients of the linear predictor divided by a reference, ln(μ it ) = β t + α i=2,⋯,S , and , where e.g.R = e β 2000 .As for the log-normal model above, the index I cP t , is invariant to including or excluding the species in the linear predictor.
3) Negative binomial MSI: The negative binomial distribution is a discrete distribution allowing for non-negative integer values and is suited for the description of count data that show a higher variance than expected by the Poisson distribution (i.e., overdispersion).We fit a negative binomial model to the population sizes to obtain an MSI: N it ∼ negative − binomial(μ it , θ) with μ it being the expected value and θ being the variance.As for the compound Poisson and lognormal model, we use the natural logarithm to link the linear pre- The MSI is, as for 1) and 2), obtained as the exponential of the year-specific coefficients divided by a reference I nb t = e β t R .If population sizes are estimated, simulated values from the (discrete) predictive distributions of N it are used and the uncertainty of the population estimate propagated to the MSI via Monte Carlo simulation as explored in the discussion.
When species is not included as a predictor variable in the negative binomial model, i.e.,ln(μ t ) = β t , the resulting exponentials of the year-specific coefficients are similar to the arithmetic mean of the absolute population sizes (path 3b2 in Fig. 1, Figure S1).Therefore, species must be included as a predictor variable in the model to ensure that the MSI measures average relative population size.
The three distributions differ in their shape (Fig. 2).On the logarithmic scale, the log-normal distribution is symmetric whereas the compound Poisson and the negative binomial distributions are leftskewed.Of those two distributions, the compound Poisson distribution shows a stronger degree of skewness.The degree of skewness for the negative binomial model depends on the expected population size: For large population sizes, the distribution is less skewed compared to low population sizes.

Comparison between the three models 2.3.1. Exploration of MSI performance using virtual data
To assess how the three MSIs perform in relation to changes in evenness, we used virtual data of a community of five species over four years with a constant total population size of 10 ′ 000.All five species had a population size of 2000 in the first year (maximal possible evenness).Over the following three years, we gradually decreased the population sizes of two species, and we increased the population sizes of two other species, while the sum of all population sizes was held constant, so that population sizes in the 4th year were 250, 500, 2000, 3500 and 3750, respectively.Thus, the virtual data showed a decrease in evenness over the four years.We applied the three MSI methods to these virtual data.
To assess how sensitive the three MSI methods are with respect to species showing strong trends, we simulated nine different sets of data.Each data set contained population sizes for 20 species over ten years.For this simulation, we assume that the true population size is known (not estimated or counted at a sample of sites).In each data set the 20 species were divided into two groups.The first group contained species showing a strong annual (multiplicative) trend of λ = 1.1, whereas the trend of the species in the second group was a moderate decrease (λ = 0.97).In all data sets, the species showing a moderate trend had an average population size of N i1 = 1000 in the first year.The first-year population size of the species showing a strong trend was 10, 1000, and 100 ′ 000 in each of three data sets, respectively.We further varied the number of species showing a strong trend, so that three data sets each included 1, 5 or 10 species with a strong trend.Thus, we had nine data sets that varied in abundance and number of species showing a strong trend.The population sizes N it were simulated as random numbers from an overdispersed Poisson distribution with means equal to N i1 *λ (t-1) *e δ , where δ is a random number from a normal distribution with a standard deviation of 0.1.Data simulation was done so that all population sizes were non-zero and positive, so that the comparison between the two new MSIs and the geometric mean MSI was not hampered by the value used to replace the zero values for calculating the geometric mean.

Application to real data and assessment of species influences
We used four MSIs from the Swiss Bird Index (Zbinden et al., 2005) to explore and illustrate the three models.These four MSIs were chosen because they differ in the number of species, the number of disappearing and newly occurring species and the range of population sizes (Table 1).The Swiss Bird Index includes all regular breeding birds, i.e., those who bred in 9 years within a 10-year period at least once during the time Fig. 2. The log-normal (corresponding to the geometric mean) and compound Poisson distributions fitted to the relative, and negative binomial (NB) distribution fitted to the absolute population sizes of an example group of wetland bird species in Switzerland (see 2.3.2).The histogram shows the distribution of the natural logarithm of the relative population sizes y it of the wetland bird species.For all distributions the x-axis was transformed to the natural logarithm of the relative population sizes.For the histogram and the log-normal distribution, prior to taking the natural logarithm, values lower than 0.01 were replaced by 0.01 (ln(0.01)= -4.6).The shape of the negative binomial distribution depends on the expected value (i.e., the species' population size).Dots correspond to the integer values of the untransformed negative binomial distributions.We illustrate the negative binomial distributions for three different population sizes (25, 100, 10 ′ 000).

series.
We used estimated annual relative population sizes y it that exist since 1990 based on a common breeding bird monitoring scheme for all species including a standard error (Knaus et al., 2021).To convert the estimated relative population sizes into estimated absolute population sizes, N it , we used population estimates for all species for the years 2013-2016 (Knaus et al., 2018;Strebel et al., 2019).We then assumed that the average population index for the years 2013 -2016, (y i2013 + y i2014 + y i2015 + y i2016 )/4, corresponds to the population estimates of the breeding bird census.We subsequently transposed all population indices y it into population estimates N it including a standard error using the laws of error propagation (Hosmer et al., 2008).For calculating the geometric mean MSI, we treated zero values by truncating the time series at 1% of the reference population size as described in Soldaat et al. (2017).
We used the Monte Carlo method to propagate the errors from the population indices or estimates to the MSI (Soldaat et al., 2017).To do so, we draw 1000 values from a log-normal distribution with mean y it and N it respectively and the corresponding standard error as scale parameter.Then, we fitted the models to each of the simulated sets of population indices and estimates respectively.Before fitting the negative binomial model, the simulated values were rounded to integers to meet discreteness.Note that we would have preferred to fit the negative binomial model to simulated values of the predictive distributions of population sizes, which would have been discrete distributions.However, then the errors would not have been comparable to the ones of the geometric mean and compound Poisson MSI.For the y it , there are no predictive distributions available (see also discussion on the errors).
As the reference R i , we use the long-term average of the population sizes R i = 1 T ∑ T t=1 N it to aid a graphical comparison between different indices.The choice of R i does not affect the resulting geometric mean based MSI and the reference of an MSI can always be changed post-hoc.
We calculated the MSI for all four groups using each of the three different models and compared the outputs graphically and to the median of the species-specific relative population sizes, a robust measure of the centre of these values.We calculated the median including the zero values.
To measure the influence of population indices of single species on the MSI, we recalculated each MSI as many times as there were species in the group, excluding each species in turn (Gregory et al., 2019).

Exploration of MSI performance using virtual data
When the total abundance (sum of population sizes of five virtual species) is kept constant and evenness decreased, the decrease in the geometric mean MSI and the negative binomial MSI was much stronger than in the compound Poisson MSI (Fig. 3).
In the simulated scenario with only one species out of 20 species

Application to real data and assessment of species influences
In the case of the forest species, a group including only abundant species, the resulting MSI looks similar for all three methods.However, when the species group includes rare, newly occurring or disappearing species, the MSI based on the compound Poisson or the negative binomial model follow more closely the median of the species-specific relative population sizes compared to the geometric mean MSI (Alpine habitats, wetlands, agricultural target species, Fig. 5).The geometric mean MSI shows larger among-year variance than the other indices, particularly when the species selection includes rare, newly appearing or disappearing species.
The influence of single species was stronger in the geometric mean MSI compared to the other two MSIs (Fig. 6).The geometric mean MSI changed substantially (up to 24%, in average 1.4%) when one species was left out from the data.In the compound Poisson MSI, only species with strong relative population changes such as the Bearded vulture markedly influenced the MSI in groups with a low number of species such as the alpine birds (maximum change was 21%, average 1.0%).In the negative binomial model, rare species had lower influence on the MSI than in the others (maximum 9%, average 0.8%).Abundant species with strong population changes, such as the Great Cormorant in the wetland index, influenced the negative binomial MSI to a similar extent as in the compound Poisson MSI.In all panels, the number of species was 20.Indices are fixed to start at 1.In the first column, there was one species showing a strongly increasing annual trend of 1.1, in the second there were five and in the third ten species with the same increasing trend.In the first row, the strongly increasing species start with a low population size of 10, in the second with 1000 and in the third row with 100'000 (number indicated on the right outer axis).All other species with a moderate negative trend of 0.97 started with a population size of 1000 in all panels.The transparent grey lines are the species-specific population indices.

Influence of the distributional assumption on the weighting of species and individuals
There is a general agreement that MSIs representing average relative population sizes should weight population trends of all species equally as is done by the geometric mean (Gregory and van Strien, 2010).The data distribution assumed in the model determines this weighting.For example, in a Poisson model, every individual is seen as an independent observation and therefore weighted equally.Fitted values in a Poisson model of population sizes are equal to the arithmetic mean of the (absolute) population sizes (e.g.Evans et al., 2000).Therefore, MSIs derived from Poisson models weight species proportional to their abundances, which is not desirable for biodiversity indices with the purpose to measure average relative population sizes.
In contrast, the negative binomial model fitted to absolute population sizes, the compound Poisson model and the log-normal model fitted to relative population sizes, in principle, give equal weight to all species.When all species have the same mean abundance and among-year variance, the three models produce the exact same MSI values.However, when species differ either in their abundances or in their amongyear variance, the three methods produce different MSI values, because 1) abundance and variance have different effects on the weight of the species in the different models, and 2) if the data contain zero values, the weight of such species depends on the value imputed for the zero in the geometric mean MSI.
Abundant species showing strong trends (simulated data example, Fig. 4) or abundant species with rather constant populations (Forest species MSI) had a large influence on the compound Poisson MSI, probably due to the narrow shape of that distribution when the among-species variance is low.Maybe for the same reason, this MSI has weak sensitivity towards changes in evenness.The negative binomial distribution has, among the three distributions studied here, the heaviest lower tail.The mass of the lower tail increases with decreasing abundance because variance due to stochasticity increases with decreasing abundances.Therefore, the negative binomial MSI gives less weight to extreme values of rare species than of abundant species.For this reason, the negative binomial MSI of the alpine species index is not overly influenced by the Bearded vulture, in contrast to the other two MSIs.

Interpretation of the MSIs, and how and when to use them
The type of MSIs discussed in this study measure average relative population sizes and they allow for comparison of average relative population sizes across different years.The geometric mean is sensitive to changes in the proportional distribution, i.e. evenness, even when total abundance is constant (Buckland et al., 2011).Similarly, the MSI based on the negative binomial model is sensitive to changes in evenness whereas the compound Poisson index is not.Several researchers have stressed that single indices always omit important parts of biodiversity (Studeny et al., 2011;Buckland et al., 2012;Leinster and Cobbold, 2012;Sattler et al., 2014;McDonald and Dimmick, 2016).The negative binomial MSI has weak sensitivity towards changes in populations of rare species and with the geometric mean MSI only species present in all (or most) years can be considered (Table 2).Therefore, we recommend reporting complementary measures, such as an evenness index (Shannon and Weaver, 1949;Buckland et al., 2005) particularly when using the compound Poisson MSI and when the aim is to measure the overall state of biodiversity of which evenness is an important aspect (Wittebolle et al., 2009).Turn-over rates (Harrison et al., 2016;Yuan et al., 2016) are valuable additional information to all three MSIs.Indicating the number of disappearing species is particularly important if the negative binomial MSI is used because this MSI gives a low weight to rare species (Knaus et al., 2021).In contrast, the compound Poisson MSI seems to weight rare species similarly to the geometric mean MSI except when zero values are present.Because zero values need to be replaced by non-zero values, which then can have a large influence on the geometric mean MSI (Buckland et al., 2005;Lamb et al., 2009), we recommend using one of the other two MSIs when the set of species includes rare, disappearing or newly occurring species and sensitivity analyses shows large influences of these species on the geometric mean MSI.When using the compound Poisson MSI, care must be taken when the selected species include abundant species with strong relative trends.
Excluding single or clusters of ecologically similar species from the MSI elucidates sensitivity of the MSI towards species with particular characteristics (Fraixedas et al., 2020).Differentiated analyses of common trends of clusters of species give a much stronger insight into the patterns of biodiversity change compared with a single index.For example, splitting the Living Planet Index into clusters showed that average trends differ across regions and taxa, with strong declines in larger animals, the Indo-Pacific region and the reptile and amphibian group (Leung et al., 2020).
To conclude, there may not be one single technique to obtain an MSI that is ideal in all situations and reporting one index only can result in hiding important information.

Quantifying uncertainty, error propagation
We will discuss three different sources of variance that we think are the most important ones an MSI is subjected to: 1) sampling error (sampling of species and years), 2) measurement error (uncertainty in estimated species-specific annual population indices or sizes; determined, e.g., by monitoring method or sampling of sites) and 3) amongyear variance in species populations.
1) The species included in the MSI could be seen as a random sample of a bigger population of species and years subject to an ecological process driving their population trends.Sometimes the set of species of interest is finite and can be sampled as a whole, e.g.all regular breeding birds of Switzerland (Zbinden et al., 2005).If all species are included in the sample, then the sampling error may be seen to be zero as done in the examples of this study.Some studies aim at understanding general ecological processes leading to specific trends in populations and therefore see the  monitored species and years as a random sample of all possible species and years.In that case the sampling error is derived from the residual variance of the model.However, the models we presented here have a too simple structure to reliably measure residual variance (Freeman et al., 2021) and, therefore, most likely are not suited to obtain an appropriate error if the MSI is interpreted as an estimate of the underlying ecological process.Theoretically, bootstrapping the sample of species and describing the variance in the resulting point estimates will reflect uncertainty measures, if the selected species are a random sample from an infinite population of species for which the MSI should measure the state (Carpenter and Bithell, 2000).However, in population monitoring programs random sampling of species is very rarely met (Buckland et al., 2005;Buckland et al., 2011).2) Independent of whether a sampling error is included or not, MSIs have a measurement error.This error is propagated to the MSI from the errors of the species-specific population indices or estimated or predicted population sizes.It is determined by the sampling of the sites where counts took place, the kind and intensity of the survey method and by the spatial variance in animal density.Most monitoring programs produce an estimate of the absolute or relative population size including a standard error (e.g.Pannekoek and van Strien, 2005).These errors are then propagated to the MSI using Monte Carlo simulations (Soldaat et al., 2017;this study).However, for the future, it may be valuable to discuss using predictive distributions of N it instead of estimates with SE as the error that is propagated to the MSI.Unfortunately, for getting the predictive distribution, solid data are needed that allow for predicting absolute population sizes.Many monitoring programs allow for robust estimation of population indices but not absolute population sizes (Pannekoek and van Strien, 2005).The negative binomial model could be fitted to the raw counts instead of the estimated population sizes and the proportion of area surveyed included as an offset.In such a model, the measurement error that is due to the sampling of the sites is automatically propagated to the error of the fitted values.The error of the fitted value will, in such a model, also include the error introduced due to the sampling of the species and years.
3) The natural among-year variance in relative population sizes differs among species.Some species (e.g., slow pace of life species) do not vary a lot, whereas others (e.g., fast pace of life) naturally show peaks and lows in their populations (Stearns, 1992;Tieleman, 2009).If the population of a latter species decreases by 10% from one year to the next, we may not be alarmed.However, if a slow pace of life species shows a similar decrease, ecologists should be alarmed.Therefore, we may wish to weigh relative population changes of species differently according to their life-history characteristics.However, it may be more feasible to report MSIs for groups of species with different life-histories rather than to weight species differently within one MSI.

Possible future directions
Using a statistical model for obtaining the MSI opens several possibilities worth exploring.Accounting for phylogeny via the correlation structure will result in an MSI that better represents genetic diversity.Fraixedas et al. (2020) introduce a bird indicator that includes phylogeny.
Potentially, the models presented in this study could also be used to create supra-national or across-taxon MSIs (Gregory et al., 2005;Gregory and van Strien, 2010;Eaton et al., 2015;van Strien et al., 2016;Burns et al., 2018).The negative binomial model only can be applied to numbers of discrete units (e.g.individuals, breeding pairs), whereas the geometric mean and the compound Poisson model can also be applied to biodiversity indices derived from other type of measurements, e.g.biomass or cover (%).How the models would perform for supra-national or across-taxon purpose remains to be explored.
The state-space models recently explored for calculating MSIs (Simmons et al., 2015;Freeman et al., 2021) may be seen as further development of the chain methods that is used for the Living Planet Index (Loh et al., 2005).Their auto-regressive nature recognises that population size depends on its size in the previous year.Instead of the log-normal models usually used for these state-space models, it may be worth exploring the compound Poisson or negative binomial models when rare, newly occurring or disappearing species are included in the sample.
Finally, complex hierarchical multi-species models allow for detection bias and sampling effort within the same model.Such models also allow estimation of general population trends averaged across species that can be used as an MSI (Iknayan et al., 2014).Although in our study we focussed on simple models, we add to the research of complex hierarchical modelling the finding that the interpretation of the resulting MSI depends on the distribution assumed in the model.Particularly, if the Poisson distribution is used, the resulting MSI may reflect the trend in the absolute number of individuals rather than a trend in the average relative population sizes.However, how the different distributions affect the interpretation of the resulting MSI in more complex models remains to be studied.

Conclusions
We presented two statistical models that naturally allow for zero 1) The influence of abundant species on the compound Poisson MSI cannot be due to the abundance per se because this MSI is based on the relative population sizes that do not contain information about the absolute abundance.However, abundant species often have low among-year variance in relative population size, and the general insensitivity towards changes in evenness implies that the compound Poisson model weights larger relative population sizes relatively stronger compared to the other two methods.
2) The large influence of single species normally is assessed by sensitivity analyses and, if present, reduced by truncation of small relative population sizes of single species, e.g.Soldaat et al. (2017).
values in population sizes and thus allow for the inclusion of rare, disappearing and newly occurring species in a relative population MSI.One advantage of the proposed methods over the geometric mean MSI is that they are no longer subject to the arbitrary influence of the values used to replace the zero values.This improvement comes at the cost of low sensitivity towards population changes in rare species (negative binomial MSI) or a high sensitivity towards trends of species with a low among-year variance and low sensitivity towards changes in evenness (compound Poisson).Although we found methods that improve the techniques for calculating MSIs, we did not find a single method that addresses all needs.If the set of considered species contains rare, newly occurring or disappearing species, we either recommend using the negative binomial MSI and indicating the number of disappearing and newly occurring species alongside or using the compound Poisson MSI together with an evenness index.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Different possible paths leading from the species-specific population sizes to multi-species indices (MSIs) either giving equal weight to species or individuals.The gaps between arrow and box indicate methods for which the weighting of species and individuals depends on factors such as population sizes and variances.Paths highlighted in orange (1a, 1c, 3b1) are explored in this study.References: (1a) e.g.Loh et al., 2005, and many others, (1b) Freeman et al., 2021, (2) Collen et al., 2009, Rosenberg et al., 2019, (3a) Freeman et al., 2020, (4) Fewster et al., 2000.Dark blue boxes indicate steps that are not clearly defined such as replacing zero values with a positive value or smoothing.
showing a strong annual trend of 1.1 (Fig.4, a, d, g) all three MSIs follow the decreasing moderate trend of the 19 species, i.e., no MSIs was influenced by the one species with the strong increase.When five out of 20 (25%) species showed a strong trend and they are rare relative to the other species (Fig.4 b), the MSI based on the negative binomial model lay within the population indices of the 15 common species, whereas the MSIs based on the geometric mean and the compound Poisson MSI lay in-between the average of the species with the moderate decrease and the ones with strong increasing trends.The negative binomial model was sensitive to strong trends only when the abundances of the strong species were equal or larger to the abundances of the species with a moderate trend of 0.97 (Fig.4 e, h).In the scenario with half of the species showing a strong increase and the other half a moderate decrease (Fig.4 c, f, i) the MSIs based on the geometric mean and on the compound Poisson model lay in-between the average of the two groups independent of the abundances.The MSI based on the negative binomial model lay in-between the average of the two groups only when the species with a strong increase had equal or higher abundances compared to the moderate trend species (Fig.4 f, i).In cases when the species with the

Fig. 3 .
Fig. 3. Change in evenness in a hypothetical community of 5 species and corresponding MSIs.Upper panel: population sizes of 5 species in a community with decreasing evenness and constant total abundance (10000) over 4 years.Lower panel: MSIs based on the population sizes of the 5 species in the upper panel.

Fig. 4 .
Fig. 4. Results of simulations that shows the influence of species with strong trends for the three MSI methods (black = geometric mean, khaki = compound Poisson model, pink = negative binomial model).In all panels, the number of species was 20.Indices are fixed to start at 1.In the first column, there was one species showing a strongly increasing annual trend of 1.1, in the second there were five and in the third ten species with the same increasing trend.In the first row, the strongly increasing species start with a low population size of 10, in the second with 1000 and in the third row with 100'000 (number indicated on the right outer axis).All other species with a moderate negative trend of 0.97 started with a population size of 1000 in all panels.The transparent grey lines are the species-specific population indices.

Fig. 5 .
Fig. 5.The geometric mean, compound Poisson and negative binomial MSI (lines) and the median (points) of the real data examples.The long-term average is set to 1 for all MSI and the median to ease comparison.

Fig. 6 .
Fig.6.Single species influence on the geometric mean MSI, the compound Poisson MSI and the negative binomial MSI.Black bold lines are MSIs based on all species in the group (number of species: forest 56, alpine habitats 13, wetlands 39, agricultural target species 28), thin red lines are obtained by leaving out one species at a time.

Table 1
Characteristics of the species sets for the four MSIs from the Swiss Bird Index.
′ 560 (274-998 ′ 290) Most species are abundant, no zero values; similar to other forest bird indicators in Europe, e.g.Gregory et al., 2019  ′ 225)Few species, one newly occurring species (Bearded vulture Gypaetus barbatus) showed a strong increase in relative population size after 17 years ′ 683) Rather high proportion of disappearing species, sequences of zero up to 14 years long at the end of the time series; selection of species of conservation interest (target species).Official name of this MSI is "Target Species Environmental Objectives Agriculture"

Table 2
Properties of the three MSI methods as evaluated in our study.