Replication Package - Does Cumulative Thermal Discomfort Increase Heating Demand?

Replication Package Does Cumulative Thermal Discomfort Reshape Heating Priorities? Evidence of Constrained Resignation among French Households Elisabeth Bourgeois (2026) Contents This repository contains the R scripts and documentation for all analyses reported in the article. Household survey data (PEPSI) are not included due to confidentiality; available upon reasonable request from the author. Analysis Scripts Script Description SCRIPT_MNL_baseline.R Multinomial logit estimation, marginal effects, benchmark comparison across 11 climatic indices SCRIPT_Brant_OL.R Brant test, ordered logit estimation SCRIPT_GOL_robustesse.R Generalized ordered logit (fully unconstrained) and partial proportional odds via VGAM (Williams, 2006; Yee, 2010) SCRIPT_heterogeneity_ztest.R Sub-sample stratification by regional TDI, z-test of coefficient differences (Paternoster et al., 1998) SCRIPT_spatial.R Moran's I, SAR, GWR, LISA diagnostics (regional NUTS-2 level) SCRIPT_LISA_departemental.R Departmental (NUTS-3) local Moran's I robustness on the logit residuals, aggregated to the 92 metropolitan departments under Queen contiguity, with Benjamini-Hochberg correction across the simultaneous tests Data Engineering Pipeline Script Description traitement_lot.R Hourly NetCDF extraction with spatial interpolation traitement_journalier_lot.R Daily aggregation from hourly CSV run_traitement.bat Batch launcher for serial processing Departmental LISA robustness The residual LISA reported at the regional level (12 NUTS-2 units) is re-estimated at the finer departmental resolution (92 metropolitan departments, Queen contiguity) on the residuals of the binary logit of stated heating priority on the cumulative TDI. Six departments reach nominal significance at p <= 0.05, a count close to the number expected by chance across 92 simultaneous tests, and none survives a Benjamini-Hochberg correction (lowest adjusted p = 0.165). The global departmental Moran's I is 0.0475 (p = 0.186), non-significant. The dispersion of constrained resignation therefore holds at fine spatial resolution. The script writes the classified map (fig_ch2_lisa_departemental.png and .pdf) and the table of local indicators with raw and corrected p-values per department (LISA_departemental_Ii_pvalues.csv). The departmental boundaries are read from the public france-geojson repository (G. David, https://github.com/gregoiredavid/france-geojson). Data Engineering Pipeline Overview Respondent postal codes were matched to the centroid of their corresponding INSEE commune. Hourly ERA5-Land data (0.1 x 0.1 degree grid, approx. 9 km resolution) were retrieved via API for the period October-March, 1991-2020, covering five meteorological variables. Data were extracted per respondent location using spatial interpolation, then aggregated to daily values in batches of 100 respondents (39 lots). Data source: Copernicus Climate Data Store, dataset reanalysis-era5-land (https://cds.climate.copernicus.eu). Registration required. Variables Extracted Variable Description t2m Temperature at 2m (K) d2m Dewpoint temperature at 2m (K) u10, v10 Horizontal wind components at 10m (m/s) ssrd Surface Solar Radiation Downwards (J/m2) Spatial Interpolation Hierarchy For each respondent, hourly values were extracted from the ERA5-Land grid cell containing their commune centroid, with interpolation from adjacent cells when available. The following hierarchy was applied: Bilinear — All four surrounding grid corners valid (standard method) Linear — Only two corners valid along one axis (unidirectional interpolation) None — Single valid corner only; raw cell value used without interpolation Substitution — No surrounding corners available; nearest valid land cell used Traceability files are generated per lot: interpolation_type_lotXX.csv — interpolation method applied per respondent ID substitution_cell_lotXX.csv — original and substituted cell coordinates (for coastal respondents) Note on grid orientation: ERA5-Land latitude is stored in decreasing order (north to south). The extraction matches grid cells to commune centroids by coordinate value, not by array index, so the ordering does not affect the assignment. Daily Aggregation For each respondent and each day, the following daily aggregates are computed from hourly values: Mean, minimum, and maximum temperature (converted from K to degrees C) Mean wind speed: sqrt(u10 squared + v10 squared), converted to km/h Mean relative humidity (derived from t2m and d2m) Solar radiation at 23h UTC (ssrd) Directory Structure traitement_clim/ lot/ # Coordinates per lot (coordonnees_lot01.csv ... lot39.csv) netcdf_files/ # Hourly ERA5-Land NetCDF files (.nc) export_traitement/ # Hourly CSV extractions per lot export_lotXX/ daily_traitement/ # Daily aggregated CSV files (1 per lot) traitement_lot.R # Hourly NetCDF extraction with interpolation traitement_journalier_lot.R # Daily aggregation from hourly CSV run_traitement.bat # Batch launcher for serial processing Execution Step 1: Hourly extraction (NetCDF to hourly CSV) run_traitement.bat Step 2: Daily aggregation (hourly CSV to daily aggregates) Rscript traitement_journalier_lot.R 01 Requirements R (version 4.3+) with packages: ncdf4, data.table, dplyr, lubridate, nnet, mlogit, spdep, sf, ggplot2, VGAM Rscript must be in the system PATH Sufficient disk space for NetCDF files (approx. 45 GB total) Notes Coastal respondents may require cell substitution due to land-sea grid boundaries. Daily CSV files contain up to 100 respondents per lot (approx. 450,000 rows each). The full extraction pipeline processes approximately 5 trillion hourly records from binary NetCDF files, reduced to 39 daily aggregate files suitable for analysis on standard hardware. Version notes This version adds the departmental (NUTS-3) LISA robustness check (SCRIPT_LISA_departemental.R) with its map and per-department table of local indicators, and updates the article title to its accepted form. The input microdata remain excluded and available on request.