Running rhap: Calculate household air pollution (HAP) impacts

Source: vignettes/run_rhap.Rmd

The primary function provided by this package is calc_hap_impacts, which is designed to estimate health impacts attributable to household air pollution (HAP) under various alternative scenarios simulated using the Global Change Analysis Model (GCAM). This function streamlines the process of assessing scenario-specific health outcomes by extracting relevant data from GCAM databases or project files. These relevant data includes key socioeconomic and environmental parameters, such as per capita GDP, emissions of primary PM2.5 and NOx, and per capita floorspace. For more details on the methodology and assumptions underlying the econometric model used in this calculation, refer to the guide on “Fitting the Econometric Model”.

Extraction of the scenario-specific data

The function calc_hap_impacts connects to a GCAM databases or project file using the rgcam package, designed to extract and import GCAM results.

First, the function extracts direct air pollutant emissions from the residential sector for each simulated GCAM scenario, GCAM region, period, and pollutant. To estimate country-specific impacts, these emissions are processed and downscaled to the country level using a static downscaling approach. This approach leverages gas-specific country-level emissions data, specifically using 2017 emissions data from the Community Emissions Data System (CEDS), updated for the Global Burden of Disease - Major Air Pollution Sources projectCEDS-GBD. Additionally, the function applies transformations such as converting black carbon and organic carbon emissions into primary PM2.5 to ensure comprehensive pollutant accounting.

The function also extracts per capita floorspace, which is calculated endogenously in GCAM by region and period. This calculation relies on a Gompertz-type function calibrated with empirical data (Sampedro et al 2022), providing a income-driven estimate of future residential space per person.

For per capita GDP, the package loads projections from the SSP database version 3.0.1,updated in 2024. These GDP projections are scenario-specific, corresponding to the Shared Socioeconomic Pathways (SSPs). By default, the package utilizes SSP2 (“middle of the road”) projections. However, if the GCAM scenario name includes a different SSP designation (e.g., “Reference_SSP3”), the package automatically selects the appropriate SSP-specific GDP projections.

Calculations

The function estimates health impacts attributable to household air pollution (HAP) using an econometric model that incorporates direct emissions from the residential sector (primary PM2.5 and NOx), per capita GDP, per capita floorspace, and a linear year trend as covariates. The model is a fixed-effects model, estimated using cross-regional and multi-year panel data compiled from various different sources. Detailed information on the econometric model, including its formulation, variable selection, and diagnostic testing, is available in the dedicated vignette: “Fitting the Econometric Model”.

To improve accuracy, the health impact predictions incorporate a country-specific, additive bias-correction term, computed by fit_model() as the mean gap between observed and model-predicted values over the most recent years of the training panel, and added to every prediction. This corrects for the country fixed effect that a fitted plm model cannot otherwise apply to new (scenario) data. An additive correction was chosen over the more theoretically obvious multiplicative one after validation against real GBD data showed the multiplicative version performing far worse (see the “Prediction and the bias adder” section of the econometric-model vignette) – the trade-off is that a small floor is applied to prevent a negative corrected rate.

The package provides flexibility to calculate three distinct health impact metrics associated with HAP:

  • Premature Mortalities
  • Years of Life Lost (YLLs)
  • Disability-Adjusted Life Years (DALYs)

In addition to country-level estimates, setting by_gr = TRUE reports within-country income-decile group shares of population, PM2.5+NOx exposure, and GDP, as a CSV table and a pie-chart snapshot (output/by_gr/). This is descriptive only, not a group-level health impact estimate: the regression model is fit on between-country variation, and applying it to within-country decile-level inputs would extrapolate it far outside the variation it was estimated on, with no way to validate the result. For a more detailed, per-pollutant breakdown of group shares, see calc_ResidEm_grp().

Outputs

The function calculates health impacts—such as deaths, Years of Life Lost (YLLs), or Disability-Adjusted Life Years (DALYs)—attributable to household air pollution (HAP) for each country, time period, and GCAM scenario. By using the normalized parameter, users can choose between producing results in absolute terms or as normalized values (per 100,000 inhabitants), offering flexibility for different analytical needs.

If saveOutput is set to TRUE, the function writes the following csv files in the output/ sub-directory: [scenario]_HAP_HIAvar.csv

In addition, by setting map to TRUE, the function generates damage maps, using the rmap package documented in the following page. The function also generates animations and individual figures.

Step-by-Step Example

This example demonstrates how to use the rhap package to estimate health impacts from GCAM outputs.

1. Install the Package

Follow the installation guide if you haven’t already.

2. Load the Required Libraries

  • If you are developing the package (e.g., using RStudio or working within the package directory), use:
devtools::load_all()
  • Otherwise, if you’ve installed the package, load it normally:

3. Set Project and Scenario Names, and Run the Model 

# Set your project and scenario names
my_prj_name  <- rhap::test_rhap  # Name of the .dat project file (avoid spaces). Note that rhap::test_rhap loads the example project file included in the package

my_scen_name <- "Reference"       # Name of the GCAM scenario to process (or vector of scenarios)

# Run the RHAP model
hap_damages <- calc_hap_impacts(
  prj_name     = my_prj_name,
  scen_name    = my_scen_name,
  final_db_year = 2050,
  saveOutput   = TRUE,
  map          = FALSE,
  anim         = TRUE,
  HIA_var      = "deaths",
  normalized   = TRUE,
  by_gr        = FALSE
)

# Preview the results
head(hap_damages)

This code will generate a dataset of health impacts and, also produce a figure similar to the one shown below:

Premature deaths per 100.000 inhabitants attributable to household air pollution in 2050