Abstract

Predictive accuracy is informative for infrastructure planning only when the target represents the quantity decision-makers intend to assess. We audit this requirement using an open inventory of 47,682 logistics sites in England and Wales whose heavy-goods-vehicle loading-bay counts were modelled from floor area and facility type. A reconstructed deterministic rule matches 47,678 records (99.992%). We reproduce an initial five-model prediction pipeline and examine feature ablation, spatially blocked validation, attribution, and residual-based flags. The reproduced pipeline achieves R² = 0.9871 and up to 98.7% of predictions within one bay. The decision tree using only floor area and category achieves R² = 0.99985. Spatial-regional features alone perform poorly under blocked validation, and four engineered features directly or indirectly incorporate the target. In the deduplicated audit configuration, features outside the reconstructed rule receive up to 33.7% of mean absolute SHAP attribution. This measures departure from the rule, not necessarily an error by the explainer. Out of the 685 baseline flags, 684 match the reconstructed rule exactly, while cross-model flag agreement is low (mean Jaccard 0.25). For simulated one-bay reductions at eligible sites, Random Forest recall is 60.0% with the full feature set and 98.8% with rule inputs alone. These perturbations assess detection of constructed deviations, not actual operational inadequacy. We propose a Target Provenance and Validity Protocol that separates target validity, predictive performance, and decision validity before model outputs are used in planning.

Keywords

Loading Bays, Freight Infrastructure Planning, Machine-Learning Benchmarking, Explainable AI, Target Validity, Target Leakage, Data Provenance, Spatial Validation

Share and Cite:

1. Introduction

Loading bays connect freight vehicles with the facilities they serve. Inadequate loading arrangements have a tendency to contribute to vehicle waiting. Also, constrained kerbside access can increase parking-search time, double parking, and disruption to other road users [1]-[4]. Clearly, studies involving kerbside loading show how access to loading space affects delivery performance. However, determining how many loading bays a warehouse or other logistics facility requires also depends on its operational demand regarding vehicle arrivals, loading times, and operating hours. In the United Kingdom, the Future of Freight strategy identifies freight as a strategic policy concern [5]. For facility-level assessment, credible benchmarks must account for site characteristics and loading demand.

Predictive modelling offers one possible source of comparative evidence. Open geocoded inventories provide information on logistics-site locations, sizes, and categories [6]. From a computational angle, studies in transport demonstrate the potential of tree ensembles for structured prediction tasks [7]. A model trained on independently observed provision has capabilities to estimate a conditional peer norm for a facility of a given size and type. Sites below that norm need to be selected for further investigation. Such a norm would describe prevailing provision, as it would not by itself establish adequate capacity. Post hoc methods like SHAP can explain how a fitted model uses its inputs [8] [9], supporting the scrutiny of operational decision tools [10] [11]. Critically, neither prediction nor explanation, however, validates the underlying measurement.

This study examines that measurement precondition. If recorded bay counts were produced by computational means rather than observed at each site, a model trained on those counts will mainly approximate the computational method. Its explanations describe the fitted surrogate, and its residuals compare two model outputs. Another problem that can be triggered is when the target signal enters the features or model-selection process, or when validation does not match the intended deployment setting [12]-[14]. The major concerns here are clear. They include how removing feature leakage or improving validation does not convert a generated target into an observed outcome. High predictive accuracy and plausible explanations can consequently coexist with an invalid interpretation of what the target measures.

We investigate this problem using an open inventory of logistics sites in England and Wales [6] [15]. The documentation of this inventory clearly identifies the bay field as modelled heavy-goods-vehicle (HGV) loading-bay provision. Floor area and category were used to extrapolate from satellite-image counts for a small sample. The field is not an independently observed count for every site. This transparency makes the inventory useful for a methodological audit. The documented construction, supplemented by reconstruction of the recorded values, supplies a reference for target generation. It does not supply ground truth for actual provision or operational adequacy.

The first author initially implemented a reference prediction-and-benchmarking pipeline, named B0. The pipeline compares Poisson and negative binomial regressions with Random Forest, XGBoost, and LightGBM. It adds spatial and target-derived features, uses SHAP, and then flags positive prediction residuals as candidate under-provision. The audit pipeline was implemented to reproduce the feature specification and tree-model results in the initial pipeline B0. Also, the features, validation, and model configuration were varied based on the objectives of this study. The reproduced pipeline reaches R² = 0.9871, up to 98.7% of predictions within one bay, and 685 residual flags. These figures motivate the main question: does benchmarking establish evidence about actual loading-bay provision when the recorded target was constructed from the available predictors? Based on the question and other claims from the introduction, the aim of this study is to distinguish rule approximation from substantive evidence and to formulate checks that should precede the planning use of such benchmarks.

To achieve this aim, the study pursues three objectives:

1) To establish the construction of the recorded bay field and assess how predictive performance changes across feature families and validation designs.

2) To assess how attribution and prediction sensitivity vary for different model configurations, distinguishing fidelity to the fitted model from fidelity to the reconstructed target rule.

3) To assess rule compliance and the stability of flagged sites, and to measure the detection of simulated downward deviations from the rule at eligible sites.

The study makes four important contributions. First, it reconstructs a deterministic approximation to the target-generating rule that matches all but four of 47,682 records. Second, it contrasts rule-input models with spatial and target-derived feature sets, separating predictive fit from the validity of the quantity being predicted. Third, it evaluates feature attribution against a known construction mechanism on real site covariates, while recognizing that attribution to a feature outside that mechanism can faithfully describe the fitted model. Fourth, it examines residual-flag stability and detection of controlled perturbations. The entire analysis carried out in the study and the findings motivate the proposition of the eight-step Target Provenance and Validity Protocol (TPVP).

2. Background and Related Work

2.1. Loading Infrastructure as a Freight-Planning Problem

Recent urban-freight research shows that loading-space allocation and parking time influence delivery operations [16] [17]. Research involving driver behaviour identifies trade-offs between actions like walking distance, parking-search time, cost, and enforcement considerations [4] [18]. Ghizzawi et al. [3] synthesise the literature around parking location, duration, inappropriate parking, and cruising. Modelling studies examine the implications of kerbside access and delivery strategies for traffic performance [2] [19], while facility-level evidence identifies operational and shipper characteristics that are associated with driver detention [1]. Clearly, provision is more than a count. It involves location, management, allocation, and competition for space [16] [20] [21].

Considering an individual logistics facility, adequacy concerns whether loading capacity meets operational demand. Floor area may be a useful proxy, but throughput, vehicle mix, operating hours, fulfillment arrangements, and network role influence the requirement for loading positions [22]. Reviews of freight planning and data infrastructure highlight the importance of better information for such assessments [23]-[25]. Open inventories can help address these information needs, provided analysts preserve the distinction between administrative attributes, estimated quantities, and independently verified infrastructure.

2.2. Statistical and Machine-Learning Approaches to Infrastructure Counts

Loading-bay counts are non-negative integers. Based on their empirical distribution, candidate models include Poisson regression, negative binomial regression, and zero-inflated or hurdle specifications [26] [27]. Transport studies illustrate the relevance of spatial dependence, measurement periods, and excess zeros to count modelling [28] [29]. These considerations concern model specification but require appropriate verification of what the count represents. A threshold-generated zero count does not establish a distinct behavioural process at logistics sites. Similarly, a large marginal variance-to-mean ratio does not establish conditional overdispersion.

Based on computational capabilities, machine learning models like Random Forests [30], XGBoost [31], and LightGBM [32] can accommodate nonlinear relationships and interactions without specifying a count-distribution likelihood. Their comparative performance depends on the data, tuning, and evaluation design. Methodological benchmarking also requires transparent design choices and sensitivity analysis [33]. Geographic features like density, DBSCAN clustering [34], and Local Moran statistics [35] must be assessed according to their information source and targeted prediction setting.

2.3. Leakage, Validation, and the Reproducibility of ML-Based Findings

Fundamentally, leakage occurs when a model uses information unavailable in the targeted prediction setting or when held-out outcomes influence model development [12] [14]. Direct outcome-derived predictors and test-set model selection are clear examples. Spatial dependence presents a related validation challenge. Random splits may overstate performance for deployment in geographically distinct areas, making spatially blocked evaluation useful when that deployment is the objective [13] [36] [37]. Unsupervised preprocessing on a full inventory is a separate issue. It can be a better option for a fixed-inventory, transductive task but does not establish performance for previously unseen locations. A generated target introduces another problem. Prediction may recover the generating model rather than the intended real-world construct. Work on generated variables and imputation addresses uncertainty and inferential compatibility [38] [39], while data-centric research emphasizes provenance, labels, and data quality [40]-[42]. Shortcut learning offers a complementary account of high accuracy that comes from an unintended input-label relationship [36].

2.4. Explainable AI and the Problem of Faithfulness

Explainable AI supports scrutiny of models used in operational decisions [10] [11]. SHAP assigns feature contributions using a Shapley-value formulation, and TreeSHAP provides efficient computation for tree models [8] [9]. Notably, interpretation depends on the value function and treatment of feature dependence [43]-[45]. Correlated predictors can redistribute attribution, and unconditional permutation can produce unrealistic input combinations [46] [47]. Works from Slack et al. [48] and Adebayo et al. [49] demonstrate manipulation of explanations and the need for sanity checks. Rudin [50] advocates interpretable models for high-stakes decisions. From all claims, three properties should be separated. These include an explanation’s faithfulness to the fitted model, the model’s fidelity to the target-generating mechanism, and the target’s validity for the domain construct. An explanation can satisfy the first while the other two are unsupported.

A documented deterministic target offers a useful testbed for this distinction. Its construction can be known even when actual real-world outcomes are not. These kinds of data combine realistic covariate distributions and geographic structure with an auditable label-generating rule. They can reveal which fitted dependencies extend beyond that rule, but they cannot establish the causal determinants of actual infrastructure or validate a planning standard.

2.5. Benchmarking as Decision Support

Benchmarking compares a unit with a defined reference and analyzes any divergence from the reference [51] [52]. A model-based peer benchmark estimates provision conditional on selected characteristics. Then, residuals identify units below or above that conditional norm. The interpretation depends on both the provenance of the outcome and the relevance of the comparison. For facility planning, even a benchmark trained on observed provision describes prevailing practice unless linked to independent evidence about operational adequacy. Data-driven loading-zone design provides related decision-support applications, but a network-design model should not be equated with a validated facility-level adequacy standard [16].

Two checks are very important in this discourse. First, analysts must distinguish deviations from observed peer provision, deviations from a constructed engineering reference, and evidence of inadequate loading capacity. These are different decision quantities. Second, flag lists depend on features, model family, tuning, rounding, and thresholds. Stability and sensitivity analyses are therefore needed before residuals are used to prioritise investigation. Controlled perturbations can assess whether a pipeline detects known departures from its reference [33].

2.6. Research Gap

The importance of data provenance, leakage control, appropriate validation, and careful interpretation of feature attribution is critically established in the literature. The unresolved question is how these checks interact when a planning benchmark is trained on a documented, nearly deterministic generated target. We examine three connected issues. The first issue is whether predictive success reflects the recovery of the target rule. The second is whether attribution and prediction sensitivity track that rule. The last is whether residual flags are stable and detect constructed downward deviations. The contribution lies in evaluating these issues together on a logistics-site inventory and translating the findings into a practical audit protocol.

3. Data and Target Provenance

3.1. The Logistics-Site Inventory

The data selected for this study are openly licensed logistics-site inventory compiled by de Saxe et al. [15]. Administrative property records were processed to obtain site locations, floor areas, and categories, while HGV loading-bay counts were subsequently modelled. The analysis file contains category, latitude, longitude, floor area in square metres (sqm), and recorded modelled bay count (bays). The audit code coerces numeric fields, excludes missing critical values, applies bounds of 100 - 500,000 m² and 0 - 200 bays, and restricts coordinates to 49˚ - 61˚N and −8˚ - 2˚ longitude. Only one of 47,683 records is excluded because latitude cannot be parsed, leaving 47,682 sites: 27,690 warehouses, 13,551 factories, and 6441 retail premises. The bounds are screening rules, not validated physical limits. Table 1 summarises the audit, and Figure 1 presents the geographic and marginal distributions.

Table 1. Data audit of the logistics-site inventory.

Figure 1. Data overview, regenerated from the source CSV. (a) Locations of 47,682 analyzed sites by category. (b) Distribution of recorded modeled HGV loading-bay counts on a logarithmic frequency scale; 43.1% record zero, and warehouses have a cap at 120. (c) Floor-area distributions by category on a logarithmic area scale; the dashed line marks 2000 m².

Four properties are necessary for interpretation. First, the minimum recorded floor area is 1000 m², so coverage does not represent all premises or all logistics activity. Second, coordinates rounded to six decimal places are shared by 64.49% of sites. The largest coordinate group contains 41 records. Co-location may reflect real shared sites or geocoding resolution, so these records should not automatically be treated as independent fine-scale locations. Third, the maximum latitude is 55.79˚N. The legacy feature labeled “scotland” for latitude ≥ 55˚N therefore denotes a northern latitude band within this England-and-Wales inventory. Fourth, the proportion at or below 2000 m² and the proportion recording zero bays both round to 43.06%. This distributional correspondence is a diagnostic clue.

3.2. Reconstructing the Data-Generating Process of the Target

de Saxe et al. [6] describe the estimation of HGV loading-bay counts from satellite imagery for 36 warehouses, 39 factories, and 22 retail premises. The documentation reports category-specific modelling, threshold assumptions, and upper limits. We reconstruct the operational mapping from the recorded values, distinguishing this empirical reconstruction from the documented model description. For warehouses, the floor-area-per-bay parameter k is searched over 700 - 1000 m² in 0.25 m² increments. For factories, it is searched over 4000 - 9000 m² in 1 m² increments. For each category, candidate caps range from three below the recorded maximum to the maximum itself. Parameters maximise exact agreement on the complete analysed sample. The reconstructed mapping is:

bays_i = 0, if sqm_i ≤ 2000;

bays_i = 1, if retail and sqm_i > 2000;

bays_i = min{cap_c, 1 + round[(sqm_i − 2000)/k_c)}, otherwise. (1)

Equation (1) is a threshold-anchored, rounded, and capped mapping. The recovered parameters k = 823.5 m² per additional warehouse bay, with cap 120, and k = 6532 m² per additional factory bay are presented in Table 2. They are shown with the scoring grid value at 6533 and cap 15. The equation reproduces 47,678 of 47,682 records (99.992%). The four discrepancies are two warehouses at exactly 2000 m² recording one bay, a warehouse at 32,058.49 m² recording 37 instead of 38 bays, and a factory at 96,742.68 m² recording 16 instead of the reconstructed cap of 15. These discrepancies reflect boundary, parameter, or implementation conventions. Figure 2 overlays the reconstructed mapping on the recorded target.

Table 2. Recovered imputation rule (Equation (1)) and exact-match rates.

Note: Parameters maximize full-sample exact agreement over the stated grids and cap candidates. The factory interval denotes two equally scoring grid values, not a statistical confidence interval. NumPy rounding uses ties to the nearest even integer. The fitted rule is a retrospective reconstruction reference, not a model selected independently of held-out labels. The separate floor-area-and-category decision tree achieves R² = 0.99985 and 99.80% exact accuracy on a random 20% hold-out (unpruned scikit-learn DecisionTreeRegressor with default settings and random_state = 42; hold-out split seed 42; full specification in Section 3.2).

The reconstruction establishes that almost all variation in the recorded target is recoverable from floor area and category. A separate decision tree trained on those two original attributes achieves R² = 0.99985, exact-integer accuracy of 99.80%, and within-one-bay accuracy of 99.98% on a random 20% hold-out. The tree is a scikit-learn DecisionTreeRegressor with default settings and random_state = 42: squared-error criterion, best-split search, no maximum depth, a minimum of two samples to split a node, a minimum of one sample per leaf, all features considered at each split, and no cost-complexity pruning. Its two inputs are raw floor area in m² and category as an integer code (factory = 0, retail = 1, warehouse = 2), without scaling. The hold-out is drawn with train_test_split (test_size = 0.2, random_state = 42), giving 38,145 training and 9537 test sites; this is the same partition as the baseline split. The fully grown tree has depth 9 and 131 leaves. Exact accuracy compares NumPy-rounded predictions with recorded counts, whereas within-one-bay accuracy uses unrounded predictions. The tree’s accuracy confirms the ease of approximating the generated target, not the accuracy of the inventory against field counts. The concentration of zeros and large marginal variance also reflect the threshold, category mix, and floor-area distribution.

Figure 2. Recorded loading bays against floor area, by category, with the recovered imputation rule (Equation (1)) overlaid. Slopes k are expressed in m² of floor area per additional bay; “match” is the percentage of sites whose recorded count equals the rule’s output exactly. The dotted line marks the 2000 m² threshold.

4. Audit Design

4.1. Overview

The audit pipeline is set up in four preparatory stages and three experiments (Figure 3). Stage 0 establishes provenance, while reconstruction of the recorded target is implemented in Stage 1. The baseline pipeline is reproduced in Stage 2. Stage 3 maps feature sources and roles. Predictive fit is examined in E1, while attribution and sensitivity relative to the rule are checked in E2. E3 examines residual flags and controlled perturbations.

4.2. The Baseline Specification (B0)

B0 is an author-developed pipeline, not an independently sourced workflow. The first author implemented it in the initial experimentation phase of this study, as a separate notebook (Loading Bay Prediction Pipeline) that ingests and cleans the inventory, engineers the spatial and target-derived features, fits the five models, computes TreeSHAP explanations, and flags under-provisioned sites from positive residuals. It follows a workflow commonly used to benchmark infrastructure counts with tree ensembles, spatial features, and SHAP, but it was not taken from any published study or deployed planning application. The audit therefore demonstrates a methodological risk that such a workflow carries, using the authors’ own initial pipeline as the case, rather than evaluating an existing application. The audit notebook reproduces the feature specification and main settings of that pipeline. Spatial and target-derived features are calculated on the full inventory before splitting. LightGBM uses the test partition for early stopping, and the reproduced baseline uses a random 80/20 split. These choices are retained deliberately so that the audit examines the pipeline as originally specified. Floor areas are used in their raw form and also as log1p(sqm), category is used as an ordinal code and two dummy indicators. Six unsupervised spatial features describe 5 km site density, DBSCAN membership and cluster size, and Gaussian kernel-density intensity. Five latitude bands are used for regional encodings. Four features are created for bay counts: Local Moran’s I, its quadrant, its permutation-significance indicator, and the regional median. Category-by-latitude-band site counts and mean floor areas complete the inputs. The original 24-column matrix repeats each of the two category-region aggregates, leaving 22 unique features. B0 retains those duplicate columns. F1 deduplicates them. B0 and F1 are consequently related but distinct configurations.

Figure 3. Revised audit framework. The preparatory stages establish documentary provenance, reconstruct the recorded rule, reproduce the author-built 24-column baseline (B0), and map feature sources. E1 - E3 examine predictive fit, rule-relative attribution and sensitivity, and residual-benchmark behaviour. Findings inform the Target Provenance and Validity Protocol.

The spatial and regional features are specified as follows. All distance-based features use haversine distances with an Earth radius of 6371 km, computed on the full inventory with a ball-tree index. density_5 km counts other sites within 5 km of each site, excluding the site itself, and log_density_5 km is its log1p transform. DBSCAN uses a neighbourhood radius ε = 3 km (3/6371 radians on the haversine metric) and a minimum of 10 points per core neighbourhood; it identifies 450 clusters and labels 6.78% of sites as noise. in_cluster indicates non-noise membership, cluster_size is the number of sites in the site’s cluster (zero for noise points), and log_cluster_size is its log1p transform. Kernel-density intensity uses SciPy’s Gaussian KDE on longitude and latitude in decimal degrees with a scalar bandwidth factor of 0.05 applied to the sample covariance. This gives kernel standard deviations of approximately 0.060˚ of longitude and 0.054˚ of latitude (about 4.1 km east-west and 6.0 km north-south at 52.5˚N); the feature is log1p (10⁶ × estimated density). Local Moran’s I uses an eight-nearest-neighbour spatial weights matrix built on longitude-latitude coordinates in decimal degrees and row-standardised. Its significance indicator uses 99 conditional permutations (seed 42) and p < 0.05, and the quadrant follows the PySAL esda coding (1 = high-high, 2 = low-high, 3 = low-low, 4 = high-low). Because 64.5% of sites share a coordinate, ties among co-located neighbours are resolved by the libpysal nearest-neighbour implementation rather than by an explicit rule. The five latitude bands, with inclusive lower bounds, are ≥ 55˚N (labelled “Scotland”; 387 sites), 53.5˚ - 55˚N (“north_england”; 9870), 52.5˚ - 53.5˚N (“midlands”; 13,359), 51.5˚ - 52.5˚N (“south_east”; 14,682), and < 51.5˚N (“south_west”; 9384). These labels are legacy names for latitude ranges, not administrative regions, and “midlands” is the omitted reference category in the dummy encoding. The regional median and the two category-band aggregates are computed over these bands on the full inventory.

Table 3. Features in B0, grouped by source and role relative to the reconstructed target rule.

Note: Moran quadrant against the site’s own bay count: quadrant 1 (n = 4395) and 4 (n = 5731) contain only sites with 3 - 120 bays; quadrants 2 (n = 7766) and 3 (n = 29,790) contain only sites with 0 - 2 bays.

Local Moran’s I embeds the site’s own label. For site i, it is proportional to zᵢ Σ_j w_ij z_j, where z is the standardized bay count. The quadrant indicator also uses the sign of the site’s centered outcome. In this dataset, it separates counts of 0 - 2 from counts of 3 - 120 exactly (Table 3). These inputs are inappropriate for predicting an unknown site count and can obscure anomalies when used for residual benchmarking. The full-inventory regional median also includes held-out outcomes, although its five-valued summary carries much less direct site information. The audit therefore distinguishes direct target encoding from aggregate target information instead of assuming all four features have equal effects.

The Poisson and negative binomial GLMs use an intercept, training-standardised features, and up to 300 iteratively reweighted least-squares iterations. The negative binomial dispersion parameter α is fixed at 1.0. Random Forest uses 300 trees, maximum depth 20, minimum leaf size 4, and max_features = √p. XGBoost uses 300 trees, depth 7, learning rate 0.05, row and column subsampling 0.8, γ = 0.1, α = 0.5, and λ = 1.0. LightGBM uses 300 trees, 63 leaves, learning rate 0.05, column subsampling 0.8, α = 0.3, and λ = 0.5. Its configured row-subsampling fraction is inactive with subsample_freq = 0. B0 early stopping has patience 30. The audit’s GLM implementation clips predictions to 0 - 200. Both the baseline pipeline and the audit standardise all columns before fitting the GLMs. The fixed-α negative binomial fit is nevertheless sensitive to the software environment: identical code gave R² = 0.7631 in the baseline run and 0.7505 in the audit run, and XGBoost differs slightly for the same reason (0.9831 versus 0.9856). The reproduction therefore preserves the main tree results and flag count but is not numerically identical for every model.

4.3. Experiment 1: Accuracy Decomposition (O1)

E1 is in six feature sets (Table 4), five model families, and three validation designs. F0 contains the rule inputs and their encodings. F1 contains the 22 unique baseline features. F2 excludes the three Moran features, and F3 excludes the regional median. F4 isolates spatial and regional inputs, while F5 isolates target-derived inputs. This design assesses family-level predictive contributions.

The designs are a random 80/20 hold-out, shuffled five-fold cross-validation, and five-fold GroupKFold validation using 420 occupied spatial blocks. Blocks are approximately 20 km × 20 km. Longitude is converted using a fixed 52.5˚ reference latitude, and latitude is calculated using a spherical Earth approximation. Entire blocks are assigned to folds. This prevents co-located records from crossing a fold boundary. E1 LightGBM early stopping uses a random 10% subset of the training fold without refitting on all training observations. Other model settings are fixed. Geographic density, clustering, and category-region summaries are full-inventory features. F0 and F3 exclude target-derived predictors, but F3 is still a transductive, fixed-inventory evaluation. F1, F2, and F5 are set up to retain held-out target information. The full-sample reconstructed rule is scored on each fold only as a reconstruction reference. Metrics are R², MAE, RMSE, exact-integer accuracy, within-one-bay accuracy, and zero-bay accuracy.

Table 4. Feature sets used in Experiment 1 (E1).

4.4. Experiment 2: Attribution and Prediction Sensitivity Relative to the Rule (O2)

The reference mechanism in Equation (1) uses floor area and category. These are designated rule-input families, not empirically identified causes of actual loading-bay provision. Category-region counts and mean floor areas are contextual aggregates that can proxy category information. E2 fits seven configurations to each of F1 and F3 on the baseline random 80/20 partition: Random Forest with max_features ∈ {√p, 0.5, 1.0}, and XGBoost and LightGBM with column subsampling ∈ {0.8, 1.0}. E2 uses 300 trees without early stopping. TreeSHAP is computed for the first 1000 test rows; permutation importance uses the first 3000, five repeats, and the estimator’s default R² score. SHAP values are computed with shap. TreeExplainer (SHAP 0.52.0), called on the fitted model alone with no background dataset and the default raw model output. Without background data, the explainer uses feature_perturbation = “tree_path_dependent”, the path-dependent TreeSHAP algorithm. The expectation over absent features is taken along each tree’s decision paths and weighted by the training-data cover stored in each node, so the implicit reference distribution is the 38,145-site training partition as recorded in the fitted trees. No separate interventional background sample is drawn, and the additivity check is disabled (check_additivity = False). The reported attribution shares are therefore path-dependent attributions. Interventional TreeSHAP with an explicit background sample could allocate attribution differently among correlated inputs, such as the several floor-area and category encodings. For continuity with the notebook, we retain the label spurious attribution share (SAS):

$\text{SAS}=\left[{\displaystyle {\sum }_{j\notin C}{I}_j}\right]/\left[{\displaystyle {\sum }_j{I}_j}\right]$ (2)

where C denotes the encoded rule-input features.

For SHAP, I_j is the mean absolute attribution for feature j. For permutation importance, negative mean score reductions are set to zero before normalising. SAS based on permutation is therefore a share of positive score reductions, not a directly equivalent quantity to SHAP SAS. Headline SAS includes category-region proxies among features outside the rule. Strict SAS retains the same total denominator but counts only spatial, regional, and target-derived attribution in its numerator. Proxy shares are reported separately. SAS measures attribution outside the reconstructed rule. It does not establish inaccurate explanation of the fitted model or causal irrelevance in the real world. Top-five precision counts encoded rule-input columns among the five highest-ranked SHAP features. Because several columns encode the same underlying attribute, it is a representation-dependent diagnostic.

We also evaluate prediction sensitivity by independently permuting each feature outside the rule-input families within test strata sharing the same reconstructed rule output. This keeps each permuted feature’s marginal distribution within an output stratum but breaks dependencies among features and with the fixed inputs. It can produce inconsistent combinations, such as a density value paired with an unrelated log-density. Consequently, the procedure is a diagnostic stress test. We report the mean absolute prediction change and the proportion whose rounded prediction changes. A predictor implementing Equation (1) exactly would be invariant, but sensitivity of a fitted model may reflect both off-rule dependence and extrapolation induced by the perturbation. Joint permutation of coherent feature groups or regeneration from valid inputs would be needed to separate these effects.

A supplementary rank-correlation diagnostic asks whether a fitted model relies more on inputs outside the rule near the rule’s rounding steps, where floor area alone is least decisive about the recorded count. It uses the F3 LightGBM configuration with column subsampling 0.8 and its 1000 SHAP-explained test rows. Eligible observations are the warehouse and factory rows above 2000 m² in that sample (n = 508: 326 warehouses and 182 factories). Retail rows and rows at or below the threshold are excluded because Equation (1) has no rounding step there. The first variable is off-rule attribution: the sum of absolute SHAP values over all features outside the floor-area and category families. The second is the distance to the nearest rounding step, |u − ⌊u⌋ − 0.5|, where u = (sqm − 2000)/k_c is floor area above the threshold in normalised bay units. The distance ranges from 0 at a rounding step to 0.5 midway between steps. Spearman’s ρ between the two variables is reported with its two-sided p-value. Four capped warehouse rows are retained, although the cap makes their step distance immaterial.

4.5. Experiment 3: Validity, Stability, and Power of Residual Benchmarking (O3)

E3a reproduces the baseline Random Forest flags using the original 24-column matrix. This includes duplicate aggregates. The model is trained on the baseline 80% partition and predicts all 47,682 sites. The resulting 685 flags combine training and test predictions. A site is flagged when its rounded predicted count exceeds its recorded count. We then assess how many flagged records match Equation (1). For a rule-compliant record, a positive residual represents a mismatch between the fitted surrogate and the reconstructed reference. It supplies no independent evidence of actual under-provision.

E3b evaluates flag stability using out-of-fold predictions from the fixed spatially blocked five-fold design. Three feature sets (F0, F1, and F3), three tree models, and five seeds (42, 7, 123, 2024, and 99) yield 45 flagged sets. Models use the E2 fitting functions without early stopping. Pairwise Jaccard similarity is $\left|A\cap B\right|/\left|A\cup B\right|$ , with two empty sets assigned a similarity of 1. Comparisons across models hold the feature set and seed fixed. Also, comparisons across feature sets hold the model and seed fixed. Seed comparisons hold the feature set, model, and fold assignment fixed. These summaries distinguish algorithmic randomness from specification sensitivity.

E3c tests the detection of simulated downward deviations from the reconstructed rule. For each decrement d ∈ {1, 2, 5}, 953 sites are sampled without replacement, corresponding to 2% of the full inventory, from rule-compliant records with at least d + 1 recorded bays. Thus, the eligibility thresholds are 2, 3, and 6 bays, respectively. The eligible populations contain 15,656, 10,124, and 4836 records for d = 1, 2, and 5, so the 953 perturbed records are 6.1%, 9.4%, and 19.7% of the respective eligible populations, while remaining 2.0% of the full inventory in every case. Because residual flags concentrate among higher-count sites (Section 5.4), the rising contamination rate within the eligible population can raise precision at larger decrements independently of detection ability. Precision is therefore not directly comparable across decrement sizes. The zero-bay and one-bay records are never perturbed. Subtracting d leaves a positive integer count, and no clipping is applied to the contaminated target. Three replicates use seeds 42, 7, and 123. Because eligible populations differ by d, results across decrement sizes are not a same-site dose-response comparison. For F1, Moran features and regional medians are recomputed from the entire contaminated target, including held-out labels. This was carried out to replicate the outcome-dependent baseline. Random Forest and LightGBM are trained on contaminated training-fold outcomes and evaluated out of fold against the known perturbation indicators. A flag requires the rounded prediction to exceed the contaminated count by at least one bay. Recall and precision quantify detection of injected deviations. They do not validate the rule as an adequacy standard or establish detection of real operational shortages.

4.6. Reporting and Reproducibility

The audit pipeline uses fixed seeds for partitions, models, and perturbation replicates. The reported analyses were run on Google Colab with Python 3.13.15, scikit-learn 1.6.1, XGBoost 3.4.1, LightGBM 4.6.0, and SHAP 0.52.0, as recorded in the results summary file that the notebook writes for each run. The GLMs use statsmodels, and the Local Moran statistics use libpysal and esda. E1 reports means across folds and exports fold-level metrics and R² standard deviations. The stability analysis holds fold assignment fixed across seeds. The summaries describe variation in the executed designs. GLM predictions are clipped to [0, 200], a numerical bound inherited from the cleaning specification rather than a verified physical maximum. All metrics use non-negative predictions; exact accuracy and flags use NumPy rounding. Within-one-bay accuracy uses the unrounded prediction. Negative binomial instability is reported as a limitation of this fixed specification, not evidence against the model family.

5. Results

5.1. The Baseline Specification: The Wrong Kind of Success

The reproduced baseline performance is shown in Table 5. Random Forest and LightGBM attain R² = 0.9871; within-one-bay accuracy is 98.1% and 98.7%, respectively. XGBoost reaches R² = 0.9856, Poisson regression 0.9648, and the fixed-α negative binomial specification 0.7505. The three tree models achieve 96.4% - 97.0% exact-integer accuracy and 99.86% - 100% zero-bay accuracy. These are strong results for reproducing the recorded generated target. They do not assess agreement with independently observed facility counts. The baseline Random Forest subsequently flags 685 sites for the full inventory.

Table 5. Test-set performance of the baseline specification (B0; random 80/20 split, n = 9537).

Note: ±1 bay = share of predictions within one bay of the recorded count; Exact = share of rounded predictions equal to the recorded count; Zero-bay acc. = exact accuracy on sites recording zero bays. LightGBM is early-stopped on the test set, as specified in B0.

5.2. Experiment 1: Accuracy Is Rule Recovery, and Spatial Context Adds Nothing

The summary of the outcome of E1 is presented in Table 6 and Figure 4. These comprise 330 scheduled learned-model fits: six feature sets × five models × (one hold-out split + five random folds + five blocked folds). From the setup, four findings emerged.

Figure 4. Predictive-fit audit (E1). Mean R² (a) and exact-integer accuracy (b) under random five-fold and spatially blocked five-fold evaluation. The plotted “F1 Baseline B0 (22)” label denotes the deduplicated audit feature set, not the original 24-column B0 matrix. The reconstructed rule, identified using all recorded outcomes, has approximately unit R² and 99.99% exact agreement; these are reconstruction diagnostics rather than independently validated predictive scores.

Table 6. Accuracy decomposition (E1): mean across folds under random (R) and spatially blocked (B) 5-fold cross-validation.

Note: RF = Random Forest (max_features = √p); LGBM = LightGBM (early-stopped on a validation split of the training fold). XGBoost results follow LightGBM closely, and the negative binomial GLM is unstable (mean R² 0.21 - 0.75 on F0 - F3); both are reported in the supplementary tables. Blocked design: 20 km grid cells.

First, rule-input features provide the strongest Random Forest fit. With F0, R² is 0.9999 under random cross-validation and 0.9998 under blocked validation. The exact accuracy is 99.65% and 99.63%, respectively. Adding the other 17 unique baseline features lowers random-fold exact accuracy to 97.18% and increases MAE from 0.0077 to 0.0944 bays. With max_features = √p, each F1 split considers four of the 22 columns after integer truncation. The probability of considering informative rule-input columns therefore changes with the feature representation. The result demonstrates specification sensitivity in this forest.

Second, spatial and regional features alone have weak predictive performance. F4 Random Forest reaches R² = 0.0986 on the random hold-out, but the three tree models have negative R² under blocked validation (approximately −0.07 to −0.12). Thus, their squared errors exceed the test-fold mean-reference error. This pattern is compatible with spatial composition supplying limited local predictive information that fails under geographic blocking. It is not evidence that geography is irrelevant to actual loading demand. Moreover, geographic features were constructed on the full inventory, so this is not a fully inductive test at unseen locations.

Third, target-derived inputs are predictive even without floor area or category. F5 Random Forest achieves R² ≈ 0.80 and approximately 89% exact accuracy. This is expected because its inputs contain the outcome, including the site’s own count through Moran statistics. When rule inputs are also present, leakage need not increase accuracy. Removing all four target-derived features improves F1 Random Forest random-fold R² from 0.9903 to 0.9946. Predictive ablation alone, therefore, cannot justify the absence of leakage. Its presence is established by feature construction, while its measured effect depends on the other available inputs.

Fourth, validation designs produce similarly high R² for models supplied with the rule inputs, but not identical scores. The E1 output includes an XGBoost F0 difference of approximately 0.0033 between random and blocked five-fold R² (0.9871 versus 0.9838). All the tree configurations nevertheless remain highly predictive of the constructed target. Blocking removes co-location across folds but cannot remove direct target leakage retained in F1, F2, and F5, or establish target validity.

The Poisson specification attains R² around 0.97 but approximately 48% exact-integer accuracy in the reproduced hold-out. A log-linear conditional mean does not readily recover the rounded, capped mapping. The fixed-α negative binomial model is unstable across feature sets and validation designs and performs particularly poorly on F5. These comparisons concern the specifications, which include redundant encodings and clipping. They do not establish the general superiority of tree models over appropriately tuned count regressions. Reporting R² with metrics like MAE and integer-based metrics makes the practical consequences of these differences clearer.

5.3. Experiment 2: Model Attribution Extends beyond the Reconstructed Rule

The attribution and prediction sensitivity for seven configurations on each of F1 and F3 are shown in Figure 5 and Table 7, respectively. Rule-input attribution and attribution to other feature families are distinguished without interpreting SHAP as a causal estimator.

Figure 5. Shares of mean absolute SHAP attribution by feature family for seven configurations on F1 and F3. The figure’s legacy “F1 Baseline B0 (22)” label denotes the deduplicated audit variant. Floor area and category are inputs to the reconstructed target rule; other families may be used by the fitted model without appearing in that rule.

Table 7. Attribution relative to the reconstructed rule and prediction sensitivity for seven configurations on F1 and F3.

Note: SAS is the share of the mean absolute SHAP attribution outside encoded rule inputs; strict SAS omits category-region proxies from the numerator but retains the full denominator. Permutation SAS normalises positive mean R² reductions after clipping negative importances to zero. Top-five precision is representation-dependent because several columns encode the same rule input. Prediction change is measured under independent within-rule-output-stratum permutations, which may disrupt feature dependencies. The first row is the deduplicated F1 audit forest, not the original 24-column B0 fit.

In the F1 Random Forest configuration with max_features = √p, 33.7% of mean absolute SHAP attribution falls outside the encoded rule-input families. Target-derived features receive 22.1%, and category-region proxies receive approximately 9.7%; the quadrant indicator ranks third among individual columns (Figure 6(b)). Three of the five highest-ranked columns encode rule inputs. These numbers describe the fitted 22-feature model. They do not establish that SHAP misrepresents that model, nor do they estimate the contribution of geography to actual loading-bay provision.

Purely spatial features receive approximately 1.1% of attribution in the F1 forest and at most 2.1% across the 14 fits. The largest off-rule share comes from outcome-derived features and category-region proxies, rather than unsupervised geographic context. The audit therefore identifies the information path behind apparently plausible explanations instead of assuming that every engineered spatial feature materially drives the model.

Figure 6. (a) Combined floor-area SHAP attribution for LightGBM on F3, coloured by category, broadly follows the threshold and category-specific slopes of the reconstructed rule. (b) Twelve highest mean absolute attributions for the deduplicated F1 Random Forest with max_features = √p. Features marked ✗ are outside the encoded rule-input families; this marker does not imply inaccurate attribution to the fitted model or demonstrated causal irrelevance in actual provision.

Three patterns qualify the attribution result. First, increasing Random Forest max_features from √p to 1.0 raises hold-out R² from 0.9895 to 0.9999 and reduces SHAP SAS from 33.7% to 4.4%; strict SAS falls from 24.0% to 0.3%. Model configuration therefore changes how closely fitted dependencies follow the rule. Second, SHAP SAS and permutation SAS differ. The √p forest has 33.7% versus 12.5%, while some boosted configurations have larger permutation shares. Because the measures use different scales, perturbations, and normalisations, their disagreement is not a direct ranking of explainer correctness. Third, the floor-area attribution broadly reflects the rule’s threshold and category-specific shape (Figure 6(a)). The rank correlation between off-rule attribution and distance to the nearest rounding step (Section 4.4) is ρ = −0.08 (p = 0.07, n = 508). Its sign is consistent with slightly more off-rule attribution near rounding steps, but the association is negligible and not distinguishable from zero at the 5% level. It therefore provides no evidence that off-rule attribution concentrates at the rule’s breakpoints.

Figure 7. Prediction sensitivity under independent permutation of features outside the rule-input families within strata sharing the same reconstructed rule output. Bars report the proportion of rounded predictions that change. F1 is the 22-feature audit variant. An exact rule predictor would be invariant; changes in fitted predictions may also reflect inconsistent feature combinations introduced by permutation.

The stress test changes rounded predictions at 14.0% of sites for the F1 √p forest and 13.5% for its F3 counterpart (Figure 7). Boosted-model values are approximately 4.6% - 8.5%, while Random Forest with max_features = 1.0 is close to 1%. These results are consistent with configuration-dependent reliance on inputs outside Equation (1). Because independent permutations can move records away from valid joint feature combinations, they do not isolate that reliance from extrapolation effects. With attribution and ablation, the test provides a diagnostic of rule fidelity. It does not establish a failure of SHAP or identify causal drivers in the real world.

5.4. Experiment 3: The Benchmark Flags Model Error, Not Under-Provision

The summary of E3 is shown in Table 8 and Figure 8, respectively. Out of the 685 baseline flags, 684 (99.85%) have recorded counts equal to Equation (1). Their positive residuals therefore represent surrogate error relative to the reconstructed rule. The discrepancy in the remaining record from the rule is not evidence of real under-provision either. One-bay residuals account for 522 flags (76.2%). Flagged sites have a median recorded count of 15 bays versus one across the inventory, with 623 warehouses, 57 factories, and five retail sites. The flag rate is 1.92% for the test partition and 1.32% for the training partition, consistent with their different fitted-prediction status.

Figure 8. Benchmark validity (E3). (a) Sites flagged as under-provisioned under out-of-fold, spatially blocked prediction, by feature set and model; bars are means over five seeds, whiskers indicate the range; the dashed line marks the 685 flags of the B0 in-sample run. (b) Jaccard similarity between all 45 flagged sets (blocks of five seeds). (c) Size of the predicted deficit among the 685 B0 flags. (d) Recall (filled markers) and precision (open markers) for placebo deficits of 1, 2, and 5 bays injected into 2% of sites.

Table 8. Validity, stability, and detection power of residual-based benchmarking (E3).

Note: Each replicate subtracts exactly d bays from 953 rule-compliant records with at least d + 1 bays, leaving positive counts. Eligible populations contain 15,656, 10,124, and 4836 records for d = 1, 2, and 5, so the perturbed share is 6.1%, 9.4%, and 19.7% of each eligible population (2.0% of the full inventory). Recall and precision are means over three seeds, with no population uncertainty interval implied. F1 target-derived features are recomputed on the entire contaminated target, deliberately retaining held-out outcome information.

Out-of-fold flag counts range from 46 to 1124 across the 45 configurations. Mean Jaccard similarity is 0.79 across seeds for a fixed configuration, 0.25 across models at a fixed feature set and seed, and 0.45 across feature sets at a fixed model and seed. Overall, 3036 sites are flagged at least once and six are flagged by every configuration. Between 97.8% and 99.9% of flagged records comply with Equation (1). These results establish substantial specification sensitivity of the residual flags of the audit.

The controlled perturbations assess detection of deviations from the rule. F0 Random Forest recall is 98.8% - 99.9% and precision is 95.9% - 98.4% across the three decrement sizes. F0 LightGBM recall is 94.9% - 97.7%, with precision of 59.0% - 76.0%. For one-bay decrements, F1 recall falls to 60.0% for Random Forest and 59.7% for LightGBM, with precision of 44.4% and 47.8%, respectively. This deterioration is consistent with contaminated outcome-derived inputs allowing models to reproduce the injected deviations. However, F0 and F1 also differ in other features and effective feature sampling, so the comparison does not isolate a causal effect of leakage alone. The decrement-specific populations and three-replicate design further limit comparisons across deficit sizes. The test demonstrates detection behaviour under the specified contamination scenario, not performance against field-verified shortages.

5.5. Consistency Checks and Scope of Robustness

The results support several consistency checks. Rule-input tree models retain high R² under random and blocked validation, but spatial-only tree performance deteriorates under blocking. Removing the regional median from F2 produces the same rounded LightGBM and Poisson metrics in Table 6. This is consistent with limited incremental information, but does not prove that those inputs are mathematically interchangeable. Off-rule SHAP attribution is evident under the strict definition excluding category-region proxies. Flag agreement is higher across seeds than across model families, and high rule compliance holds for out-of-fold flags as well as baseline flags. These checks strengthen the core provenance conclusion within this case.

6. Discussion

6.1. What the Models Learned, and What They Did Not

The main finding is that strong predictive performance can coexist with a mismatch between the modelled target and the phenomenon decision-makers intend to assess. The baseline achieves R² near 0.99 and up to 98.7% within-one-bay accuracy (Table 5), but the recorded outcome is almost entirely reproduced by a floor-area-and-category rule. The results establish a successful approximation of that constructed field. They provide no independent evidence that the models recover actual loading-bay counts. The failure occurs when this performance is interpreted as evidence of observed provision.

The distinction concerns analytical use and not flaws in the source inventory. de Saxe et al. [6] [15] explicitly describe the bay field as modelled and designed the inventory to support system-level freight analysis. Such a field can be appropriate for scenario modelling when its assumptions and uncertainty remain explicit. Using it to infer actual counts or site-level inadequacy creates a different evidential demand. This case illustrates why provenance documentation must be incorporated into task definition and interpretation, rather than cited only as a data-source description [40].

The geographic results are likewise conditional on the constructed target. Equation (1) has no geographic input, and the selected spatial feature family provides little transferable predictive fit under the executed blocked design. Contextual additions also change Random Forest feature sampling and residual behavior. These findings justify testing feature families against a simple reference rather than assuming that more spatial detail improves a benchmark. They do not imply that location is unimportant for actual freight operations, or that every geographic feature requires a causal interpretation to be useful for prediction.

6.2. Why Explanation Plausibility Does Not Establish Domain Validity

The attribution analysis clarifies the role of XAI in operational governance [10] [11]. SHAP explains a fitted model under a chosen attribution formulation. It does not independently verify the target or recover causal effects. When a model uses outcome-derived or contextual inputs, attribution to those inputs can be faithful to the model while extending beyond the rule that generated its training labels. The 33.7% share therefore quantifies off-rule attribution in one configuration, not a general SHAP failure rate. The distinction between model behaviour and domain interpretation is consistent with the dependence-sensitive attribution literature [43] [44].

Changing the forest’s feature-sampling setting substantially improves both predictive fit and concentration of attribution on rule inputs. This demonstrates that interpretation depends on model design as well as on the explanation method. In this dataset, the explicit reconstructed mapping provides the clearest account of target generation. It should not be treated as a validated physical model because it is transparent or highly predictive.

Differences between SHAP and permutation importance reinforce the need to specify what each measure evaluates [45]-[47]. They should not be presented as interchangeable confirmations of causality. For planning applications, a useful reporting standard is to disclose feature construction, attribution assumptions, sample selection, and sensitivity to coherent perturbations. Agreement between explanation methods may support consistency, but agreement is neither necessary for their distinct estimands to be valid nor sufficient to establish the target’s substantive validity.

Documented generated targets can nevertheless support methodological research. They combine real covariates with an auditable construction reference. This allows researchers to examine rule approximation and attribution under realistic correlations. Such testbeds should be explicitly labeled as target-construction audits. Evaluation against the known rule and validation against actual field outcomes answer different questions and should remain separate.

6.3. Implications for Loading-Bay Planning and Benchmarking

For planning practice, there are three implications. First, the inventory is a useful source of processed site attributes and modelled HGV bay estimates. However, its bay field should not be treated as independently observed provision. Second, residuals from an ML surrogate of that field identify discrepancies from a constructed reference rather than operational shortages. Even with independently observed counts, a peer norm would describe typical provision unless linked to demand and service requirements. Third, any residual-based screening tool should disclose threshold sensitivity, out-of-fold flag stability, and the consequences of false positives and missed cases. Controlled perturbations can assess detection of specified deviations, but field verification is needed before flags support site-level planning decisions.

The reconstructed parameters also show how engineering assumptions can be inherited by downstream models. Above 2000 m², warehouse increments occur at approximately 823.5 m² per additional bay and factory increments at approximately 6532 m², with caps; larger retail records receive one bay. The imagery sample comprises 97 sites in total, with modeling based on selected observations and category-specific assumptions [6]. These choices were made for the inventory’s modeling purpose. An ML surrogate cannot turn them into representative empirical evidence of adequate capacity. A defensible adequacy benchmark would require verified bay counts linked to loading demand, throughput, vehicle mix, and operating hours, with uncertainty and sampling coverage reported.

6.4. Theoretical Implications: Targets as Instruments

The conceptual contribution is to treat target construction as a measurement procedure that defines the learning task. Here, the targeted construct is physical HGV loading-bay provision, but the available label is an estimate produced from floor area and category. Training on that estimate evaluates approximation of the measurement procedure. Literature on generated variables and imputation highlights the need to account for construction assumptions and uncertainty [38] [39]. It does not establish that every generated label necessarily biases every predictive analysis. Suitability depends on whether the intended task is surrogate modelling, scenario analysis, estimation of actual provision, or operational decision support.

This interpretation connects target validity with shortcut learning and leakage audits while keeping the concepts distinct [12] [53]. Rule inputs can predict a generated label without leaking an unknown outcome, but the result may still be misinterpreted as evidence about actual provision. Conversely, target-derived features constitute an additional information problem even when rule inputs already determine the label. Target provenance, feature availability, and validation independence should therefore be audited separately. The proposed protocol makes these dependencies explicit rather than treating high accuracy as a substitute for measurement validity.

6.5. A Target Provenance and Validity Protocol

Table 9 presents the Target Provenance and Validity Protocol for planning benchmarks. It is an evidence-informed checklist derived from this case, not a prospectively validated standard. Provenance and distributional inspection should precede predictive modelling; a simple surrogate test can then reveal whether recorded outcomes are nearly recoverable from a documented rule. The remaining steps examine feature availability, deployment-aligned validation, incremental predictive value, interpretation, and decision consequences. A failed check requires clarification of the permissible analytical use instead of automatic rejection of every modelled target.

6.6. Limitations

There are a few limitations that define the scope of the findings. F3 excludes target-derived features but still uses transductive preprocessing. With this, the blocked folds are approximate and unbuffered. The independent within-stratum permutations disrupt joint feature relationships, so their sensitivity estimates cannot isolate off-rule dependence from extrapolation. Also, TreeSHAP settings rely on the installed implementation, additivity checking is disabled, and attribution samples were the first selected test rows rather than repeated samples. Lastly, no independently verified site counts, demand data, or operational adequacy assessments are available.

6.7. Future Research

Future work should prioritise independent measurement. A representative sample of verified HGV loading-bay counts, linked to throughput, vehicle mix, and service requirements, would permit evaluation of both the inventory’s imputation rule and any operational benchmark. Methodological extensions should reconstruct parameters within training folds when estimating generalisation. Inductive geographic features should be derived from training data. Repeated attribution sampling, explicit SHAP background settings, and more deficit-injection replicates would improve uncertainty assessment. Importantly, prospective applications of TPVP should test whether it changes modelling choices and decisions.

Table 9. Target provenance and validity protocol (TPVP) for ML benchmarking in planning.

7. Conclusions

This audit shows how a highly accurate loading-bay predictor can reproduce a constructed target without establishing actual site provision. A thresholded, rounded, and capped rule matches 99.992% of recorded counts, and a simple model using its inputs outperforms the reproduced spatially engineered baseline. Outcome-derived inputs create direct information leakage. Other contextual inputs influence fitted dependencies and residual flags without appearing in the reconstructed rule. Out of the 685 baseline flags, 684 comply exactly with that rule. Flag lists vary across specifications, and detection of simulated downward deviations is weaker with the full feature set than with rule inputs alone. These findings concern target approximation and constructed perturbations, not field-verified loading shortages.

The central implication is that target validity must precede predictive benchmarking, and decision validity must follow it. Accuracy evaluates agreement with the supplied label, while attribution describes a fitted model. Neither independently validates what the label measures or whether a residual warrants action. The Target Provenance and Validity Protocol provides a structured way to make those distinctions explicit. Used within its stated limits, it can help planning analysts distinguish useful surrogate and scenario models from claims that require independently observed infrastructure and operational evidence.

Data Availability

The logistics-site inventory is openly available from the University of Cambridge Apollo repository [15]. The reported analyses use its five-field final_dataset.csv file. No independent site-level observations of actual bay counts were collected for this audit.

Code Availability

The analysis is documented in LoadingBay_GroundTruth_Audit_RUN. ipynb, with saved execution outputs. Loading_Bay_Prediction_Pipeline. ipynb supplies the baseline specification. The audit notebook exports tables, figures, and a results summary when run with the source CSV. The saved outputs support the numerical results reported here; exact regeneration also requires the data file and a compatible software environment. Public archival availability of the complete reproducibility package should be established before publication.

Author Contributions

Conceptualization, Tewogbade, S.A.; Methodology, Oraegbunam, L.; Software, Tewogbade, S.A.; Validation, Omowaiye, R.; Formal analysis, Ayodeji, B.; Investigation, Tewogbade, S.A.; Data curation, Oraegbunam, L.; Writing—Omowaiye, R.; Writing—review and editing, Tewogbade, S.A.; Visualization, Oraegbunam, L.; Supervision, Tewogbade, S.A.

All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References