The analysis in CA2013 is based on comparisons of species accumulation curves (Colwell et al. 2012) between periods (Figure 1). These curves express the dependence of the number of species in a sample on a variable representing sampling effort, and the horizontal asymptote of the curve is species richness. The comparisons in CA2013 are between three 20-year periods (1950–1969, 1970–1989 and 1990–2009). Per pair of successive periods, richness change was estimated as the difference between the log transformed numbers of species in the second minus the first period, where these numbers of species were predicted at a sampling effort (number of records) specific to each difference. Differences were calculated per group of species, country or per grid cell at smaller spatial scales. The differences between logged numbers of species were predicted at three times the numbers of records of the least sampled period, either by extrapolation when numbers of records in both periods were smaller than that, or by inter- and extrapolation combined (Fig. 1). The standard deviation of each difference was estimated using a bootstrap approximation. For spatial scales with multiple grid cells, random effects models that use squared standard deviations as the known error variances (Viechtbauer 2010) were fitted to the estimated changes per grid cell. The average effect across grid cells was used as a measure of richness change. As a check of robustness of results, the analysis was repeated with only interpolation (rarefaction) or only extrapolation and also with standard deviations of the logged difference estimated with an analytical expression (Colwell et al. 2012, CA2013). CA2013 inferred that rates of decline have decreased from observing that estimates became less accentuated, or that significant species number increases occur between the most recent periods when there had been a decrease before. Per taxon and country, Table 1 lists the statements from the text that the authors used to support their conclusion of a decelerating decline in species richness in several taxa and countries. All statements in the table can be found in CA2013’s Results section on changes since 1990. When a spatial scale is mentioned in CA2013, I list the spatial scales to which the statement applies. Eight out of fifteen taxon/country combinations have statements in support of a decelerating decline. CA2013 state that this is independent of the way in which they carried out the analysis, hence robust.

Figure 1.

Species accumulation curves. Species richness is the asymptote of a species accumulation curve, which expresses the dependence on sampling effort of the number of species sampled from an assemblage. In CA2013, sampling effort is given by the number of records from which the number of species is calculated. For illustrative purposes, an example with three arbitrary samples (for 10000, 5000 and 2000 records, labeled from one to three) is drawn. For sample one, a predicted species accumulation curve is added that gradually increases from one species sampled to the predicted species richness for that assemblage (full line). Such curves are constructed on the basis of interpolation and extrapolation. For samples two and three only segments of extrapolated curves are drawn (dotted lines). For sample two, a curve that crosses the species accumulation curve of sample one is sketched. For samples one and three species accumulation curves are more or less proportional. The way in which the species richness differences between samples are assessed in CA2013 is illustrated by indicating on the species accumulation curves at which numbers of records pairwise comparisons would be made between two sample pairs (1 vs. 2 and 1 vs. 3). The number of species of the sample with the smallest number of records is extrapolated to the number expected at three times the number of records. When the number of records of the other sample is still larger than that, the number of species of the second sample is interpolated (rarefied), otherwise it is extrapolated as well.

Note that species richness was nowhere estimated in CA2013, rather numbers of species at particular finite sampling efforts were used as proxies for species richness.

As Table 1 shows, CA2013 does not contain any test for differences in rates of richness change between the pairs of twenty-year time periods considered. They did not aim to reject the null hypothesis that species numbers change at a constant rate. One cannot infer a change in species number decrease by checking which confidence intervals overlap with zero change and which ones not. For example if confidence intervals for declines would have been [-2, -0.2] between the first pair of twenty-year periods and [-1, 0.2] between the second pair, these intervals should not be interpreted as evidence of a change in decline as they do overlap, while the estimates would indeed have become less accentuated.

While appropriate tests for a slowing down of richness decline are lacking in CA2013, we can still check ourselves whether confidence intervals for changes overlap between interval pairs, and conclude on the significance of decelerations in species richness from the absence of overlaps.

In part of their calculations, CA2013 have used a bootstrap estimate of the variance of predicted species number. After inspection of the R code used, it turns out that their estimate is neither bootstrap variance nor bootstrap accuracy (e.g., Walther and Moore 2005), so not a regular bootstrap estimate of variance. In their paper as published before correction, they sum the absolute value of the bootstrap estimated bias and the squared bootstrap average of Colwell et al.’s (2012) expression for the standard deviation of predicted species richness. CA2013 intended the bootstrap to account for additional uncertainty caused by potential non-random sampling, thus variances “corrected” in this manner should rarely become smaller than the analytical expression based on multinomial sampling. However, simulations using samples from the EIS bee and bumblebee dataset used in CA2013 indicate that they often do.

CA2013’s script contained a calculation error: the bias on species number in the earlier period is used in calculations where it should have been that of the later period. In a modified script distributed with their corrigendum (Carvalheiro et al. 2013b) and used to produce a modified Figure 1, this error has been corrected. At the same time, however, a second change has been implemented: the bootstrap standard deviation of species number is now calculated as the product of the bootstrap average of the estimate of standard deviation based on multinomial sampling, times a scaling factor which is equal to one plus the absolute value of the ratio of the bootstrap estimated bias in predicted species number divided by the original estimate of species number. In my simulations, this quantity is smaller than the analytical expression in over 90% of samples simulated using random (multinomial) draws from the EIS bee data species distribution. As the CA2013 bootstrap confidence intervals are obtained using non-standard approaches, we know little about their performance in frequentist inference, except for the simulations I mention here which suggest that they have undesirable properties that lead to anticonservative inference.

The unweighted tests in Table S5 of CA2013, where we expect standard deviations to be calculated automatically from the data variances, all use averages of the bootstrap standard deviations as weights. They are therefore all weighted in an unexpected manner and should not be considered.

In Table S5 of CA2013, the three listed tests per taxon/country for effects at the national scale are in fact always the same test pasted in three times, namely the test based on a bootstrapped variance estimate. The R script of the authors does not contain any calculations for non-bootstrap weighted tests for national data, such that there is no robustness assessment possible. An unweighted test is impossible to carry out at the national level, as there are too few values to calculate data variances per species group/country. I am forced to rely solely on CA2013’s bootstrap statistics for the country-level comparisons. Statistics in Tables S2 and S5 of CA2013 and the supplementary figures have not been adapted to the new heuristic to estimate bootstrap standard deviation, we therefore need to inspect the corrected Figure 1 of CA2013 to assess comparisons based on bootstrapped statistics and the uncorrected figures of the supplement (by comparing plotted limits of confidence intervals with a ruler). Table 2 summarizes the different estimates CA2013 provides and problems arising when one wants to use them further. A proper assessment of robustness is difficult, as in fact nearly any category of estimates either has uncorrected errors, is not provided, or not available at the national level. Nevertheless, I will use rarefaction and extrapolation results in my assessment, and will proceed as if their standard deviations had been correctly calculated. When such estimates will affect conclusions, I will note whether the uncorrected standard deviations play an important role in that or whether just the estimates of species number change would lead to the same conclusion.

CA2013 assume that their response variable (log-ratio of predicted species numbers at a particular sampling effort) is directly comparable between different sampling efforts. This is only expected when non-linear species accumulation curves for different periods are proportional and will often not hold good when these curves for example cross (Fig. 1).

CA2013 call it a bias that their response variable often becomes larger with larger differences in data sample sizes between periods. To correct for this presumed bias, they included the difference of the logged numbers of records in the two time periods as covariate in the random effects models. That the difference in species number often becomes larger with larger differences in sample sizes is an absolutely normal pattern if species accumulation curves are not proportional and sample sizes where species numbers are interpolated are not randomly distributed (for example consistently smaller in one period). It is not difficult to sketch a pair of species accumulation curves where unbiased estimation and a non-random set of sample sizes per period would produce this pattern. Applying the proposed regression correction here would distort the pattern of real differences. On the other hand, we also need to check whether it removes estimation bias when differences are calculated between different samples from the same assemblage.

Simulating data from a single species accumulation curve based on the EIS bee data, I found that the “regression correction” approach CA2013 used in their analysis of species richness change patterns at spatial resolutions smaller than the national scale does not remove bias on the intercept. When the ratio of sample sizes between periods is on average different from one in these simulations, tests on the intercept or on the partial residuals used in this regression approach too often conclude that species accumulation curves have changed. These tests have anticonservative amounts of type I errors and also show estimation bias. When the sample in the second period is on average larger (smaller), the intercept is negative (resp. positive).

Additionally, the variances of partial residuals used for “corrected” per-grid tests in CA2013 were not taking variances of estimated regression parameters into account neither the random effect variance of the model fitted (Viechtbauer 2010). Also here, simulations indicate that the variance used is often smaller than that of the partial residual when calculated correctly. Thus tests on small spatial scale differences can be expected to be anticonservative due to a biased estimator and underestimated errors, even if the regression correction were appropriate. This completes my arguments to give sub-national analyses much less importance. I will focus on the tests at the national level as much as possible, even while anticonservative inference is expected there too.

There are further irregularities in the analysis of CA2013. In the functions used to predict species accumulation curves, which I originally wrote, zero values have been replaced by positive numbers, and where functions should produce missing values or zeroes, shortcuts have been inserted that return other values. The threshold for the rma.uni() function used for weighted regression is not set at the default value of 10^-5 but at 0.01, which makes convergence of the algorithm on a decent estimate of unexplained heterogeneity in the data uncertain. In one instance, a 0.06 tail probability is used to conclude on significance of a test and not 0.05 as would be standard. Numbers of cells significantly declining or increasing in Tables S2 and S5 of CA2013: The authors have applied a two-sided test, not one-sided ones. They have counted as significant declines cells that had tail probabilities below 0.05 for that two-sided test and with a negative log-ratio value. These are thus not one-sided tests at 5% level as suggested by the verbal statement, but one-sided tests at the 2.5% significance level. The numbers of significant increasing cells are wrong altogether. Authors have counted number of grid cells with tail probabilities above 0.05 and a positive coefficient as significant, instead of the cells with a probability below 0.05. This error is made for all weighting methods, inter- and extrapolation. Thus entire columns “significant increases” are wrong.