Triple-stack geostatistical modeling for urban injury: Integrating grids,built environment,and Poisson kriging for pedestrian fatalities in Cali,Colombia

Abstract

We introduce a novel geostatistical framework to map and predict pedestrian fatalities in Cali, Colombia (2008–2010), addressing critical gaps in research on built environment risk factors particularly in low- and middle-income countries. We triple stack: (1) a grid-based global moving window reducing the modifiable areal unit problem, redistributing events and population; (2) a Poisson-based land use regression incorporating built environment data; and (3) Poisson kriging addressing small number issues. Nine built environment variables were significant predictors of pedestrian fatalities. Health centers had the highest risk increase (54%), while parks and population density were protective. The triple-stack model markedly improved risk prediction, sharpening hotspot delineation and reducing background noise; at high population densities, predictive errors were as low as 0.26 deaths. The framework is a transferable tool for other urban contexts; thus, we provide method decision-making support helping users select between a simplified, codeless, GUI based platform or an advanced, custom-coded approach. Our results demonstrate gains in model predictive accuracy may be primarily attributable to the first two layers, the grid and LUR. Although, Poisson kriging is traditionally used to address the “small number problem,” the first two layers may substantially stabilize rates, resulting in more realistic, interpretable risk maps.

1. Introduction

Traffic-related injuries and fatalities among pedestrians remain a significant global concern.^1–5 Approximately 1.19 million people die globally each year from road traffic crashes, making them the leading cause of death for people aged 5–29 years.^4,5 Although road safety improvements have been implemented, challenges persist, even in many high-income countries. For example, in the United States, pedestrian deaths due to motor vehicle crashes rose by 59% between 2009 and 2020. ² Globally, more than half of all fatalities occur among vulnerable road users, including pedestrians, cyclists, and motorcyclists.^2–5 These persistent trends have complicated international road safety targets. The United Nations’ goal to halve road traffic injuries and fatalities between 2015 and 2020 was not met. This shortfall prompted renewed global commitments to achieve a 50% reduction by 2030, with responsibility assigned across governmental, corporate, and civil sectors.^6–10

In the Latin America and Caribbean (LAC) region, traffic deaths have remained persistently high despite progress in other health domains. Unlike in many high-income countries where road safety improvements have contributed to declining fatality rates, most low- and middle-income countries have seen little progress.^1,3,11,12 Road traffic crashes claimed approximately 110,000 lives in the region in 2016, ranking as the fifth leading cause of disability-adjusted life years lost. ¹ Furthermore, among individuals aged 15–49 years, traffic injuries were the second leading cause of death after interpersonal violence.^1,12 By 2021, countries in the LAC region reported approximately 107,000 annual traffic-related deaths. These deaths represented roughly 9% of global fatalities, indicating a persistently high regional burden.^3,11 Despite modest global declines in traffic-related fatalities, pedestrian safety remains a persistent crisis in the LAC region. This underscores the need for evidence-based strategies tailored to its urban contexts.^1,12

Built environment factors including roadway design, traffic calming measures, and land use are crucial determinants of pedestrian safety.^3,5,9,10 Even small modifications can substantially reduce risk.^13–21 However, pedestrian safety research often lacks comprehensive longitudinal exposure data, particularly in low- and middle-income countries, prompting calls for localized studies to fill critical evidence gaps.^{1,2,12,20,22,23} Additionally, most pedestrian fatality research focuses primarily on individual-level behaviors, neglecting geo-temporal dynamics despite their proven efficacy in public health. By identifying fine-scale spatial clusters and temporal shifts, geostatistical models can help policymakers better understand localized risk patterns and inform planning decisions.^24–36

A prior study examined pedestrian fatalities in Cali, Colombia, where pedestrian deaths represent a substantial road safety burden. The study used Bayesian Maximum Entropy (BME), an advanced geostatistical technique which yields kriging as its linear limiting case, and incorporated population smoothing via Poisson kriging. ³⁰ Although the analysis was informative, it relied on administrative boundaries and did not integrate built-environment variables. These limitations reduced its ability to capture fine-scale spatial variation and made the results vulnerable to the small number problem, boundary sensitivity, discontinuities, and the Modifiable Areal Unit Problem (MAUP).^{6,28,30,33,37}

In this study, we analyze pedestrian fatalities in Cali, Colombia, from 2008 to 2010. We build on the prior study by improving fine-scale estimation of pedestrian risk and evaluating how built-environment characteristics contribute to spatial variation. In doing so, this study addresses both study-specific limitations and a broader methodological gap in spatial epidemiology. We introduce the Triple Stack (TS) approach, a unified, high-resolution analytic framework that integrates three geostatistical components. The three layers are: Layer One, overlapping static grid-based aggregation; Layer Two, Land Use Regression (LUR); and Layer Three, population-adjusted Poisson kriging (PK).

Each layer enriches the model with additional information to improve its estimation ability. Layer One replaces administrative units with overlapping grid-based incidence areas, allowing both pedestrian fatality counts and population denominators to be estimated over a consistent spatial support. This structure may improve alignment with underlying residential patterns and reduce instability caused by sparse event data. Layer Two incorporates built-environment predictors to explain systematic variation in risk.^{25,27,29,31,33,38–40} Layer Three applies population-adjusted geostatistical smoothing to stabilize rare-event rates and reduce misleading hotspots caused by small population denominators. This challenge is common in rare-event modeling contexts such as pedestrian fatalities, syphilis, HIV, and Guillain-Barré Syndrome.^{6,28,30,31,33,37,41,42}

To evaluate the added value of this framework, we compare the Triple Stack approach with a simpler Poisson kriging-only model. This comparison allows the assessment of improvements in model accuracy when integrating grid-based aggregation with overlapping rates areas (Layer 1) and built-environment predictors (Layer 2). The study also provides a step-by-step walkthrough of the Triple Stack geostatistical workflow. This walkthrough aims to lower methodological barriers and support broader application of advanced spatial modeling in urban health and safety research. It also offers practical guidance on spatial model selection.

A major barrier to broader use of BME and related geostatistical methods is the workflow, terminology, and computational steps can appear opaque to interdisciplinary users. This can limit uptake even when the underlying modeling decisions are conceptually interpretable. To address this barrier, we include a conceptual schematic of the modeling process. Additionally, we clarify key BME and geostatistical terminology for readers who may wish to interpret, evaluate, or apply similar methods.

For users seeking a more accessible entry point into geostatistical modeling, BMEGUI provides a free, standalone platform. It supports basic space-time BME analyses and map generation without advanced programming. ³¹ However, BMEGUI does not currently implement the full Triple Stack framework, including grid-based rate construction and LUR-integrated mean-trend modeling. We therefore provide a comparative guide that contrasts analyses feasible in BMEGUI with the Triple Stack approach across analytic complexity, data needs, implementation effort, and predictive value. This guide supports the selection of appropriate tools for equity-sensitive planning.

Although the data originate from the early 2010s, the study remains highly relevant given the recent surge in pedestrian mortality.^3,4,18 Local interventions, development projects, and infrastructure changes may have occurred since the study period, particularly given ongoing urban-environment and health-disparity concerns documented in Cali. ⁴³ However, the spatial patterns estimated by the Triple Stack framework are likely to remain relevant because key pedestrian-safety-related built-environment features, including roadway configuration, land-use structure, density, and destination accessibility, typically change gradually rather than abruptly. Prior research has therefore measured built environment change across multi-year or decadal periods. Broader built-environment scholarship also emphasizes the long-term evolution of urban form, infrastructure, and neighborhood systems across extended temporal scales.^44–46

2. Methods

2.1. Triple stack framework

The analytical framework consists of three-layered methods that each add information to the model, together forming the Triple Stack approach (Figure 1). In Layer 1, pedestrian fatalities, population data, and built-environment variables are integrated on a regular spatial grid using overlapping circular rate areas. This common spatial support reduces sensitivity to administrative boundaries, MAUP effects, and edge effects. It also reduces kriging islands and may reduce small-number instability.^6,28,31,33 The resulting grid structure generates population-normalized incidence rates and standardized spatial predictors for subsequent modeling.

Figure 1.

Conceptual workflow of the Triple Stack (TS) framework. The framework integrates three sequential layers. Layer 1 constructs grid-based incidence rates and built-environment predictors using a regular grid, overlapping circular rate areas, redistributed population denominators, and GIS-derived covariates. Layer 2 estimates the global offset of the process. The global offset consists of the sum of a spatial mean trend and a temporal mean trend. These mean components are then removed to generate a residual spatiotemporal field. Layer 3 models the residual field using Bayesian Maximum Entropy (BME)/Poisson kriging. In this layer, the covariance structure defines a globally informed prior probability distribution, which is updated using nearby observed data to produce local posterior estimates of the residual field. In the final step, the spatial and temporal mean trends are then added back to the BME residual estimates to generate the final smoothed risk surfaces and support cross-validation-based model comparison.

In Layer 2, the broad background pattern in the data was estimated. A Poisson land use regression (LUR) model was used to estimate the spatial mean trend of pedestrian fatalities, with built-environment characteristics included as predictors. An exponential kernel-smoothed temporal mean trend was estimated separately to capture the dataset’s overall temporal structure. The spatial and temporal mean trends were then combined to define the dataset’s global offset.^{24,25,27,35,36,38–40}

In Layer 3, the global offset was removed from the observed incidence rates to create a residual field. A space-time covariance model was then estimated from this residual field and used to define the prior dependence structure for Bayesian Maximum Entropy (BME)/Poisson kriging (PK). BME/PK then uses the covariance model and observed residual data to condition the prior structure and generate local residual estimates. This process uses Bayesian knowledge blending to reflect both the general knowledge base, represented by the global offset and covariance model, and the influence of nearby residual observed data. The global offset was added back to the BME residual estimates to generate the final BME estimates of incidence rates. These estimates produced the final maps, or smoothed risk surfaces.^{24,25,27,29,31,35–40,42,47}

Conceptually, the baseline model can be summarized as PK = census-sector support + Poisson kriging only, whereas the Triple Stack model can be summarized as TS = grid-based incidence areas + LUR mean trend + Poisson kriging. Specifically, the baseline PK model only used Colombian census sectors as the spatial support. Thus, the comparison isolated the added value of replacing administrative boundaries with grid-based incidence areas (Layer 1) and incorporating built-environment predictors before geostatistical estimation (Layer 2).

2.2. Data sources

To implement the Triple Stack framework, three datasets were used: 1) pedestrian mortality data, 2) census population and administrative areas, and 3) urban infrastructure. Pedestrian mortality data were obtained from the Cali Injury Mortality Surveillance System.^48,49 This dataset includes the location and date of all transport-related fatal injuries in the urban area of Cali from January 2008 through December 2010. A total of 356 pedestrian deaths occurred during this 36-month period, and 97% (345 of 356) were successfully geocoded to an address using GEOBIS EZU software. ⁵⁰ To study the built environment, we used Colombian Census data provided by the Administrative Department of Municipal Planning of Cali. This dataset consists of GIS shapefiles containing administrative boundaries, population data, and geographic infrastructure information. Infrastructure variables included roads, libraries, hospitals, schools, community centers, police and fire stations, sidewalks, and fire hydrants. A full list of all available infrastructure variables examined in this study is provided in Supplemental Materials Table 1.

All computational analyses were conducted on UNC’s LongLeaf compute cluster using MATLAB, ⁵¹ MATLAB’s parallel computing toolbox ⁵² and with the BMElib 2.0c package, a freely available MATLAB library.^52,53 A simplified BME analysis can also be conducted in BMEGUI. This free standalone platform supports basic space-time BME modeling and map generation without advanced programming. ³¹ Both BMElib and BMEGUI can be freely downloaded at https://mserre.sph.unc.edu/BMElab_web/. However, BMEGUI does not currently implement the full Triple Stack framework nor support high performance computing environments. ⁵³

2.3. Layer One: Grid construction, incidence rate calculation, and built-environment predictors

As the first layer of the Triple Stack framework, we constructed a regular grid across the urban area of Cali. Parameter choices for the circular rate-area radius, grid spacing, and SEDC influence distance were guided by prior Cali-specific spatiotemporal analyses, sensitivity testing, and prior grid-based spatial epidemiologic work. Previous BME analyses showed pedestrian fatality hotspots in Cali were highly localized, with limited spatial influence beyond 500–750 m. Thus, circular rate areas were defined using a 500 m radius.^28,30 Grid centroids, indexed as i, were spaced 650 m apart, equal to 0.65 times the diameter of each rate area. Each grid point served as the centroid of a circular rate area (A _i ). Incidence rates (IR _i ), pedestrian fatality counts (Y _i ), population denominators (n _i ), and land use characteristics were then calculated within each rate area (Eqs. 1–2). Unless noted otherwise the fatality counts are summed over the area of interest and over the study duration, and the population is estimated at the midpoint of the study duration.

This spacing intentionally produced overlapping windows, allowing spatial transitions to be smoothed while preserving fine-scale variation. ²⁸ In contrast to discrete areal aggregation approaches, overlap among circular rate areas does not constitute double counting. Rather, each grid point defines a distinct spatial support over which a local rate is estimated. Dependence among overlapping supports is directly accounted for and explicitly modeled within the BME covariance structure.^{25,27,34,36,38} The overlap is therefore a deliberate smoothing and support-design feature. It reflects the continuous nature of the underlying spatial process, which is a central assumption of geostatistical modeling. Unlike traditional statistical approaches that often treat observations as independent, geostatistical methods explicitly model spatial dependence among nearby observations.^{25,27,28,34,36,38}

Next, the individual built-environment predictors, j, such as libraries, fire hydrants, roads, or sidewalks were analyzed and linked to each rate area (A _i ) for the Poisson LUR. All 18 infrastructure variables initially analyzed are listed in Supplemental Table 1, and the significant predictors used in our final model are shown in Table 1. Built-environment data were provided as ESRI shapefiles and categorized as points, lines, or polygons. For the LUR model, each predictor value (X _ij ) represents the value of built-environment predictor j at rate area i. The predictors have different feature geometries, thus X _ij is calculated differently for point, line, and polygon features.

Table 1.

Statistically significant built environment predictors from the Poisson Land Use Regression (Layer 2), with back-transformed incidence rate ratios (IRRs), log-scale coefficients, and corresponding standard errors. For point-based predictors, exposure was quantified using SEDC, a dimensionless, spatially weighted index reflecting the proximity and density of nearby built-environment features.

Variable	P-value	β: Back transformed (log)	Directionality	SE: Back transformed (log)
Intercept	<.0001	1.9e-5 (-10.873)	Protection	1.24 (0.2177)
Libraries	<.0001	1.01 (0.005)	Escalation	1.06 (0.054)
Fire Hydrants	<.0001	1.06 (0.0541)	Escalation	1.01 (0.0108)
Health Centers	<.0001	1.54 (0.432)	Escalation	1.09 (0.0886)
Educational Centers	0.0225	1.05 (0.052)	Escalation	1.01 (0.0092)
Open Areas/Parks	0.0072	0.92 (-0.086)	Protection	1.03 (0.031)
Sidewalk Density	0.0044	1.0 (0.00016)	None	1.00 (9.486E-05)
Roads (Primary & Secondary)	<.0001	1.03 (0.0296)	Escalation	1.01 (0.0074)
Population Density	<.0001	1.0 (-8.51E-10)	None	1.00 (1.665E-10)

For predictors represented as point data, such as libraries and fire hydrants, we calculated the Sum of Exponentially Decaying Contribution (SEDC _i,j ). SEDC is a dimensionless, spatially weighted exposure index. It allows for built-environment context to be incorporated beyond the 500 m rate area. Features closer to the centroid receive greater weight, with influence decreasing exponentially to approximately 5% at the 750 m radius. Each individual point feature (k) within the point predictor category (j) is at a distance D _i , _k , _j from grid centroid i. For predictors (j) represented as line features, such as roads and sidewalks, segment lengths were calculated within each rate area (A _i ). These segment lengths (l) were then summed to estimate total line-feature density (L _ij ) within each rate area (A _i ).

2.4. Calculation details for key variables

To improve accessibility for interdisciplinary readers, Supplemental Materials Table 2 summarizes the notation used in Equations 1–5. Incidence rates were calculated differently for the LUR and BME modeling stages. For the LUR model, we calculated full-period incidence rates, while for the BME analysis, we calculated 9-month rolling-window incidence rates. These differences are described in the geostatistical modeling section (2.5).

Parameter selection was informed by both prior published work and additional sensitivity analyses conducted for the current study. Earlier published work on Cali pedestrian fatalities compared multiple spatial and temporal aggregation levels using exploratory maps, mean-trend patterns, and covariance diagnostics.^28,30 For the current analysis, additional sensitivity analyses examined alternative radii for key built-environment covariates, including roads, hydrants, and health centers. Across these checks, the main hotspot structure and dominant spatial patterns were broadly stable across a reasonable range of parameter values. These analyses supported the use of fine-scale spatial supports and 9-month temporal windows as a balance between retaining spatial resolution and limiting excess variability.

These parameters should therefore be interpreted as empirically tuned to the Cali context rather than arbitrarily selected. The general framework remains transferable to other cities and public health outcomes, but parameter values should be locally calibrated rather than directly reused. Smaller radii and tighter grid spacing preserve local variation but increase sensitivity to sparse-event noise. Larger radii and broader spacing increase smoothing and stability but can reduce fine-scale hotspot delineation. The selected parameters were intended to balance this bias-variance trade-off for the sparse Cali pedestrian fatality data.

For the LUR model, a full-period incidence rate (IR _i ) was calculated for each grid point i (Eq 1). Pedestrian deaths were aggregated within the circular rate area centered at i using a radius of r = 500 m. This count, Y _i , represents the number of unique pedestrian deaths within the rate area across the 36-month study period. The population denominator for each rate area is represented by n _i and and T represents the duration of the observation window. For the full-period LUR model, T = 36months, equivalent to 3 person-years of exposure.

${IR}_i=\frac{{\mathrm{Y}}_i}{{\mathrm{n}}_{i\hspace{0.17em}}T}$

1

The population denominator, n _i , was calculated using areal redistribution of census population data (Eq. 2). This approach estimated the portion of each census sector’s population falling within the circular rate area, A _i , centered at grid point i. The radius, r, was 500 m. Therefore, A _i = πr² = 0.785 km². A_a denotes the area of census sector, a that overlaps A _i , and P _Aa denotes the population of that sector. Population was assumed to be uniformly distributed within each census sector.

$n_i={\displaystyle \sum }_{a\in A_i}\left(\frac{\Vert A_{a\hspace{0.17em}}\cap A_i\Vert }{\Vert A_{a\hspace{0.17em}}\Vert }\right)*P_{A_a}$

2

Equation 3 shows the Poisson LUR model used to estimate the spatial mean trend in pedestrian fatality risk from built-environment predictors. The broader Triple Stack framework can incorporate different regression models for the mean-trend component. In this study, a Poisson LUR was used because pedestrian fatalities are sparse count data modeled as population-normalized incidence rates. The log-link structure of the Poisson model ensures predicted rates remain non-negative. Equation 3 estimates the expected incidence rate for each rate area based on its built-environment predictors. This predicted value is the LUR-based spatial mean trend (LUR_i) used in the subsequent geostatistical model. In Equation 3, β ₀ is the intercept, X _ij is the value of built-environment predictor j at rate area i, β_j is the coefficient for predictor j, and p is the total number of predictors included in the model.

${LUR}_i=\exp \left({\beta }_{0\hspace{0.17em}}+{\displaystyle \sum }_{j=1}^p{\beta }_{j\hspace{0.17em}}{X}_{i,j\hspace{0.17em}}\right)$

3

For point-based predictors in the LUR model, the predictor value X _ij in Equation 3 was calculated using the Sum of Exponentially Decaying Contribution (SEDC_ij, Eq. 4). SEDC_ij was calculated for each grid point i and point-based predictor j. It sums the exponentially decaying distance contributions of all individual features belonging to that predictor type. For example, if j represents health centers, then k indexes each individual health center included in the calculation. D_{i ,k} ,_j represents the distance from the centroid, i to individual feature, k of predictor type, j. Features closer to the centroid received greater weight, with influence decreasing exponentially to approximately 5% at the influence range, α, set here to 750 m.

${SEDC}_{i,j}=\underset{k}{\sum e^{-3\frac{D_{i,k,j}}{\alpha }}}$

4

For line-based predictors in the LUR model, the predictor value X _ij in Equation 3 was calculated as the total length of line features within each rate area (Eq. 5). Because each rate area had the same area (A _i = 0.785 km²), this total length represented line-feature density for predictor j at rate area i. In Equation 5, L _ij is the line-feature density for predictor j at rate area i, and l _k ,_i is the length of individual line segment k of predictor j. For polygon-derived predictors used in the LUR model, such as population density, the predictor value X _ij was calculated as n _i /A _i .

$L_{i,j}={\displaystyle \sum }_k{l}_{k,j}\cap A_i$

5

Further details on the mathematics and methodology behind the grid construction can be found in Fox et al. ²⁸ To initiate Layer 2 of the Triple Stack framework, the independent grid-level built-environment predictor values (X _ij ) were entered into the Poisson LUR model shown in Equation 3. The full-period incidence rate (IR _i ) was used as the dependent variable. The regression was conducted in Stata. ⁵⁴ Backward selection and likelihood ratio tests were used to identify the most parsimonious model. We began with a saturated model and removed variables that did not improve model fit or act as confounders. The final Poisson LUR model was then used to generate the spatial mean trend surface, which formed one component of the subsequent geostatistical modeling framework.

2.5. Geostatistical mapping

In geostatistical terminology, the global component represents information estimated from the full study domain, while the local component represents observed data near a specific prediction point. To separate broad mean structure from local space-time dependence, we removed an additive global offset composed of spatial and temporal mean-trend components. The residual covariance model then quantified the remaining local space-time dependence after these broad patterns were removed.^{25,27,34,36,38} BME is “maximum entropy” because the prior PDF is constructed as the widest, least-assumptive distribution consistent with the available global information, including the mean trend and covariance structure. It is “Bayesian” because this prior maximum entropy PDF is updated using nearby observed data.^27,34,38

Thus, the final estimate reflects both the global structure of the process and local observed evidence. After estimation, the spatial and temporal mean trends were added back to the smoothed residual estimates to produce the final rate surface.^{24,27,34,39,40} For the LUR model, incidence rates were calculated across the full study period. For BME/Poisson kriging, incidence rates were recalculated using rolling 9-month time windows so each observation was indexed by both grid location i and time t. The time-indexed incidence rate was defined as IR_i, _t = Y_i, _t /(n_i, _t T), where Y_i, _t is the number of pedestrian fatalities within rate area i during time window of duration T=9 months, and n_i, _t is the corresponding population denominator at time t.

In standard BME applications, both the spatial and temporal mean trends are estimated using exponential kernel-smoothing functions available in the BME library.^{24–27,29,35,36,38–40} However, the BME framework allows mean-trend components to be estimated using models other than kernel smoothing.^28,38,39 In this study, the spatial component of the global offset, go _s (i), was estimated using the Poisson LUR model and corresponds to LUR_i (Eq. 3). This enriched the BME knowledge base by incorporating built-environment predictors directly into the spatial component of the global offset. The temporal component, go _t (t), was estimated separately across the rolling 9-month windows using exponential kernel smoothing. Together, these components defined the global offset, go(i,t) = go _s (i) + go _t (t). The spatial component, go _s (i), varies across location but not time because it represents the full-period spatial pattern. The temporal component, go _t (t), varies across time but not location because it represents the study-wide temporal pattern.^{24,27,29,39,40} Therefore, for BME modeling, the incidence rate was re-estimated as IR_i, _t so each observation was indexed by both location i and time t. The resulting residual field used for BME modeling was calculated as R(i,t) = IR_i, _t − go(i,t).

Covariance modeling was performed on the residual field, R(i,t) to quantify the degree to which pedestrian fatality rates were related across space and time. In this context, the covariance function identifies which nearby observations are informative for estimating risk at a prediction location. Stronger covariance indicates greater similarity across nearby rate areas or adjacent time windows, giving those observations greater weight in the estimate. As covariance decays, that influence weakens.

The residual covariance structure is then used in the Layer 3 BME/Poisson kriging estimation to estimate the residual risk R̂(i,t) at unsampled space-time locations. Poisson kriging was selected because the outcome was a population-normalized incidence rate based on sparse pedestrian fatality counts. When the population denominator (n_i, _t ) is small, a single fatality can produce an unstable or artificially high incidence rate. Unlike simple kriging, Poisson kriging incorporates the population denominator into the estimation process as a measure of rate uncertainty. Rates based on smaller denominators are treated as less reliable and receive less weight in the local estimate, while rates based on larger denominators are treated as more stable. This allows the model to rely more heavily on surrounding space-time observations when estimating risk in areas with sparse counts or small population support.^{6,28,30,32,33,37,41,42}

Before final maps were created and cross-validation was conducted, the global offset was added back to the smoothed residual estimates to generate the final incidence-rate estimates, IR̂(i,t) = R̂(i,t) + go(i,t). To compare the implemented methods and identify the most accurate approach, a mean squared error (MSE) based cross-validation analysis was conducted. This leave-one-out cross-validation method was applied at the grid-point level (i). Each observed incidence rate was sequentially removed and re-estimated using the remaining data. Estimation errors were computed as the squared differences between observed and predicted incidence rates. These squared errors were then averaged across all grid points and observation times to obtain the MSE, as shown in Equation 6. m is the method of choice (e.g. Triple Stack (TS) vs Poisson kriging (PK)), $\widehat{I{R}_{i,t}}$ is the estimated incidence rate via the model at location i and time t, the summation is across all i and t, and n is the number of space/time values ${IR}_{i,t}$ . The method with the lowest MSE was considered the best-performing predictor.

Cross-validation was restricted to observed locations and time periods available under the traditional administrative-boundary PK framework. This ensures the baseline PK and Triple Stack models were evaluated against the same validation set rather than against interpolated prediction-only points. In addition to MSE, root mean squared error (RMSE) and mean absolute error (MAE) were calculated from the same leave-one-out prediction errors to provide complementary measures of predictive accuracy. RMSE and MAE were rescaled to deaths per 10,000 person-years for interpretability, whereas MSE was retained on the squared incidence-rate scale.

${\text{MSE}}_m=\frac{1}{n}{\displaystyle \sum }_{i,t\hspace{0.17em}}{\left({\widehat{I{R}_{i,t}}}^m-I{R}_{i,t}\right)}^2$

6

In disease mapping contexts where outcomes are rare, overall error metrics can be misleading because observed rates in low-population areas are often dominated by stochastic variability. Latent rate theory posits that as population size increases, observed incidence rates converge toward the underlying latent disease rate. Consistent with this framework, validation metrics were calculated separately across population-size groups. An effective method would be expected to show stable or improving predictive performance as population size increases. Operationally, validation support locations were grouped into population percentiles b, and MSE (Equation 7), RMSE, and MAE were calculated separately within each percentile .^27,33

${\text{MSE}}_m\left(b\right)=\frac{1}{n_b}{\displaystyle \sum }_{i,t\in b}{\left(\widehat{{I{R}_{i,t}}^m}-I{R}_{i,t}\right)}^2$

7

Relative performance across methods was further quantified using the percent change in mean squared error (PCMSE; Equation 8). PCMSE was calculated as the relative difference in MSE between the baseline Poisson kriging model and the Triple Stack model, which incorporated Layers 1 and 2. Negative PCMSE values indicate improved predictive performance. Increasingly negative values at higher population percentiles reflect enhanced estimation of the latent incidence rate.^27,33 Because PCMSE is a deterministic, descriptive measure derived from cross-validation errors, confidence intervals are not reported. Instead, robustness was assessed by examining the stability and monotonic behavior of PCMSE across increasing population percentiles. Leave-one-out cross-validation was preferred over K-fold cross-validation to maximize use of the sparse pedestrian fatality data. This approach allowed each observation to contribute to both model fitting and validation while preserving a consistent basis for method comparison.

${\text{PCMSE}\left(\mathrm{b}\right)}_{TS}=100*\left(\frac{{MSE}_{TS}\left(b\right)-{MSE}_{PK}\left(b\right)}{{MSE}_{PK}\left(b\right)}\right)$

8

3. Results

Of the 18 built environment variables assessed in the Poisson regression, nine were retained as statistically significant predictors of pedestrian fatalities in Cali (Table 1). Built-environment predictors were modeled using continuous, spatially weighted exposure metrics rather than simple counts. For point-based features, coefficients represent the change in pedestrian fatality risk associated with a one-unit increase in SEDC. SEDC is a dimensionless index that captures both feature proximity and density relative to each grid centroid. For line-based features, coefficients correspond to a one-unit increase in total feature length within the incidence area. Population effects reflect changes per unit increase in population density. Effect sizes should therefore be interpreted as relative changes in risk per unit increase in the modeled exposure metric rather than per individual infrastructure element.

Variables linked to increased risk included fire hydrants (6%), health centers (54%), educational centers (5%), libraries (1%), and roads (3%). These associations suggest pedestrian fatalities tend to cluster around these features, potentially due to increased exposure or inadequate pedestrian infrastructure. The particularly high effect size near health centers stands out, exceeding that of all other variables. The low standard errors for these estimates suggest the results are precise, though the larger SE for health centers (1.09) indicates greater variability in how fatalities cluster around these locations. Conversely, open areas/parks and population density were protective, each associated with an approximate 8% reduction in fatalities. These environments may promote safety by providing pedestrian-dedicated space, slowing vehicle speeds, and enhancing driver vigilance. The low standard error indicates a stable relationship across the study area. Sidewalk density was statistically significant (p = 0.0044) but had a negligible effect (β ≈ 1.00), indicating limited practical influence on fatality risk. The intercept (1.9e-5) reflects the baseline fatality rate, with a larger SE indicating expected variation in fatality distribution across the city.

The temporal mean trend, shown in the upper portion of Figure 2, compares the raw temporal trend with the smoothed mean trend. The raw trend fluctuates with short-term variation, while the smoothed mean trend captures the overall trajectory. The results show an initial increase in pedestrian fatalities. Fatalities peaked between 10 and 15 months into the study period, followed by a steady decline that stabilized toward the end of the observation window. The lower portion of Figure 2 depicts the spatial mean trend through two maps, displaying the distribution of pedestrian fatality rates across Cali. The TS method showed that pedestrian fatalities were highly concentrated rather than uniformly distributed across space.

Figure 2.

Raw and smoothed mean temporal and spatial trends used for mean-trend removal in the Triple Stack (TS) framework. The top panel shows the raw temporal trend (dashed line) and the smoothed temporal mean trend (solid line), estimated using a global exponential kernel. The bottom panels show the raw spatial pattern (left) and the smoothed spatial mean trend derived from the Poisson land use regression (LUR; right). Together, these panels illustrate the large-scale temporal and spatial structure removed prior to covariance modeling, with the LUR-smoothed surface capturing the dominant spatial concentration of pedestrian mortality in central Cali while reducing local noise.

Figure 3 presents the spatial and temporal covariance analysis of pedestrian mortality in Cali, Colombia, from 2008 to 2010. The top panel illustrates the spatial covariance, showing how mortality rates are correlated across increasing distances. The observed spatial experimental covariance decreased sharply within 0.25–0.5 km. This indicates strong spatial clustering of fatalities in close proximity, with correlation weakening beyond 0.75–1 km. In the bottom panel, the covariance decreased gradually as temporal lag increased. Autocorrelation persisted for approximately 12–24 months before losing 95% of its initial autocorrelation. Overall, these results suggest pedestrian mortality in Cali clusters spatially within small geographic areas and maintains temporal dependence over years. In both panels, the fitted covariance model closely followed the experimental covariance values, indicating good agreement between the observed covariance structure and the model used for BME/Poisson kriging estimation.

Figure 3.

Spatial and temporal covariance functions for pedestrian mortality in Cali, Colombia (2008–2010) using the Triple Stack (TS) framework. The top panel shows the spatial covariance as a function of spatial lag distance, and the bottom panel shows the temporal covariance as a function of temporal lag in months. In both panels, the experimental covariance is shown in red and the fitted model is shown in blue. The figure indicates strong short-range spatial dependence, with covariance declining rapidly within approximately 500–750 m, and longer temporal persistence, with dependence extending across roughly 12–24 months. These patterns support the interpretation that pedestrian fatality risk is highly localized in space but more persistent over time.

Figure 4 presents a comparative visualization of pedestrian fatality hotspots derived from PK versus the Triple Stack approach. To ensure a clean methodological comparison, both maps omit administrative boundaries and built-environment overlays. This allows emergent hotspot patterns to be evaluated without visual interference. In the standard Poisson kriging map on the left, risk is overgeneralized. Hotspots cluster around centroids, likely reflecting artifacts of data structure and spatial smoothing. In contrast, the TS model (right) reveals more spatially discrete and realistic clusters. These clusters highlight high-risk intersections aligned with built-environment features such as schools, health centers, and arterial corridors. This refinement likely stems from the TS model’s ability to remove background mortality effects caused by kriging islands,^28,31 which often exaggerate centrality-based risk.

Figure 4.

Comparison of pedestrian fatality hotspot patterns estimated using Poisson kriging only (left) and the Triple Stack (TS) approach (right) for a representative rolling 9-month window (July 2009–March 2010). Both panels show predicted pedestrian mortality rates in Cali for the same time period using a shared color scale. The Poisson kriging-only map produces broader, more diffuse hotspot patterns, whereas the Triple Stack approach yields more spatially discrete and localized hotspots. This comparison highlights the improved hotspot delineation and reduced background smoothing achieved by integrating grid-based aggregation and LUR-derived mean trends prior to Poisson kriging.

Figure 4 captures only a single time point. To assess comparative model performance over time, we also generated an animated time series covering the full study period. This animation demonstrated large differences in spatial and temporal fidelity between the models. The animation and analysis code are available at: https://mserre.sph.unc.edu/BMElab_web/mappingStudies/LF_PedMort_Cali_BMELUR/index.htm.

Table 2 presents the cross-validation results comparing the traditional PK baseline with the full Triple Stack framework. The traditional PK model represents the Layer 3-only approach, whereas the Triple Stack model includes the added grid-based spatial support and LUR-derived spatial mean trend before BME/Poisson kriging smoothing. Both models showed lower MSE values as population percentile increased, consistent with greater rate stability in areas with larger population denominators. However, the TS model produced lower MSE values than the PK model at every population percentile. The relative improvement was consistently large across the population range, with MSE reductions of approximately 84.9%–89.5%. To better contextualize these improvements, RMSE and MAE were expressed as deaths per 10,000 person-years.

Table 2.

Cross-validation comparison of traditional Poisson kriging and the Triple Stack model across population percentile strata. Traditional Poisson kriging represents the Layer 3-only baseline, whereas the Triple Stack model integrates grid-based spatial support, LUR-derived mean-trend structure, and Poisson kriging/BME smoothing across Layers 1–3. MSE is reported on the squared incidence-rate scale. RMSE and MAE are expressed as deaths per 10,000 person-years. Percent change in MSE (PCMSE) quantifies the relative reduction in prediction error for the Triple Stack model compared with traditional Poisson kriging, with negative values indicating improved predictive performance. The average row represents the arithmetic mean across population-percentile strata.

Population percentile	MSE PK	RMSE PK	MAE PK	MSE TS	RMSE TS	MAE TS	PCMSE
0	1.93E-07	4.40	1.28	2.68E-08	1.64	0.43	-86.12
10	4.73E-08	2.18	0.95	7.14E-09	0.84	0.32	-84.92
20	3.36E-08	1.83	0.85	4.71E-09	0.69	0.28	-85.99
30	3.06E-08	1.75	0.80	4.18E-09	0.65	0.27	-86.35
40	2.58E-08	1.61	0.76	3.88E-09	0.62	0.26	-84.99
50	2.30E-08	1.52	0.71	2.96E-09	0.54	0.23	-87.12
60	2.16E-08	1.47	0.68	2.63E-09	0.51	0.21	-87.84
70	2.18E-08	1.48	0.68	2.50E-09	0.50	0.21	-88.55
80	1.49E-08	1.22	0.61	1.77E-09	0.42	0.19	-88.16
90	7.56E-09	0.87	0.46	7.93E-10	0.28	0.14	-89.52
Average	4.20E-08	1.83	0.78	5.74E-09	0.67	0.25	-86.96

The largest practical implication was observed at the lowest population percentile, where the small-number problem is expected to be most pronounced. In this group, the traditional PK model had an RMSE of 4.40 deaths per 10,000, whereas the TS model reduced RMSE to 1.64 deaths per 10,000. MAE showed a similar reduction, decreasing from 1.28 deaths per 10,000 for PK to 0.43 deaths per 10,000 for TS. As population density increased, both models improved, but the TS model continued to show lower error across the full population range. At the highest population percentile, the TS model had an RMSE of 0.28 deaths per 10,000 compared with 0.87 deaths per 10,000 for the traditional PK model. MAE at the highest population percentile was also lower for TS, decreasing to 0.14 deaths per 10,000 compared with 0.46 deaths per 10,000 for PK. Across population-percentile strata, average RMSE decreased from 1.83 deaths per 10,000 for PK to 0.67 deaths per 10,000 for TS, while average MAE decreased from 0.78 to 0.25 deaths per 10,000.

Figure 5 depicts the spatial distribution of pedestrian mortality in Cali, Colombia, over three consecutive years (2008, 2009, and 2010) for the period between July and March. The mortality rate is represented using a color gradient, with darker shades indicating higher fatality rates. Across all three years, pedestrian fatalities exhibit clear spatial clustering, with persistent high-risk areas concentrated in the central and northern regions of the city. In 2008, mortality hotspots are primarily concentrated in central urban areas, with several smaller clusters emerging in peripheral regions. By 2009, these central hotspots remain, but additional clusters appear along neighborhood boundaries, particularly in the northern sector. In 2010, the distribution of fatalities shifts slightly, with some previously dense clusters dispersing while others intensify, particularly in the northernmost areas. The temporal comparison highlights both stability and evolution in pedestrian mortality patterns. An animated time series illustrating the full evolution of pedestrian mortality throughout the study period can be found online at: https://mserre.sph.unc.edu/BMElab_web/mappingStudies/LF_PedMort_Cali_BMELUR/index.htm. The visualization provides a dynamic representation of how pedestrian fatality patterns shifted over time.

Figure 5.

Predicted pedestrian mortality rates in Cali, Colombia, for representative July–March rolling 9-month windows in 2008, 2009, and 2010. Each panel shows the spatial distribution of predicted mortality rates (deaths per 10,000 person-years) for the same seasonal window across three consecutive years, using a common color scale. The maps illustrate both persistence and temporal evolution in hotspot structure, with recurrent high-risk areas in central and northern Cali alongside year-to-year shifts in the location and intensity of localized clusters.

These differences underscore the ability of the TS approach to enhance hotspot detection beyond what can be inferred from visual overlays alone. To systematically assess the two approaches and aid interpretation, Table 3 compares them across key spatial modeling dimensions and includes illustrative public health use-case examples across epidemiological and planning contexts. Although the TS method entails greater setup and computational cost, it improves hotspot delineation, enhances alignment between predictors and outcomes, and mitigates bias associated with MAUP. These strengths make it better suited for high-resolution urban safety analysis (Table 3).

Table 3.

Comparison of single-layer areal aggregation approaches (e.g. Simple kriging) implemented in BMEGUI (an open-source, no-code platform) and the Triple Stack (TS) framework across key spatial modeling dimensions and illustrative use cases.

Dimension	BMEGUI single layer implementation	Best use example (kriging only)	Triple stack implementation	Best use example (TS)
Setup/Coding Complexity	Low; predefined zones	Routine surveillance	High; custom grid + LUR	Urban injury or vector-borne disease mapping
Computational Load	Low	Rapid assessments or rural health monitoring	High	Urban air pollution or tsetse fly risk
Cluster Accuracy	Moderate; coarse	Chronic disease trends	High; fine-grained	STI clusters or pedestrian risks
Temporal Consistency	Inconsistent	Seasonal outbreak summaries	Stable across time	Long-term tracking of injury or HIV
Emerging Risk Detection	Weak	Flu or COVID-19 trend detection	Strong	Zika/mosquito or gun violence emergence
Predictor Alignment	Moderate	Basic environmental correlations	High; reflects BE factors	Built environment risks: schools, roads
MAUP Bias	Present	County cancer rates	Mitigated	Neighborhood maternal mortality or overdose spikes
Policy/Intervention Utility	Generalized	Mobile clinic routing	Precise interventions	Intersection upgrades or equity-based policies
Edge Detection	Weak in peripheries	Core city crash mapping	Sensitive to fringe clusters	Peripheral injury detection in underserved or informal settlements
Stakeholder Clarity	Intuitive zones	Community-focused reporting	Clear with context	Equity planning with technical input
Portability	Easy across cities	Nationwide chronic disease review	Adaptable with setup	City-specific road injury analysis
Spatiotemporal Modeling Constraints	Moderate; resolution-limited	Malaria or flu surveillance	Minimal; grid enhances detection	Time-resolved outbreaks (STIs, mosquito)

4. Discussion

This study demonstrates the methodological and practical value of the Triple Stack framework for pedestrian fatality modeling. By integrating grid-based aggregation, built-environment-informed mean-trend modeling, and population-adjusted Poisson kriging, the framework produced more discrete, localized areas of elevated risk than single-layer Poisson kriging. It also reduced broader smoothing artifacts and background noise. Although each layer has been applied independently in prior work, this study is, to our knowledge, the first to combine all three in this context. One important methodological insight is that the observed improvement may not arise solely from the final geostatistical smoothing step (PK). Instead, the first two layers, particularly grid-based rate construction and LUR-based mean-trend modeling, may stabilize sparse rates before Poisson kriging is applied. Beyond pedestrian safety, this framework offers a scalable template for spatially clustered health outcomes where rare events, unstable denominators, population structure, and local environmental context shape apparent risk patterns.

The cross-validation results support this broader methodological relevance by showing improved prediction across the full range of population densities. Across population percentile strata, average MSE decreased from 4.20E−08 with traditional Poisson kriging to 5.74E−09 with the Triple Stack model, corresponding to an average percent change in MSE of −86.96%. RMSE decreased from 1.83 to 0.67 deaths per 10,000 person-years, while MAE decreased from 0.78 to 0.25 deaths per 10,000 person-years. Because MSE and RMSE are more sensitive to larger errors, these improvements suggest the Triple Stack framework reduced both average prediction error and larger localized discrepancies. In the lowest population percentile, where sparse denominators would be expected to produce the greatest instability, MSE decreased from 1.93E−07 to 2.68E−08, RMSE decreased from 4.40 to 1.64 deaths per 10,000 person-years, and MAE decreased from 1.28 to 0.43 deaths per 10,000 person-years.

Given that observed pedestrian mortality rates rarely exceed 15 deaths per 10,000, reductions of this magnitude are important for identifying high-risk urban areas. This is especially true in low-population areas, where small-number instability can create or exaggerate apparent hotspots. In Figure 3, hotspot regions become visually distinguishable at approximately 4–6 deaths per 10,000; therefore, a PK RMSE of 4.40 deaths per 10,000 is large enough to materially affect hotspot interpretation. However, because leave-one-out cross-validation was conducted at observed support locations and time periods shared by both models, predictive performance may remain somewhat optimistic under strong spatial autocorrelation. These results should therefore be interpreted as descriptive, comparative measures of model performance rather than formal tests of statistically significant improvement or absolute out-of-sample prediction accuracy.

An additional methodological insight is that although the grid-based framework was originally designed to better represent the spatial distribution of health outcomes, it may also improve rate estimation by implicitly redistributing population denominators across adjacent areal units. This unintended benefit may yield estimates that better reflect local population structure, particularly in transitional areas where administrative boundaries impose abrupt denominator changes. ²⁸ This is likely most pronounced near administrative edges and in transition zones between high-density urban or residential areas and lower-density commercial, industrial, or rural areas. By redistributing both cases and population denominators across overlapping local supports, the moving-window approach may reduce overestimation in low-population areas adjacent to denser neighborhoods and reduce underestimation or spatial displacement in nearby urban areas. Similar improvements have been reported in prior work, including syphilis incidence in Forsyth County, NC. ³² Formal evaluation of this population-redistribution mechanism should be addressed in future methodological work.

Beyond improvements in rate estimation, the spatiotemporal covariance analysis clarifies the scale and persistence of pedestrian fatality risk. Spatial correlation declined sharply within approximately 500–750m, indicating risk is highly localized within small urban catchments. In contrast, temporal correlation persisted over roughly 12–24 months, suggesting high-risk areas are not merely short-term fluctuations but may reflect stable structural conditions, such as inadequate infrastructure, consistently high pedestrian activity, or persistent traffic patterns. The short spatial range supports geographically targeted interventions, such as intersection redesign, crossing improvements, traffic calming, and sidewalk or visibility enhancements, rather than diffuse citywide treatment. The longer temporal dependence suggests these locations may warrant sustained monitoring and medium-term intervention planning.^{15,19,21,55–59} The covariance ranges also support the modeling choices described in the methods, including the 500 m circular rate areas, 650 m grid spacing, and 750 m SEDC influence distance, which align with the observed short-range spatial structure of pedestrian fatality risk. ³⁰

Although complete homogeneity and stationarity cannot be proven empirically for a complex urban injury process, several diagnostic steps were used to assess whether the residual structure was sufficiently stabilized for BME modeling. Following removal of the spatial mean trend via LUR and the temporal mean trend via kernel smoothing, residual behavior was evaluated through comparison of the spatial mean trend, temporal mean trend, covariance structure, and cross-validation performance before and after mean-trend removal. Additional sensitivity work from the earlier Cali analysis informed the choice of temporal and spatial aggregation levels by comparing exploratory patterns and autocorrelation behavior across multiple spatial and temporal resolutions. This supports the use of 9-month temporal windows and fine-scale spatial structures for the current application. ³⁰ Taken together, these checks do not establish perfect homogeneity or stationarity, but they provide empirical support the trend-removal strategy substantially reduced large-scale structure and improved suitability for subsequent covariance modeling.

Our analysis identified several built-environment features associated with increased pedestrian fatalities, highlighting persistent infrastructure challenges in urban settings. For example, fire hydrants are typically considered safety-related infrastructure because they can reduce on-street parking and improve visibility, yet poorly placed hydrants may obstruct sidewalks or create barriers that force pedestrians into the street.^60,61 This interpretation is especially relevant in Cali, where sidewalks are not ubiquitous or uniformly maintained, complicating efforts to directly link sidewalk presence or quality to pedestrian fatality risk.

Elevated risk near educational institutions is consistent with prior global studies and likely reflects areas where high pedestrian activity near schools intersects with insufficient crossing infrastructure and discontinuous pedestrian facilities.^59,62 In contrast, library proximity was associated with reduced risk, potentially due to differences in pedestrian behavior, such as lower alcohol consumption near libraries. ⁵⁹ The concentration of fatalities near health centers is particularly striking, as it exceeds the impact of all other built environment variables. It may reflect concentrated pedestrian activity among vulnerable road users, including older adults, people with disabilities, patients, and transit users, as well as increased vehicle movements from drop-offs, ambulances, parking access, and large parking areas with multiple entry points. However, the high standard error for health centers suggests that risk may vary substantially by location, indicating that some facilities may be much more hazardous than others.

Road characteristics further reinforce these risk patterns, as primary and secondary roads may increase pedestrian fatality risk through higher speeds, wider lanes, and limited pedestrian crossings. Conversely, the protective association with population density is consistent with prior research showing that denser, more connected neighborhoods often experience fewer crashes because lower speeds, narrower roads, smaller block sizes, and improved pedestrian infrastructure can reduce injury risk. For example, a 2019 study in the Philadelphia region found denser neighborhoods had fewer collisions, injuries, and fatalities, underscoring the protective effect of compact urban form observed in our own analysis.^63,64

From a public health perspective, these findings show how high-resolution spatial modeling can support more targeted pedestrian safety surveillance and planning.^21,59 In Cali, the Triple Stack results could help public health and transportation agencies prioritize smaller areas for field audits, infrastructure review, traffic-calming assessment, and ongoing monitoring.^{14,15,19,21,55,59} The built-environment findings can further guide what to examine within those locations, including pedestrian access near health centers and schools, crossing availability, sidewalk continuity, visibility, pedestrian volume, traffic speed, enforcement, maintenance conditions, and road design along primary and secondary corridors.^{14–17,19–21,55,56,58,59,62–65}

These findings should be interpreted as predictive spatial associations rather than causal effects. Features such as health centers, educational institutions, primary and secondary roads, sidewalks, fire hydrants, parks, libraries, and other services may identify locations where pedestrian activity, roadway complexity, infrastructure quality, or historical risk patterns coincide. However, the analysis cannot determine whether any single feature directly increases or reduces fatality risk. Some associations may also reflect spatial confounding or reverse causality, whereby infrastructure and services are concentrated in areas that were already high risk. This limitation reflects both the observational design and the available covariate data, which captured the presence and spatial configuration of infrastructure but not detailed measures of condition, maintenance, pedestrian volume, crossing quality, enforcement, or traffic speed. Accordingly, the results should be used to identify places and contextual features warranting closer investigation, not as direct evidence that modifying any single feature would necessarily reduce pedestrian fatalities.

One limitation is that the analysis was based on pedestrian fatality data from a single city, Cali, Colombia, during 2008–2010. This may limit direct generalizability to other urban contexts or time periods. Nevertheless, rigorously analyzed historical data remain valuable in data-constrained settings. Pedestrian fatality risk is strongly shaped by structural built-environment characteristics, including roadway design, traffic speeds, land-use patterns, pedestrian infrastructure, road hierarchy, intersection layout, land-use configuration, and access to major destinations. Many of these features change slowly over time, meaning historical spatial patterns can remain informative for understanding persistent urban safety conditions.^44–46 Recent evidence also suggests infrastructure inequalities may persist alongside urban development rather than resolving quickly over time. This supports the use of historical built-environment data to identify enduring spatial risk patterns, even when not all local conditions remain static. ⁶⁶ Although recent studies have examined emerging behavioral risk factors, such as pedestrian mobile phone distraction, available evidence suggests these factors play a comparatively minor role in fatality risk. Enduring structural features appear to play a larger role in shaping pedestrian fatality risk. ⁶⁷

Our analysis is also constrained by the availability and resolution of covariate data. Built-environment variables were limited to the presence and spatial configuration of infrastructure rather than its quality, condition, or maintenance. Detailed information on sidewalk integrity, crosswalk design, signage, enforcement intensity, or traffic calming measures typically requires field-based audits or high-frequency administrative data, which were not available at the citywide scale for the study period. These attributes may also vary rapidly over time, making them more appropriate for targeted, site-specific investigations rather than retrospective spatial modeling at the metropolitan scale.

The scope of the current study was further restricted to pedestrian fatalities. While extending this framework to other vulnerable road users, such as cyclists and motorcyclists, would provide a more comprehensive assessment of urban traffic safety, these groups should not simply be included within a single pooled model. Different transport modes involve distinct exposure mechanisms, infrastructure interactions, and risk profiles, and therefore warrant dedicated analyses. This distinction is also consistent with broader Vision Zero and vulnerable road user research, which emphasizes pedestrian and bicyclist risks are shaped by overlapping but not identical built-environment, roadway, and socioeconomic factors,^56,58 and adoption of safety measures can vary substantially across urban settings and user groups. ⁶⁸ Future studies should therefore leverage longitudinal and multimodal data to evaluate how targeted infrastructure improvements, safety policies, and urban planning interventions alter risk over time as transportation systems and travel behaviors evolve.

Although the Triple Stack framework improves spatial resolution and rate stability, it requires geocoded events, population denominators, GIS-derived covariates, and greater technical capacity than simpler boundary-based approaches. BMEGUI partially reduces this accessibility barrier by supporting basic BME analyses without advanced programming, although it does not implement the full Triple Stack workflow. As summarized in Table 3, BMEGUI-capable single-layer approaches are most appropriate for routine surveillance, rapid assessments, broad chronic disease trends, and settings where intuitive administrative zones are useful for stakeholder communication. In contrast, the added complexity of the Triple Stack approach is more warranted when the research or policy question depends on fine-scale spatial support, locally calibrated covariates, and localized risk detection beyond what administrative boundaries can provide. Relevant applications may include pedestrian injury prevention, STI clustering, vector-borne disease emergence, overdose spikes, or equity-focused infrastructure planning.

Although the full Triple Stack framework is computationally more demanding than administrative-boundary-based Poisson kriging, its first layer remains conceptually straightforward. This is especially important within a BME workflow, where computational demands can be substantial. Grid-based rate construction requires four main steps: defining regular grid centroids, creating circular local rate areas, aggregating events and population within each area, and linking GIS-derived covariates to those standardized supports. Even when the full Triple Stack approach is not feasible, this first layer may still be transferable.

The broader value of the Triple Stack lies in the structure of the approach rather than in any city-specific parameterization. Recent work in higher-income settings has similarly shown pedestrian and non-motorized road user risk is shaped by income, race, land use, and roadway conditions, indicating these relationships are not unique to Cali or Latin America.^56,58,68 Cali-specific research on neighborhood environmental resources and health disparities further shows urban conditions are socially patterned, while evidence from Santiago, Chile, demonstrates pedestrian injury severity varies with both socioeconomic and built-environment conditions.^43,65 Together, these findings reinforce the need for locally calibrated models rather than assuming that built-environment associations transfer uniformly across cities.

Accordingly, the framework should be viewed as adaptable rather than fixed. Its implementation in other cities would require local calibration of spatial parameters, covariate selection, smoothing choices, aggregation scale, and model inputs rather than direct reuse of the Cali specification, particularly given cross-city differences in infrastructure inequality, urban development, and data availability. ⁶⁶ More broadly, the Triple Stack framework illustrates BME mean-trend components can be adapted to the research question rather than fixed to a single default specification. In this study, the spatial component was enriched using LUR because built-environment context was central to pedestrian fatality risk. In other applications, the temporal component could be enriched to capture seasonal, climatic, intervention-related, or other time-varying processes.^29,39,40

Supplemental material

Supplemental material - Triple-stack geostatistical modeling for urban injury: Integrating grids, built environment, and poisson kriging for pedestrian fatalities in Cali, Colombia

Supplemental material by Triple-stack geostatistical modeling for urban injury: Integrating grids, built environment, and poisson kriging for pedestrian fatalities in Cali, Colombia by Lani Fox, Marc Serre, Daniel A. Rodríguez, Shrikant I. Bangdiwala and Andrés Villaveces in Health Informatics Journal.