Spatial INLA Modelling Workflow

Data Cleaning And Processing Overview

The modelling workflow in this folder starts from a mix of raw downloads in raw_data/ and a few already-curated inputs in clean_data/. The cleaned tables used by the INLA scripts were produced in two stages: first, each raw source was reduced to a planning-area-level or coordinate-level analysis file; second, prepare_inla_data.R merged those cleaned files into a single spatial modelling dataset.

Upstream cleaning scripts

• extract_planning_area_resident_count.py reads the June 2025 resident Excel file, keeps only planning-area totals where subzone, age group, and floor area are all reported as totals, and writes clean_data/planning_area_total_residents_2025.csv.
• extract_planning_area _worker_counts.py identifies planning areas whose 2025 resident count is missing, then attaches 2020 workplace counts for those areas and writes clean_data/planning_areas_worker_counts_2020.csv.
• extract_planning_area_visitor_counts.py cleans attraction names, joins them to a place-to-planning-area lookup table, sums annual visitor counts by planning area, divides those totals by 2 to approximate the second half of 2024, and writes clean_data/planning_area_visitors_2024_halved.csv.
• aggregate_weekly_rainfall.py reads the extracted rainfall CSV, groups rows by station coordinates, aggregates rainfall into consecutive 7-day bins, and writes a complete location-by-week table to clean_data/rainfall_weekly_by_location.csv, including zero-rainfall weeks when a station has no positive readings in a week.
• extract_flood_alert_locations.py reads the saved flood-alert JSON, keeps only records where msgType == "Alert", extracts the alert-circle coordinates, removes duplicates, and writes clean_data/flood_alert_locations.csv.
• filter_landuse_exposure_areas.py filters the Master Plan 2025 land-use layer to the exposure-related classes used in this analysis (BEACH AREA, WATERBODY, AGRICULTURE, PARK, OPEN SPACE, CEMETERY) and writes clean_data/HighExposureLandUse2025.geojson.

Not every modelling input is generated by a script in this folder. The case coordinates in clean_data/Case.txt are used as the observed outcome input, and clean_data/Non-resident_population.csv is a small curated fallback table used for planning areas that do not have resident counts but still need an exposure proxy.

Final processing before modelling

prepare_inla_data.R performs the last integration step before model fitting. It standardizes planning-area names, spatially joins case points and flood alerts to planning areas, averages weekly rainfall to a station-level mean and then to a planning-area mean, intersects the filtered land-use polygons with planning-area boundaries, and builds a single population_proxy using the resident, worker, non-resident, and visitor sources in that order. It also applies nearest-feature fallback rules when points or rainfall stations do not fall cleanly inside a planning area, computes expected counts E, and creates the adjacency graph used by the BYM2 model.

This folder contains R scripts for planning-area-level spatial modelling using INLA.

prepare_inla_data.R is the shared data-preparation module. It reads the cleaned inputs, aggregates cases to planning areas, builds the population proxy and covariates, computes the expected counts E, creates the adjacency graph, and returns the prepared objects used by the modelling scripts. It also contains the shared fit_poisson_inla() helper used to fit Poisson INLA models with the common settings for this project.

compare_inla_models.R is the model-comparison script. It sources prepare_inla_data.R, then fits a greedy sequence of candidate Poisson INLA models, starting from an intercept-only model, adding fixed covariates when WAIC improves, and adding BYM2 only when the spatial term further improves WAIC. It saves model-level comparison metrics to results/model_comparison.csv and combined fixed-effect summaries to results/model_fixed_effects_comparison.csv.

build_inla_model.R is the main analysis script. It sources prepare_inla_data.R, fits the Poisson spatial model with the best WAIC in results/model_comparison.csv, which is with rainfall, high-exposure land-use share, and a BYM2 random effect, then saves the detailed outputs, tables, and maps used for the main analysis.

Together, these files are used to:

• read the cleaned inputs from clean_data/
• aggregate case points to planning areas
• create planning-area covariates
• build a proxy population offset
• fit the full Poisson BYM2 model
• compare simpler candidate models against the full spatial specification
• save model outputs, tables, and maps into results/

Main Files

• prepare_inla_data.R: shared data preparation and common INLA fitting helper
• build_inla_model.R
• compare_inla_models.R

R Packages Used

The modelling script uses these R packages:

• sf
• dplyr
• spdep
• INLA
• ggplot2

If INLA is not installed yet, install it from the INLA repository before running the script.

Input Files

The workflow expects these files inside clean_data/:

• 2025PlanningAreaBoundary.geojson
• Case.txt
• rainfall_weekly_by_location.csv
• HighExposureLandUse2025.geojson
• flood_alert_locations.csv
• planning_area_total_residents_2025.csv
• planning_areas_worker_counts_2020.csv
• planning_area_visitors_2024_halved.csv
• Non-resident_population.csv

If any of these files are missing, the shared preparation step stops with an error before fitting any model.

What The Shared Workflow Does

1. Load spatial boundaries

The planning-area boundary file is read as an sf object and standardized to a single planning-area key using uppercase names. A projected version is also created for area-based calculations.

2. Aggregate cases to planning areas

The case coordinates in Case.txt are converted into points and spatially joined to planning areas.

If a case point falls on a boundary and is not matched by st_within, the script assigns it to the nearest planning area using st_nearest_feature.

The final response variable is:

• cases: number of observed case points per planning area

3. Build the population proxy

Because some planning areas do not have resident population counts, the script constructs a single population_proxy variable using this order:

1. total_residents
2. workplace_count_2020
3. non_resident_population
4. annual_visitors_2024_halved
5. 0

The script also records a proxy_source column so it is clear which source was used for each planning area.

Planning areas with:

• population_proxy = 0
• and cases = 0

are kept in the exported tables and maps, but they are excluded from the model fit itself because a zero exposure cannot be used as a stable Poisson offset.

If any planning area has:

• population_proxy = 0
• and cases > 0

the script stops and reports that problem.

4. Create covariates

The script derives three planning-area covariates.

Rainfall

rainfall_weekly_by_location.csv contains weekly rainfall totals by station location. The script:

• averages weekly rainfall over time for each station
• assigns stations to planning areas
• computes the mean rainfall for each planning area from stations inside it
• uses the nearest station as a fallback when an area has no station inside it

The final rainfall covariate is:

• rainfall_mm

High-exposure land-use share

HighExposureLandUse2025.geojson contains polygons representing filtered land uses such as waterbody, park, open space, agriculture, beach area, and cemetery.

The script intersects those polygons with planning areas and computes:

• high_exposure_area_m2
• high_exposure_share

where high_exposure_share is:

high_exposure_area_m2 / planning_area_area_m2

Flood-alert count

flood_alert_locations.csv is converted into spatial points and joined to planning areas.

The script counts how many flood-alert locations fall inside each planning area. If a point does not fall cleanly inside an area, it is assigned to the nearest planning area.

The final flood covariate is:

• flood_alert_count

5. Compute the expected count E

For planning areas included in the model, the script computes:

E = total_cases * population_proxy / sum(population_proxy)

This is the expected number of cases under a baseline where cases are distributed proportionally to the proxy population.

The model then estimates whether each planning area has higher or lower risk than expected after accounting for the covariates and spatial structure.

6. Fit the INLA model

The fitted model is the selected rainfall_landuse_bym2 Poisson BYM2 model:

cases ~ rainfall_z + high_exposure_share_z +
  f(region_id, model = "bym2", graph = graph)

Where:

• rainfall_z and high_exposure_share_z are standardized covariates
• region_id is the planning-area index used by INLA
• graph is the neighborhood graph built from shared planning-area borders

The response is cases, and the exposure term is E.

The script also requests:

• fitted values
• DIC
• WAIC
• CPO

7. Attach fitted quantities

After fitting, the script adds:

• fitted_mean
• fitted_lower
• fitted_upper
• relative_risk
• rr_lower
• rr_upper
• pp_rr_gt_1

relative_risk is computed as:

relative_risk = fitted_mean / E

pp_rr_gt_1 is the posterior probability that relative risk is greater than 1.

Output Files

When the script runs successfully, it writes these files into results/:

• model_data.csv
• fixed_effects.csv
• hyperparameters.csv
• fit_metrics.csv
• risk_ranking.csv
• inla_model.rds
• adj.graph
• planning_area_risk.geojson
• cases_map.png
• relative_risk_map.png
• posterior_prob_map.png

Output Description

model_data.csv

Planning-area table containing:

• case counts
• population proxy information
• covariates
• expected counts
• fitted summaries
• relative risk summaries

fixed_effects.csv

Posterior summaries of the fixed-effect coefficients.

hyperparameters.csv

Posterior summaries of the BYM2 hyperparameters.

fit_metrics.csv

Model fit statistics such as:

• DIC
• WAIC
• mean negative log CPO

risk_ranking.csv

Planning areas sorted from highest to lowest estimated relative risk.

planning_area_risk.geojson

Spatial output for mapping in GIS tools or other spatial software.

Map PNGs

• cases_map.png: observed case locations over planning-area boundaries
• relative_risk_map.png: modelled relative risk by planning area
• posterior_prob_map.png: posterior probability that relative risk exceeds 1

Important Assumptions

• The case file represents the observed outcome being modeled.
• Rainfall is treated as a spatial covariate by averaging the available weekly station records in the cleaned file.
• Flood-alert locations are treated as a static environmental covariate, not as the response.
• Planning areas with no positive proxy population and no cases remain in the exports, but not in the fitted model.

Notes

• The script is designed for the cleaned inputs currently present in this project.
• If you change column names or file names in clean_data/, the script will need to be updated.
• The script does not modify the input data; it only reads from clean_data/ and writes new outputs into results/.