Child Survival and Early Lifetime Exposures to Ambient Fine Particulate Matter in India: A Retrospective Cohort Study.

BACKGROUND: Ambient fine particulate matter [PM ≤2.5μm in aerodynamic diameter (PM2.5)] is a major health risk for children, particularly in South Asia, which currently experiences the highest PM2.5 levels globally. Nevertheless, there is comparatively little epidemiological evidence from this region to quantify the effects of PM2.5 on child survival.
OBJECTIVES: We estimated the association between PM2.5 exposure and child survival in India.
METHODS: We constructed a large, retrospective, and nationally representative cohort of children <5 years of age, born between 2009-2016, from the publicly available, cross-sectional 2015-2016 Demographic Health Surveys in India. In utero and post-delivery lifetime average ambient PM2.5 exposures were estimated with data from satellite remote sensing, meteorology, and land use information (model R2= 0.82). We used Cox proportional hazards regression to estimate the association between both average in utero and post-delivery lifetime PM2.5 and all-cause child mortality, controlling for individual- and household-level covariates, seasonality, location, and meteorology.
RESULTS: Over 7,447,724 child-months of follow-up, there were 11,559 deaths at <5 years of age reported by the children's mothers. The mean concentrations of 9-month in utero and post-delivery lifetime average ambient PM2.5 exposure were 71.1 μg/m3 (range: 20.9-153.5 μg/m3) and 73.7 μg/m3 (range: 14.0-247.3 μg/m3), respectively. Estimated child mortality adjusted hazard ratios were 1.023 [95% confidence interval (CI): 1.008, 1.038] and 1.013 (95% CI: 1.001, 1.026) per 10-μg/m3 increase of in utero and post-delivery lifetime PM2.5, with both exposures in the model. DISCUSSION: This study adds to the growing body of evidence about the adverse health effects of PM2.5 by demonstrating the association between exposure, both in utero and post-delivery, on child survival at the national level in India. Strategies to reduce ambient air pollution levels, including steps to minimize in utero and early life exposures, are urgently needed in India and other countries where exposures are above recommended guideline values. https://doi.org/10.1289/EHP8910.

Introduction

Fine particulate matter [airborne particles with an aerodynamic diameter of ()] has been shown to be associated with adverse health effects in children, including preterm birth (Balakrishnan et al. 2018; Huynh et al. 2006; Li et al. 2018; Stieb et al. 2012), low birth weight (LBW) (Bell et al. 2007; Pedersen et al. 2013), low height-for-age -score (Spears et al. 2019), respiratory infections (Gurley et al. 2013; Smith et al. 2011), and mortality (Heft-Neal et al. 2018; Son et al. 2017; Woodruff et al. 2006). In India, ambient levels are among the highest in the world, with an estimated population-weighted annual average (India State-Level Disease Burden Initiative Air Pollution Collaborators 2021) in 2019. The Global Burden of Disease (GBD) Project identified ambient as the second largest risk factor for disease burden in India (GBD 2019 Diseases and Injuries Collaborators 2020) and estimated that premature deaths every year are associated with ambient air pollution (India State-Level Disease Burden Initiative Air Pollution Collaborators 2021), of which are premature deaths in children at (India State-Level Disease Burden Initiative Air Pollution Collaborators 2019, 2021). To date, however, most evidence on child mortality due to exposure to ambient air pollution during pregnancy and early childhood was from regions with lower ambient levels (Glinianaia et al. 2004; Hans et al. 2011; Son et al. 2017; Woodruff et al. 1997, 2006). Although prospective pregnant mother–child cohort studies have recently been initiated in India to examine the association between air pollution exposure and indicators of child health, such as birthweight (Balakrishnan et al. 2015, 2018; Clasen et al. 2020), to our knowledge, estimates of child mortality have not been reported. Furthermore, owing to the low density of ambient air pollution monitoring stations in India and other low- and middle-income countries (LMICs) (Brauer et al. 2019; Landrigan et al. 2018), current studies of air pollution and infant mortality in India and LMICs are mainly based on annual average ambient prediction models (Shaddick et al. 2018; van Donkelaar et al. 2016). Such exposure estimates may be prone to error or misclassification, especially for outcomes such as child mortality, which has been associated with monthly or daily changes in air pollution levels (Hans et al. 2011; Son et al. 2017). The recent expansion in the number of ground monitoring stations in India and the development of spatial prediction models enable more spatiotemporally resolved estimation of ambient exposure throughout India.

The primary objective of our study was to evaluate the association between child mortality and ambient exposure during in utero and post-delivery lifetime periods among a nationally representative, retrospective cohort based on maternal-reported data in the Demographic and Health Surveys [DHS, also known as the National Family Health Survey 2015–2016 (NFHS-4)] in India. To estimate time-dependent ambient exposure for the cohort, we developed a monthly ambient model with 0.1° resolution between 2009–2018 over the entire Indian subcontinent. Our ground station–validated model was developed using machine learning, with multiple inputs from satellite remote sensing, meteorological data, and land use information.

Methods

Study Design and Population

We conducted a retrospective cohort study of a nationally representative sample of children 0–59 months of age born throughout India (as of 2019, 29 states and 5 Union Territories, excluding the Andaman and Nicobar Islands and Lakshadweep Union Territories), whose mothers participated in the 2015–2016 DHS/NFHS-4 between January 2015 and November 2016. The design of the DHS/NFHS-4 has been described elsewhere (IIPS and ICF 2017). Briefly, DHS/NFHS-4 is cross-sectional survey using a two-stage stratified sampling framework to collect information from a nationally representative sample of households with women 15–49 years of age and all of their children. Each respondent answered survey questions separately for all children no matter whether the child was still alive at the time of interview. We included children born to women no more than 5 y preceding the survey. Three tiers of all-cause mortality—including neonatal mortality (death at of age), infant mortality (death at of age), and child mortality (death at of age)—were outcome variables in our study. Child birth year and month, child death year and month (if the child died before interview), and potential confounders were reported by mothers during the DHS/NFHS-4 interview during the 2015–2016 period. The child-, maternal-, and household-level covariates were also reported by mothers at the end of the follow-up period. We calculated follow-up time in months, starting from child birth until child death as the failure event and the month of the DHS/NFHS-4 survey or passing the at-risk age as the censoring event.

Because this study analyzed only publicly available data without personal identifiers, the institutional review board (IRB) of Emory University determined this study not to be human subject research and did not require IRB review. Data analyses were conducted between March and September 2019.

Ambient Exposure and Meteorological Variables

In this study, we developed a random forest (RF) model to predict monthly ambient concentrations at a 0.1° spatial resolution ( at the equator) over the entire Indian subcontinent, from 2009 to 2018 (10 y). This approach has been shown to have reliable historical prediction capabilities, to allow a variety of input variables, and to have good agreement with ground measurements from previous studies for modeling ambient concentrations in East Asia and Latin America (Vu et al. 2019; Xiao et al. 2018). Other advantages of RFs are that these models can work with correlated predictors, can provide variable importance vectors to compare the influence of variables, and are not likely to overfit (Xiao et al. 2018).

First, we collected daily ambient concentration between 2015 and 2018 from air monitoring stations across India from the Central Pollution Control Board (https://app.cpcbccr.com/ccr/#/caaqm-dashboard-all/caaqm-landing) and hourly data at the U.S. Embassy and Consulates in India from AirNow (https://airnow.gov/). We excluded daily concentrations that were outside of 3 standard deviations of the mean for the log-transformed data (0.05%; exclusions were values outside of ) because we believed these measurements may not be plausible. For hourly concentration data from U.S. embassies or consulates, we removed data that lacked valid quality checks and daily average concentrations with at hourly measurements (2.2%).

Daily average measurements from stations within the same grid were averaged, resulting in 56,834 grid-day measurements (Table S1). We averaged gridded daily measurements into gridded monthly averages if more than 15 d had valid daily measurements for a given grid cell in a given month. We ended up with 1,446 grid-months of measurements during 2017–2018 from 134 ground air pollution monitoring stations for model development and 465 grid-months of measurements during 2015–2016 from 57 stations for model hindcasting evaluation.

We used multiple satellite-retrieved aerosol, nitrogen dioxide (), fire spot data, and vegetation index products; global meteorological and aerosol reanalysis models; land use information; and population density data for the period between 2009 and 2018 as inputs for our model. Briefly, we used Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 6 level-2 aerosol optical depth (AOD) products (MOD04_L2 and MYD04_L2) at a 10-km resolution collected by the Aqua and Terra satellites, from the Distributed Active Archive Center (https://ladsweb.modaps.eosdis.nasa.gov/) (Levy and Hsu 2015). We used Deep Blue (DB) and Deep Blue/Dark Target combined parameters (Combined) from MODIS retrials and assigned and averaged the centroid of each retrial to create grid cells and calculate monthly average AOD value for DB and Combined algorithms. We then conducted gap-filling of the missing AOD values for the period similar to previous study (Xiao et al. 2018). Table S2 summarizes the AOD data missingness patterns, and Figure S1 shows the Combined AOD parameter and DB AOD parameter spatial distribution after gap-filling. The aerosol absorbing index in visible and ultraviolet light and tropospheric density data were obtained from the Ozone Monitoring Instrument, downloaded from the Goddard Earth Sciences Data and Information Service Center (https://mirador.gsfc.nasa.gov/). We also collected active fire data obtained from the Information for Resource Management System (https://earthdata.nasa.gov/earth-observation-data/near-real-time/firms) and the Normalized Difference Vegetation Index from the MODIS Vegetation Indices in the Distributed Active Archive Center (https://ladsweb.modaps.eosdis.nasa.gov/).

Daily meteorological data and aerosol diagnostics were obtained from the Goddard Earth Observing System Data Assimilation System 5/Modern-Era Retrospective analysis for Research and Applications (MERRA-2; https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/) (Bosilovich et al. 2015). Elevation data was collected from the National Aeronautics and Space Administration’s Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (version 2). Highway, primary road, and other road density parameters were obtained from the Global Roads Inventory Project (Meijer et al. 2018). Yearly varying population numbers for 2009–2018 were obtained from the LandScan Global Population Database. All of the above variables were assigned into grid cells, and their monthly averages were used as predictors in the model. Details of these variables are listed in Table S3.

We trained the RF model on gridded monthly mean concentrations from the 2017–2018 period () and evaluated model performance using standard, spatial, and temporal 10-fold cross-validation (CV). Spatial CV refers to using data from 90% of grids for developing the model and then testing the model on the remaining 10% of grids; temporal CV relies on using 90% of monthly data to develop the model and then testing the model based on the remaining 10% of monthly data. We also evaluated the hindcasting performance of the model using monthly data from the 2015–2016 period (). Root-mean-squared error (RMSE) and values from 10-fold cross-validation and from hindcasting were reported to evaluate the historical prediction performance of the model. We then estimated population-weighted annual mean concentrations over India during 2009–2018, evaluated the change of levels over those 10 y, and compared our estimates with other published databases (India State-Level Disease Burden Initiative Air Pollution Collaborators 2019; Shaddick et al. 2018; van Donkelaar et al. 2016). Predicted ambient monthly levels were matched to each child with geo-coordinates at the cluster level [a primary sampling unit (PSU) or a segment of a PSU, of around 100–150 households] provided by DHS/NFHS-4. For each child, we assigned the 9-month average concentrations prior to birth as the in utero exposure. We treated post-delivery average exposures as time dependent and defined them as the average concentration from the calendar month of child birth through each month the child was at risk, until the month of death or censoring.

Meteorological variables including monthly air temperatures at and monthly rainfall precipitation were included as well, and they were also matched with the geo-coordinates of households (in the same way as performed for exposure). Monthly air temperatures at values were extracted from MERRA-2 data (Bosilovich et al. 2015), and monthly rainfall precipitation in mm/month were extracted from Climate Hazards Group InfraRed Precipitation with Station data (Funk et al. 2015). R (version 3.5; R Development Core Team) with the randomForestSRC (version 4.6) package was used for modeling, prediction, and assignment of and meteorological variables for each child.

Statistical Analysis

We used Cox proportional hazards regression models with follow-up in calendar months as the time variable and child death as the outcome. We included 9-month in utero and time-dependent, post-delivery average ambient exposure to estimate adjusted hazard ratios (HRs) relating child mortality risk with increases of exposure. In this model, the estimated effects of exposure during in utero and post-delivery periods were adjusted for each other. We also used single-exposure models with exposures from either the in utero or post-delivery periods for comparison. We used a set of confounders from the previous literature and a directed acyclic graph (Figure S2) to identify covariates to be included in the Cox regression model (Heft-Neal et al. 2018; Spears et al. 2019; Subramanian et al. 2009). The prespecified time-independent covariates included child sex, birth month and year, birth order, and location of birth (institutional birth or not); multibirths; and mother’s age at child birth, height, marital status, education (whether above secondary level), maternal smoking (yes/no), secondhand smoke exposure in the home (yes/no), wealth index, urban or rural location of households, geographical region/zonal council (Northern, Central, North Eastern, Eastern, Western, and Southern, based on https://www.mha.gov.in/zonal-council), primary cooking fuel, and toilet facilities. The household wealth index was calculated and provided by DHS/NFHS-4 using principle component analysis of household ownership of assets (Staveteig and Mallick 2014). We dichotomized toilet facilities as whether they were safely managed based on the World Health Organization (WHO) and United Nations Children’s Fund Joint Monitoring Program (Croft et al. 2018; WHO/UNICEF JMP for Water Supply, Sanitation and Hygiene 2018), and dichotomized cooking fuel as clean (electricity, liquefied petroleum gas, natural gas, and biogas) or dirty (kerosene, coal, wood, straw, agriculture crop residue, and dung). Time-dependent covariates included monthly rainfall precipitation and monthly temperature. Child birth year and month, sex, and the geographical region of households were modeled as strata. Unique baseline hazard functions were given for all 927 strata of combinations of these strata. Robust variance estimation for coefficients was used. For tied outcomes, Breslow’s method was used. Table S4 provides the details on the variables used in the Cox regression models. Although we did not include children born with LBW because it has been found to be a potential intermediate variable between air pollution and childhood death (Hernández-Díaz et al. 2006), we conducted sensitivity analysis adjusting for LBW in our study population, as described below. We also tested the association between LBW and in utero exposure, and between LBW and child mortality. We stratified the Cox regression by child birth month and year, child sex, and geographical regions (927 strata) to account for spatiotemporal and sex differences in child survival. In addition, we included the sample probability weight from DHS/NFHS-4 in the Cox regression (Equation 1).

is the hazard function (probability of death) at month t, and refers to the (1–927) strata of child birth month, sex, and geographical location. is the vector containing time-independent variables, including in utero exposure, and are the time-dependent variables including post-delivery exposure. The association between and child mortality was reported as the HR per increase of , and we also calculated the HR per interquartile range (IQR) increase of () to compare other studies. To evaluate whether effects of differed across covariates, we conducted subgroup stratified analysis of the two-exposure models. The proportional hazards assumption for in utero and post-delivery average exposure was tested on the basis of Schoenfeld residuals on time, and we found it did not violate the proportional hazards assumption (). The correlation matrix of the predictor estimation from the Cox regression model was calculated to check multicollinearity, and the Akaike information criterion (AIC) was calculated to compare single- and two-exposure models. We conducted subgroup analysis by conducting separate models by stratifying variables including sex, LBW, delivery mode, whether multibirth, maternal age, maternal height, education, maternal smoking, secondhand smoking, household wealth index, cook fuel type, toilet type, urbanity, marital status, and geographical location, and we tested the hypothesis that there were no interaction effects between in utero and post-delivery effects with each of the stratifying variables.

To examine the robustness of our results, we conducted additional sensitivity analyses. We added LBW back to the model (modify 1); conducted analysis using stratified variables as time-independent covariates (modify 2); adding additional cluster-level covariates, including the percentage of household with improved toilets and clean cooking fuels (modify 3); performing unweighted analysis without sample weight (modify 4); using generalized estimating equations to account for the additional correlation within household (modify 5) and cluster level (modify 6); replacing post-delivery lifetime exposure with monthly exposure (modify 7); and performing crude analysis without adjusting covariates (modify 8). In addition, we also conducted an analysis treating in utero and post-delivery as categorical variables using quartiles of exposure in the model and tested the linear trend of the categorized exposure.

Last, we calculated the population attributable fraction (PAF) as the burden of child mortality from ambient exposure. We estimated the PAF from the modeled in utero exposure distribution by every 10th percentile, and the risk of mortality (derived from the single-exposure model only) at a given exposure level compared baseline exposure levels. Stata (version 14; StataCorp) and R (version 3.6; R Development Core Team) were used in the statistical analyses.

Results

The RF prediction model characterized monthly levels well, with CV 0.82 at the monthly level (RMSE of ) in 2017–2018. The spatial and temporal CV (RMSE) of the RF model was 0.77 () and 0.77 () respectively, indicating a stable and good fit of monthly ambient measurements from ground stations. Compared with historical measurements in 2015–2016, our model provided good predictions of monthly with (RMSE) of 0.82 (). Figure S4 shows the performance of the CV and hindcasting of the RF model. Figure S3 shows changes in model performance across different numbers of decision trees and the variable importance matrix for the model. The historical prediction shows that our model underestimated levels during periods with extreme levels, when monthly levels were . Regardless, our model still had relatively good hindcasting capacity to predict historical monthly levels. The annual within-sample RMSE from the model CV in 2017–2018 was , and the annual RMSE from historical hindcasting in 2015–2016 was , with better performance compared with previous ambient air pollution models used in the GBD study in the South Asia Region (Shaddick et al. 2018), which had within-sample RMSE of .

Figure 1 shows the predicted annual in India in 2018 and the annualized change of over a 10-y period. The predicted annual population-weighted mean in India in 2018 was . The highest annual concentrations () were in the Indo-Gangetic Plain Region, covering the states of Haryana, Uttar Pradesh, and Bihar and the National Capital Territory of Delhi. These values are three times , the recommended limit set by the National Ambient Air Quality Standards of India (Gautam 2009). We observed increased ambient concentrations over a 10-y period in most areas, with the strongest increase of ambient concentrations at per year in the Indo-Gangetic Plain Region and along the west coast of the Western Region, whereas the Northeastern Region and the north part of the Northern Region showed decreases in ambient levels. Details of predicted in India (Figures S5 and S6), major metropolitan areas (Delhi, Mumbai, Kolkata, Chennai, and Bangalore) (Table S5), and population-weighted exposures from our model vs. other global models are listed in Table S6. The predicted ambient air pollution can be downloaded at https://doi.org/10.15139/S3/RFILYH.

Figure 1.

Spatial patterns of annual concentrations. (A) Model-predicted annual concentrations with measurements from 133 ground-based monitoring stations in subcontinental India in 2018. (B) Annualized change rate of concentrations between 2009–2018. Circles represent ground station with an annual mean measurement of . Note: , particulate matter in aerodynamic diameter (fine particulate matter).

The present study included data from 259,627 live-birth children born 5 y preceding the survey from 175,865 women in the DHS/NFHS-4 data. We excluded 6,839 (2.6%) children in the analysis due to a) living in the Andaman and Nicobar Islands and Lakshadweep Union Territories, which are islands more than from the India subcontinent (, 0.9%); b) not living with their mothers (, 0.7%); or c) missing variables (, 1.1%). A total of 252,788 children born between January 2010 and November 2016 were included from 27,853 clusters (Table 1). In the follow-up period of 7,447,724 child-months, 11,559 deaths were reported in children by their mothers. Among these deaths, 7,520 deaths occurred before the first month of life, and 10,862 deaths occurred during the first year of life (including 7,520 neonatal deaths). For the entire cohort, the mean (SD) of 9-month in utero and post-delivery lifetime average exposure until failure or censoring was 71.1 (28.2) and 73.3 , respectively. The characteristics of the cohort are shown in Table 1, stratified by 9-month in utero average exposure, at levels , 49.7 to , 62.4 to , and . In general, compared with children in the lowest quartile, children in the highest quartile of exposure were more likely to reside in rural areas (80.1% vs. 70.8%), were born to mothers with shorter stature (47.4% vs. 34.2% of mothers were in height), lower educational levels (44.3% vs. 19.4% were below the secondary education level), and lived in households with a lower wealth index (35.6% vs. 15.6% were in the lowest quartile of the wealth index). Children residing in central (Chhattisgarh, Madhya Pradesh, Uttarakhand, and Uttar Pradesh) and eastern (Bihar, Jharkhand, Odisha, and West Bengal) India states had higher estimated ambient levels compared with other states. Most children included in the analysis were born between 2012 and 2014, with the least number of children born in 2016 (Table 2). The correlation between average in utero and post-delivery lifetime exposure for the same child was 0.74.

Table 1

Unweighted characteristics of children years of age by four quartiles of 9-month in utero exposure who were born between 2010–2016 in the Demographic and Health Surveys [DHS, also known as the National Family Health Survey 2015–2016 (NFHS-4)] in India.

Characteristics	Entire cohort	Average 9-month in utero PM_2.5 exposure quartile (μg/m^3) (Q)
		Q1 (20.9–<49.7)	Q2 (49.7–<62.4)	Q3 (62.4–<90.3)	Q4 (90.3–153.5)
Average 9-month in utero PM_2.5 exposure (μg/m^3)
Children (n)	252,788	63,197	63,197	63,197	63,197
Total follow-up in month (n)	7,447,724	1,860,527	1,891,613	1,873,198	1,822,386
Child death at <5 years of age (n)	11,559	2,187	2,688	3,001	3,683
All-cause child mortality rate (n/1,000 child-years)	18.62	14.11	17.05	19.22	24.25
Average 9-month in utero PM_2.5 {μg/m^3 [mean (SD)]}	71.1 (28.2)	41.4 (6.0)	55.4 (3.6)	75.0 (8.1)	112.6 (13.9)
Average post-delivery lifetime PM_2.5 exposure until death or censoring {μg/m^3 [mean (SD)]}	73.7 (29.8)	46.6 (10.6)	57.1 (11.5)	83.6 (25.6)	107.6 (20.1)
Child covariatesa
Sex (%)
Female	47.9	48.4	47.8	47.6	47.3
Male	52.1	51.6	52.1	52.4	52.6
LBW (%)b	13.1	14.2	13.4	13.9	11.2
Birth order (%)
First	36.9	40.8	38.0	37.1	31.5
Second	30.7	34.2	31.1	30.0	27.5
Third	16.1	14.2	15.7	16.4	18.1
≥Fourth	16.3	10.8	15.2	16.5	22.9
Multiple birth (%)	1.65	1.55	1.53	1.81	1.71
Birth at institution (%)	75.3	81.8	75.7	74.8	68.7
Maternal covariatesa
Age at child birth [y (%)]
<19	12.2	12.6	11.3	12.9	11.8
20–24	43.4	43.4	42.2	44.5	43.4
25–29	28.4	28.6	29.2	27.8	28.3
≥30	16.0	15.4	17.5	14.8	16.5
Mothers height [cm (%)]
<149.9	38.3	34.2	34.0	37.5	47.4
150–154.9	33.8	35.3	34.6	33.3	32.2
155–159.9	19.7	21.5	21.9	20.4	15.0
>160	8.2	9.0	9.5	8.8	5.4
Marital status (%)
Not married	1.4	1.9	1.8	1.1	0.9
Married	98.6	98.1	98.2	98.9	99.1
Education (%)
Below secondary level	31.5	19.4	27.3	34.8	44.3
Above secondary level	68.5	81.6	72.7	66.2	55.7
Current smoker (%)	10.0	13.4	14.0	7.5	5.2
Secondhand smoke exposure (%)	56.6	51.8	60.5	56.5	58.3
Household covariatesa
Cook with clean fuel (%)	27.8	38.5	28.4	23.1	21.1
Access to improved toilet (%)	41.5	48.0	44.7	39.3	33.9
Urban/Rural (%)
Urban area	23.6	29.2	23.3	22.1	19.9
Rural area	76.4	70.8	76.7	77.9	80.1
Zone area (%)
Northern	16.6	3.5	25.7	29.4	7.9
North Eastern	13.8	20.9	23.9	9.8	0.5
Central	31.8	18.6	31.2	27.2	50.3
Eastern	21.1	5.0	9.6	28.5	41.3
Western	7.1	14.5	8.7	5.1	0
Southern	9.6	37.5	0.9	0	0
Wealth quartiles (%)c
First (lowest)	25	15.6	21.7	27.1	35.6
Second	25	23.4	26.3	24.5	25.9
Third	25	31.3	26.4	22.4	19.8
Forth	25	29.7	25.6	26.0	18.7
Meteorological covariatesd
Temperature {°C [mean (SD)]}	27.1 (6.4)	26.7 (4.9)	24.7 (7.8)	27.6 (6.0)	29.4 (5.5)
Total rainfall {cm [mean (SD)]}	103 (144)	115 (144)	98 (134)	90 (145)	110 (151)

Note: LBW, low birth weight; , particulate matter in aerodynamic diameter (fine particulate matter); SD, standard deviation.

Child, maternal and household-level covariates collected at the end of the follow-up during the DHS/NFHS-4 interview between 2015–2016.

Assuming children with missing LBW values are not LBW infants.

The household wealth index developed by the DHS/NHFS-4 was calculated using principle component analysis of household ownership of assets.

Meteorological covariates were averaged post-delivery monthly estimates.

Table 2

Number of children and clusters and their estimated annual ambient in utero and post-delivery lifetime exposure by child birth year for children in the Demographic and Health Survey [DHS, also known as the National Family Health Survey 2015–2016 (NFHS-4)] in India.

	2010	2011	2012	2013	2014	2015	2016	All years
Children/cluster (n)
Children	15,206	40,727	51,960	50,769	51,668	33,951	8,507	252,788
Clustera	8,445	19,766	22,548	22,353	22,419	17,217	5,590	27,835
In utero PM_2.5 exposure (μg/m^3)b
Mean (SD)	72.5 (31.7)	68.3 (27.0)	72.3 (28.4)	72.1 (29.1)	71.2 (28.0)	70.8 (26.3)	72.0 (27.8)	71.1 (28.2)
Min	21.2	20.9	22.6	22.4	23.1	24.7	23.7	20.9
Q1	45.7	47.9	51.1	49.9	49.7	50.6	50.2	49.7
Median	61.3	60.1	63.6	63.2	62.4	62.5	62.0	62.4
Q3	96.2	86.9	91.8	91.1	90.4	88.2	95.6	90.3
Max	145.1	142.1	153.5	153.0	144.1	144.2	141.4	153.5
Post-delivery lifetime PM_2.5 exposure (μg/m^3)c
Mean (SD)	73.7 (30.7)	73.0 (28.9)	72.7 (28.3)	73.1 (28.8)	76.0 (31.7)	75.2 (31.3)	67.9 (28.5)	73.7 (29.8)
Min	15.7	14.4	14.0	14.7	15.6	15.1	17.4	14.0
Q1	48.2	50.5	50.5	50.5	51.0	52.6	48.6	50.6
Median	61.5	63.0	62.8	63.3	65.1	67.2	61.6	64.0
Q3	106.4	101.5	98.8	98.5	101.3	90.9	81.6	98.3
Max	211.2	238.3	247.3	234.6	241.8	232.2	228.6	247.3

Note: Max, maximum; min, minimum; , particulate matter in aerodynamic diameter (fine particulate matter); Q, quartile; SD, standard deviation.

Cluster is a primary sampling unit (PSU) or a segment of PSU, of around 100–150 households, in the two-stage clustering sampling design of the DHS/NFHS-4.

9-month average of ambient exposure prior to child birth month.

Average of post-delivery lifetime ambient exposure from month of child birth until the month of death, or month of DHS/NFHS-4 interview of the child’s mother (censoring).

In the two-exposure Cox regression models, both in utero and post-delivery lifetime average exposures were significantly () associated with increases in child mortality, infant mortality, and neonatal mortality, after controlling for covariates (Table 3). A increase of in utero exposure was associated with an increase in neonatal mortality [; 95% confidence interval (CI): 1.001, 1.035], infant mortality (; 95% CI: 1.006, 1.037), and child mortality (; 95% CI: 1.008, 1.038). A increase of post-delivery average exposure was associated with an increase in neonatal mortality (; 95% CI: 1.003, 1.030), infant mortality (; 95% CI: 1.003, 1.027), and child mortality (; 95% CI: 1.001, 1.026). In the single-exposure Cox regression models, with either in utero exposure or post-delivery average exposure, we found similar but relatively stronger HRs compared with HRs in the two-exposure models (Table 3). In the parameter correlation estimates in the two-exposure models, we found moderate-to-high correlation between exposure in utero and post-delivery (Table S7). However, the two-exposure models had lower AICs compared with the single-exposure models (Table S8). Table S9 lists the HRs and 95% CIs of covariates in the two-exposure models of child mortality, and Table S10 lists the HRs and 95% CIs of mortality per IQR increase of exposure in the single- and two-exposure models.

Table 3

Hazard ratios (95% confidence intervals) of neonatal (), infant (), and child () mortality per increase of ambient exposure for children in the Demographic and Health Surveys [DHS, also known as the National Family Health Survey 2015–2016 (NFHS-4)] in India.

Models		All-cause neonatal mortality (<1 month of age)	All-cause infant mortality (<12 months of age)	All-cause child mortality (<60 month of age)
Deaths (n)		7,520	10,862	11,559
Two-exposure modelsa	In utero PM_2.5	1.018 (1.001, 1.035)	1.021 (1.006, 1.037)	1.023 (1.008, 1.038)
	Post-delivery average PM_2.5	1.017 (1.003, 1.030)	1.015 (1.003, 1.027)	1.013 (1.001, 1.026)
Single-exposure modelsb	In utero PM_2.5	1.032 (1.019, 1.045)	1.033 (1.021, 1.044)	1.033 (1.022, 1.044)
	Post-delivery average PM_2.5	1.025 (1.015, 1.035)	1.025 (1.016, 1.034)	1.025 (1.016, 1.033)

Note: Adjusted hazard ratios and 95% confidence intervals for all-cause neonatal mortality, infant mortality and child mortality are shown for increase of ambient during 9-month in utero period before child birth and post-delivery periods. Models stratified on child sex, birth month and year, geographical zone (927 strata), adjusted for birth order, multibirth, birth location, mother’s age, height, marital, education, maternal smoking, household wealth, secondhand smoking, cooking fuel, improved toilet, urban or rural location of households, monthly temperature, monthly precipitation. , particulate matter in aerodynamic diameter (fine particulate matter).

Two-exposure models include both ambient exposure during 9-month in utero period before child birth and post-delivery lifetime average until death, or censoring.

Single-exposure models include either one of exposure during 9-month in utero period before child birth or post-delivery average until death, or censoring.

We found little evidence of effect modification by most covariates on exposure on mortality in subgroup analysis. Figure 2 shows the estimated in utero and post-delivery lifetime average ambient effects on child mortality across covariates, and the results were similar for infant mortality and neonatal mortality (Figures S7 and S8). We observed some heterogeneity of post-delivery HRs across different maternal height, maternal age, and LBW, indicating some interactions between these covariates and post-delivery exposure. Numerical values for these results are shown in Table S11. For children with LBW, we observed in utero HRs of . In addition, we found that LBW was associated with both in utero exposure and child mortality, as shown in additional analyses (Tables S12 and S13). Our results were robust to different variable selections (modify 1–3 in Figure S9, Table S14), unweighted analysis (modify 4 in Figure S9, Table S14), and Cox regressions with generalized estimating equations to account for the correlation of outcomes within households and clusters (modifies 5 and 6 in Figure S9, Table S14), inclusion of time-dependent post-delivery monthly exposure to replace post-delivery lifetime exposure (modify 7 in Figure S9, Table S14), and crude models with single or two exposures (modify 8 in Figure S9, Table S14). The results for the categorized exposure models are listed in Table S15, showing a consistent linear trend for both categorized in utero and post-delivery exposure on mortality.

Figure 2.

Hazard ratios (HRs) and 95% confidence intervals (CIs) of child mortality per increase of ambient exposure, according to subgroup analysis for the two-exposure models for children in the Demographic and Health Surveys [DHS, also known as the National Family Health Survey 2015–2016 (NFHS-4)] in India. Adjusted HRs and 95% CIs for all-cause child mortality are shown for increases in ambient exposure during both (A) in utero and (B) post-delivery lifetime periods, stratified on subgroups by individual and household-level covariates. Some error bars and points are off the limits outside the image. All numerical values can be found in Table S8. All model specifications are the same as for the main analysis except for the subgroup analysis of low birth weight (LBW), which is not included in the main analysis. Red circles (A) represent effects of in utero exposure, and blue triangles (B) represent the effects of post-delivery lifetime average exposure. -Values for testing differences among strata for each stratifying variable are listed. Note: , particulate matter in aerodynamic diameter (fine particulate matter); Q, quartile; SHS, secondhand smoking.

As for the child mortality burden from air pollution, we estimated that 17.6% (95% CI: 12.3%, 22.4%) of mortality could be reduced if exposures were reduced to and that 9.3% (95% CI: 6.5%, 11.9%) of mortality could be reduced if exposures were reduced to the National Ambient Air Quality Standard of India of (Table S16).

Discussion

In a large, nationally representative retrospective cohort of children in India, we found a consistent association between exposure to ambient during in utero and post-delivery periods and elevated child mortality. The estimated associations were stronger in single-exposure models compared with two-exposure models, perhaps due to the strong correlation between in utero and post-delivery exposures. Similar results were observed for infant and neonatal mortality. Even though in utero and post-delivery lifetime exposure showed some multicollinearity in the two-exposure models, two-exposure models had smaller AICs compared with single-exposure models, indicating the proper fit of the two-exposure models. To our knowledge, this is the first cohort study conducted in India to estimate the effects of ambient exposure on child survival in a nationally representative cohort of children.

To estimate monthly concentrations, we developed a model based on multiple inputs from satellite observations, meteorological data sets, and land use information from between 2009 and 2018. We identified disproportionally high levels of in the Indo-Gangetic Plain Region, where the population density is high. In addition, we found increasing levels of ambient over the past 10 y, emphasizing an urgent need to control ambient air pollution in India. Children included in our study had high exposure levels to ambient at nearly twice the National Ambient Air Quality Standard of India () (Gautam 2009) and more than seven times the WHO Air Quality Guidelines () (WHO 2005). Given the high levels of ambient throughout India, reducing ambient air pollution levels could provide a relatively large reduction in child mortality and substantial health benefits.

The child mortality-HR estimates per of exposure in our study were in the range of 1–3%, lower than previous estimates in LMICs. In sub-Saharan Africa, Heft-Neal et al. (2018) found a 9% increase in risk of infant mortality associated with a increase in annual concentration, and Loomis et al. (1999) found a 7% increase of infant mortality associated with a increase of concentration during the 3–5 d before death in Mexico City. In the present study, given that the exposure range was three times higher than those in previous studies, the estimated increases of mortality per IQR (one ) increase of were in the range of 5–12%, more similar to earlier findings (Heft-Neal et al. 2018; Loomis et al. 1999). One study, which assessed early life ambient exposure and child mortality in 43 LMICs, did not find a significant association between air pollution and child mortality (Goyal et al. 2019). This could be due to their use of annual exposure as a proxy for exposure at year of age, which thus introduced measurement error and which may have biased results toward the null (Zeger et al. 2000). Another study conducted in Beijing using a time-series design also did not find a significant association between current-month exposure and infant mortality (Wang et al. 2019).

The present study has several strengths. First, we used machine learning to develop a high-quality model of ambient concentrations that had monthly temporal resolution and high spatial resolution that could reflect seasonal differences in ambient air pollution. These estimates of ambient levels could be potentially applied in future health impact assessment and epidemiological analysis in India. Second, we used a large, retrospective, and nationally representative cohort of children in India based on the well-documented and -conducted DHS/NHFS-4 survey. Third, we applied Cox regression with time-dependent air pollution, temperature, and rainfall precipitation on child mortality using the calendar month as the time variable. This approach allowed us to examine the longitudinal association between air pollution and child mortality and compare post-delivery exposure for children with the same follow-up time. This approach could largely reduce exposure misclassification compared with previous studies (Goyal et al. 2019; Heft-Neal et al. 2018), which were based on yearly averages to assess a binary outcome of child death at the end of 1 year of age. Our analysis controlled for potential bias from seasonality and different weather patterns, in addition to using individual-level, household-level, and cluster-level covariates. Notably, we did not include LBW as a covariate in the model due to the potential selection bias of including LBW as a covariate, and for children with LBW, we found protective effects of in utero exposure. This is similar to the protective effects of maternal smoking on infant mortality for children with LBW observed in the United States as a result of selection bias (Hernández-Díaz et al. 2006). Finally, the associations were consistently positive in single- and two-exposure models, and these associations were robust in additional sensitivity analyses.

Our study also has some limitations. First, no national monitoring data were available during the study period of 2009–2014, so we lacked ground-level measurement data before 2015 with which to validate our prediction model. The monitoring stations were more heavily concentrated in urban areas, which could have led to differential errors across the country as well. Similar to previous prediction models (Xiao et al. 2018), our model tended to underestimate some very high values, when monthly concentrations were (Figure S3). Second, we used only total ambient mass concentration at the cluster level as our exposure for participants at the household level, without estimation of different components, and did not include the impact of other ambient air pollutants, such as PM with an aerodynamic diameter of , , and ozone, or direct household air pollution exposures from cooking, which are significant in many parts of India. Third, our retrospective cohort was based on a cross-sectional survey of the DHS/NFHS-4. We focused only on all-cause mortality reported in the survey, whereas it is likely that ambient is associated with specific cause of death. Child birth year, birth month, vital information, and covariates were collected at the end of the study and recalled by mothers and may have differed from values at the time the children were born. In addition, we did not have the children’s gestational age at birth and assigned in utero exposure levels as a 9-month average prior to the birth, which may not be the case for preterm infants who have a shorter in utero exposure period. Given this data collection limitation, our study could have suffered from recall bias and misclassification. However, we believe most of the maternal and household-level covariates will be nondifferential misclassification of outcomes and could most likely bias our results toward the null. Last, as with all observational studies, our study may have potential unmeasured confounders that were inadequately controlled for in this analysis. Our results indicate a need for future prospective cohort studies (Balakrishnan et al. 2015), including intervention trials of reducing air pollution levels (Clasen et al. 2020), to investigate air pollution effects on child health.

In conclusion, this study found increased mortality risk associated with ambient during the 9-month in utero and post-delivery periods for children . This study, based on a nationally representative retrospective cohort in India, provides evidence suggesting that in utero and post-delivery exposure contribute to child mortality in developing countries. Given the high levels of ambient throughout India, expanding air pollution monitoring stations (both in rural and urban areas), adding more epidemiological research, and making a substantial effort to reduce ambient air pollution and early life exposures are needed.

Click here for additional data file.

Click here for additional data file.