1 Purpose
2 1. Setup and data loading
- 2.1 1.1 Packages, configuration, and modules
- 2.2 1.2 Source-data loader
3 2. Variogram modeling
4 3. Kriging and cross-validation
5 4. Synthetic groundwater-level generation
6 5. Monte Carlo Experiment I
- 6.1 5.1 Construct the Experiment I weights
- 6.2 5.2 Run the Monte Carlo experiment
7 6. Experiments II-V
8 7. Model comparison and coefficient exports
9 8. Formal SLX temperature tests
10 9. Figures and archival record
11 Generate manuscript outputs
12 Interpretation safeguards

1 Purpose

This document executes and explains the complete replication flow. It begins with the source-data loader, estimates groundwater and rainfall variograms, performs regression-kriging cross-validation, generates synthetic groundwater fields, runs the Monte Carlo Moran’s I experiment, and then estimates the full set of spatial-regression models.

The code is modular. The same functions are called by run_replication.R, so the walkthrough and the one-command driver implement the same statistical workflow.

2 1. Setup and data loading

2.1 1.1 Packages, configuration, and modules

source(file.path("R", "00_utils.R"))

check_packages(c(
  "sf",
  "spdep",
  "spatialreg",
  "gstat",
  "Matrix",
  "sandwich",
  "ggplot2",
  "dplyr",
  "forcats",
  "readr",
  "terra",
  "knitr"
))

source("config.R")

module_files <- c(
  "01_data.R",
  "02_variograms.R",
  "03_kriging.R",
  "04_simulation.R",
  "05_weights.R",
  "06_monte_carlo.R",
  "07_models.R",
  "08_diagnostics.R",
  "09_temperature_tests.R",
  "10_plots.R"
)

for (module in module_files) {
  source(file.path("R", module))
}

set.seed(config$seed)
ensure_dir(config$output_dir)
ensure_dir(file.path(config$output_dir, "tables"))
ensure_dir(file.path(config$output_dir, "figures"))
ensure_dir(file.path(config$output_dir, "rasters"))
ensure_dir(file.path(config$output_dir, "spatial"))

The configuration file contains every analysis choice that should be changed by a replicator: input paths, year and region filters, formulas, CRS, distance bands, inverse-distance power, variogram settings, kriging neighbourhood size, and the Monte Carlo iteration count.

2.2 1.2 Source-data loader

The archive retains the original function name load.data(). The rewritten loader uses cross-platform paths, validates and projects the prepared groundwater observations, loads the study-area polygon, and loads the covariate-complete kriging grid. The prediction grid is stored in the input CSV as WGS84 longitude/latitude (longitude, latitude), converted to an sf point object, and transformed to the same projected CRS as the groundwater observations. The original longitude/latitude columns are retained for raster export, while projected coordinates are used for kriging and distance calculations.

spatial_data <- load.data(config)

## Removed 652980 duplicated longitude/latitude rows from the interim prediction grid.

samples <- prepare_analysis_samples(spatial_data, config)

write_csv(
  samples$audit,
  file.path(config$output_dir, "tables", "sample_audit.csv")
)

knitr::kable(samples$audit, digits = 0)

sample	observations	rows_removed
Loaded prepared spatial data	1310	0
Model sample	1310	0
Variogram sample	1310	0

Two complete-case observation samples are retained. The variogram sample uses the fuller trend specification required by the manuscript’s variogram analysis. The model sample uses the common dependent variable and three regressors used in Experiments II-V. The returned samples object also contains samples$study_area, samples$full_grid, and samples$full_grid_sf, so all downstream prediction surfaces use the same validated study boundary and covariate-complete grid.

3 2. Variogram modeling

The replication now retains two groundwater residual variograms for two different purposes. The full environmental-trend variogram is the original manuscript variogram and remains the covariance model used by the spatial-dependence analysis, variogram-informed weights, simulations, and the spatial-regression experiments. A second prediction-specific variogram is estimated after removing the reduced longitude–latitude mean used by the final universal-kriging surface. Keeping the two objects separate prevents the final prediction uncertainty from being calculated with a covariance model fitted under a different mean specification.

3.1 2.1 Full manuscript groundwater variogram

The full groundwater semivariogram is estimated from residual dependence after conditioning on the environmental trend in config$variogram_regressors:

full_variogram_formula <- stats::reformulate(
  config$variogram_regressors,
  response = config$dependent
)
full_variogram_formula

## PREMONSOON ~ L1RF + L0T2M.FULL + L1RO.SUM + PL1DTC + PL1RABIZAID + 
##     PL1BUILTUP + DNRV + DNRD + SUMPHB5K

The empirical omnidirectional variogram is calculated with the configured 5 km bins and 150 km cutoff. Directional empirical variograms are also calculated at 0, 45, 90, and 135 degrees and a spherical model is fitted to each. These are the variograms used by the directionality analysis and by all later VIW-based work.

3.2 2.2 Prediction-specific groundwater variogram

The final universal-kriging prediction uses the reduced deterministic mean in config$kriging$regressors, currently longitude and latitude. Its residual variogram is re-estimated separately rather than reusing the full environmental variogram.

kriging_variogram_formula <- stats::reformulate(
  config$kriging$regressors,
  response = config$dependent
)
kriging_variogram_formula

## PREMONSOON ~ LONGITUDE + LATITUDE

The spatial coordinates used to calculate variogram separation distances remain in the projected analysis CRS. LONGITUDE and LATITUDE enter only as covariates in the deterministic mean.

3.3 2.3 Rainfall variogram and variogram fitting

Rainfall receives a separate empirical variogram and spherical fit. The default rainfall formula contains a coordinate trend, matching the supplied source analysis. This rainfall covariance model is later used only where a rainfall-specific spatial lag is required.

variogram_results <- fit_variogram_models(
  spatial_data = samples$variogram,
  config = config
)

export_variogram_results(
  variogram_results,
  config$output_dir
)

knitr::kable(
  variogram_results$parameter_table,
  digits = 3,
  caption = "Estimated variogram parameters"
)

Estimated variogram parameters
specification	process	direction	model	partial_sill	total_sill	range_map_units	range_km
Full environmental trend	Groundwater residual	Omnidirectional	Sph	20.462	20.462	21967.47	21.967
Full environmental trend	Groundwater residual	0	Sph	23.084	23.084	39100.43	39.100
Full environmental trend	Groundwater residual	45	Sph	20.978	20.978	21707.84	21.708
Full environmental trend	Groundwater residual	90	Sph	19.724	19.724	16324.10	16.324
Full environmental trend	Groundwater residual	135	Sph	21.003	21.003	23961.15	23.961
Kriging trend	Groundwater residual	Omnidirectional	Sph	21.969	21.969	24312.03	24.312
Kriging trend	Groundwater residual	0	Sph	24.210	24.210	43583.13	43.583
Kriging trend	Groundwater residual	45	Sph	22.918	22.918	24096.44	24.096
Kriging trend	Groundwater residual	90	Sph	21.746	21.746	18241.38	18.241
Kriging trend	Groundwater residual	135	Sph	21.999	21.999	25178.55	25.179
Rainfall trend	Rainfall	Omnidirectional	Sph	0.026	0.026	136680.82	136.681

plot(
  variogram_results$groundwater$empirical,
  variogram_results$groundwater$fitted
)

Full environmental-trend groundwater variogram and fitted spherical model

plot(
  variogram_results$groundwater_kriging$empirical,
  variogram_results$groundwater_kriging$fitted
)

Prediction-specific longitude-latitude residual variogram and fitted spherical model

plot(
  variogram_results$rainfall$empirical,
  variogram_results$rainfall$fitted
)

Rainfall variogram and fitted spherical model

For backwards compatibility, variogram_results$groundwater and variogram_results$formula continue to refer to the full environmental trend. The prediction-only covariance model is stored separately in variogram_results$groundwater_kriging, with its formula in variogram_results$kriging_formula.

4 3. Kriging and cross-validation

Two distinct cross-validation exercises are retained and should not be conflated.

The directionality cross-validation holds the full environmental trend fixed and compares the full-trend omnidirectional and directional residual variograms. This is the original exercise used to assess covariance directionality.
The prediction-model cross-validation uses the reduced deterministic mean selected for the final universal-kriging map. Its residual variogram is fitted under the same reduced mean. The separate trend-selection robustness exercise is implemented in R/13_tests.R.

4.1 3.1 Full-trend directionality cross-validation

directionality_cv <- run_full_trend_directionality_cross_validation(
  spatial_data = samples$variogram,
  variogram_results = variogram_results,
  config = config
)

## Full-trend directionality CV: omnidirectional

## Full-trend directionality CV: direction_0

## Full-trend directionality CV: direction_45

## Full-trend directionality CV: direction_90

## Full-trend directionality CV: direction_135

export_directionality_cv_results(
  directionality_cv,
  config$output_dir
)

knitr::kable(
  directionality_cv$metrics,
  digits = 3,
  caption = "Full environmental-trend variogram directionality diagnostics"
)

Full environmental-trend variogram directionality diagnostics
	variogram_model	observations	mean_error	MAE	RMSE	correlation	fitted_R2	calibration_intercept	calibration_slope	mean_kriging_variance	mean_kriging_se	mean_standardized_error	RMS_standardized_error
omnidirectional	omnidirectional	1310	-0.075	2.540	3.680	0.785	0.616	2.390	0.760	8.085	2.646	-0.014	1.467
direction_0	direction_0	1310	-0.072	2.554	3.728	0.780	0.608	2.418	0.757	5.148	2.110	-0.013	1.841
direction_45	direction_45	1310	-0.075	2.538	3.676	0.785	0.617	2.389	0.760	8.387	2.695	-0.014	1.439
direction_90	direction_90	1310	-0.089	2.532	3.647	0.788	0.621	2.410	0.759	10.498	3.018	-0.019	1.287
direction_135	direction_135	1310	-0.076	2.546	3.696	0.783	0.613	2.395	0.759	7.619	2.567	-0.015	1.514

This table belongs to the original directionality assessment. It uses variogram_results$groundwater, not the prediction-specific groundwater_kriging object.

4.2 3.2 Final longitude–latitude universal kriging

The final prediction surface uses the reduced mean in config$kriging$regressors and the matching residual covariance model in variogram_results$groundwater_kriging. The default neighbourhood contains at most config$kriging$nmax observations.

kriging_results <- run_complete_kriging_analysis(
  spatial_data = samples$variogram,
  variogram_results = variogram_results,
  config = config,
  boundary = samples$study_area,
  prediction_grid = samples$full_grid_sf,
  clip_to_boundary = config$kriging$clip_to_study_area
)

## Kriging cross-validation: omnidirectional

## Generating kriging surface: omnidirectional

raster_resolution_m <- if (!is.null(config$kriging$raster_resolution_m)) {
  config$kriging$raster_resolution_m
} else {
  2000
}

export_kriging_results(
  results = kriging_results,
  output_dir = config$output_dir,
  spatial_data = samples$variogram,
  raster_resolution_m = raster_resolution_m
)

## `geom_smooth()` using formula = 'y ~ x'

## Exporting kriging surfaces: omnidirectional

## Wrote GeoTIFF: C:/Users/hp/Dropbox/root/work/iiitd/Research/analyses/spauto_ecol_econ_replication_aug_1_2026/output/rasters/kriging_prediction_omnidirectional.tif

## Wrote GeoTIFF: C:/Users/hp/Dropbox/root/work/iiitd/Research/analyses/spauto_ecol_econ_replication_aug_1_2026/output/rasters/kriging_variance_omnidirectional.tif

## Wrote GeoTIFF: C:/Users/hp/Dropbox/root/work/iiitd/Research/analyses/spauto_ecol_econ_replication_aug_1_2026/output/rasters/kriging_standard_error_omnidirectional.tif

## Kriging rasters available for omnidirectional: kriging_prediction_omnidirectional.tif, kriging_variance_omnidirectional.tif, kriging_standard_error_omnidirectional.tif

## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
## Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.

knitr::kable(
  kriging_results$metrics,
  digits = 3,
  caption = "Leave-one-out diagnostics for the final prediction specification"
)

Leave-one-out diagnostics for the final prediction specification
	variogram_model	observations	mean_error	MAE	RMSE	correlation	fitted_R2	calibration_intercept	calibration_slope	mean_kriging_variance	mean_kriging_se	mean_standardized_error	RMS_standardized_error
omnidirectional	omnidirectional	1310	-0.006	2.139	3.007	0.851	0.724	2.233	0.769	5.236	2.215	-0.001	1.809

The exporter writes the point surface as CSV and GeoPackage for overlay and, critically, writes the three analysis rasters before attempting any map plots: prediction, prediction variance, and prediction standard error. This ordering ensures that a plotting problem cannot leave only the point GeoPackage behind.

selected_surface_model <- if (!is.null(config$kriging$surface_models)) {
  config$kriging$surface_models[[1]]
} else {
  config$kriging$models[[1]]
}

expected_rasters <- data.frame(
  layer = c("prediction", "variance", "standard_error"),
  file = file.path(
    config$output_dir,
    "rasters",
    c(
      paste0("kriging_prediction_", selected_surface_model, ".tif"),
      paste0("kriging_variance_", selected_surface_model, ".tif"),
      paste0("kriging_standard_error_", selected_surface_model, ".tif")
    )
  ),
  stringsAsFactors = FALSE
)
expected_rasters$exists <- file.exists(expected_rasters$file)

knitr::kable(
  expected_rasters,
  caption = "Kriging raster export check"
)

Kriging raster export check
layer	file	exists
prediction	output/rasters/kriging_prediction_omnidirectional.tif	TRUE
variance	output/rasters/kriging_variance_omnidirectional.tif	TRUE
standard_error	output/rasters/kriging_standard_error_omnidirectional.tif	TRUE

if (!all(expected_rasters$exists)) {
  stop(
    "Kriging raster export is incomplete: ",
    paste(expected_rasters$file[!expected_rasters$exists], collapse = ", ")
  )
}

Prediction standard error is computed as sqrt(var1.var) and is expressed in the same units as groundwater depth. The kriging variance and standard-error surfaces are conditional on the fitted longitude–latitude mean and its matching residual variogram.

4.3 3.3 Optional trend-selection robustness exercise

R/13_tests.R contains the second LOOCV exercise used to choose the kriging mean specification: ordinary kriging, longitude–latitude universal kriging, and the full environmental trend. It includes both the primary comparison with a separately re-fitted residual variogram for each mean specification and the fixed-variogram sensitivity test. This exercise is optional in the walkthrough because it is computationally redundant once the prediction specification has been selected.

source(file.path("R", "13_tests.R"))

knitr::kable(
  trend_test$metrics,
  digits = 3,
  caption = "Trend-selection LOOCV with model-specific residual variograms"
)

knitr::kable(
  trend_test_fixed$metrics,
  digits = 3,
  caption = "Trend-selection sensitivity test with the variogram held fixed"
)

5 4. Synthetic groundwater-level generation

Synthetic GWL is generated at the actual monitoring locations. The generator combines:

fitted values from the regression trend;
Gaussian residual fields with the fitted groundwater nugget, partial sill, and spherical range;
optional truncation at zero.

The covariance matrix and Cholesky factor are calculated once and reused across all Monte Carlo iterations.

synthetic_generator <- prepare_synthetic_generator(
  spatial_data = samples$variogram,
  variogram_results = variogram_results,
  config = config
)

synthetic_first <- simulateDataReal(
  generator = synthetic_generator,
  seed = config$seed,
  truncate.at.zero = config$simulation$truncate_at_zero
)

export_synthetic_realisation(
  synthetic_first,
  config$output_dir,
  dependent = config$dependent
)

summary(synthetic_first$simulated_GWL)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   6.337   9.631   9.765  12.953  28.263

hist(
  synthetic_first$simulated_GWL,
  breaks = 35,
  main = "First synthetic GWL realization",
  xlab = "Synthetic groundwater level"
)

Distribution of the first synthetic groundwater realization

6 5. Monte Carlo Experiment I

Experiment I isolates the sensitivity of Moran’s I to the weighting scheme. Each Monte Carlo iteration uses the same fitted trend and covariance process but a new Gaussian realization. Moran’s I is then evaluated under:

binary distance bands of 22, 27, 37, 47, and 57 km;
unrestricted inverse-distance weights with exponent 1;
groundwater variogram-informed weights.

Rainfall VIW is not part of Experiment I because the simulated process is specifically a groundwater process.

6.1 5.1 Construct the Experiment I weights

experiment_i_weights_object <- create_all_weights(
  spatial_data = samples$variogram,
  groundwater_variogram_fit = variogram_results$groundwater$fitted,
  rainfall_variogram_fit = variogram_results$rainfall$fitted,
  distance_bands_km = config$distance_bands_km,
  inverse_power = config$inverse_power,
  style = config$weight_style,
  zero_policy = config$zero_policy
)

## Warning in spdep::mat2listw(weights, style = style, zero.policy = zero_policy):
## neighbour object has 2 sub-graphs
## Warning in spdep::mat2listw(weights, style = style, zero.policy = zero_policy):
## neighbour object has 2 sub-graphs

experiment_i_weights <- experiment_i_weights_object$listw[
  setdiff(
    names(experiment_i_weights_object$listw),
    "rainfall_variogram"
  )
]

experiment_i_connectivity <- weight_connectivity(experiment_i_weights)
write_csv(
  experiment_i_connectivity,
  file.path(
    config$output_dir,
    "tables",
    "experiment_i_weight_connectivity.csv"
  )
)

knitr::kable(experiment_i_connectivity, digits = 2)

weight_name	connected_components	minimum_neighbours	mean_neighbours	maximum_neighbours
d22km	2	1	37.20	77
d27km	1	4	53.57	107
d37km	1	13	93.42	167
d47km	1	20	142.00	250
d57km	1	27	196.59	323
inverse_distance	1	1309	1309.00	1309
groundwater_variogram	2	1	37.10	77

6.2 5.2 Run the Monte Carlo experiment

if (!isTRUE(params$run_full_monte_carlo)) {
  config$monte_carlo$iterations <- min(
    config$monte_carlo$iterations,
    25L
  )
}

experiment_i <- run_monte_carlo_experiment_i(
  generator = synthetic_generator,
  weights = experiment_i_weights,
  config = config
)

## Experiment I: generating 500 synthetic groundwater fields.

## Warning in ifelse(grepl("^d[0-9.]+km$", weight_names),
## as.numeric(sub("d([0-9.]+)km", : NAs introduced by coercion

export_monte_carlo_results(
  experiment_i,
  config$output_dir
)

knitr::kable(
  experiment_i$summary[
    order(experiment_i$summary$mean_moran_I, decreasing = TRUE),
  ],
  digits = 4,
  caption = "Monte Carlo distribution of Moran's I"
)

Monte Carlo distribution of Moran’s I
	weight_name	simulations	mean_moran_I	sd_moran_I	median_moran_I	confidence_low	confidence_high	mean_difference_from_oracle	RMSE_from_oracle	weighting_scheme	distance_km
6	groundwater_variogram	500	0.5071	0.0348	0.5076	0.4436	0.5762	0.0000	0.0000	Groundwater VIW	NA
1	d22km	500	0.3222	0.0443	0.3226	0.2420	0.4108	-0.1850	0.1857	Binary distance band	22
2	d27km	500	0.2683	0.0452	0.2683	0.1782	0.3573	-0.2389	0.2398	Binary distance band	27
3	d37km	500	0.2084	0.0447	0.2077	0.1213	0.2963	-0.2987	0.2999	Binary distance band	37
4	d47km	500	0.1739	0.0434	0.1736	0.0923	0.2624	-0.3332	0.3346	Binary distance band	47
5	d57km	500	0.1501	0.0416	0.1487	0.0749	0.2358	-0.3570	0.3584	Binary distance band	57
7	inverse_distance	500	0.1238	0.0188	0.1232	0.0898	0.1625	-0.3833	0.3840	Inverse distance	NA

The groundwater VIW result is the process-aligned reference for within-realization comparisons. It should not be described as a weights-free “true Moran’s I,” because Moran’s I is itself conditional on the selected weights matrix.

7 6. Experiments II-V

7.1 6.1 Model-sample weights

The weights are reconstructed on the model sample to guarantee identical row counts and ordering across the outcome, regressors, and every spatial matrix. Both GWL-VIW and rainfall-VIW use the variograms fitted in Section 2.

model_weights_object <- create_all_weights(
  spatial_data = samples$model,
  groundwater_variogram_fit = variogram_results$groundwater$fitted,
  rainfall_variogram_fit = variogram_results$rainfall$fitted,
  distance_bands_km = config$distance_bands_km,
  inverse_power = config$inverse_power,
  style = config$weight_style,
  zero_policy = config$zero_policy
)

## Warning in spdep::mat2listw(weights, style = style, zero.policy = zero_policy):
## neighbour object has 2 sub-graphs
## Warning in spdep::mat2listw(weights, style = style, zero.policy = zero_policy):
## neighbour object has 2 sub-graphs

model_weights <- model_weights_object$listw

model_connectivity <- weight_connectivity(model_weights)
write_csv(
  model_connectivity,
  file.path(config$output_dir, "tables", "weight_connectivity.csv")
)

knitr::kable(model_connectivity, digits = 2)

weight_name	connected_components	minimum_neighbours	mean_neighbours	maximum_neighbours
d22km	2	1	37.20	77
d27km	1	4	53.57	107
d37km	1	13	93.42	167
d47km	1	20	142.00	250
d57km	1	27	196.59	323
inverse_distance	1	1309	1309.00	1309
groundwater_variogram	2	1	37.10	77
rainfall_variogram	1	316	695.82	1069

7.2 6.2 OLS, SAR, and SEM

ols <- fit_ols(
  samples$model,
  config$dependent,
  config$regressors
)

sar <- fit_spatial_family(
  samples$model,
  model_weights,
  family = "SAR",
  dependent = config$dependent,
  regressors = config$regressors,
  method = config$spatial_model_method,
  zero_policy = config$zero_policy
)

## SAR: d22km

## SAR: d27km

## SAR: d37km

## SAR: d47km

## SAR: d57km

## SAR: inverse_distance

## SAR: groundwater_variogram

## SAR: rainfall_variogram

sem <- fit_spatial_family(
  samples$model,
  model_weights,
  family = "SEM",
  dependent = config$dependent,
  regressors = config$regressors,
  method = config$spatial_model_method,
  zero_policy = config$zero_policy
)

## SEM: d22km

## SEM: d27km

## SEM: d37km

## SEM: d47km

## SEM: d57km

## SEM: inverse_distance

## SEM: groundwater_variogram

## SEM: rainfall_variogram

7.3 6.3 SLX specifications

The full SLX spatially lags rainfall, temperature, and multiple cropping. The reduced SLX spatially lags rainfall only. Both are estimated under every common weights matrix.

slx_full <- fit_spatial_family(
  samples$model,
  model_weights,
  family = "SLX",
  dependent = config$dependent,
  regressors = config$regressors,
  durbin = TRUE,
  zero_policy = config$zero_policy,
  model_label = "SLX_all_lags"
)

## SLX_all_lags: d22km

## SLX_all_lags: d27km

## SLX_all_lags: d37km

## SLX_all_lags: d47km

## SLX_all_lags: d57km

## SLX_all_lags: inverse_distance

## SLX_all_lags: groundwater_variogram

## SLX_all_lags: rainfall_variogram

slx_rain <- fit_spatial_family(
  samples$model,
  model_weights,
  family = "SLX",
  dependent = config$dependent,
  regressors = config$regressors,
  durbin = stats::reformulate(config$rainfall_variable),
  zero_policy = config$zero_policy,
  model_label = "SLX_rainfall_lag"
)

## SLX_rainfall_lag: d22km

## SLX_rainfall_lag: d27km

## SLX_rainfall_lag: d37km

## SLX_rainfall_lag: d47km

## SLX_rainfall_lag: d57km

## SLX_rainfall_lag: inverse_distance

## SLX_rainfall_lag: groundwater_variogram

## SLX_rainfall_lag: rainfall_variogram

7.4 6.4 SDM specifications

The archive estimates both the full SDM and a rainfall-lag-only SDM under every common weights matrix.

sdm <- fit_spatial_family(
  samples$model,
  model_weights,
  family = "SDM",
  dependent = config$dependent,
  regressors = config$regressors,
  durbin = TRUE,
  method = config$spatial_model_method,
  zero_policy = config$zero_policy,
  model_label = "SDM_all_lags"
)

## SDM_all_lags: d22km

## SDM_all_lags: d27km

## SDM_all_lags: d37km

## SDM_all_lags: d47km

## SDM_all_lags: d57km

## SDM_all_lags: inverse_distance

## SDM_all_lags: groundwater_variogram

## SDM_all_lags: rainfall_variogram

sdm_rain <- fit_spatial_family(
  samples$model,
  model_weights,
  family = "SDM",
  dependent = config$dependent,
  regressors = config$regressors,
  durbin = stats::reformulate(config$rainfall_variable),
  method = config$spatial_model_method,
  zero_policy = config$zero_policy,
  model_label = "SDM_rainfall_lag"
)

## SDM_rainfall_lag: d22km

## SDM_rainfall_lag: d27km

## SDM_rainfall_lag: d37km

## SDM_rainfall_lag: d47km

## SDM_rainfall_lag: d57km

## SDM_rainfall_lag: inverse_distance

## SDM_rainfall_lag: groundwater_variogram

## SDM_rainfall_lag: rainfall_variogram

7.5 6.5 Mixed-weights robustness models

The reduced hybrid model uses GWL-VIW for the endogenous groundwater lag and rainfall-VIW for lagged rainfall. The full hybrid SDM additionally includes GWL-VIW lags of temperature and multiple cropping.

hybrid_sdm_full <- fit_hybrid_sdm(
  spatial_data = samples$model,
  groundwater_weights = model_weights$groundwater_variogram,
  rainfall_weights = model_weights$rainfall_variogram,
  dependent = config$dependent,
  rainfall = config$rainfall_variable,
  other_regressors = setdiff(
    config$regressors,
    config$rainfall_variable
  ),
  method = config$spatial_model_method,
  zero_policy = config$zero_policy
)

hybrid_sar_rainfall <- fit_hybrid_sar_rainfall_result(
  spatial_data = samples$model,
  groundwater_weights = model_weights$groundwater_variogram,
  rainfall_weights = model_weights$rainfall_variogram,
  dependent = config$dependent,
  rainfall = config$rainfall_variable,
  other_regressors = setdiff(
    config$regressors,
    config$rainfall_variable
  ),
  method = config$spatial_model_method,
  zero_policy = config$zero_policy
)

8 7. Model comparison and coefficient exports

model_results <- list(
  sar,
  sem,
  slx_full,
  slx_rain,
  sdm,
  sdm_rain,
  hybrid_sdm_full,
  hybrid_sar_rainfall
)

model_table <- extract_model_results(
  ols_model = ols,
  model_results = model_results,
  weights = model_weights,
  zero_policy = config$zero_policy
)

## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page
## This method assumes the response is known - see manual page

all_coefficients <- extract_all_coefficients(
  ols_model = ols,
  model_results = model_results
)

write_csv(
  model_table,
  file.path(config$output_dir, "tables", "model_comparison.csv")
)
write_csv(
  all_coefficients,
  file.path(config$output_dir, "tables", "all_model_coefficients.csv")
)

model_errors <- collect_model_errors(model_results)
write_csv(
  model_errors,
  file.path(config$output_dir, "tables", "model_fit_errors.csv")
)

knitr::kable(
  head(model_table, 20),
  digits = 4,
  caption = "Top-ranked model specifications"
)

Top-ranked model specifications
model	weight_name	observations	rho	lambda	log_likelihood	AIC	BIC	RMSE	MAE	fitted_R2	residual_moran_I	residual_moran_p	elapsed_seconds	delta_AIC	residual_moran_significance
SAR	groundwater_variogram	1310	0.8937	NA	-3502.333	7016.665	7047.732	3.3308	2.4545	0.6659	0.0306	0.0017	21.37	0.0000	**
SDM_rainfall_lag	groundwater_variogram	1310	0.8937	NA	-3502.332	7018.664	7054.909	3.3308	2.4545	0.6659	0.0306	0.0017	22.82	1.9994	**
Hybrid_SAR_rainfall	hybrid_rainfall_only	1310	0.8936	NA	-3502.332	7018.665	7054.909	3.3309	2.4546	0.6659	0.0307	0.0017	19.89	1.9997	**
SEM	groundwater_variogram	1310	NA	0.9017	-3505.196	7022.391	7053.458	3.3322	2.4537	0.6681	0.0273	0.0050	14.47	5.7263	**
Hybrid_SDM_full	hybrid_full	1310	0.8937	NA	-3502.240	7022.481	7069.081	3.3306	2.4535	0.6660	0.0298	0.0022	21.54	5.8157	**
SDM_all_lags	groundwater_variogram	1310	0.8944	NA	-3502.260	7022.519	7069.119	3.3301	2.4531	0.6662	0.0297	0.0023	21.23	5.8541	**
SAR	d22km	1310	0.9073	NA	-3645.306	7302.611	7333.678	3.8101	2.8459	0.5568	0.0060	0.3276	23.33	285.9464
SEM	d22km	1310	NA	0.9173	-3645.571	7303.142	7334.209	3.8063	2.8356	0.5592	0.0034	0.5409	14.31	286.4771
SDM_rainfall_lag	d22km	1310	0.9089	NA	-3645.011	7304.022	7340.266	3.8085	2.8430	0.5573	0.0052	0.3889	22.00	287.3569
SDM_all_lags	d22km	1310	0.9115	NA	-3644.598	7307.196	7353.796	3.8061	2.8400	0.5581	0.0040	0.4933	22.66	290.5306
SEM	d27km	1310	NA	0.9271	-3699.926	7411.851	7442.918	3.9958	2.9648	0.5124	-0.0027	0.7330	14.38	395.1863
SAR	d27km	1310	0.9160	NA	-3699.942	7411.883	7442.950	4.0000	2.9754	0.5100	-0.0006	0.9769	20.47	395.2181
SDM_rainfall_lag	d27km	1310	0.9195	NA	-3699.237	7412.473	7448.718	3.9966	2.9698	0.5111	-0.0016	0.8770	21.21	395.8080
SDM_all_lags	d27km	1310	0.9218	NA	-3698.954	7415.909	7462.509	3.9949	2.9681	0.5117	-0.0026	0.7525	22.33	399.2437
SEM	d37km	1310	NA	0.9511	-3744.555	7501.109	7532.176	4.1636	3.1150	0.4688	-0.0043	0.4080	15.03	484.4441
SDM_rainfall_lag	d37km	1310	0.9448	NA	-3744.168	7502.335	7538.580	4.1642	3.1172	0.4677	-0.0037	0.4900	19.86	485.6702
SAR	d37km	1310	0.9364	NA	-3746.297	7504.594	7535.661	4.1732	3.1298	0.4649	-0.0004	0.9228	25.08	487.9289
SDM_all_lags	d37km	1310	0.9470	NA	-3743.947	7505.894	7552.494	4.1629	3.1162	0.4683	-0.0043	0.4072	22.37	489.2288
SEM	d47km	1310	NA	0.9799	-3777.264	7566.528	7597.594	4.2805	3.2306	0.4407	0.0126	0.0001	14.70	549.8626	***
SDM_rainfall_lag	d47km	1310	0.9785	NA	-3776.875	7567.751	7603.995	4.2798	3.2282	0.4405	0.0124	0.0001	20.04	551.0856	***

The comparison table reports likelihood, AIC, BIC, RMSE, MAE, fitted correlation-squared, the spatial parameter, residual Moran’s I, its p-value, and elapsed estimation time. The coefficient file contains every local and lagged coefficient from every fitted specification.

9 8. Formal SLX temperature tests

The formal temperature analysis is performed in the full rainfall-VIW SLX model. Classical and HC3-robust covariance estimates are used to test:

the local temperature coefficient;
the lagged temperature coefficient;
their joint significance;
the sum of the local and lagged coefficients.

temperature_tests <- run_slx_temperature_tests(
  full_slx_model = slx_full$models$rainfall_variogram,
  temperature_variable = config$temperature_variable,
  reduced_slx_model = slx_rain$models$rainfall_variogram,
  listw = model_weights$rainfall_variogram,
  zero_policy = config$zero_policy
)

export_temperature_tests(
  temperature_tests,
  file.path(config$output_dir, "tables")
)

knitr::kable(temperature_tests$local, digits = 4)

test	covariance	estimate	standard_error	statistic	degrees_freedom	p_value_two_sided	confidence_low	confidence_high
L0T2M.FULL	Classical	5.2603	1.4462	3.6373	1303	3e-04	2.4232	8.0975
L0T2M.FULL	HC3 robust	5.2603	1.5802	3.3288	1303	9e-04	2.1603	8.3604

knitr::kable(temperature_tests$lagged, digits = 4)

test	covariance	estimate	standard_error	statistic	degrees_freedom	p_value_two_sided	confidence_low	confidence_high
lag.L0T2M.FULL	Classical	-10.1211	3.3488	-3.0223	1303	0.0026	-16.6907	-3.5514
lag.L0T2M.FULL	HC3 robust	-10.1211	3.1847	-3.1781	1303	0.0015	-16.3687	-3.8735

knitr::kable(temperature_tests$joint, digits = 4)

test	covariance	restrictions	wald_chisq	chisq_df	chisq_p_value	F_statistic	numerator_df	denominator_df	F_p_value
L0T2M.FULL = 0 and lag.L0T2M.FULL = 0	Classical	2	13.8507	2	0.0010	6.9254	2	1303	0.0010
L0T2M.FULL = 0 and lag.L0T2M.FULL = 0	HC3 robust	2	13.2442	2	0.0013	6.6221	2	1303	0.0014

knitr::kable(temperature_tests$total, digits = 4)

test	covariance	estimate	standard_error	statistic	degrees_freedom	p_value_two_sided	confidence_low	confidence_high
1L0T2M.FULL + 1lag.L0T2M.FULL	Classical	-4.8607	2.6105	-1.8620	1303	0.0628	-9.9820	0.2605
1L0T2M.FULL + 1lag.L0T2M.FULL	HC3 robust	-4.8607	2.5654	-1.8948	1303	0.0583	-9.8934	0.1720

10 9. Figures and archival record

plots <- plot_model_comparison(model_table)
save_replication_plots(
  plots,
  file.path(config$output_dir, "figures")
)

## Warning: The shape palette can deal with a maximum of 6 discrete values because more
## than 6 becomes difficult to discriminate
## ℹ you have requested 8 values. Consider specifying shapes manually if you need
##   that many of them.

## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_point()`).

plots$delta_aic

plots$fit_vs_moran

## Warning: The shape palette can deal with a maximum of 6 discrete values because more
## than 6 becomes difficult to discriminate
## ℹ you have requested 8 values. Consider specifying shapes manually if you need
##   that many of them.
## Removed 16 rows containing missing values or values outside the scale range
## (`geom_point()`).

writeLines(
  capture.output(sessionInfo()),
  file.path(config$output_dir, "sessionInfo.txt")
)

saveRDS(
  list(
    config = config,
    samples = samples,
    variograms = variogram_results,
    kriging = kriging_results,
    synthetic_first = synthetic_first,
    experiment_i = experiment_i,
    experiment_i_weights = experiment_i_weights_object,
    model_weights = model_weights_object,
    models = list(
      ols = ols,
      sar = sar,
      sem = sem,
      slx_full = slx_full,
      slx_rain = slx_rain,
      sdm = sdm,
      sdm_rain = sdm_rain,
      hybrid_sdm_full = hybrid_sdm_full,
      hybrid_sar_rainfall = hybrid_sar_rainfall
    ),
    temperature_tests = temperature_tests,
    model_table = model_table,
    all_coefficients = all_coefficients,
    model_errors = model_errors
  ),
  file.path(config$output_dir, "replication_objects.rds")
)

11 Generate manuscript outputs

source(file.path("R", "11_manuscript_outputs.R"))

manuscript_outputs <- generate_manuscript_outputs(
  samples = samples,
  variogram_results = variogram_results,
  kriging_results = kriging_results,
  model_table = model_table,
  all_coefficients = all_coefficients,
  temperature_tests = temperature_tests,
  experiment_i = experiment_i,
  config = config,
  study_region = samples$study_area,
  prediction_grid = samples$full_grid_sf,

  # Figure 6 reuses the already-computed final universal-kriging surface,
  # ensuring the manuscript figure matches the exported prediction rasters.
  generate_figure_6 = TRUE,
  
  output_dir = file.path(
    config$output_dir,
    "manuscript"
  )
)

source(file.path("R", "12_model_specific_tables.R"))

model_specific_tables <- export_model_specific_tables(
  all_coefficients = all_coefficients,
  model_table = model_table,
  model_results = model_results,
  output_dir = file.path(
    config$output_dir,
    "manuscript",
    "tables",
    "models"
  )
)

12 Interpretation safeguards

Moran’s I is conditional on the selected weights matrix. The Monte Carlo experiment therefore compares weighting-induced variation rather than estimating a weights-free true value.
The rainfall variogram is used for rainfall lags, not for temperature or cropping lags.
In the full hybrid SDM, temperature and cropping lags use groundwater VIW; no temperature or cropping variograms are estimated.
A small residual Moran’s I can coexist with poor likelihood and predictive performance. Residual whitening should not be used as the sole model-selection criterion.
The local-plus-lagged SLX temperature sum is an association under the fitted spatial specification, not a causal effect.

Spatial Dependence and Groundwater Governance: A Geostatistical-Econometric Approach (Replication Document - R)

Saif Ali

2026-09-20