Analysis overview

We analysed data from a cluster randomized trial (NCT04864444) evaluating MDA for accelerating malaria elimination in low to moderate transmission settings in Senegal. We summarize key trial information below; additional trial details were previously published⁴.

Study site

The trial was conducted in the Tambacounda Health District of southeastern Senegal, where malaria transmission is low to moderate and highly seasonal, with most cases occurring from July to December. At this site, the National Malaria Programme (NMP) implements standard malaria control interventions, including mass distribution of insecticide-treated nets, health facility-based malaria case management and SMC for children aged 3–120 months. In remote areas with limited healthcare access, Senegal uses malaria case management through the Prise en Charge à Domicile (PECADOM) model, whereby village-level community health workers (dispensateurs de soins à domicile (DSDOMs)) are trained to test and treat febrile malaria cases using rapid diagnostic tests (RDTs) and first-line antimalarials. Some villages use PECADOM+, in which DSDOMs conduct weekly proactive household visits to improve early detection and treatment.

In other regions, malaria endemicity is heterogeneous, driven largely by geography, climate, population density and age distribution. To accommodate varying transmission intensities, the NMP adopts a tailored set of interventions. In the northern and central regions, where transmission is very low due to the semi-arid Sahelian climate, the NMP prioritizes pre-elimination activities consisting of case investigations and proactive response measures. In the central and southern regions, where transmission ranges from low to high due to denser vegetation and higher rainfall, the NMP prioritizes malaria control interventions, including routine net distribution, prompt case management and SMC in highly seasonal areas.

Study design

The original study was a two-arm, open-label, cluster randomized controlled trial. Sixty villages were randomly selected based on the following eligibility criteria: (1) population size between 200 and 800; (2) location within a health facility catchment area with an annual incidence of 60–160 cases per 1,000 population; and (3) the PECADOM+ model was established or eligible for roll-out. To minimize intervention contamination, a buffer zone of 2.5 km was maintained between village centroids.

Ethical approval

The study protocol was approved by the Comité National d’Ethique pour la Recherche en Santé (Dakar, Senegal) and the University of California, San Francisco Human Research Protection Program (San Francisco, CA, USA). Stanford investigators had access only to de-identified, cluster-aggregated monthly data.

Randomization and blinding

Villages were randomized 1:1 using a stratified, constrained randomization approach based on the following covariates: baseline DSDOM presence, health facility, distance to the nearest health facility, baseline microscopy-confirmed malaria prevalence, population in 2019, and the population of children

Interventions

During the intervention year (2021), intervention villages received three rounds of MDA with dihydroartemisinin–piperaquine plus a single low dose of primaquine every 6 weeks. Control villages received three rounds of SMC with sulfadoxine–pyrimethamine plus amodiaquine administered to children aged 3–120 months every 4 weeks, according to the standard of care. MDA was administered every 6 weeks based on discussions with the Tambacounda District Medical Office to ensure adequate spacing of drug courses, minimizing potential side effects. The schedule was supported by modelling evidence⁷ suggesting that this interval would provide adequate protection, given piperaquine’s long half-life²⁸. The timing also ensured comparable coverage of the transmission season relative to SMC. After the three rounds of MDA, the study villages were followed for an additional year (2022), during which all villages resumed receiving SMC (Fig. 1).

Procedures

Community sensitization was conducted before the intervention delivery. Before the intervention implementation, all study villages received a mass distribution of pyrethroid–PBO bed nets and year-round PECADOM+. A baseline survey was conducted to assess pyrethroid–PBO net coverage and malaria prevalence (Fig. 1). SMC and MDA were delivered door-to-door using an age-based dosing strategy. All three doses of SMC and MDA were directly observed. During the campaign, suspected cases of malaria were confirmed by RDTs and treated with artemether–lumefantrine. In positive cases, chemoprevention was deferred until the next cycle.

Throughout the study period, data on RDT-confirmed malaria cases were collected from health facilities and PECADOM+ registries. Average village population size was estimated by calculating the mean of the two censuses conducted before and after the intervention implementation. Additional details of the procedures are provided in Supplementary Appendix 13.

Outcomes

The primary outcome of the trial was village-level P. falciparum incidence during the postintervention transmission season (July–December 2022), defined as the number of RDT-confirmed, symptomatic malaria cases divided by the mean village population size (Fig. 1). Malaria incidence during the transmission season of the intervention year was a secondary outcome.

Effect measure modifiers

We obtained high-resolution data for Senegal on the following potential effect measure modifiers: P. falciparum malaria prevalence, precipitation, temperature, vegetation density, percentage of the population under 10 years old and travel time to the nearest health facility. These variables were hypothesized to influence intervention effectiveness by affecting transmission intensity, vector ecology and access to care. For example, moderate rainfall can increase vector breeding sites, while prolonged and heavy rainfall may flush mosquito larvae from breeding sites, thereby reducing malaria transmission²⁹. Vegetation density, captured by EVI, may enhance mosquito survival and breeding^30,31. Temperature has been shown to have a non-linear association with the Anopheles biting rate, vector competence, survival and P. falciparum development³². Population density has a non-monotonic relationship with malaria; transmission models have linked higher population densities to persistent malaria transmission outside of the peak season³³. The proportion of children is a predictor of transmission intensity, given their disproportionate malaria burden³⁴. We ultimately excluded surface water, predicted malaria incidence, predicted malaria prevalence, night-time light radiance, travel time to the nearest health facility and population density from the transportability analyses due to a lack of variation or other reasons detailed in Supplementary Appendix 14.

For estimates of the percentage of the population aged 35. Daily precipitation data at 0.05° resolution were sourced from the Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) dataset³⁶. Using Google Earth Engine, we aggregated daily data into monthly totals. Mean monthly temperature data at 0.05° resolution were obtained from the MODIS Monthly CMG Land Surface Temperature and Emissivity (MOD21C3)³⁷. Because our intent was to transport estimates to small administrative areas, we used a dataset with higher spatial resolution (0.05°) but with monthly, rather than daily, data. Thus, we included only the mean monthly temperature.

While prior studies have shown that daily temperature variation may influence malaria transmission³⁸, high-resolution daily remote sensing datasets exhibited missingness rates of 10–20% during high-rainfall months, probably due to cloud cover. Monthly data on the presence of surface water at 30-m resolution were sourced from the Joint Research Centre Monthly Water History, v1.4 dataset³⁹. Monthly EVI data at 1-km resolution were obtained from the MODIS/Terra Vegetation Indices Monthly L3 Global 1 km SIN Grid⁴⁰. We extracted the mean values of the percentage of the population under 10 years old and the monthly temperature for the geocoordinates of each trial village centroid and for all communes in Senegal. For precipitation, we calculated the monthly minimum, mean and maximum levels. To account for delayed effects of environmental factors on malaria incidence, we used 1–2-month lags for precipitation, 0–2-month lags for temperature and 0–3-month lags for EVI, informed by prior studies. Justifications for these lags are provided in Supplementary Appendix 14.

Statistical analysis

First, we assessed the effect modification of the MDA intervention by geographical, environmental and demographic variables. Analyses were conducted at the village-month level across the preintervention, intervention and postintervention transmission seasons (that is, July–December of 2020, 2021 and 2022). Indicator variables were used to test for effect modification by commune. For continuous modifiers, we generated indicators of whether values were above or below the median within each follow-up year. To compare incidence rate ratios for MDA versus control within different levels of potential modifiers, we fitted mixed-effects Poisson regression models restricted to each level of the potential effect modifier, allowing for potential interactions between modifiers and baseline covariates. Models included village-level random intercepts, a log link, an offset for the mean village population size during follow-up and robust standard errors. Consistent with the original trial’s intention-to-treat analysis⁴, malaria case counts were modelled as the dependent variable, with covariates including the following: trial-year fixed effects; an MDA indicator equal to 1 for randomization to MDA in the intervention year (2021) and 0 otherwise; a post-MDA indicator equal to 1 for randomization to MDA in the postintervention year (2022) and 0 otherwise; an indicator equal to 1 for periods and villages with an existing PECADOM+ model and 0 otherwise to account for differential capture of malaria cases at baseline; and variables included in the constrained randomization (health facility, distance to a health facility, baseline microscopy-confirmed malaria prevalence, village population size and population size under 10 years old).

To test for effect modification on the ratio and additive scales¹⁹, we fitted models with two-way interaction terms between each modifier and year, as well as for the MDA indicator in intervention year analyses and the post-MDA indicator in postintervention year analyses. We defined R_xz as the incidence rate under treatment x and modifier z. The ratio scale measure of interaction was defined as (R₁₁R₀₀)/(R₁₀R₀₁), which can be expressed in terms of incidence ratios as IR₁₁/(IR₁₀IR₀₁) (ref. ¹⁹). We assessed the additive scale interaction using the relative excess risk due to interaction (RERI), defined as IR₁₁ − IR₁₀ − IR₀₁ + 1 (ref. ¹⁹). Because MDA had a protective effect on malaria incidence, we recoded the reference levels for the MDA indicator variable and each modifier to the lowest risk level⁴¹. We computed standard errors and CIs using the delta method. Effect measure modification analyses were conducted using Stata version 16.

To transport trial estimates to non-trial areas, we applied a doubly robust transportability modelling approach from Dahabreh et al.²⁰. We used a non-nested design in which the trial data were concatenated with data from an external target population that partially overlapped with the trial population⁴². When the difference between the trial and external target populations can be characterized by baseline covariates, transportability is conceptually similar to direct standardization across multiple covariates. These methods rely on measuring variables that both modify the intervention effect and differ in distribution between trial and target populations¹⁶. A flow diagram illustrating the step-by-step procedure of our transportability analyses is provided in Fig. 4.

We restricted analyses to communes where SMC was routinely offered and population density was similar to that of trial sites (≤152 people per 100 m²). This was done to minimize potential violations of positivity—that is, to ensure that trial populations could reasonably represent non-trial areas¹⁵. We combined trial and non-trial data into a single dataset: trial observations included treatment assignment, outcome data and effect modifiers, while non-trial observations included only effect modifiers. An indicator variable for trial participation (S_i) was included in the combined dataset for each village i (for trial data) and commune i (for non-trial data), where S_i = 1 for trial observations and S_i = 0 for non-trial communes. Estimates were transported to the commune level, the smallest standard administrative unit used for public health planning in Senegal. We aggregated malaria incidence data by summing the total number of cases and the population size, and we aggregated effect modifier data by calculating the minimum, mean, maximum or total of covariates as described in the prior section.

The doubly robust approach fits a model for the probability of trial participation and a model for the expectation of the outcome. If either is correctly specified, the estimator is consistent²⁰. First, for each commune, we fitted a model for the probability of trial participation Pr(S_i = 1|X_i,t), where X_i,t is a matrix of effect measure modifiers for commune i at month t that includes the variables listed in the previous section. Models were fitted using data from trial sites combined with data from non-trial communes. We then calculated weights by intervention arm (A_i) in each commune i at month t normalized to the sum of the number of non-trial areas as follows:

where A_i is an indicator for random assignment to MDA versus control, \(p({X}_{i,t};\hat{\beta })\) is an estimator of the probability that a commune would be included in the trial at month t (Pr(S_i = 1|X_i,t), \({e}_{a}({X}_{i,t};\hat{\gamma })\) is an estimator of the probability of assignment to treatment arm a among trial clusters (Pr(A = a|X_i,t,S_i = 1)), and I(S_i = 1|A_i = a) represents an indicator function equal to 1 for communes that were included in the trial and assigned to arm a and 0 otherwise.

Second, using trial data, we estimated separate outcome models within each intervention arm, modelling monthly malaria cases as a function of effect modifiers. \({g}_{a}({X}_{i,t};\hat{\theta })\) is an estimator for the conditional expectation of the number of malaria cases (E(Y_i,t|X_i,t, S_i = 1, A_i = a)), where Y_i,t is the number of malaria cases for each cluster i at month t. These models were then applied to covariate patterns in non-trial communes to generate arm-specific counterfactual predictions. Outcome models did not adjust for covariates used in the constrained randomization or include a log population offset, as in the original trial⁴, because these data were unavailable in non-trial areas. However, trial estimates were similar with or without these terms (Supplementary Appendix 10). Additionally, we did not adjust for baseline malaria incidence as was done in the original trial due to the lack of such data in non-trial areas. A re-analysis comparing the original trial’s estimates restricted to the intervention year (2021) or the postintervention year (2022) found minimal differences in point estimates and 95% CIs for the intervention year and moderate differences for the postintervention year (Supplementary Appendix 9). However, the interpretation of findings was consistent between the two models, with overlapping 95% CIs.

Given the relatively small number of communes and the potentially large number of covariates, we fitted both trial participation and outcome models using elastic net models with ∝ = 0.5, which use an equal mix of LASSO (least absolute shrinkage and selection operator) and ridge regularization to balance feature selection and the handling of correlated predictors. We fitted models for \({g}_{a}({X}_{i,t};\hat{\theta })\) using a quasi-Poisson family and \(p({X}_{i,t};\hat{\beta })\) and \({e}_{a}\left({X}_{i,t};\hat{\gamma }\right)\) using a binomial family. The relax parameter was set to true so that models were refitted without penalty after the initial variable selection, potentially improving prediction accuracy. To identify the optimal level of regularization, we used tenfold cross-validation and a lambda sequence ranging from 0.01 to 100, with a length of 100. We used elastic net model fits with the lambda value that minimized the cross-validated error. If only one covariate was included in the model, we used a generalized linear model with a binomial family instead. Transportability analyses were performed using R version 4.1.0; elastic net models were fitted with the glmnet package (version 4.1.2).

Before fitting the elastic net models, we performed several forms of covariate screening to mitigate multicollinearity, streamline the covariate set and improve model stability. First, we excluded covariates with no overlap between S = 1 and S = 0 to minimize empirical violations of the positivity of trial participation assumption. Second, we screened covariates for multicollinearity within predefined variable groups (precipitation, temperature and EVI) across lag specifications. For each group of covariates, we calculated pairwise Pearson correlations between variables and estimated their correlations with the outcome. If any correlations were >0.7 within a group of covariates, we retained the covariate with the strongest correlation to the outcome. Third, we screened covariates for inclusion in outcome models using Poisson likelihood ratio tests, retaining those with a P value of 15.

We then implemented a doubly robust estimator to estimate the intervention effects for each commune. The estimator uses data from trial clusters as well as from non-trial communes. For trial clusters, it uses outcome model predictions as bias correction terms, along with observed outcomes, weighted by the inverse probability of selection into the trial. For non-trial communes in the target population, it relies solely on outcome model predictions. The double robust estimator \(\hat{\mu }(a)\) is defined as follows:

We separately estimated weighted averages for the intervention and control arms and calculated percentage effectiveness as (1 − (incident cases in intervention/incident cases in control)) × 100% to estimate the transported effect. To obtain 95% CIs, we performed non-parametric bootstrapping by resampling communes with replacement, using 1,000 iterations.

To validate our transportability model, we transported trial estimates from the intervention year (2021) and the postintervention year (2022) to the geographical area surrounding the trial villages by creating a convex hull of the trial village centroids (Supplementary Appendix 11). Given that effect modifier values should be highly similar to those in geographical areas approximate to the trial site, we expected this analysis to produce effect estimates close to those of the original trial. We extracted effect modifier values from remote sensing data within the trial site, set S = 0, and then refitted the transportability models. We obtained 95% CIs using the bias-corrected and accelerated bootstrap method, as described above.

To make causal inferences from our transported estimates, several assumptions are required (Supplementary Appendix 16). First, we assumed conditional exchangeability between clusters assigned to MDA and those assigned to control. Second, we assumed positivity of intervention assignment, meaning that there was a non-zero probability of intervention assignment across covariate strata. We expect the first two assumptions to hold by randomization. Third, we assumed conditional exchangeability over trial participation, meaning that trial clusters were exchangeable with target communes conditional on covariates and study arm. Fourth, we assumed positivity of trial participation, meaning that there was a non-zero probability of trial participation in any covariate stratum needed to ensure conditional exchangeability. To minimize possible violations of this positivity assumption, we restricted analyses to communes where SMC was offered during the trial period and proceeded with estimation only in communes where \(p({X}_{i,t};\hat{\beta })\) was ≥0.75. Fifth, we assumed consistency, such that each commune’s potential outcome under an intervention equalled its observed outcome when the intervention was implemented similarly in both trial and target settings. We expect this to hold because the trial was designed as a pragmatic trial with MDA delivered as close to routine programmatic conditions as possible. Consistency also requires no interference between clusters; we expect this assumption to hold based on the ≥2.5-km buffer zone between trial clusters⁴. Sixth, we assumed correct model specification. To mitigate violations, we used a doubly robust estimator, which yields consistent estimates if either the trial participation model or the outcome model is correctly specified.

As a sensitivity analysis, we evaluated correlations between MDA coverage and effect modifiers to assess whether the variables used in the transport model could explain heterogeneity in coverage across trial clusters. We calculated Spearman correlation coefficients between coverage and each effect modifier across intervention clusters.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.