Productivity, biodiversity, and pathogens influence the global hunter-gatherer population density

Data.

All data manipulation and analyses were conducted using R statistical software (41). The R script and the datasets are freely available at the Zenodo repository (https://doi.org/10.5281/zenodo.1069786), which makes it possible to reproduce this study. SI Materials and Methods gives more details regarding datasets and methods.

We extracted the annual mean temperature and annual precipitation values from the WorldClim.org database (42), and based on these values, we calculated climatic NPP using the empirical Miami model (43) (Fig. S1A). The Miami model, although relatively simple, has been shown to give reasonable estimates of NPP in current climatic and vegetation conditions (44) (SI Materials and Methods). Instead of using satellite-based values, we decided to use climate-based NPP estimates to avoid the potential problems caused by recent human appropriations of NPP and to reach an estimate that reflects more natural levels of NPP.

We used combined mammal, bird, and vascular plant richness as a surrogate of overall biodiversity. Global mammal and bird richness data were extracted from the BiodiversityMapping.org database (45). Global native vascular plant richness data are model based, because native plant species’ ranges are insufficiently documented in most regions of the world (46, 47). We interpolated plant richness data to cover the same global grid with the same spatial resolution as the NPP and mammal and bird richness data using inverse distance weighting interpolation with the “gstat” package in R (48⇓–50) (SI Materials and Methods and Fig. S5 A and B). The biodiversity variable was calculated as an average of the scaled animal (mammal + bird) and plant richness (Fig. S1B).

The estimations of environmental pathogen stress are based on observations of the local pathogen conditions prevailing among 186 nonindustrial societies around the globe (51⇓–53). These societies form the Standard Cross-Cultural Sample (www.worldcultures.org/) (54). We used the pathogen prevalence index, which is an average of the z scores of the prevalence of 10 pathogens (malaria, dengue, filariae, typhus, trypanosomes, leishmaniasis, schistosomes, plague, leprosy, and spirochetes). We interpolated the pathogen data to cover the same global grid and the same spatial resolution as the NPP and biodiversity data using kriging interpolation and an exponential variogram model with the gstat package in R (48⇓–50) (SI Materials and Methods and Figs. S1C and S5 C and D). We assume that the global interpolation of the pathogen data reflects the relative pathogen stress levels in the local environments of the study populations.

The hunter-gatherer population density data were combined from two ethnographic datasets compiled by Binford (10) and Kelly (11), with the addition of Saami in northern Europe from another source (55). To reduce the uncertainty related to ethnographic population density estimates, if a population was present in both datasets, we averaged the population density values from both sources. When a hunter-gatherer population was present in only one of the datasets, we used the value given by that dataset (SI Materials and Methods and Fig. S6). For each population in the data, we extracted the values of the environmental explanatory variables from the global environmental layers of NPP, biodiversity, and environmental pathogen stress.

To analyze biomewise differences in population density, we used the Food and Agriculture Organization of the United Nations (FAO) Global Ecological Zones data on the distribution of 20 natural vegetation zones (28, 29) available at www.fao.org/geonetwork/srv/en/main.home?uuid=2fb209d0-fd34-4e5e-a3d8-a13c241eb61b.

Methods.

Visual inspections of the relationships between NPP and the log_e population density of hunter-gatherers and between NPP and environmental pathogen stress indicated shifts in the effects of NPP. This impression was confirmed using breakpoint analysis with the segmented package in R, which allows for the fitting of piecewise linear regression models and provides estimates of breakpoints (56, 57). This package uses an iterative model-fitting algorithm. We did not provide any starting values for the breakpoint parameters, which means that the algorithm used the median of the predictor variable as a starting value for breakpoint estimation. The Davies test (58) included in the segmented package was performed to test whether the difference between the slopes before and after the estimated breakpoint was different from zero. We calculated the mean of the two NPP breakpoints (for population density and pathogen stress) and used this value to distinguish between low- and high-productivity environments.

Structural equation modeling is an approach that interprets information about the observed correlations between system components to evaluate potentially complex causal relationships (23). The approach begins with hypothesizing the underlying structure of causal pathways in the system (a conceptual model or metamodel). Then, this structure is translated into regression equations, and finally, these equations are evaluated against data to support or refute the hypothesized structure. Through this process, structural equation modeling provides an understanding of the nature and magnitude of direct and indirect effects in a system (23), which in our case, is composed of three environmental factors and human population density.

Here, we performed structural equation modeling as a piecewise estimation of local relationships, in which the parameters are estimated separately for each regression equation (59). The fit of the structure of the piecewise model is evaluated using a test of directed separation, with Fisher’s C as the test statistic (59, 60). This test indicates whether missing paths between variables in a system should be included in the model. We performed modeling using the R package piecewiseSEM, which allows for the automated fitting of lists of structured equations and goodness-of-fit tests for full and component models (59). For the sake of interpretability, we fit separate models for low- and high-productivity environments.

We began with the models (for low- and high-productivity environments) that omitted the direct effects of NPP on log_e population density. For low- and high-productivity structural equation models, the goodness-of-fit test strongly indicated that the direct effect of NPP should be included. In such saturated models, however, the effect of pathogen stress on population density (low productivity) and the effect of biodiversity on both pathogen stress and population density (high productivity) were nonsignificant. After omitting these effects, the goodness-of-fit test indicated that the structural equation models fit the data well in both low- and high-productivity environments.

However, after adjusting the P values for spatial autocorrelation (61) (SI Materials and Methods), the effect of NPP on pathogen stress turned out to not be statistically significant in the low-productivity environments. Therefore, we omitted this effect from the final structural equation model, although Fisher’s C indicated a lower fit than for the model that included the effect of NPP on pathogen stress.

For global projections of hunter-gatherer population density, we used two regression methods: ordinary least squares linear regression model (LM) (Fig. S3 and Table S2) and support vector regression (SVR) (62⇓–64). Both of these models included interaction terms between biodiversity and NPP and between pathogen stress and NPP. LM and SVR produce highly similar predictions (Fig. S7). However, because fitting algorithms are different (62⇓–64), their predictions are not identical, and therefore, these techniques are complementary (Fig. S7). Neither of the techniques seem to extrapolate, and the predicted values are within the boundaries of the original hunter-gatherer data (Fig. S7). In the context of global prediction and extrapolation, we nevertheless feel that it safer to use a combination of two techniques instead of single technique.

We fit global models using NPP, biodiversity, and pathogen stress as predictors of log_e population density. SVR was implemented with the e1071 package in R (65). The LM and SVR (with linear kernel) models were trained without splitting the ethnographic data into low- and high-productivity subsets. After the training, we used model algorithms to predict the log_e population density using global environmental data. We calculated the mean of the predictions of the LM and SVR models and used this average as our global projection of hunter-gatherer population densities as shown in Fig. 3B. The performance of the averaged (LM and SVR) model was evaluated using h-block cross-validation (SI Materials and Methods) (66, 67); h-block cross-validated (h = 500 km) prediction root mean squared error is 1.26, and R² is 0.46. These metrics indicate that the model is performing well with the unseen data.

The limiting effect of a predictor variable is calculated as the difference between the population density of the prediction, in which the given predictor is set to its optimum value, and the “true” prediction (see above), in which all of the predictor variables are allowed to vary. The larger the difference, the stronger the limiting effect of the predictor. The map of the limiting effects of NPP, biodiversity, and pathogen stress was constructed by determining the optimum values of each predictor in the averaged (LM and SVR) model. The optimum value of a predictor is the value where the log_e population density reaches its highest value in the ethnographic data, as seen in the response shapes of the averaged (LM and SVR) model. We then created three global datasets, in which one predictor variable at a time was set to its optimum and the other two predictors were allowed to vary. Next, we produced three global projections of population density using these datasets. From these projections, we subtracted the true prediction. The resulting differences were scaled to vary between 0 and 255. The final map was produced using the plotRGB function in the “raster” package in R (68). As humans have had profound impact on biodiversity, the limiting effect of biodiversity may be artificially high in certain regions. This is because biodiversity in these regions is further away from its optimum value than it would be without the recent human influence.

For biomewise comparisons, we rasterized the biome dataset (see above) and aligned it with the global projected limiting effects and the population density data. This allowed us to extract information about the distributions of the limiting effects and the log_e population density in different biomes. For the biomewise comparisons, the limiting effects were scaled to sum to unity.