We extracted and post-processed routine surveillance data for malaria from the DHIS2 database in collaboration with the NMCP of Côte d’Ivoire. In the country, the reporting follows the hierarchy of referral in the sanitary pyramid, i.e., health facilities in the periphery with recently added community health workers (CHW) and local general hospitals (Hôpital Général, HG) refer to regional hospitals (Centre Hospitalier Régional, CHR) who refer to university hospitals (Centre Hospitalier Universitaire, CHU). The nested spatial organization consists of 113 health districts within 26 health regions that report to the Ministry of Public Hygiene and Universal Health Coverage at the national level. We considered clinical malaria cases of all ages reported at the level of health districts. These monthly data were entered into DHIS2 from paper registries in health facilities, HG, and CHR by local health workers and validated by the NMCP. As for CHW and CHU, inconsistency in reporting across districts and time, prohibited further consideration for this analysis. Since presumed cases (i.e., without confirmed tests) have been reported into the DHIS2 database only since 2021, we considered the difference between suspected and tested cases as presumed cases whenever the latter were not directly available in the database.

In preparation of the impact evaluation for IRS in the districts Nassian and Sakassou with control districts Dabakala and Beoumi, PMI decided to collect data directly from health facility registries in the four study districts due to insufficient data quality [15]. Patient data were abstracted from consultation registries (“registres de consultations curatives”) for the period September 2018 to April 2022, we refer to details of inclusion criteria to recently published study results [9].

In order to control for factors unrelated to the intervention, we also included environmental, entomological and surveillance covariates into the causal impact analysis (see Supplementary Table S1). Remote-sensed data (rainfall, average daily temperature, leaf cover, average daily humidity) was extracted from the ERA5 database [16], aggregated monthly at the health-district level. Entomological covariates (monthly indoor and outdoor biting rates and monthly entomological inoculation rates) were extracted from the PMI VectorLink annual entomological reports for the years 2019–2021 [17–19] for the four districts of interest. To take into account changes in the health surveillance we also derived from the DHIS2 database the following district-level monthly quantities: health facility reporting rate (short: hf_reporting), monthly testing rate (short: testing), test positivity rate (short: tpr_total), the number of active health facilities (short: N_hf_active) and coverage of intermittent preventative therapy for pregnant women for the first trimester (IPT1) (short: ipt1_coverage). Since the 2021 mass ITN campaign was rolled out in the four districts from May until June 2021, we assume that residual effects of the 2017 ITN campaign are negligible during the IRS evaluation period from July 2020 to July 2021. Although IPT1 concerns only a small proportion of the population at risk, we consider it as a general proxy for access to healthcare. We normalized and scaled all covariates prior to further analysis.

To adjust incidence (to both DHIS2 and PMI crude data), we used the following epidemiological indicators. First, we extracted monthly confirmed malaria cases reported and validated at the level of health districts in the DHIS2 system or by PMI respectively. Second, we used presumed cases (short: pr) from the DHIS2 database which are defined as cases with malaria-like symptoms without having been tested. Third, we calculated the test positivity rate (short: tpr) as the ratio of confirmed over tested cases. Here, we assumed that cases were tested with either blood spot or rapid diagnostic test (RDT). For the reporting rate, we considered for each health district the ratio of months within a year, where health facilities reported into the DHIS2 system. Finally, the healthcare seeking rate (short: s) was extracted from the 2016 Malaria Indicator Survey [20] and the 2021 Demographic and Health Survey (DHS) report [21] estimates of healthcare seeking for children with fever made at the regional level. Sakassou, Beoumi and Dabakala are located in the Center-North region and Nassian in the North East region. We interpolated these values linearly between the years. We then calculated adjustments to incidence rates for both data sources DHIS2 and PMI using the same quantities as detailed in Supplementary Table S2. District-level population data were extracted from the 2021 census data with 3% annual growth rate extrapolation. We calculated incidence rate as the ratio of adjusted incidence over district population.

In the period of interest, we have detected 6.7% of missing data among all contemporaneous covariates across the four health districts (see Supplementary Figure S1). We used the R package Amelia [22] to perform multiple imputations of missing data for a subset of contemporaneous variables (biting_rate_hcl_outdoor, biting_rate_hcl_indoor, EIR_ibpn), with health district as cross-sectional variable and the month during the year as temporal variable to account for seasonality. The variable EIR_ibpn refers to the entomological inoculation rate (EIR), i.e., the number of infectious bites received by a host within a year. Based on the rule-of-thumb we performed 7 imputations for 7% of missing data. We averaged imputed data points and used them to replace missing data in the original data set. In addition, we also imputed two incidence points for the Beoumi district, which had suspicious outliers.

In order to identify associations between difference in incidence rates recorded in the DHIS2 and PMI data sets, we performed a mixed-effect linear model analysis [23] on the difference of log-transformed incidence data with a random intercept for the health district. We used a backward stepwise model selection approach [24] to eliminate sequentially factors from two models (see Supplementary Tables S3, S4): model1 had only health system covariates (IRS intervention, hf_reporting_rate, N_hf_active, tpr_total, ipt1_coverage), model2 had in addition also environmental variables (rainfall, temperature, leaf, humidity, biting rates, EIR) without and with 1 month lag obtained from cross-correlation analyses. For the resulting parsimonious model, we determined the coefficient estimators and two-sided p-values.

We performed a Pearson correlation analysis between monthly incidence rates and contemporaneous variables for lags of up to 1 month in the past. This analysis was performed across health districts and per incidence data source.

Bayesian structural times series analysis [25, 26] is an extension of classical time series analysis techniques such as auto-regression moving average models. The method integrates elements of state-space models (i.e., observation and hidden state parameters are estimated concurrently) and Bayesian feature selection (e.g., through spike-and-slab priors, acting like a probabilistic “on-off switch” for every variable in a regression model) and provides forecasting of time series with high-dimensional time-dependent covariates. We utilize this method to construct counterfactuals for intervention impact quantification. For technical details we refer to the Supplementary Figure S2. For monthly outcome measures (malaria incidence rates or EIR) we implement local linear and seasonal effects, and use contemporaneous regressors for lagged time series that are independent from the incidence rates (e.g., rainfall, temperature). In order to train the pre-intervention model we use either the incidence data in the IRS districts (control: “none”), incidence data in both IRS and non-IRS districts (control: “incidence), or all contemporaneous covariates including incidence from all districts (control: “all covariates) but excluding entomological covariates, since those are supposed to be statistically dependent on incidence. Similarly, for the impact analysis with EIR, we use EIR data in IRS districts only (control: “none”), EIR from both IRS and non-IRS districts (control: “EIR”), or all contemporaneous covariates except incidence covariates, since those are supposed to be statistically dependent on EIR.

From the Bayesian feature selection for models with contemporaneous covariates, we extract marginalized posterior distribution of regression parameters to quantify how often a regressor has been included in the model, but also its magnitude and whether 95% of the probability mass was either positive or negative (“consistent regressors”). Causal impact is determined by the difference of predicted incidence in the absence of the intervention (i.e., from the counterfactual) and the actually recorded data in presence of the intervention, averaged over one-year after the intervention deployment. We consider causality as established if in less than 5% of 10,000 simulated counterfactuals we observe incidence inferior to the observed incidence with IRS.

After adjusting crude incidence of malaria cases from DHIS2 and PMI data sources (Supplementary Table 2, adjustment 2), we observed a range between 15 and 100 monthly incidence points per 1000 population between September 2018 and April 2022 with peaks occurring generally between July and October (Figure 1). Two consecutive months for the control district of Sakassou had zero incidence recorded in the DHIS2 data. We considered these as missing data and imputed them using multiple imputation algorithms (see Methods). PMI incidence rates were generally higher than those based on DHIS2 data, and Sakassou followed more consistently the rainfall seasonality trends than Nassian. The impact analysis focused on incidence data until August 2021, date of the start of the second deployment of IRS. Analyzing the difference between incidence rates from the two data sources including model selection for a comprehensive set of contemporaneous and categorical variables (see Supplementary Tables S2, S3) suggested that there were several statistically significant factors explaining the difference. We identified the presence of the first IRS intervention, increased IPT1 coverage and test positivity rates to be statistically associated with a decrease in difference (P < 0.05) Table 1.

For the correlation analysis between monthly incidence rates and contemporaneous variables, we considered both data from the same month and with a one-month lag across all four districts before and after the intervention (Figure 2). We noticed that entomological variables were highly correlated amongst each other (e.g., 0.57 and 0.53 Pearson correlation coefficient between EIR and indoor or outdoor biting rates). Biting rates also positively correlated with malaria incidence for the PMI data (0.53 and 0.57 correlation coefficient). Strikingly, none of the entomological variables had significant correlations with malaria incidence data from DHIS2. On the other hand, DHIS2 incidence data significantly correlated with rainfall and humidity, albeit with smaller coefficients than PMI data. A similar pattern was observed for correlations between incidence and health surveillance yield and quality covariates (number of health facilities active, health facility reporting rates) which had significant correlations with PMI incidence data (0.27, and −0.32 resp.), whereas these same covariates did not correlate with DHIS2 incidence data. Only test positivity rates positively correlated with incidence data from both sources.

Based on our findings from the correlation analysis, we included all except the entomological covariates into our Bayesian structural time series analysis (see Methods) as these were causally dependent on the intervention. We trained a model on incidence data prior to the intervention including contemporaneous variables and utilized the resulting model estimates to predict a counterfactual for the Nassian and Sakassou district during the intervention period. Finally, to measure impact, we compared the counterfactual incidence with actually recorded incidence data (see Supplementary Figure S3). Using the incidence in the IRS districts alone during the model training (Figure 3, control = none), we obtained an estimated average decrease in incidence during the first year post intervention of 30% (95%CI: 11%–43%) for Nassian and 36% (95%CI: 20%–48%) for Sakassou using DHIS2 data. Incidence decrease according to PMI data was not significant. Including incidence data from the control districts (Dabakala for Nassian and Beoumi for Sakassou) to the pre-intervention model training resulted in significant impact of 12% (95%CI: 3%–19%) resp. 19% (95%CI: 10%–27%) for PMI data and 30% (95%CI: 19%–38%) resp. 43% (95%CI: 34%–50%) for DHIS2 data. Considering all covariates besides those related to entomology yielded an estimated impact for PMI data of 9% (not significant) for Nassian and 23% (95%CI: 12%–31%) for Sakassou. No significant impact for DHIS2 data was determined in that case. The large error estimates from the credence intervals suggest a high variability in posterior variable selection.

Utilizing the Bayesian structural time series model, we obtained regression coefficients for each of the contemporaneous covariates. Posterior distributions from the Bayesian updating stem from spike and slab priors (see Methods), such that not all coefficients were necessarily selected in each model run. We focused on the marginalized regression coefficients that consistently either had negative or positive values for 95% of the model runs in order to quantify covariates influencing the counterfactual, and consequently the evaluated impact. E.g., a consistently positive regression coefficients would mean increased incidence for the counterfactual for an increase in the covariates, such that the overall intervention impact would be have been overestimated without the covariate, and vice versa for negative regression coefficients.

Considering empirical distributions of posterior regression coefficients (Figure 4), we observed that incidence from the control district had positive estimates across datasets and health districts with high probability of inclusion for the model selection process. Health facility reporting rates were consistently associated with negative coefficients for PMI data, suggesting that for periods with increased health facility reporting during the intervention, one would have estimated a decrease in the counterfactual and consequently a decreased impact of IRS. For DHIS2 on the other hand, test positivity rate, number of health facilities active and rainfall had positive coefficients, such that an increase in these factors might be the main driver of the overestimated impact of IRS (Figure 3). We stress that the posterior distribution of covariate regression coefficients used for the impact quantification are high-dimensional and not independent from each other such that the underlying causality might rather be a combination of many factors than just those mentioned above.

Utilizing the same approach as for the causal impact analysis towards incidence, we also investigated causal impact towards EIR. This time we excluded malaria incidence as covariates, as these were not independent from EIR rates. We observed a strong and statistically significant causal impact of IRS on the decrease of EIR at the order of 70% averaged during 1 year post intervention (Figure 5). Among consistent regressors, only EIR values from control districts were chosen with high inclusion probability (Figure not shown).