Fitting the econometric model

Panel data

To build the econometric model to estimate health impacts attributable to household air pollution, we built a cross-regional and multi-year panel dataset, that gathers almost all the countries in the world for a 30-year period (1990–2019). The dependent variable is “health impacts” which could be measured as premature deaths, years of life lost, or disability-adjusted life years. In order to explain this variables we collect different key socieconomic and environmental variables.

First, the panel gathers data on different gases emitted in the residential sector, sourced from the Community Emissions Data System (CEDS) database. Specifically, the panel includes the following pollutants: black carbon (BC), methane (CH4), carbon monoxide (CO), nitrous oxide (N2O), ammonia (NH3), non-methane volatile organic compounds (NMVOC), nitrogen oxides (NOx), organic carbon (OC), and sulfur dioxide (SO2).

Another included variable is outdoor air pollution (AAP) for each country, measured in PM2.5. This is particularly relevant in urban areas, where it can influence HAP.

A range of socioeconomic variables are also incorporated to explore if/how greater development affects the health impacts. These include:

• GDP per capita in Purchasing Power Parity (PPP), sourced from the Our World in Data website.
• Urban and rural population sizes, obtained from the World Bank.
• Access to clean cooking and heating technologies and fuels, provided by the World Health Organization.
• Average household size per person, generated using data from various sources, such as national household surveys (e.g., the China Statistical Yearbook, the Residential Energy Consumption Survey [RECS], the India NSSO Household Consumer Expenditure Surveys, and the European Union Statistics on Income and Living Conditions [EU-SILC]). Note that for regions where detailed microdata is not directly accessible, global data sources such as Odyssee and the International Energy Agency (IEA) were used.

Finally, the panel includes a set of climatic variables, such as average maximum, minimum, and mean temperatures; average precipitation for each country per year (see the climate knowledge portal); and indicators for heating and cooling needs. These are:

• Heating Degree Days (HDD): calculated as the difference between a base temperature (typically 18C) and the average daily temperature.
• Cooling Degree Days (CDD): representing days when cooling is required, calculated similarly to HDD but using outdoor high temperatures instead of low ones. Data for these indicators was sourced from the Community Climate System Model.

Econometric regression

While the panel provides a wide range of variables for the regression model, not all are statistically significant. Moreover, some are highly correlated, making it essential to select covariates carefully to ensure independence and avoid multicollinearity. After evaluating multiple models with different variable combinations, the final regression includes the following socioeconomic and environmental variables: per capita GDP, emissions of primary PM2.5 (comprising BC and OC), nitrogen oxides, and per capita floor space, plus a linear year trend. Non-methane volatile organic compounds (NMVOCs) were tested and dropped: their coefficient is statistically indistinguishable from zero under every specification and standard-error estimator tried, and fails a joint significance test both alone and jointly with PM2.5.

Once the covariates are selected, the variables are first normalized by converting them to units per 100,000 inhabitants. This step ensures that the scales of the variables are consistent, making them more comparable across different regions or populations. Subsequently, a logarithmic transformation is applied to these normalized variables. The logarithmic scale helps linearize any nonlinear relationships, addresses issues related to skewness, and further stabilizes variance, ultimately improving the reliability and interpretability of the regression model.

To determine whether to use a fixed or random effects model, we apply the Hausman test. Since the null hypothesis is rejected at conventional significance levels, this suggests that unobservable individual-specific characteristics are correlated with the explanatory variables. As a result, we opt for a fixed effects model over a random effects model. Note that the estimation is performed using the plm package.

Time trend

Mortality attributable to HAP has an underlying downward trend over 1990-2019 that is largely common across countries (driven by factors outside the model, such as general improvements in medicine and public health infrastructure). Leaving this trend unmodeled risks confounding it with the genuinely declining trend in residential pollutant emissions over the same period. We compared four ways of controlling for it: no trend control, a linear year trend, a flexible spline (4 degrees of freedom), and full year fixed effects (year dummies). The linear trend, spline, and year-dummy specifications all converge to similar pollutant coefficients; only the no-trend-control specification disagrees, indicating that the “reasonable-looking” coefficients from that specification were actually the biased ones (omitted common-trend bias). Full year dummies over-absorb the common trend to the point of washing out the pollutant coefficients almost entirely. We use the linear trend: it is the most parsimonious option that agrees with the more flexible alternatives, and unlike year dummies it extrapolates sensibly to future (post-2019) scenario years.

Inference

Country-year panels of this kind routinely violate the classical OLS/within-estimator assumptions used for standard errors: residuals from the same country in adjacent years are correlated (serial correlation), and global shocks (oil prices, recessions, multilateral pollution treaties, pandemics) hit many countries in the same year (cross-sectional dependence). We tested for both directly: a Breusch-Godfrey/Wooldridge test rejects the no-serial-correlation null (p < 0.001), and a Pesaran CD test rejects the no-cross-sectional-dependence null (p < 0.001) – so both are actually present in this panel, not just theoretical concerns. Naive within-estimator standard errors, and even standard errors clustered only by country (Arellano), are not valid under cross-sectional dependence. We therefore report Driscoll-Kraay standard errors (plm::vcovSCC()), which are robust to serial correlation, cross-sectional dependence, and heteroskedasticity (also confirmed present via a Breusch-Pagan test) simultaneously.

In summary, we estimate health impacts attributable to household air pollution using the following fixed effects model, with country and year specified as index variables:

$$log(Deathsper100K)_{i,t} = \alpha_{i} + \beta_{1}log(PrimPM2.5per100K)_{i,t} + \beta_{2}log(NOXper100K)_{i,t} + \beta_{3}log(GDPpc)_{i,t} + \beta_{4}log(FLSPpc)_{i,t} + \gamma \cdot t + e_{i,t}$$ Where:

• $log(Deathsper100K)_{i,t}$ is the dependent variable for region $i$ at period $t$.
• $\alpha_{i}$ is the fixed effect for region $i$.
• $\beta_{1}$, … $\beta_{4}$ are the coefficients for the corresponding explanatory variables.
• $\gamma$ is the linear year-trend coefficient, with $t$ the calendar year.
• $e_{i,t}$ is the error term.

** Note that the dependent variable can also be set to be YLLs or DALYs

Prediction and the bias adder

predict() on a fitted plm “within” model returns fitted values from the coefficients alone, without the country fixed effect $\alpha_i$ – and this is unavoidable for scenario prediction, since GCAM scenario data has no fixed effect to draw on. To correct for this, fit_model() also computes a per-country additive bias adder: the mean of (observed $-$ coefficient-only-predicted) values, in per-100k units, over the most recent years of the training panel. This is added to (not multiplied into) every future per-100k prediction for that country.

This choice was not obvious going in. In a log-linear model, $log(Y) = X\beta + \alpha_i$ implies $Y = \exp(X\beta) \cdot \exp(\alpha_i)$ – the fixed effect is technically a multiplicative scalar in levels, not an additive one, which argues for a multiplicative (geometric-mean ratio) correction instead. We implemented that version and validated it against real GBD data before adopting it – and it performed far worse: cross-country correlation with observed 2023 GBD rates collapsed from ~0.86 (additive) to ~0.18 (multiplicative), even though the multiplicative version looked more accurate on the model’s own training data. The reason is that the model’s single common linear year trend doesn’t match every country’s own trajectory (mortality in a country like Uganda declines far more slowly than the trend assumes), so the gap between the trend-implied prediction and reality keeps growing even within the training years – for about 22% of countries the implied fixed effect exceeds a factor of 10x, and for some it exceeds 1000x. A multiplicative correction amplifies that already-large, still-growing gap catastrophically the moment it’s applied to a slightly different (scenario) year; a correction bounded by the actual scale of the outcome variable does not. We therefore use the additive form despite its weaker theoretical justification, because it is dramatically more reliable in the use this package is actually built for.

The one real failure mode of the additive form is that it can push a country’s corrected rate below zero (this happened for ~12% of country-year observations on the model’s own training data). Rather than clamping that to exactly 0 – which silently reports “zero risk” for a country the model has no real basis to call zero – calc_hap_impacts() floors the corrected value at a small fraction (1%) of the naive, uncorrected prediction, so a low but non-degenerate estimate is reported instead.

Known limitations

A RESET test on the pooled specification rejects correct functional form (p < 0.001), and individual country series for deaths, GDP, and floorspace show unit-root-like (non-stationary) behavior under an intercept-only panel unit-root test – both expected for 30-year macro/demographic series, but neither is fully resolved by entity fixed effects plus a single common linear trend if countries have meaningfully different trends of their own. This is a standard trade-off in applied fixed-effects panel work with many countries and a short time dimension (T=30), and the model’s cross-country ranking has been validated against real GBD data (~0.86 correlation on out-of-sample deaths/YLL/DALY rates); but formal cointegration testing or country-specific trends were not pursued and would be a reasonable next step for anyone using the model’s point coefficients for causal interpretation rather than for scenario comparison.