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

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.


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][2][3][4][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][10][11][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][22][23][24] confirmed the empirical validity of TL using largescale 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][31][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][34][35][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?

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 2 ).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.5km 2 ).The smallest county (in area) is Winchester City of VA in the 1900 census (0.28 km 2 ).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 2 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 10 .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: where the TL parameters a and b are respectively the intercept and slope of TL, and � is the normal error with zero mean and constant variance.Following [21] (p.388, their Eq 14), we also fit the same log(mean)-log(variance) pairs using the ols quadratic regression to examine the curvilinearity between log(mean) and log(variance): Here e is the quadratic coefficient.If e > 0, the nonlinearity between log(mean) and log(variance) is convex.If e < 0, the nonlinearity is concave.We call the mean-variance relationship described by Eq 2 the quadratic TL (QTL).If b of Eq 1 is significantly nonzero and e of Eq 2 is not significantly different from zero, then we do not reject the hypothesis that TL is an appropriate model for the mean-variance relationship.If e of Eq 2 is significantly different from zero, then we prefer QTL over TL as a description of the mean-variance relationship.Neither TL nor QTL is suitable for describing the relationship when the linear and quadratic regressions are not significant.Here "significance" means that the p-value of the corresponding parameter or regression is less than 0.05 or, after Bonferroni correction, the p-value is less than 0.05/number of tests.We obtain the point estimate, standard error and p-value of each coefficient parameter of TL (a and b) and QTL (c, d and e), and the p-value of the overall regression of QTL.We also record the adjusted coefficient of determination (adj.R 2 ) for each ols regression.We perform regression diagnostics for the linear and quadratic models and present the results in the S1 File.
For each type of TL, we use one minus adj.R 2 (1-adj.R 2 ) 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 2 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 2 correspond to larger deviations, indicating a greater degree of dissimilarity.Specifically, for the spatial hierarchical (Q)TL and spatial (Q)TL, 1-adj.R 2 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 2 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 meanvariance 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   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  and d).Sign of the quadratic coefficient is determined by the coefficient significance after Bonferroni correction.Specifically, "zero" is defined when the p-value of the quadratic coefficient (e of Eq 2) is greater than 0.05/number of tests (55 for spatial (Q)TL and 64 for temporal (Q)TL), meaning that the coefficient is not significantly different from zero."minus" is defined when the p-value is smaller than 0.05/number of tests and the point estimate of e is less than zero."plus" is defined when the p-value is smaller than 0.05/number of tests and the point estimate of e is larger than zero.This US map is made with the package fiftystater [43] in R. fiftystater is free software that can be redistributed and/or modified under the terms of the GNU General Public License as published by the Free Software Foundation, version 3. https://doi.org/10.1371/journal.pone.0226096.g003linearly 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 2 = 0.34), followed by the 1810 and 1800 censuses (with 1-adj.R 2 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 acrossstate dissimilarity (deviation from TL) among all censuses with 1-adj.R 2 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 2 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 2 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 2 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 2 equaling 0.01, 0.02 and 0.04 respectively.

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 ( P density i � area i area ¼ area ), 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][47][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][50][51][52][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 multiyear 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.
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,S3and 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,S5and S6 Tables).

Fig 3 .
Fig 3. US maps showing the slope (a and b) and the sign of quadratic coefficient (c and d) of the spatial (Q)TL (a and c) and temporal (Q)TL (band d).Sign of the quadratic coefficient is determined by the coefficient significance after Bonferroni correction.Specifically, "zero" is defined when the p-value of the quadratic coefficient (e of Eq 2) is greater than 0.05/number of tests (55 for spatial (Q)TL and 64 for temporal (Q)TL), meaning that the coefficient is not significantly different from zero."minus" is defined when the p-value is smaller than 0.05/number of tests and the point estimate of e is less than zero."plus" is defined when the p-value is smaller than 0.05/number of tests and the point estimate of e is larger than zero.This US map is made with the package fiftystater[43] in R. fiftystater is free software that can be redistributed and/or modified under the terms of the GNU General Public License as published by the Free Software Foundation, version 3.

Fig 5 .
Fig 5. Temporal mean-temporal variance relationships on log-log scale.These illustrative examples are selected from a continuum of variation; they do not represent distinct categories.(a) States showing sublinear (superlinear) increase in log(spatial variance) as a function of log(spatial mean) are defined as decentralizing (centralizing) states, compared to the balancing states where log(spatial variance) grows linearly with log(spatial mean).(b) States are classified by the temporal fluctuation of county population count of their most populous (on average) counties: hibernating county (county with the least temporal fluctuation), hybrid county (county with intermediate temporal fluctuation), less vibrant county (county with large temporal fluctuation) and vibrant county (county with the largest temporal fluctuation).https://doi.org/10.1371/journal.pone.0226096.g005

S40Fig.
Spatial variance of county count against spatial mean of county count across censuses in each state that occurs in at least 5 censuses.(PDF) S41 Fig. Spatial variance of county density against spatial mean of county density across censuses in each state that occurs in at least 5 censuses.(PDF) S42 Fig. Point estimates of the slope ((a) and (b)) of ols linear regressions and the sign of the quadratic coefficient ((c) and (d)) of ols quadratic regressions for spatial TL for each state that occurs in at least 5 censuses, using count ((a) and (c)) and density ((b) and (d)).(PDF) S43 Fig. Temporal variance of county count against temporal mean of county count across counties in each of 36 states.(PDF) S44 Fig. Temporal variance of county density against temporal mean of county density across counties in each of 36 states.(PDF) S45 Fig. Point estimates of the slope ((a) and (b)) of ols linear regressions and the sign of the quadratic coefficient ((c) and (d)) of ols quadratic regressions for spatial TL for each state, using count ((a) and (c)) and density ((b) and (d)).(PDF) S46 Fig. Scatterplots of the point estimate of b of the spatial hierarchical TL (shtl) and e of the spatial hierarchical quadratic TL (shqtl) using count, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement (n � 15).(PDF) S47 Fig. Scatterplots of the point estimate of b of the spatial hierarchical TL (shtl) and e of the spatial hierarchical quadratic TL (shqtl) using density, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement (n � 15).(PDF) S48 Fig. Scatterplots of the point estimate of b of the spatial TL (stl) and e of the spatial quadratic TL (sqtl) using count, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement (n � 15).(PDF) S49 Fig. Scatterplots of the point estimate of b of the spatial TL (stl) and e of the spatial quadratic TL (sqtl) using density, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement (n � 15).(PDF) S50 Fig. Scatterplots of the point estimate of b of the temporal TL (ttl) and e of the temporal quadratic TL (tqtl) using count, with (vertical axis) and without (horizontal axis) the minimum number of sampling units requirement (n � 15).(PDF) S51 Fig. Scatterplots of the point estimate of b of the temporal TL (ttl) and e of the temporal quadratic TL (tqtl) using density, with (vertical axis) and without (horizontal axis) the

.05/number of tests n N adj. R 2 of Eq 1 adj. R 2 of Eq 2 b in Eq 1 e in Eq 2 b in Eq 1 e in Eq 2
b denotes the slope estimate of TL and e denotes the quadratic coefficient estimate of QTL.n represents the number of spatial or temporal units used to calculate a mean-variance pair.N represents the number of mean-variance pairs in a single testing of TL or QTL.The proportion of significantly nonzero regression coefficients was calculated without (significance level equals 0.05) and with (significance level equals 0.05/number of tests) Bonferroni correction separately.https://doi.org/10.1371/journal.pone.0226096.t001