Soil Sample Assay Uncertainty and the Geographic Distribution of Contaminants: Error Impacts on Syracuse Trace Metal Soil Loading Analysis Results

A research team collected 3609 useful soil samples across the city of Syracuse, NY; this data collection fieldwork occurred during the two consecutive summers (mid-May to mid-August) of 2003 and 2004. Each soil sample had fifteen heavy metals (As, Cr, Cu, Co, Fe, Hg, Mo, Mn, Ni, Pb, Rb, Se, Sr, Zn, and Zr), measured during its assaying; errors for these measurements are analyzed in this paper, with an objective of contributing to the geography of error literature. Geochemistry measurements are in milligrams of heavy metal per kilogram of soil, or ppm, together with accompanying analytical measurement errors. The purpose of this paper is to summarize and portray the geographic distribution of these selected heavy metals measurement errors across the city of Syracuse. Doing so both illustrates the value of the SAAR software’s uncertainty mapping module and uncovers heavy metal characteristics in the geographic distribution of Syracuse’s soil. In addition to uncertainty visualization portraying and indicating reliability information about heavy metal levels and their geographic patterns, SAAR also provides optimized map classifications of heavy metal levels based upon their uncertainty (utilizing the Sun-Wong separability criterion) as well as an optimality criterion that simultaneously accounts for heavy metal levels and their affiliated uncertainty. One major outcome is a summary and portrayal of the geographic distribution of As, Cr, Cu, Co, Fe, Hg, Mo, Mn, Ni, Pb, Rb, Se, Sr, Zn, and Zr measurement error across the city of Syracuse.


Introduction
The geography of error (arising from, e.g., calculation, sampling, measurement, specification, and stochastic sources) dates back many decades, with Goodchild and Gopal's [1] book, followed by the first International Symposia on Spatial Accuracy Assessment in Natural Resources and Environmental Sciences convened in Williamsburg, VA, in 1994 (http://www.spatial-accuracy.org/History, accessed on 8 February 2021), signifying more concerted geospatial research efforts to address this theme. The history of this topic parallels that of popular applied statistics in general, which evolved from reporting only central tendency values to also reporting their associated uncertainties (e.g., margins of error), which became a common practice only after 1928. In parallel, maps began including error statements, mostly after 1941, although they were global map-wide ones. However, reporting a global error measure on an increasing number of maps furnishes little knowledge and understanding about the geography of error in general (even though a meta-analysis might provide some insights). Meanwhile, as geostatistical analysis became more prevalent (the software packages GEO-EAS [2] and GSLIB [3] appeared at the beginning of the 1990s), prediction error maps began accompanying krigged surface (i.e., mean response) maps, a practice that continues today. The principal disclosure of this geography of error exception is that its map pattern relates to the underlying sampling network, which almost always is visible in a prediction error map. Hexagonal-tessellation stratified random spatial sampling 2 of 28 (e.g., [4]) error constitutes another exception, with a primary concern of its geography being the relationship between spatial landscape coverage and sampling error. Current georeferenced data releases, such as those from the United States (USA) American Community Survey (ACS), include sampling errors, furnishing an increasing number of databases to study the geography of sampling error in a more comprehensive way.
The objective of this paper is to contribute to this geography of error literature, focusing on combinations of measurement (i.e., soil geochemistry content analysis accuracy/precision), sampling (i.e., the selection of only a tiny fraction of soil in a geographic landscape), and specification (i.e., the correctness of assumptions and/or functional forms of equations for an analysis) errors in the context of spatial autocorrelation, a fundamental feature of geospatial data. Its empirical analyses and simulation experiments exploit a unique database containing analytical measurement errors (i.e., uncertainty introduced by chemical assay procedures) for trace metals contaminating soil in the city of Syracuse, NY, USA.

Background
Errors in georeferenced data can occur in both of their components: location and attribute. The geography and GIScience literature contains extensive investigations about location errors in various contexts, including geocoding (e.g., [5,6]) and raster modeling (e.g., [7]). This literature also includes studies about georeferenced data attribute errors. Griffith et al. [8] discuss four major sources of georeferenced data errors, namely sampling error, measurement error, specification error, and analytical assaying error. The literature addresses these first three sources for attribute errors. For example, Griffith et al. [9] argue that sampling error of estimates is heterogeneous across spatial units, specifically census tracts, and tends to correlate with both the size of tracts and their socio-economic characteristics. Wang et al. [10] discuss sampling as a major source of error, and present a sampling approach to reduce error variances. In addition, Leung et al. [11] discuss a general error analysis framework concerning georeferenced data measurements, whereas many literature entries (e.g., [12][13][14]) recognize model specification error. In contrast, assaying errors have not been a popular topic in spatial analysis and/or spatial statistics, although they generally occur in field surveys of, for example, soil [15,16]. Kriging is one spatial analysis topic that does not avoid referencing assaying errors; it links these errors to the nugget effect [17]. Otherwise, assaying errors are rarely investigated in spatial analysis.
Griffith et al. [8] discuss ways to enhance research about uncertainty. They focus on the following four research themes: visualizing error; spatial patterns and spatial modelling acknowledging error; spatial data aggregation; and data quality. Visualizing error, the first theme, can improve the current understanding about spatial patterns of information uncertainty. Research dealing with mapping uncertainty presents various approaches focusing on mapping observations together with their data quality information [18,19]. Recently, Koo et al. [20] discussed a framework for uncertainty mapping. Nevertheless, the common protocol has two side-by-side maps, displaying a geographic distribution of values juxtaposed with a geographic distribution of individual value uncertainties, because concurrently exhibiting the two sets of information in a single map tends to be overwhelmed by the portrayal of a large amount of data in a limited two-dimensional space.
Spatial patterns and spatial data modelling constitute a second research theme that can benefit from incorporating uncertainty into spatial analysis. Sun et al. [21] discuss a map classification technique that seeks to avoid adverse impacts by the presence of uncertainty, one of which is the possibility of a different map pattern outcome. They present a new map classification approach incorporating uncertainty, demonstrating it with American Community Survey data, which is accompanied by sampling error measures. Meanwhile, Koo et al. [22] and Mu and Tong [23] present alternative approaches to incorporating uncertainty into map classifications. New approaches to spatial modelling recognizing uncertainty, going beyond widely recognized model specification issues such as omitted variables [24], also need development. Hu et al. [14] furnish one promising effort that models the existence of a mixture of positive and negative spatial autocorrelation in spatial data, showing how successfully accounting for such a mixture pattern can enhance a spatial analysis.
A third research theme, spatial data aggregation, can also raise uncertainty issues. The modifiable areal unit problem (MAUP) describes a convolution of spatial analysis, with observations based on polygon area units [25]. Lee et al. [26] present simulation experimental results demonstrating how spatial autocorrelation levels may be affected by areal unit aggregation (both geographic resolution and zonation). Despite these spatial data aggregation issues being recognized in the literature for decades now, a universally acceptable solution to the MAUP remains elusive for geographic analysis. Finally, Griffith et al. [8] discuss the importance of spatial metadata for data quality so that users understand the accuracy level of data and, when necessary, more cautiously assess spatial analysis results.
Given this preceding context, this paper investigates assaying measurement errors of heavy metal amounts in soil sample observations. It focuses mainly on the geographic patterns that materialize from mapping individual soil samples with point locations and then further aggregates these data by census tracts. In addition, this paper explores both individual and a mixture of different error types with simulation experiments. By doing so, its principal knowledge contribution is to the geography of error literature.

About the Data
A total of 3628 useful surface soil samples were collected across the city of Syracuse, New York, USA, mostly during the summers of 2003 and 2004, under the auspices of the National Science Foundation (BCS-0552588), with Institutional Review Board (IRB) oversight by that organization for Syracuse University, then the University of Miami, and finally the University of Texas at Dallas. Of these, 3324 unique location samples-multiple samples for the same locations were averaged by location-had useful (i.e., positive) analytical assay error values ( Figure 1); two of these soil samples had locations outside of Syracuse, resulting in their being removed from census tract resolution spatial analyses. This dataset contained a rich set of observations, with extensive sample coverage across the city of Syracuse, including measurements for 15 heavy metals and corresponding assay errors with accurate point location tags. Accordingly, it furnished a rare opportunity to investigate geographic patterns of assay error at the individual point level as well as to compare these individual results with ones at an aggregated geographic resolution, specifically the census tract level. Currently, this type of dataset is scarce. A NITON XL-700-series x-ray fluorescence (XRF) instrument (NITON Corporation 900 Middlesex Turnpike, Billerica, MA 01821)was used in a chemistry laboratory to measure 15 trace metals (As, Cr, Co, Cu, Fe, Hg, Mn, Mo, Ni, Pb, Rb, Se, Sr, Zn, and Zr; see Table 1) in these soil samples, based on 120 s testing time and NIST 2711 standard reference materials (SRM). Assaying also supplied standard deviations for each of these quan-   Table 1) in these soil samples, based on 120 s testing time and NIST 2711 standard reference materials (SRM). Assaying also supplied standard deviations for each of these quantities (i.e., analytical assay errors) by soil sample. Calculated quantities are in milligrams of trace metal per kilogram of soil, or parts per million (ppm). Table 1 summarizes detection and natural background thresholds for these metals as well as the number of useful samples (i.e., with at least a non-negative assay trace metal measurement quantity) exceeding either a published maximum permissible level (MPL; [27]) or the worldwide average contamination/risk threshold [28]; all soil samples have an analytical error calculation, but not all have a valid (i.e., non-negative) trace metal quantity calculation. A number of alreadypublished studies ( [29][30][31][32][33][34]) summarize targeted analyses of the geographic distribution of this sample of trace metal measurement quantities and/or their location error. In contrast, assaying-furnished error standard deviations for all soil samples constitute a dataset yet to be analyzed, until now. A useful expectation is that a log-normal distribution describes the ordered distribution of soil contaminant errors (after [39]; also see [40]), especially for a cross-sectional study. Diagnostic results appearing in Table 2 imply that a log-normal approximation furnishes a suitable, albeit not perfect, description of the sample trace metal errors studied here. With so little residual variance (Table 2), any residual soil sample geographic resolution level spatial autocorrelation that is present has little impact upon variance estimates; this topic constitutes an appealing future research theme. The classical form of the log-normal random variable implies the following theoretical standard deviation (std), acknowledging that its minimum is zero (i.e., variance is non-negative): where r o is the rank of the existing/new value, y o , MSE denotes the appropriate regression mean squared error (see Appendix B), and Φ denotes the standard normal random variable's cumulative distribution function (the mean and variance of these n or n + 1 z-scores respectively are 0 and 1). This data feature suggests that theoretical standard errors may be posited for each of the 15 trace metal variances studied in this paper (based on Appendix B, Table A1). These quantities support a resampling-based uncertainty measure for each of the quantified analytical assay standard deviations. Anderson-Darling; * Overall, a fitted log-normal corresponded more closely than a fitted gamma distribution (an asterisk denotes the exceptions to this finding); outlier identification resulted from robust estimation using an M-estimator coupled with Huber's weight and scale. Table 3 summarizes selected descriptive spatial statistics about the trace metal analytical errors, recognizing that individual soil contaminants, because of their positive skewness (e.g., the minimum is zero), often are approximately log-normally distributed (e.g., [40]). In general, their latent spatial autocorrelation is positive and weak, which is sensible given that analytical error basically should be locationally independent across samples. Similarities common to nearby sample soil compositions could introduce some spatial autocorrelation; however, the assaying sequence should not.

Dimensions of Analytical Assay Error
One relevant research question asks whether or not the geographic distribution of trace metal soil sample analytical assay measurement error has distinct dimensions. These errors for the 15 heavy metals were analyzed with factor analysis. Table 4 reports factor loadings of the four prominent uncovered data dimensions (after being subject to a varimax rotation), which summarizes a simple structure (i.e., varimax rotated) version for them. Both the point and census tract averaged data render the same four dimensions, which account for more than 90% of the total generalized variance for the 15 trace metal analytical assay errors, with all of these errors positively correlated with their respective data dimensions. In accordance with their toxicity to living organisms, a subset of these heavy metals can be ordered as follows: Hg > Cu > Zn > Ni > Pb > Cr > Fe > Mn [41]. These four dimensions are as follows: (1) Co, Cu, Hg, Ni, Se, and Sr; (2) Co, Fe, Mn, and Rb; (3) As, Pb, and Zn; and, (4) Mo and Zr. This first dimension most likely embodies the most toxicity. Factor-1 essentially reflects variability patterns arising from biosolids and manure usage coupled with airborne vehicle pollution. Factor-2 signals variability patterns arising from the use of modern yard landscape organic and mineral fertilizers. Factor-3 signifies variability patterns typical of urban soil contamination arising from the use of outdoor wood and metal material treatments coupled with pesticide chemicals for landscape plant protection. Factor-4 may mirror variability patterns attributable to a geogenic dimension, given that lower amounts of Zr tend to be found in glacial drift soils (e.g., drumlins), whereas Mo contamination of soils is widespread in historically industrial places, particularly those that produced steel alloys (the 100+ years of production continues in Syracuse to this day in the form of specialty steel) (See Appendix C).   Figure 2c highlight much of the low density housing and other low population density land use sections of the city. Finally, the red census tracts in Figure 2d basically align with major drumlins across the city. Overall, these four map patterns are consistent with the preceding trace metal dimension interpretations. Furthermore, the spatial autocorrelation exhibited by these four dimensions signals the swarthiness tendencies rather than any coherent geographic clusters depicted by their map patterns. with the preceding trace metal dimension interpretations. Furthermore, the spatial autocorrelation exhibited by these four dimensions signals the swarthiness tendencies rather than any coherent geographic clusters depicted by their map patterns.  Figure 3 presents the geographic distributions of the analytic assay measurement errors for the 15 trace metals using a quantile classification scheme. These maps generally confirm the results of the dimension analysis. The maps for the metals in Factor-1 (i.e., Co, Cu, Hg, Ni, Se, and Sr) have a similar pattern, with high values in the northwest and southeast areas and low values in the northeast area. The four metals in Factor-2 (i.e., Co, Fe, Mn, and Rb) have a similar pattern, with the first three having high values in the northwest area and low values in the south area; Rb has a slightly different geographic pattern from the others in Factor-2, with high value clusters in the east as well as the west side  Figure 3 presents the geographic distributions of the analytic assay measurement errors for the 15 trace metals using a quantile classification scheme. These maps generally confirm the results of the dimension analysis. The maps for the metals in Factor-1 (i.e., Co, Cu, Hg, Ni, Se, and Sr) have a similar pattern, with high values in the northwest and southeast areas and low values in the northeast area. The four metals in Factor-2 (i.e., Co, Fe, Mn, and Rb) have a similar pattern, with the first three having high values in the northwest area and low values in the south area; Rb has a slightly different geographic pattern from the others in Factor-2, with high value clusters in the east as well as the west side and low values through the central north-south axis of the city. The three metals in Factor-3 (i.e., As, Pb, and Zn) have a similar geographic pattern. In particular, the maps for As and Ni are almost identical with the quantile taxonomy, depicting high values occupying the center and the northwest areas and low values occurring in the northwest and south areas. The maps for Mo and Zr generally show high values in the northern part and low values in the southern part of the city.  Figure 3 presents the geographic distributions of the analytic assay measurement errors for the 15 trace metals using a quantile classification scheme. These maps generally confirm the results of the dimension analysis. The maps for the metals in Factor-1 (i.e., Co, Cu, Hg, Ni, Se, and Sr) have a similar pattern, with high values in the northwest and southeast areas and low values in the northeast area. The four metals in Factor-2 (i.e., Co, Fe, Mn, and Rb) have a similar pattern, with the first three having high values in the northwest area and low values in the south area; Rb has a slightly different geographic pattern from the others in Factor-2, with high value clusters in the east as well as the west side and low values through the central north-south axis of the city. The three metals in Factor-3 (i.e., As, Pb, and Zn) have a similar geographic pattern. In particular, the maps for As and Ni are almost identical with the quantile taxonomy, depicting high values occupying the center and the northwest areas and low values occurring in the northwest and south areas. The maps for Mo and Zr generally show high values in the northern part and low values in the southern part of the city.

Mapping Analytical Assay Measurement Error and Sampling Error
Although the choropleth maps present the spatial patterns of the assay errors for the heavy metals, they do not consider the uncertainty of these assay errors. Studies discuss that a map classification result may not be robust when the uncertainty of observations is not considered [42,43]. This section presents choropleth maps for the heavy metals utiliz-

Mapping Analytical Assay Measurement Error and Sampling Error
Although the choropleth maps present the spatial patterns of the assay errors for the heavy metals, they do not consider the uncertainty of these assay errors. Studies discuss that a map classification result may not be robust when the uncertainty of observations is not considered [42,43]. This section presents choropleth maps for the heavy metals utilizing the approach presented in Koo et al. [22] and the separability index [21], which consider observed values and their uncertainties simultaneously in map classification. The choropleth mapping methods are implemented in the SAAR software package that utilizes R functions through R.NET [44]. It is available from https://thesaar.github.io/ (accessed on 30 April 2021). Figure 4 presents visual representations of classification results for the assay errors and their theoretical analytical measurement errors for three selected heavy metals (As, Cu, and Fe). These three heavy metals are highly associated with Factors 3, 1, and 2, respectively, reported in the dimensional analysis section ( Table 4). Note that plots for the other heavy metals are presented in Appendix D ( Figure A2). In these plots, blue dots represent assay errors for the census tracts, and the corresponding bars represent the 95% confidence intervals calculated with analytical measurement errors. The analytical measurement error is relatively consistent across the census tracts. This consistent pattern is conspicuous for Cu and Fe. In contrast, the analytical measurement error for As tends to increase as assay error increases. In other words, the length of the bars gets longer as the As value increases. Nevertheless, the analytical errors for two consecutive observations in the graph are similar. Because of this error consistency, at least between consecutive observations, the assay error values (blue dots) profoundly dominate, whereas the analytical error has little or no impact on the classification results. The classification break points (the vertical bars) appear where the assay error value separations are sizeable. The plots for the other heavy metals have this same pattern.     All heavy metals have a similarly larger variation in their resampling errors ( Figures 5  and A3). This resampling error variation has an impact on classification results. For Cu, although most observations fall into the first class in Figure 4b, most observations fall into the fifth class in Figure 5b. For As, most observations fall into the fifth class in Figure 5a, whereas observations are relatively equal in number across the fourth and fifth classes in Figure 4a.  Figure 6 presents the geographic distribution of the census tract aggregated simulated log-normal analytical measurement error for each metal. The magnitudes of the errors are represented with proportional circle sizes (the symbol sizes in the legends for their corresponding labels), whereas the choropleth maps represent mean values for the heavy metals. Some metals, such as Co (Figure 6c), display relatively sizeable error dispersions, whereas others, such as Hg (Figure 6h), display relatively small error dispersions, and yet others, such as As (Figure 6a), display more of a mixture of dispersion magnitudes. Regardless, all are well described by a log-normal distribution ( Table 2).  Figure 6 presents the geographic distribution of the census tract aggregated simulated log-normal analytical measurement error for each metal. The magnitudes of the errors are represented with proportional circle sizes (the symbol sizes in the legends for their corresponding labels), whereas the choropleth maps represent mean values for the heavy metals. Some metals, such as Co (Figure 6c), display relatively sizeable error dispersions, whereas others, such as Hg (Figure 6h), display relatively small error dispersions, and yet others, such as As (Figure 6a), display more of a mixture of dispersion magnitudes. Regardless, all are well described by a log-normal distribution ( Table 2). Figure 7 reveals that the resampling error dispersion tends to be substantially larger than its corresponding analytical measurement error. Sampling error across the Figure 7 maps partially reflects variation in the post-tabulated sample sizes, which is far from constant. Accordingly, many maps (e.g., Figure 7c,h, respectively portraying Co and Hg) have relatively wide ranges of dispersions, which are quite conspicuous. Some (e.g., Figure 7b,k) have relatively narrow ranges of small dispersions, whereas none have relatively constant and more substantial resampling errors (similar to Figure 6b).
The geographic error variability displays some noticeable features and evokes several implications. First, analytic measurement error generally is stable across census tracts. Most of the heavy metals have essentially the same analytical error for all of the census tracts, although some heavy metals (As, Pb, and Zn) show an increasing pattern that covaries with the assay error. Second, the resampling error varies considerably, unlike the analytical measurement error. The geographic distributions of the resampling error show pattern similarity to the factors uncovered as data analytic dimensions of the assay error. The heavy metals that are mainly associated with the same factor in Table 4 have a map pattern similar to that of the resampling error. For example, the geographic patterns for Co, Fe, Mn, and Hg, which mainly contribute to Factor-2, are alike. This outcome may imply that resampling error results can be affected by samples in an areal unit, consequently affecting its geographic pattern. Int. J. Environ. Res. Public Health 2021, 18, x 11 of 28   Figure 7 reveals that the resampling error dispersion tends to be substantially larger than its corresponding analytical measurement error. Sampling error across the Figure 7 maps partially reflects variation in the post-tabulated sample sizes, which is far from constant. Accordingly, many maps (e.g., Figure 7c,h, respectively portraying Co and Hg) have relatively wide ranges of dispersions, which are quite conspicuous. Some (e.g., Figure  7b,k) have relatively narrow ranges of small dispersions, whereas none have relatively constant and more substantial resampling errors (similar to Figure 6b).  The geographic error variability displays some noticeable features and evokes sev eral implications. First, analytic measurement error generally is stable across census tract Most of the heavy metals have essentially the same analytical error for all of the censu tracts, although some heavy metals (As, Pb, and Zn) show an increasing pattern that co varies with the assay error. Second, the resampling error varies considerably, unlike th analytical measurement error. The geographic distributions of the resampling error show

Error Mixtures in Data: Some Simulation Experiments
This section summarizes output from exploratory simulation experiments addressing the sensitivity of map patterns to analytical assay measurement error (via sampling from appropriate log-normal random variables), sampling error (via bootstrap resampling), and their combination. Table 4 reports census tract resolution summary results for the observed data. Concluding comments focus on specification error with regard to estimating the nature and degree of spatial autocorrelation in geographic dimensions of data.

Geographic Dimension Sensitivity to Analytical Assay Measurement Error
Generation of the assaying sensitivity analysis simulation experiment results summarized in Table 5 utilized random samples from appropriate heavy metal log-normal distributions for each of the 3322 individual soil samples located in the city of Syracuse. These census tract averages then were georeferenced data input to a varimax-rotated factor analysis. Results for the data analytic factor structure of the 15 Syracuse, NY, trace metal measurement errors imply that changes in assay measurement error based on theoretical frequency distributions corrupt findings in subtle ways (e.g., average Factor-2 and Factor-3 loadings decrease, but not significantly), although it preserves simple structure as well as the dimensional clustering of heavy metals. Nevertheless, the percentage of variance accounted for by each dimension remains essentially the same as that for the observed data. Table 5. Log-normal distribution simulated Syracuse trace metal analytical assay measurement error varimax-rotated factor dimensions: census tract averaged soil sample data (10,000 replications).

Trace Metal
Factor-  Table 4); s r denotes the standard deviation of a simulated assaying data correlation coefficient.

Geographic Dimension Sensitivity to Sampling Error
This resampling experiment addressing sensitivity to soil sampling error employed a bootstrap census tract tessellation stratified random sampling (with replacement) design. Table 6 summarizes sampling error resampling findings for the data analytic factor structure of the 15 Syracuse trace metal measurement errors. Sampling error corrupts results in various ways. Four heavy metals exhibit a propensity to load onto two factors, compromising simple structure. Now Zn loads onto the first dimension (which arguably occurs for the collected data), and Rb fails to load onto any of the first four dimensions. Prominent loadings for the second and third dimensions markedly decrease (although not significantly). Meanwhile, as with the assaying simulation findings, the percentage of variance accounted for by each dimension remains essentially the same as that for the observed data. NOTE: bold font denotes a prominent factor loading (see Table 4); S r denotes the standard deviation of a resample correlation coefficient.
A third experiment, addressing sensitivity to specification error, employed the gamma distribution (with LN r o −3/8 n−r o +5/8 + 0.25 as its covariate; representing minimal specification error) because it furnishes a competitive empirical distribution to the log-normal one. It also employed the beta distribution (with LN r o −3/8 n−r o +5/8 as its covariate; representing noticeable specification error), which provides a conceptually appealing alternative, given that the trace metal quantities are proportions in ppm. It additionally employed the uniform distribution (with p i = r i −3/8 n+1/4 as its covariate; representing severe specification error), whose probabilities are much smaller for the concentration portion and much larger for the right-hand tail of the positively skewed empirical frequency distribution. Based upon order statistics combined with six-sigma theory-e.g., covariate values have a variable transformed beta distribution with variance r i (n−r i +1) (n+2)(n+1/4) 2 for the ith rank-sampling draws were from the uniform distribution over the interval where p i (which also is the covariate) denotes the empirical cumulative distribution function probability for soil sample i. In each case, 10,000 analytical assay measurement error samples were drawn. Tables 7-9 summarize averaged varimax-rotated factor analysis results from these experiments. Meanwhile, Figure 8 portrays selected extreme gamma and beta distributions embraced by these simulations.
Approximating the analytical assay measurement error that is a proportion restricted to the interval [0, 1] with a log-normal distribution whose support is the interval [0 , ∞) appears to introduce little specification error, most likely because the skewness of heavy metal error tends to concentrate its values very close to zero. However, specification error introduced by replacing this log-normal with its gamma distribution competitor [45] corrupts a number of results (see Table 7): Mn loads onto the first rather than the second dimension; Pb and Zn fail to load onto any of the first four dimensions; all of the prominent loadings are lower, some substantially so, verging on statistical significance; and the percentage of variance accounted for decreases noticeably for the first three dimensions. Figure 8 portrays the scatterplot association between the log-normal and gamma individual heavy metal goodness-of-fit regression (pseudo-)R 2 values (also see Appendix E, Table A3); overall, the gamma assumption appears to be inferior to the log-normal assumption.  Table 4); S r denotes the standard deviation of a combined simulated assaying and resampling error correlation coefficient.  Approximating the analytical assay measurement error that is a proportion restricted to the interval [0, 1] with a log-normal distribution whose support is the interval 0, ∞) appears to introduce little specification error, most likely because the skewness of heavy metal error tends to concentrate its values very close to zero. However, specification error introduced by replacing this log-normal with its gamma distribution competitor [45] corrupts a number of results (see Table 7): Mn loads onto the first rather than the second dimension; Pb and Zn fail to load onto any of the first four dimensions; all of the prominent loadings are lower, some substantially so, verging on statistical significance; and the percentage of variance accounted for decreases noticeably for the first three dimensions. Figure 8 portrays the scatterplot association between the log-normal and gamma individual heavy metal goodness-of-fit regression (pseudo-)R 2 values (also see Appendix E, Table  A3); overall, the gamma assumption appears to be inferior to the log-normal assumption.
Meanwhile, although the beta distribution and analytical assay measurement error share the same support interval, the beta distribution's performance is poorer than that produced by the gamma assumption (Tables 5, 7, and 8). Here the first dimension gains Co and Mn and loses Hg and Sr. The second dimension resembles the original third dimension, whereas the third dimension lacks substance. Zn fails to load onto the fourth dimension; more specifically, Fe, Hg, Rb, Se, and Zr fail to load onto any of the first four dimensions. Again, all of the prominent loadings are lower, some substantially so, verging on statistical significance. Figure 9 also portrays the scatterplot association between the log-normal and beta individual heavy metal goodness-of-fit regression (pseudo-)R 2 values (see Appendix E, Table A3); overall the beta appears inferior to the gamma, and markedly inferior to the log-normal, assumption. Meanwhile, although the beta distribution and analytical assay measurement error share the same support interval, the beta distribution's performance is poorer than that produced by the gamma assumption (Tables 5, 7 and 8). Here the first dimension gains Co and Mn and loses Hg and Sr. The second dimension resembles the original third dimension, whereas the third dimension lacks substance. Zn fails to load onto the fourth dimension; more specifically, Fe, Hg, Rb, Se, and Zr fail to load onto any of the first four dimensions. Again, all of the prominent loadings are lower, some substantially so, verging on statistical significance. Figure 9 also portrays the scatterplot association between the log-normal and beta individual heavy metal goodness-of-fit regression (pseudo-)R 2 values (see Appendix E, Table A3); overall the beta appears inferior to the gamma, and markedly inferior to the log-normal, assumption.
Approximating the analytical assay measurement error that is a proportion to the interval [0, 1] with a log-normal distribution whose support is the inter appears to introduce little specification error, most likely because the skewnes metal error tends to concentrate its values very close to zero. However, specifica introduced by replacing this log-normal with its gamma distribution competito rupts a number of results (see Table 7): Mn loads onto the first rather than t dimension; Pb and Zn fail to load onto any of the first four dimensions; all of t nent loadings are lower, some substantially so, verging on statistical significanc percentage of variance accounted for decreases noticeably for the first three di Figure 8 portrays the scatterplot association between the log-normal and gamm ual heavy metal goodness-of-fit regression (pseudo-)R 2 values (also see Append A3); overall, the gamma assumption appears to be inferior to the log-normal as Meanwhile, although the beta distribution and analytical assay measurem share the same support interval, the beta distribution's performance is poorer produced by the gamma assumption (Tables 5, 7, and 8). Here the first dimen Co and Mn and loses Hg and Sr. The second dimension resembles the origina mension, whereas the third dimension lacks substance. Zn fails to load onto dimension; more specifically, Fe, Hg, Rb, Se, and Zr fail to load onto any of the dimensions. Again, all of the prominent loadings are lower, some substantially s on statistical significance. Figure 9 also portrays the scatterplot association be log-normal and beta individual heavy metal goodness-of-fit regression (pseud ues (see Appendix E, Table A3); overall the beta appears inferior to the gamma, a edly inferior to the log-normal, assumption. Figure 9. Goodness-of-fit (pseudo-)R 2 scatterplots labelled by heavy metal (each label pertains to a horizontal row of three scatterplot points). Denoted distribution results are: solid black circles denote gamma; solid gray circles denote beta; and, Hebrew pictograms tet (⊗) denote uniform. Figure 9 endorses the contention that the uniform distribution represents the most extreme specification error assumption studied here. Table 9 reveals that the simulation error is negligible because selection is from relatively narrow intervals (regardless of the use of the six-sigma principle). This assumption essentially preserves the multivariate statistical structure of Table 5 but with an overfitting cost (i.e., analyses over-emphasize assumption-specific information). The second and third dimensions switch (which is not surprising, given that the respective percentages of accounted variance are almost the same), and many of the loadings decrease somewhat. Nevertheless, a number of the bivariate regression R 2 values are only moderate in degree and hence relatively low (see Appendix E, Table A4), with corresponding coefficients that deviate considerably from an intercept of 0 and a slope of 1.

Selected Error Propagation Illustrations
This section reviews three instances of error propagation. The first, motivated by Koo et al. [46], assesses analytical measurement error, resampling error, and specification error on indices of spatial autocorrelation. The second, inspired by a widespread recognition of the existence of interacting error sources, evaluates a mixture of sampling and specification error latent in georeferenced data. The third illustrates variance inflation, a common geospatial data complication. Table 10 summarizes Moran Coefficient (MC) and Geary Ratio (GR) spatial autocorrelation index results by data factor dimensions. With regard to spatial autocorrelation measures, based upon a root mean squared error criterion, the presumably correct statistical distributional assumption results most closely align with the original data results. Sampling error noticeably impacts these measures, with a tendency to decrease them. Specification error markedly impacts these measures, with the beta and uniform distribution assumptions corrupting them far more than the gamma assumption. The MC indicates that the beta, whereas the GR indicates that the uniform, assumption introduces more specification error. Meanwhile, the standard errors display considerably more of an impact than the spatial autocorrelation indices themselves. Interestingly, the resampling error does not tend to exhibit the closest alignment with the theoretical standard errors reported for the original data. Specification error introduced by the uniform distribution assumption deviates the most, which is rather obvious from a visual inspection of the table entries. The asymptotic standard error (i.e., 0.85) for the MC is consistent with these results; that for the GR (i.e., 0.214) is not.

Error Propagation to Spatial Autocorrelation Indices
The general implication here is that error propagation impacts spatial autocorrelation index values and their standard error estimates. Specification error may or may not be a more serious source of this corruption. This topic merits considerably more future research.

Error Propagation from a Mixture of Error Sources
Another experiment, inspired by Gustavsson et al. [47], involved the following twostep procedure: draw a bootstrap census tract tessellation stratified random sample (with replacement) of 3322 log-normal fitted values, and then use these measurements to generate simulated (log-normal distribution assumption) assaying analytical errors. This sequence was executed 10,000 times. Table 11 summarizes output for this combined error source sensitivity analysis. As expected, deviations from the original data results increase, with sampling error markedly dominating analytical measurement error (also see Table 12).  Table 4); S r denotes the standard deviation of a combined simulated assaying and resampling error correlation coefficient.

Error and Spatial Autocorrelation Induced Variance Inflation
The final exploratory analysis concerns an additional evaluation of error, partially reminiscent of the arguments of Koo et al. [46]. A principal impact of positive spatial autocorrelation is variance inflation. Therefore, a conventional variance formula includes specification error, whereas adjusting for spatial autocorrelation renders an uninflated variance estimate; adjustment here is with the popular pure spatial simultaneous autoregressive (SAR) model specification. Table 12 shows the relevant results. The various pairs of data results demonstrate how variance inflation increases with increasing positive spatial autocorrelation; this inflation is rather modest for this example because the spatial autocorrelation parameter ρ is not vary large (e.g., approximately 0.4), generating roughly a 20% increase in variance. Random sampling error impacts variance estimates more than the other reported statistical quantities; Table 12 suggests that it also might slightly dampen variance inflation vis-à-vis the other error sources. In contrast, random measurement error has similar but less pronounced impacts on the reported statistical quantities; it appears to have a rather neutral impact on variance inflation. When mixed, random sampling error impacts tend to dominate analytical measurement error impacts. To quantify this situation based on the Table 12 information, for Factor-1.
The corresponding mixture result confirms this outcome. Results for the other three factors are more ambiguous.
Because the gamma distribution often is an assumption competitor for the log-normal distribution, it furnishes the basis for a minimum specification error assessment here. Specification error impacts are detectable in Factor-1 and considerable in the other three factors. These impacts are at least as substantial as those for sampling error, being far more extreme for Factor-4. Accordingly, in terms of trace metal analytical error, the evidence presented in this paper supports the following rank ordering: measurement < sampling < specification This ordering merits further research. Table 13 furnishes a basic summary revealing prominent similarities and differences displayed by the presence of measurement, sampling, and specification error. The random nature of these error sources partially functions as noise, deteriorating latent patterns. Regardless, Factor-4 is a robust dimension, and Factor-1 is a reasonably robust dimension, with error-created confusion primarily materializing in the second and third factors.

Discussion
Sampling error tends to dominate analytical measurement error, with specification error introducing serious corruption. Meanwhile Rb results appear to be the most sensitive, of the 15 heavy metals studied here, to sources of error.

Conclusions and Implications
This paper investigates assaying errors in surface soil samples that were collected across the city of Syracuse, NY, mostly during the summers of 2003 and 2004. These samples are composed of a total of 3628 observations, with measurements for 15 heavy metals. The geographic distributions of the 15 heavy metals were examined at the individual point level and at the aggregate census tract levels. In addition, the data analytic dimensions of these heavy metals were examined using varimax-rotated factor analysis. The results show that the assaying errors do not present a high level of spatial autocorrelation. This result may indicate that the assay errors generally are independent across geographic locations. The factor analysis results imply that the geographic distributions of the 15 heavy metals display four prominent data analytic dimensions and that their patterns are associated with objects in the physical environment (e.g., expressways), historical settlement patterns in the city, and low housing/population densities. In addition, these assay error patterns are consistent with those of the corresponding observed values of the 15 heavy metals. Empirically, this outcome shows that the size of assay errors tends to be positively associated with the sizes of observed values.
This paper also investigates impacts of sampling error, measurement error, and their mixture on assaying error for the 15 heavy metals using explanatory simulation experiments. The results show that although both sampling error and measurement error have an impact on assaying errors, measurement error has a greater impact than sampling error. Furthermore, geographically, the measurement error is stable across the census tracts, whereas the sampling error varies considerably. These findings may provide useful insights for constructing a more rigorous sampling design as well as a measurement plan. Also, the results suggest that specification error would have a severe impact; however, this conjecture needs further investigation to better understand this impact in a general context. Finally, the results of the simulation experiments confirm a geographic resolution impact.
Because this paper investigates a specific dataset, the findings may not be directly applicable to other geographic landscapes. Nevertheless, this paper shows impacts of the different error sources on assaying error and posits a potential importance ranking of different error sources.
Author Contributions: The co-authors D.A.G. and Y.C. meet the authorship criteria and are in agreement about the submission of the manuscript. They collaboratively conceived the research idea, designed the methodology, and processed the data, and then executed different parts of the statistical analyses. Both contributed to the writing of this paper, Both authors have read and agreed to the published version of the manuscript.

Appendix B. Tabulation of Important Back-Transformation Statistics
Griffith [48] discusses back-transformations, specifically commenting about transformation-introduced specification error that can be indexed by the mean of the ratioŷ i /y I appearing in Table A1. All of the Syracuse heavy metals contain less than 1% specification error based upon this criterion. Furthermore, 11 of the heavy metals have a bivariate regression slope of approximately 1, its expectation for a perfect match between the two sets of values. In particular, Zn exhibits what appears to be excessive bivariate regression parameter values, implying that its log-normal characterization may contain excessive specification error.