Analyzing and interpreting spatial and temporal variability of the United States county population distributions using Taylor's law

Summary

We study the spatial and temporal variation of the human population in the United States (US) counties from 1790 to 2010, using an ecological scaling pattern called Taylor's law (TL). TL states that the variance of population abundance is a power function of the mean population abundance. Despite extensive studies of TL for non-human populations, testing and interpreting TL using data on human populations are rare. Here we examine three types of TL that quantify the spatial and temporal variation of US county population abundance. Our results show that TL and its quadratic extension describe the mean-variance relationship of county population distribution well. The slope and statistics of TL reveal economic and demographic trends of the county populations. We propose TL as a useful statistical tool for analyzing human population variability. We suggest new ways of using TL to select and make population projections.

Keywords:

Population growth – Census – United States – Population density – Urban areas – Test statistics – Norwegian people

Introduction

Spatial and temporal distributions of human populations hold important information about human societies and environments. City planning, business development, immigration policy, and many other activities rely on knowledge about the spatial and temporal distribution of people and their mobility [1–5]. Human population distributions in different areas are related. For example, metropolitan areas with high population density require large and sustained supplies of natural resources and products, which can affect the distribution of humans in neighboring areas. Human population distribution interacts with the environment and its resources, such as desertification [6] and biodiversity loss [7,8], and is affected by external conditions such as climate change and other environmental shifts [9–12]. Studying the spatial and temporal distribution of humans is a step towards describing and understanding such complex social and ecological interactions.

Human populations are not uniformly distributed in space or time, but depend on resources (e.g., food, jobs, and weather) and social development (e.g., urban and rural). Human population densities and numbers vary spatially and temporally. Specifically, the spatial unevenness or "centralization" of the population distribution is commonly characterized using two measures from economic theory: the Gini coefficient [13] and the Hoover index [14]. Both measures quantify the dissimilarity between the population count and inhabited area, and are intrinsically affected by the spatial scales used in the calculation. On the other hand, temporal variation of the population distribution is often measured by the coefficient of variation, which may vary with the mean population abundance. While these statistics have their own caveats, they are useful tools of describing the variation patterns to understand human behaviors (e.g., migration) and their impact on the environment.

Here we propose another measure to describe human population variation in space and time: Taylor's law (TL). TL is an ecological scaling pattern stating that the variance of species abundance (measured by number of individuals, i.e., "count," or by count per unit area, i.e., "density") is a power function of the mean species abundance. TL holds when the logarithm of the variance is well approximated as a linear function of the logarithm of the mean. The slope of the log(mean)-log(variance) relationship equals the exponent of the power function from the mean to the variance. For a given unit change in the log(mean) of population abundance, TL's slope equals the change in log(variance), which measures heterogeneity or scatter in the distribution of population abundance.

TL has been confirmed for numerous biological taxa and non-biological variables [15,16]. Bliss [17], Fracker and Brischle [18], Hayman and Lowe [19] and Taylor [20] proposed and showed empirically that a power law described well the variance of population count or population density as a function of the mean population count or population density in multiple samples. Taylor and colleagues [21–24] confirmed the empirical validity of TL using large-scale sample data of aphid, moth, and bird species across the United Kingdom. Most but not all empirical estimates of TL's slope fell between one and two (see Fig 2 in [25]). TL has also been verified for non-biological variables such as crime incidents [26], stock illiquidity [27], prime numbers [28], and tornado outbreaks [29]; in this last case, the slope exceeded four. Several models have been proposed to study the biological or statistical mechanisms of TL [30–32]. A unified explanation of TL is still lacking and probably does not exist, given the diversity of empirical situations where TL applies. Tests of TL against human population data are few [16,21,33–36]. Only Taylor et al. [21] (p. 392) and Taylor [16] (p. 392–394) analyzed the United States (US) population data. Other works tested TL using human population data from Norway, Italy and the whole world.

The slope of TL has three advantages for measuring the variability of sets of human population distributions. First, TL's empirical and theoretical roots in ecology can suggest possible mechanisms of human population distribution. Second, TL's slope measures the relative change in log(variance) with respect to a unit increase in log(mean) and, when TL is valid, the slope is independent of the average population abundance. This feature overcomes the density-dependence of other measures (e.g., coefficient of variation and the relative dispersion index from the mean-variance relationship of negative binomial distribution [37]) and makes TL ideal for comparisons of multiple populations across time or space. Third, TL can be tested under various spatial and temporal structures of population distribution. Such flexibility makes TL's slope a consistent measure for spatial and temporal variation.

Our goal here is to give a detailed empirical example, based on the history of the human population of the United States, to highlight for demographers and human ecologists the usefulness of TL as a statistical tool for analyzing human population variation. While previous studies have confirmed the usefulness of TL for human and nonhuman populations, we provide the first comprehensive analysis of different types of TL using US county populations, and we interpret the results from a demographic perspective.

We use TL to characterize the spatial and temporal variation of US county population count and population density reported in the decennial censuses from 1790 to 2010. We use counties as the basic spatial units. As the subordinate administrative divisions of the states, counties allow us to study the within-state population distribution and change. Although analyses of TL can be done similarly for populations from larger or smaller spatial scales (e.g., states and census tracts), they are not included in this paper due to the length limit. We test empirically three types of TL: spatial hierarchical TL, spatial TL and temporal TL. Then we focus on the interpretation and implication of statistics derived from these tests. We address the following specific questions: 1. What type of variation is quantified by each TL? 2. Is TL valid for all censuses and all states? If yes, what demographic information can be extracted from TL and what is the usefulness of such information? If not, what causes the exception(s)? 3. What are the similarities and differences in TL patterns between the US and other countries?

Materials and methods

Data

County data downloaded from NHGIS (National Historical Geographic Information System) Data Finder on 19 July 2018 [38] provide county population count (number of individuals residing in a county) in the decennial censuses from 1790 to 2010. County area (in square meters) and coordinates (latitude and longitude of the geographic centroid of each county) come from the 2000 TIGER/Line shapefile from 1790 to 2010. We join 23 shapefiles (one for each census), which specify each county's area, with the county population count data using the identifier "GISJOIN." GISJOIN is the GIS Join Match Code unique for each county that ever existed between 1790 and 2010. A "GISJOIN" code starts with letter "G" and is followed by two digital parts, the three-digit NHGIS Integrated State Code (STATENH) and the four-digit NHGIS Integrated County Code (COUNTYNH). For example, Autauga County of Alabama has GISJOIN code "G0100010," where "010" and "0010" are its state and county codes respectively. We denote states by their initials and territory by letter "T" (e.g., ND for North Dakota, NMT for New Mexico Territory and UTT for Utah Territory). Using these joined data, we define and calculate the population density (individuals km^-2) of a county in a census as the ratio of county population count divided by the county area (in square kilometers, or km²). County population count and density are the two abundance measures used in this work.

The data include 50 states and 31 state-equivalents (29 historical territories, Puerto Rico (PR), and District of Columbia (DC)) and 4,128 counties and county-equivalents (e.g., parishes, boroughs, independent cities, and municipalities), historical or existing. We do not distinguish between a state and its equivalents, or between a county and its equivalents, and refer to them uniformly as state or county, respectively. On average, the number of counties per state is about 51 (≈ 4,128/81) during the entire census period (1790–2010). The largest county (in area) that ever existed is Northern District of AKT in the 1880 census (899,538.5 km²). The smallest county (in area) is Winchester City of VA in the 1900 census (0.28 km²). The most populous county is Los Angeles County of CA in the 2010 census (9,818,605 individuals). 72 combinations of counties and censuses report zero individuals. The densest county is New York County of NY in the 1920 census (38,279.08 individuals km^-2). Only 62 combinations of counties and censuses have zero individuals per km² because 10 of the 72 combinations with zero individuals lack area measurements. We do not correct the early censuses for their omission or misrepresentation of native American and slave populations.

Three types of Taylor's law

Testing Taylor's law (TL) requires a hierarchical structure of the spatial and temporal population distribution. First, the mean and the variance of population abundance must be calculated at the lower spatial or temporal levels. Second, the calculated means and variances are grouped by the higher spatial or temporal levels to detect their relationship. US census data contain two nested spatial levels (state and county) and one temporal level (census), allowing us to test three types of TL.

If we calculate the spatial mean and spatial variance of county population abundance over counties within each state and each census, and find a linear relationship between the logarithmic mean and logarithmic variance across states within each census, it yields a spatial hierarchical TL. If we use the same set of means and variances and observe a linear relationship between their logarithms across censuses within each state, it corresponds to a spatial TL. Last, if we calculate the temporal mean and temporal variance of county population abundance over censuses within each county and each state, and find a linear relationship between their logarithms across counties within each state, it corresponds to a temporal TL.

Both the spatial hierarchical TL and the spatial TL examine the spatial distribution of population abundance of counties within each state, across states and across censuses respectively. Specifically, the spatial hierarchical TL analyzes the spatial heterogeneity (variation across states) of the spatial mean and the spatial variance of county population abundance. A greater slope of the spatial hierarchical TL means a greater degree of change between states in the spatial variance of county population abundance with respect to its spatial mean. This comparison is useful for the study of urbanization or "centralization" pattern in human geography [39].

The spatial TL analyzes the variation across censuses of the same means and variances. A greater slope of the spatial TL represents a greater degree of change between censuses in the spatial variance of county population abundance for a given unit of change in its spatial mean. This comparison can be used to study the temporal pattern of county population distribution for each state.

The temporal TL analyzes the variation across counties of the temporal mean and temporal variance of county population abundance. A greater slope in the temporal TL means a greater degree of change between counties in the temporal variance of county population abundance with respect to its temporal mean. This comparison can help us identify counties that experience high or low fluctuation in population abundance over time within a state.

Repeated testing of the same type of TL can extract further information about variation in the county population distribution. Specifically, repeated testing of the spatial hierarchical TL for each census allows us to detect a temporal trend of urbanization or counterurbanization. Repeated testing of spatial TL and temporal TL for each state provides a new way of classifying states according to their characteristic county population temporal dynamics.

Statistical analysis of Taylor's law

We define the number of spatial or temporal units used in calculating a single mean-variance pair as the number of sampling units, denoted by n. We denote the number of pairs of finite positive mean and finite positive variance used in a single test of TL as N. In the main text we test each TL with N ≥ 5, with no constraint on the number n of sampling units used to calculate a mean-variance pair. In S1 File we repeat the analysis of all three types of TL when the number of sampling units for each mean-variance pair is at least 15 (n ≥ 15) and the number of pairs of finite positive mean and finite positive variance is at least five (N ≥ 5) (following the recommendation in [40]).

We test whether log(variance) is an approximately linear function of log(mean) using ordinary least-squares (ols) linear and quadratic regressions. Throughout log = log₁₀. For each test of TL (for a unique combination of TL type and abundance measure), we first fit log(variance) as a linear function of log(mean) using the ols linear regression:

For each type of TL, we use one minus adj. R² (1-adj. R²) of the best-fitting model (TL or QTL) to measure the dissimilarity of county distribution among spatial and temporal units relative to the standard of TL or QTL. We may think of 1-adj. R² as one measure of the degree to which a log(mean)-log(variance) pair deviates from the linear or quadratic function (Eqs 1 and 2). Higher values of 1-adj. R² correspond to larger deviations, indicating a greater degree of dissimilarity. Specifically, for the spatial hierarchical (Q)TL and spatial (Q)TL, 1-adj. R² reflects the dissimilarity (from the (Q)TL standard) in the spatial distribution of county population abundance among states within each census and among censuses within each state, respectively. For the temporal (Q)TL, 1-adj. R² measures the dissimilarity (from the (Q)TL standard) in the temporal distribution of county population abundance among counties within each state.

All figures and statistical analyses, including summary statistics and ordinary least-squares regressions, are done in R [41]. The data that support the findings of this study are openly available in IPUMS National Historical Geographic Information System: Version 13.0 [Database] at http://doi.org/10.18128/D050.V13.0.

Results

Results when county population count is used as the abundance measure are in the main text. Results regarding county population density are in the S1 File.

Empirical confirmation of TL

Visually, the spatial hierarchical TL describes well the relationship between mean and variance of county population count (S1 Fig). Spatial TL, temporal TL and their quadratic forms are adequate models for the corresponding mean-variance relationship in a majority of states (S2 and S3 Figs). In only a few states neither the linear nor the quadratic model depicts the mean-variance relationship. Regression analysis confirms these observations (Table 1). The quadratic model (Eq 2) is favored over the linear model (Eq 1) in more than half of the states for the spatial TL, and in one quarter of the states for the temporal TL. The linear model is preferred in more states after Bonferroni correction of regression significance. Variability in the slope of spatial hierarchical TL is smaller than that of the spatial and temporal TL (Fig 1). Quadratic coefficients of all three types of QTL are close to zero on average, with average values 0.045 for the spatial hierarchical QTL, -0.12 for the spatial QTL, and -0.14 for the temporal QTL, but there are extreme outliers for each form of QTL.

Specifically, regardless of Bonferroni correction, spatial hierarchical TL (Eq 1) is not rejected in favor of its quadratic form (Eq 2) in all but one of the 23 censuses (S1 and S2 Tables). In the 1790 census, neither the linear nor the quadratic model is superior to the constant model (S1 and S2 Tables), indicating no consistent relationship between the means and the variances of county population across states during the early settlement of US.

Spatial TL and QTL are tested in 55 states (S3 and S4 Tables). Without Bonferroni correction, in two of the 55 states (AKT and NMT), neither the linear nor the quadratic model describes the mean-variance relationship well (regression p-value>0.05). In the other 53 states, spatial TL is not rejected in 23 states and spatial QTL is not rejected in 30 states. With Bonferroni correction, neither the spatial TL nor the spatial QTL describes the mean-variance relationship well in five states (AKT, AZT, NMT, ND and UTT). In the other 50 states, spatial TL is not rejected in 32 states and spatial QTL is not rejected in 18 states.

Temporal TL and QTL are tested in 64 states (S5 and S6 Tables). Without Bonferroni correction, neither the temporal TL nor the temporal QTL describes adequately the mean-variance relationship in four states (AKT, ART, FLT and WYT). In the remaining 60 states, temporal TL is not rejected in 44 states and temporal QTL is not rejected in 16 states. With Bonferroni correction, neither the temporal TL nor the temporal QTL describes adequately the mean-variance relationship in 11 states (AKT, AZT, ART, DC, FLT, IDT, MIT, MTT, NMT, RI and WYT). In the remaining 53 states, temporal TL is not rejected in 46 states and temporal QTL is not rejected in seven states.

For the spatial hierarchical TL, the mean-variance pairs are scattered in the 1790 census (S1 Fig). In particular, in DE (lowest point) and VA (highest point) we see the largest deviations from the fitted regression lines. The scatter disappears in the subsequent censuses. The violation of spatial TL or QTL of county population count in AKT and NMT is due to their drastic territorial changes preceding their transition to statehood (S2 Fig, S3 and S4 Tables). Among all the 50 current states, ND shows a unique mean-variance relationship: in most censuses the average county population count changes slightly while the variance continues to rise. The reason is that, during most of the 20th century (after the large increase in population count during the earlier censuses due to the farmers migrating into ND), ND people moved from agricultural lands to large cities, especially young adults emigrating from rural areas [42]. This internal migration caused an uneven population distribution among ND counties. The temporal TL and QTL are rejected in states with either large territorial changes (AKT) or short census periods (two censuses available in ART, FLT and WYT) (S3 Fig, S5 and S6 Tables).

Population growth trend and pattern illustrated by TL and QTL

TL and QTL of all three types reveal useful information about the population dynamics of US counties. First, the slope of the spatial hierarchical TL fluctuates before the 20th century and reaches a record high in 1910 (Fig 2), reflecting the early boundary changes and the population centralization caused by urbanization in the 19th century. Since the 1910 census, the slope has been declining. This decline occurred because suburbanization (or counterurbanization) or the emergence of new metropolitan areas across US counties during the 20th century led to increasing spatial evenness in county population count as a function of increasing mean county population count. Compared to other types of TL, the spatial hierarchical TL shows smallest variability in its slope, ranging from 1.5 to 2.5. If the trend in slope over the past five decades continues, the slope will likely fall in the 2020 census.

The spatial TL exhibits great differences among states. With Bonferroni correction, TL is the superior model in 31 current states. The slope of TL is not different from two in AK, AR, CO, HI, IA, MA, MT, OK, OR and WY, and is significantly greater than two in AL, CA, FL, GA, ID, IL, IN, KY, ME, MI, MN, MO, NV, NH, NC, OH, TN, TX, VT, VI and WI (Fig 3A). The higher slopes reflect the greater degree of centralization in county population distribution.

With Bonferroni correction, the spatial QTL depicts the mean-variance relationship better than TL in 18 current states. With the exception of SD, the quadratic coefficient (e in Eq 2) is significantly smaller than zero in AZ, CT, DE, LA, MD, NJ, NM, NY, PA, RI, UT and WA, and greater than zero in KS, MS, NE, SC and WV (Fig 3C). The former group contains states with declining centralization over time (i.e., log(variance) grows slower than linearly with log(mean)), probably due to counterurbanization during the past century. The latter group includes states with increasing centralization over time (i.e., log(variance) grows faster than linearly with log(mean)), suggesting rapid metropolitan development in these states. The mean-variance relationship in some states may be better fitted by a cubic function than by a linear or quadratic function (e.g., AR, IN and MS in S7 Fig).

Temporal TL and QTL reveal differences in the temporal variability of population count among counties within a state. Since most states have counties with small population change over the censuses, counties with great population fluctuation will dictate the slope of temporal TL. With Bonferroni correction, the temporal TL is the better model compared to the temporal QTL in 41 current states. The slope is greater than two in AL, AR, CA, CO, GA, ID, IL, IN, KY, LA, MI, MO, NY, NC, OR, PA, SC, TN, UT, VI, WV and WI, and not different from two in AK, AZ, CT, HI, ME, MD, MA, MN, MS, NV, NH, NJ, NM, ND, OK, SD, VT, WA and WY (Fig 3B). A higher slope implies that some counties experience greater population change. With Bonferroni correction, the temporal QTL is the better model compared to the temporal TL in seven current states. The quadratic coefficient is smaller than zero in FL and OH, and greater than zero in IA, KS, MT, NE and TX (Fig 3D).

Spatial and temporal dissimilarity of county population distribution

According to the spatial hierarchical TL, except the 1790 census, the dissimilarity (deviation from TL) of county spatial population distribution among states is highest in the 1900 census (1-adj. R² = 0.34), followed by the 1810 and 1800 censuses (with 1-adj. R² equaling 0.29 and 0.28 respectively) (S1 and S2 Tables). This observation may be attributed to large wave of international immigration to the US at the beginning of the 20th century and the geographic change during the earlier census periods. The 1870 census corresponds to the smallest across-state dissimilarity (deviation from TL) among all censuses with 1-adj. R² of 0.13. The dissimilarity in the most recent 2010 census ranks 12th among all 23 censuses.

The spatial TL (or QTL) points to IA, OK and MT as the three states with the largest dissimilarity (deviation from TL or QTL) of county spatial population distribution among censuses (with 1-adj. R² equaling 0.26, 0.18 and 0.15 respectively) (S3 and S4 Tables), which may be explained by various economic and political events in their state history. First, IA experienced two distinct periods of population change. From 1850 to 1900, farming and coal industries brought many out-of-state immigrants to IA who contributed to rapid population growth [44]. From 1910 to 2010, IA experienced gradual change in its economic sectors, mainly from agriculture to manufacturing operations, which led to small population growth and population redistribution from rural to urban areas due to within-state migration. During 1930s-1950s, OK's population declined (decreasing mean) and lost farming populations to urban areas (increasing variance) due to the Dust Bowl [45]. In MT, the Enlarged Homestead Act in 1909 (a law passed by the US to promote dryland cultivation in MT) seems to have contributed to the spatial evenness of county population distribution, as reflected by the drop-off in variance in the 1910 census. In the other 50 states where spatial TL or QTL describes the data well (without Bonferroni correction), 45 states have a 1-adj. R² lower than 0.1.

Using the temporal TL (or QTL), the three states with the largest dissimilarities (deviation from TL or QTL) are MS, SD and OK, with 1-adj. R² equaling 0.35, 0.31 and 0.29 respectively (S5 and S6 Tables). In MS, the relatively small change in population of several counties (e.g., Jefferson Davis County, Hinds County and Benton County) increases the overall dissimilarity of temporal distribution of county population among counties. In SD, the presence of historically short-lived counties (e.g., Sterling County and Schnasse County) and two counties (Minnehaha County and Pennington County) with fast-growing cities (Sioux Falls and Rapid City) contribute to the across-county dissimilarity. In OK, the dissimilarity is driven by Oklahoma County and Tulsa County, sites of the state's two most populous cities (Oklahoma City and Tulsa). HI, CA and OR have the least dissimilarity among current states, with 1-adj. R² equaling 0.01, 0.02 and 0.04 respectively.

Discussion

Comparison of TL patterns in the US and other countries

To the authors' knowledge, all prior works on TL for human population abundance examined the spatial hierarchical or spatial (Q)TL. In addition, most studies focused on testing the spatial TL at the national level, instead of at the subnational level as done in the current work. (The exceptions in [16,33] are described below.) Nevertheless, the comparison of spatial TL across countries provides useful insights about changes in population distributions.

Taylor et al. [21] were the first to test a spatial TL against human population census data. Despite omitting details about the data and fitting method, they observed that the spatial TL described the relationship of mean and variance of male population abundance across multiple censuses, with a slope significantly greater than two (2.04±0.01). This result probably describes variation among states, not counties, and is consistent with our finding: at the county level within states, in most states the spatial TL slope was greater than two (Fig 1). Independently, using the US county-level data, Taylor [16] examined the spatial hierarchical TL (which he called "ensemble" TL) and observed the same temporal trend of its slope as in our Fig 2A (see Fig 11.6 in [16]).

Two independent works studied the spatial TL using human population density outside the US. Cohen et al. [34] tested a spatial TL for the mean and variance of Norwegian county (equivalent to the US state) population densities across 33 years. They also repeated the same analysis using Norwegian region (equivalent to US region or division) and municipality (equivalent to the US county) population densities. They found that the spatial TL pattern depends on the temporal scale and reflects the population change caused by economic development. This finding is consistent with our findings that the spatial TL can reveal different development stages of county populations within a state (see Demographic insights from observed TL patterns below).

Using the methods of Cohen et al. [34], Benassi and Naccarato [33] tested a spatial TL for the mean and variance of Italian regional (equivalent to US divisional) population densities across 41 years. Both studies used three different weighting methods (equally weighted, areally weighted and population weighted) to calculate the set of spatial means and spatial variances. Since the areally weighted means and variances of population density are equivalent to the equally weighted means and variances of population count (∑densityi⋅areaiarea=∑countiareai⋅areaiarea=∑countiarea), we compare our results to the spatial TL of population density under areal weighting [33,34]. The spatial TL slope is greater than two for Norwegian county (or region or municipality) population densities and less than two for Italian region population densities. These cross-country and cross-state comparisons suggest that the spatial TL slopes reflect the different patterns of change in the spatial population distribution.

Benassi and Naccarato [33] also studied the spatial TL subnationally for each of the three Italian macro-geographical areas (North, Centre, South and Island, equivalent to US regions). A spatial mean and a spatial variance of population density were calculated across the provinces (equivalent to US states) within each area. Using the seemingly unrelated regression method, they found that the spatial TL slope was smaller in the North area than in the South area. The authors attributed this observation to the steady influx of international immigration to the North, the less-developed economic condition of the South, and the internal migration from South to North. Comparing the state-level slopes of spatial and temporal TL among US regions shows that the South has the highest median slope (Fig 4). These observations reflect the similarity of regional characteristics in economic development and migration between the US and Italy [46–48]. They also indicate the usefulness of spatial TL in reflecting economic and demographic patterns.

Our analysis shows that the spatial TL is not the only relation between the mean and variance of spatial population distribution of US counties. The spatial QTL in multiple states reflects a different evolution of spatial population distribution. Other works that tested spatial TL using the Norwegian [34] and world [35] population density have also documented curvature in the mean-variance relationship.

Demographic insights from observed TL patterns

TL quantifies the variability of the spatial or temporal distribution of county population abundance in relation to the mean county population abundance. We show that the observed TL pattern and slope can identify economic, demographic, and political events. The insights derived from TL thus can help government and policy makers to plan economic and demographic development at local or regional levels. To illustrate this, we summarize and distill the main points from our results.

The temporal trend in the slope of the spatial hierarchical TL classifies the US county population distribution dynamics into two periods. During the 19th century, the US frequently acquired land and established states, causing fluctuation in the slope. In addition, the industrial revolution propelled population urbanization, leading to an overall increase in the slope starting 1820. The territorial boundaries of the contiguous US became stable after the state of Arizona was added as the 48th state in 1912 (before Alaska and Hawaii in 1959). During the 20th century, suburbanization and counterurbanization were the dominating factors that shaped population distribution [49–53]. In recent decades, the formation of megalopolises enhanced the connectivity between urban areas [54,55]. These factors together contribute to the decrease of the spatial hierarchical TL slopes. The spatial hierarchical TL does not indicate a "clean break" of US population distribution around 1970 compared to previous decades [56,57].

The spatial TL and spatial QTL reveal a continuum of patterns in the log(mean)-log(variance) relationship of county population distribution: convexity (or superlinearity), linearity, and concavity (or sublinearity) (Fig 5A). States showing a convex relationship are likely in their early developmental stage, experiencing fast growth in their metropolitan areas. States showing a concave relationship have probably passed the early developmental stage and entered a period when suburbanization or counterurbanization leads to a more evenly distributed county population. The cubic pattern observed in some of the early states corroborates this claim. States with a linear relationship are probably in the stage where both centralization and decentralization are occurring in multiple areas simultaneously within the state. These illustrative examples are selected from a continuum of variation; they do not represent distinct categories.

The temporal TL and temporal QTL reflect the spatial heterogeneity of counties' temporal population distribution, since the linearity or curvature of temporal variances as a function of temporal means is largely determined by the most populous county in a state (Fig 5B). When the most populous county exhibits higher temporal fluctuation in population abundance (vibrant county), the mean-variance relationship is convex. A county with large time-averaged population that experiences small population change over time (a hibernating county) will produce a concave relationship. The linearity in the relationship could describe counties with intermediate population change over time, in either less vibrant county (TL slope > 2) or hybrid county (TL slope = 2).

What demographic factors determine the values of the TL slope? Our analysis here does not provide a direct answer. Cohen et al. [34] used a stochastic multiplicative population growth model to predict the slope of a spatial hierarchical TL at large time. They found that the size of the slope is analytically related to the population growth pattern: slope > 2 for supercritical population (population abundance increases on average) and slope < 2 for subcritical population (population abundance decreases on average). Our current work shows that the slope of the spatial hierarchical TL is less than two in recent censuses, implying that the average county population is losing population under the assumption of multiplicative growth. This claim is contradicted by the finding that population gain is still the main trend for a majority of US counties in the past century [58]. A plausible explanation for such a discrepancy is that the observed TL slope cannot be predicted solely by multiplicative demographic events (e.g., birth and death), but may be substantially affected by additive processes such as internal migration and international immigration. The demographic interpretation of the values of the TL slope deserves further research.

Other uses of TL in demographic research

This work gives a comprehensive analysis of TL for US human population data at the county level and shows the multifaceted uses of TL in extracting important socioeconomic and demographic signals. In this subsection, we address two questions. First, is the widespread (though not universal) agreement of county population distribution with TL a result of intrinsic characteristics of human populations or a purely statistical consequence of random sampling from some underlying distribution? Second, can TL be used to construct and select among alternative demographic projections?

First, Cohen and Xu [32] predicted TL parameters using random samples of probability distributions. They showed that TL's slope is proportional to the skewness of a distribution, when sampled observations are identically and independently distributed. Given the strong temporal and spatial autocorrelation in the US census data (see S1 File), the simple random sampling model by Cohen and Xu alone cannot explain TL for human population data. It is therefore necessary to construct regression models capable of incorporating the spatial and temporal structure of US county populations to analyze the statistical mechanisms of TL in this context. Alternatively, it is possible to build mathematical models to study the impact of birth, death, and migration on observed TL patterns [25,59]. The interactions between demographic stochasticity and data structure, and their contributions to TL, remain an important topic for future research.

Second, Xu et al. [36] used the spatial TL to select the most plausible projections of Norwegian human population, based on the similarity of TL slopes between data and projections. They found that the short-term population projections selected by TL agreed most closely with the recent data. Their method is purely empirical and statistical, and does not rely on any assumptions regarding population growth or migration. Their approach can be used for the US census data studied here. In addition to the spatial TL used by Xu et al. [36], the current work points to the possibility of selecting the most plausible population projection(s) using multiple TLs and QTLs. This multiplicity of approaches can be useful in case one type of TL fails to distinguish among several projections (as happened in [36]). For example, when multi-year state-level projections are provided, one could test and compare the parameters of spatial and temporal TL (or QTL) between data and projections. If multi-year county-level projections are available, one could use the spatial, temporal, and spatial hierarchical TL together to select among alternative projections.

Another potential application of TL in demographic research is to project population distributions based on TL parameters estimated from historical data. Tweedie [60], Jørgensen [61], Kendal and Jørgensen [62,63] and many other papers have explored various probability distributions that follow TL either exactly or asymptotically. If TL and appropriate parametric models can be further confirmed for human population counts or densities, then one can estimate parameters of one or several parametric probability distributions from the historical data and construct (and test) the associated abundance distribution. If multiple censuses are available, one can predict the TL parameters based on their historical trend and construct corresponding population distribution projections. This method can be tested using populations at different spatial scales (e.g., states, counties, census tracts). This approach would offer a novel but still untested way of making population projections based on the temporal or spatial population distribution in historical data described by TL.

Supporting information

S1 File [docx]
Regression diagnostics and miscellaneous testing of Taylor's law models.

S1 Fig [pdf]
Spatial variance of county population count against spatial mean of county population count on log-log scale across states within each census.

S2 Fig [pdf]
Spatial variance of county population count against spatial mean of county population count on log-log scale across censuses within each state.

S3 Fig [pdf]
Temporal variance of county population count against temporal mean of county population count on log-log scale across counties within each state.

S4 Fig [pdf]
Regression coefficient estimates for TL and QTL using county population density.

S5 Fig [pdf]
Spatial variance of county population density against spatial mean of county population density on log-log scale across states within each census.

S6 Fig [a]

S7 Fig [pdf]
Spatial variance of county population density against spatial mean of county population density on log-log scale across censuses within each state.

S8 Fig [a]

S9 Fig [pdf]
Temporal variance of county population density against temporal mean of county population density on log-log scale across counties within each state.

S10 Fig [pdf]
Residuals against fitted values of the linear regression for log(spatial variance of county count) as a function of log(spatial mean of county count) across states in a given census.

S11 Fig [pdf]
Residuals against fitted values of the linear regression for log(spatial variance of county density) as a function of log(spatial mean of county density) across states in a given census.

S12 Fig [pdf]
Residuals against fitted values of the quadratic regression for log(spatial variance of county count) as a function of log(spatial mean of county count) across states in a given census.

S13 Fig [pdf]
Residuals against fitted values of the quadratic regression for log(spatial variance of county density) as a function of log(spatial mean of county density) across states in a given census.

S14 Fig [pdf]
Standardized residuals against their normal quantiles for the ols linear regression fitted to log(spatial variance of county count) as a function of log(spatial mean of county count) across states in a given census.

S15 Fig [pdf]
Standardized residuals against their normal quantiles for the ols linear regression fitted to log(spatial variance of county density) as a function of log(spatial mean of county density) across states in a given census.

S16 Fig [pdf]
Standardized residuals against their normal quantiles for the ols quadratic regression fitted to log(spatial variance of county count) as a function of log(spatial mean of county count) across states in a given census.

S17 Fig [pdf]
Standardized residuals against their normal quantiles for the ols quadratic regression fitted to log(spatial variance of county density) as a function of log(spatial mean of county density) across states in a given census.

S18 Fig [pdf]
Residuals against fitted value of the ols linear regression for log(spatial variance of county count) as a function of log(spatial mean of county count) across censuses for a given state.

S19 Fig [pdf]
Residuals against fitted value of the ols linear regression for log(spatial variance of county density) as a function of log(spatial mean of county density) across censuses for a given state.

S20 Fig [pdf]
Standardized residuals against their normal quantiles for the ols linear regression fitted to log(spatial variance of county count) as a function of log(spatial mean of county count) across censuses for a given state.

S21 Fig [pdf]
Standardized residuals against their normal quantiles for the ols linear regression fitted to log(spatial variance of county density) as a function of log(spatial mean of county density) across censuses for a given state.

S22 Fig [pdf]
Autocorrelation function against lag for the residuals of ols linear regression (up to lag 10) fitted to log(spatial variance of county count) as a function of log(spatial mean of county count) across censuses in a given state.

S23 Fig [pdf]
Autocorrelation function against lag for the residuals of ols linear regression (up to lag 10) fitted to log(spatial variance of county density) as a function of log(spatial mean of county density) across censuses in a given state.

S24 Fig [pdf]
Residuals against fitted value of the ols quadratic regression for log(spatial variance of county count) as a function of log(spatial mean of county count) across censuses for a given state.

S25 Fig [pdf]
Residuals against fitted value of the ols quadratic regression for log(spatial variance of county density) as a function of log(spatial mean of county density) across censuses for a given state.

S26 Fig [pdf]
Standardized residuals against their normal quantiles for the ols quadratic regression fitted to log(spatial variance of county count) as a function of log(spatial mean of county count) across censuses for a given state.

S27 Fig [pdf]
Standardized residuals against their normal quantiles for the ols quadratic regression fitted to log(spatial variance of county density) as a function of log(spatial mean of county density) across censuses for a given state.

S28 Fig [pdf]
Autocorrelation function against lag for the residuals of ols quadratic regression (up to lag 10) fitted to log(spatial variance of county count) as a function of log(spatial mean of county count) across censuses in a given state.

S29 Fig [pdf]
Autocorrelation function against lag for the residuals of ols quadratic regression (up to lag 10) fitted to log(spatial variance of county density) as a function of log(spatial mean of county density) across censuses in a given state.

S30 Fig [pdf]
Residuals against fitted value of the ols linear regression for log(temporal variance of county count) as a function of log(temporal mean of county count) across counties for a given state.

S31 Fig [pdf]
Residuals against fitted value of the ols quadratic regression for log(temporal variance of county count) as a function of log(temporal mean of county count) across counties for a given state.

S32 Fig [pdf]
Standardized residuals against their normal quantiles for the ols linear regression fitted to log(temporal variance of county count) as a function of log(temporal mean of county count) across counties for a given state.

S33 Fig [pdf]
Standardized residuals against their normal quantiles for the ols quadratic regression fitted to log(temporal variance of county count) as a function of log(temporal mean of county count) across counties for a given state.

S34 Fig [pdf]
Residuals against fitted value of the ols linear regression for log(temporal variance of county density) as a function of log(temporal mean of county density) across counties for a given state.

S35 Fig [pdf]
Residuals against fitted value of the ols quadratic regression for log(temporal variance of county density) as a function of log(temporal mean of county density) across counties for a given state.

S36 Fig [pdf]
Standardized residuals against their normal quantiles for the ols linear regression fitted to log(temporal variance of county density) as a function of log(temporal mean of county density) across counties for a given state.

S37 Fig [pdf]
Standardized residuals against their normal quantiles for the ols quadratic regression fitted to log(temporal variance of county density) as a function of log(temporal mean of county density) across counties for a given state.

S38 Fig [pdf]
Spatial variance of county count against spatial mean of county count across states in each census.

S39 Fig [pdf]
Spatial variance of county density against spatial mean of county density across states in each census.

S40 Fig [pdf]
Spatial variance of county count against spatial mean of county count across censuses in each state that occurs in at least 5 censuses.

S41 Fig [pdf]
Spatial variance of county density against spatial mean of county density across censuses in each state that occurs in at least 5 censuses.

S42 Fig [a]

S43 Fig [pdf]
Temporal variance of county count against temporal mean of county count across counties in each of 36 states.

S44 Fig [pdf]
Temporal variance of county density against temporal mean of county density across counties in each of 36 states.

S45 Fig [a]

S46 Fig [pdf]
Scatterplots of the point estimate of of the spatial hierarchical TL (shtl) and of the spatial hierarchical quadratic TL (shqtl) using count, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement ( ≥ 15).

S47 Fig [pdf]
Scatterplots of the point estimate of of the spatial hierarchical TL (shtl) and of the spatial hierarchical quadratic TL (shqtl) using density, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement ( ≥ 15).

S48 Fig [pdf]
Scatterplots of the point estimate of of the spatial TL (stl) and of the spatial quadratic TL (sqtl) using count, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement ( ≥ 15).

S49 Fig [pdf]
Scatterplots of the point estimate of of the spatial TL (stl) and of the spatial quadratic TL (sqtl) using density, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement ( ≥ 15).

S50 Fig [pdf]
Scatterplots of the point estimate of of the temporal TL (ttl) and of the temporal quadratic TL (tqtl) using count, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement ( ≥ 15).

S51 Fig [pdf]
Scatterplots of the point estimate of of the temporal TL (ttl) and of the temporal quadratic TL (tqtl) using density, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement ( ≥ 15).

S1 Table [xlsx]
Ordinary least-squares (ols) linear regression statistics of spatial hierarchical Taylor's law, using count and density separately.

S2 Table [xlsx]
Ordinary least-squares (ols) quadratic regression statistics of spatial hierarchical quadratic Taylor's law, using count and density separately.

S3 Table [xlsx]
Ordinary least-squares (ols) linear regression statistics of spatial Taylor's law, using count and density separately.

S4 Table [xlsx]
Ordinary least-squares (ols) quadratic regression statistics of spatial quadratic Taylor's law, using count and density separately.

S5 Table [xlsx]
Ordinary least-squares (ols) linear regression statistics of temporal Taylor's law, using count and density separately.

S6 Table [xlsx]
Ordinary least-squares (ols) quadratic regression statistics of temporal quadratic Taylor's law, using count and density separately.

S7 Table [xlsx]
Summary of the regression statistics of the Taylor's law (TL) () and quadratic Taylor's law (QTL) () models using county population density.

S8 Table [xlsx]
Summary statistics of county area, county count and county density in each state and census.

S9 Table [xlsx]
Diagnostic statistics of the ols linear regression for spatial hierarchical Taylor's law, using count and density separately.

S10 Table [xlsx]
Diagnostic statistics of the ols quadratic regression for spatial hierarchical quadratic Taylor's law, using count and density separately.

S11 Table [xlsx]
Diagnostic statistics of the ols linear regression for spatial Taylor's law, using count and density separately.

S12 Table [xlsx]
Diagnostic statistics of the ols quadratic regression for spatial quadratic Taylor's law, using count and density separately.

S13 Table [xlsx]
Diagnostic statistics of the ols linear regression for temporal Taylor's law, using count and density separately.

S14 Table [xlsx]
Diagnostic statistics of the ols quadratic regression for temporal quadratic Taylor's law, using count and density separately.

S15 Table [xlsx]
Ordinary least-squares (ols) linear regression statistics of the spatial hierarchical TL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S16 Table [xlsx]
Ordinary least-squares quadratic regression statistics of the spatial hierarchical quadratic Taylor's law (QTL) with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S17 Table [xlsx]
Ordinary least-squares linear regression statistics of the spatial TL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S18 Table [xlsx]
Ordinary least-squares quadratic regression statistics of the spatial QTL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S19 Table [xlsx]
Ordinary least-squares linear regression statistics of the temporal TL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S20 Table [xlsx]
Ordinary least-squares quadratic regression statistics of the temporal QTL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S21 Table [xlsx]
Diagnostic statistics of the ols linear regression for the spatial hierarchical TL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S22 Table [xlsx]
Diagnostic statistics of the ols quadratic regression for the spatial hierarchical QTL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S23 Table [xlsx]
Diagnostic statistics of the ols linear regression for the spatial TL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S24 Table [xlsx]
Diagnostic statistics of the ols quadratic regression for the spatial QTL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S25 Table [xlsx]
Diagnostic statistics of the ols linear regression for the temporal TL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.

S26 Table [xlsx]
Diagnostic statistics of the ols quadratic regression for the temporal QTL with the minimum number of sampling units requirement (n ≥ 15), using count and density separately.