Spatiotemporal analysis of the association between Kawasaki disease incidence and PM_2.5 exposure: a nationwide database study in Japan

Data source

In this retrospective study, we extracted clinical data from the Japanese administrative claims database named the Diagnosis Procedure Combination (DPC) database, comprising anonymised clinical and administrative claims data featuring baseline information of patients and facilities, diagnostic records, procedural data, device utilisation and prescription details. As of 2023, over 2000 hospitals had implemented DPC-based reimbursement systems. This database substantiated its reliability through prior research.25 Data were accessed on 16 August 2023. Among hospitalisation data from April 2014 to March 2022, we extracted clinical information on children under 5 diagnosed with KD, identified by the International Classification of Diseases, Tenth Revision (ICD-10) code of M30.3. To minimise bias associated with misclassification, we focused on hospital admissions where patients received KD-specific medications, namely intravenous IG, cyclosporine A, infliximab or ulinastatin.3 6 7 We considered the date of first admission with KD treatment as the onset date, excluding cases with unclear onset dates, specifically transfer cases and those not administered intravenous IG within 7 days of the first admission. To address uncertainties associated with identifying of initial hospitalisations, cases of KD that occurred in the first 3 months of the observation period were excluded, given the risk of misinterpreting the middle of a series of hospitalisations that began before the observation period as the onset. Cases from the last 3 months of the period were also excluded, as the number of onsets during this period may be underestimated due to administrative delays in medical claims processing. Then, the timeframe from July 2014 to December 2019 was defined as the period before the COVID-19 pandemic, whereas from January 2020 to December 2021, it was defined as the period during the COVID-19 pandemic.

The atmospheric environment database of the National Institute for Environmental Studies publishes pollution data from 2184 monitoring stations across 319 (95%) of the 335 secondary medical care areas in Japan.26 Each secondary medical care area, established across 1718 of the 1724 municipalities and managed by the 47 prefectural governments, ensures general inpatient treatment, including initial treatment of KD. We extracted daily exposure to PM_2.5, nitric monoxide (NO), nitrogen dioxide (NO₂) and sulphur dioxide (SO₂) for each medical care region, imputed missing values using the prefectural average and calculated monthly exposure. As a result, we obtained 22 100 and 8040 spatiotemporal units based on the exposure status in 335 secondary medical care areas over 66 months and 24 months before and after the onset of the COVID-19 pandemic, respectively.

Outcomes and variables

As an outcome measure, the monthly incidence of KD was counted for each secondary medical care area associated with facilities. The monthly or annual exposure to PM_2.5, NO, NO₂ and SO₂ in the corresponding area was incorporated in the analysis as continuous variables. The logarithm of person-days for each spatiotemporal unit based on the under 5 population in the Population Census 2020 was implicitly incorporated in all the statistical models as an offset variable.27

Statistical analysis

To capture the fundamental relationship between KD incidence and exposure to PM_2.5, NO, NO₂ and SO₂, we developed non-Bayesian Poisson regression models, both univariable and multivariable, using overall exposure levels during the two distinct periods before and after the onset of the COVID-19 pandemic. Subsequently, we performed Markov chain Monte Carlo (MCMC) simulations with the CARBayes library V.6.1 and CARBayesST library V.5.0 in R V.4.3.2 to create four types of multivariable Bayesian Poisson regression models predicting the monthly incidence of KD based on 1-month and 12-month exposure to these air pollutants: ‘GLM model’ is a Bayesian implementation of a generalised linear model that ignores spatiotemporal autocorrelations; ‘CARar(1) model’ is a first-order CAR model, where ‘first-order’ indicates that the model accounts for dependencies on the immediately previous time step; ‘CARar(2) model’ is an extension of the CARar(1) model, incorporating dependencies on the past two time steps; and ‘CARadaptive model’ is another first-order CAR model, which includes an adapted spatial weight matrix to handle spatial heterogeneity.20 21 28–31 We adopted the model with the lowest widely applicable information criterion (WAIC) among these four Bayesian models.32 Univariable models were also developed to assess the impact of individual air pollutants. The parameters were estimated from distributions derived from 40 000 MCMC samples, equating to 400 000 iterations with a thinning factor of 10 to reduce autocorrelation. This estimation followed an initial burn-in period of 100 000 iterations to stabilise the sampling process. In the sensitivity analysis, we developed comparative Bayesian models with subjects divided into three age groups: 0 years, 1 year and 2–4 years.

Patient and public involvement

Patients and/or the public were not involved in this study’s design, conduct or dissemination.

Limitations

Selection bias is a concern in observational studies. In light of the incidence rates of KD reported in previous studies, it can be estimated that approximately 70% of the domestic cases were included.23 Although the inclusion criteria were carefully constructed based on the ICD-10 code and KD-specific medications, the level of concordance between the judged and actual onset of KD is yet to be confirmed. In this real-world data study, information on symptoms and clinical findings was not available. We considered the risk of misclassification with multisystem inflammatory syndrome in children (MIS-C) to be negligible based on the rarity of MIS-C cases in Japan.8 38 The exclusion of untreated cases can be expected to be marginal, considering the ubiquity of early intravenous IG administration in Japan.2 Imputation of exposure at the prefectural level for the small amount of missing data may have biased the analyses toward the null. Although the dose–response relationship observed in this study aligns with previous research conducted in geographically close Korea, different results might be obtained in distant countries due to varying sources of PM_2.5. Unmeasured substances or microorganisms dispersing similarly to PM_2.5, rather than PM_2.5 itself, might be involved in the onset of KD.39 40 Besides, it should be noted that spatiotemporal analysis with different granularities of spatiotemporal units may yield different results.41 Given the limited geographic activity range of children under the age of five, the impact of exposure outside their secondary medical care area would be minimal. Analysis with a finer granularity would pose challenges due to boundary-crossing admissions, while extensive unit aggregation would reduce statistical power. We handled data at the spatiotemporal unit level, thereby not distinguishing between prenatal and postnatal exposures at the individual level. While the impact of annual PM_2.5 exposure in infants under 1 year may imply potential influences of prenatal exposure, these effects have not been explicitly examined.

In conclusion, we used the CAR models to address the spatiotemporal aggregation of KD, confirming the robust association between the incidence of KD and annual exposure to PM_2.5. Further investigation is required to clarify the underlying mechanism of association between the spatiotemporal distribution of KD and PM_2.5.