Malaria prevalence data
We used a recently published database of P. falciparum prevalence in sub-Saharan Africa2. This compendium, compiled by Snow et al. over more than two decades, is one of the most spatially and temporally complete publicly available databases of infectious disease burden. The database covers the period from 1900 to 2016, although sampling has increased substantially since the turn of the century (pre-2000: n = 32,533; post-2000: n = 17,892). Most prevalence surveys used microscopy for diagnostics (n = 36,805) but a substantial portion of data also derive from rapid diagnostic tests (n = 11,154). The data have been compiled from a mix of archival research through public health documents, including the records of colonial governments and elimination campaigns from different periods; national survey data; electronic records published in peer-reviewed journals and grey data sources (for example, World Health Organization technical documents); and a mix of other sources compiled by international organizations. Records were georeferenced in the original study using a standard set of protocols, with a 5-km grid uncertainty threshold for point data, and broader areas stored as administrative polygons. In total, the data include a total of 50,425 prevalence surveys at a total of 36,966 unique georeferenced locations.
The Snow et al. data cover all available prevalence surveys, including all age ranges, but were converted by the authors of the original study to a standardized estimate of prevalence in children 2–10 years of age (PfPR2−10), using a catalytic conversion Muench model. We chose to use these standardized estimates of childhood malaria prevalence because falciparum malaria has the highest mortality in children and pregnant women. The trends that we infer should generally be representative of broader transmission across age groups. In some cases, we note that declines in early-life exposure can lead to increases in incidence in adults61; however, these impacts are likely to be small, particularly given that active and passive improvements in malaria prevention, control and treatment much more directly determine trends in adult malaria risk.
Climate data
We used two sets of climate data in this study. The first is an observational dataset from the Climatic Research Unit (CRU-TS; v4.03 for model training and bias correction), which is constructed from monthly observations from extensive networks of meteorological stations from around the globe62. CRU-TS provides land-only climatic variables at a spatial resolution of 0.5° × 0.5° extending from 1901 to present (although our analysis is limited to the period 1901–2016). The second set of data is from ten global climate models (GCMs) selected from the sixth phase of the Coupled Model Intercomparison Project (CMIP6): ACCESS-CM2, ACCESS-ESM1-5, BCC-CSM2-MR, CanESM5, FGOALS-g3, GFDL-ESM4, IPSL-CM6A-LR, MIROC6, MRI-ESM2-0 and NorESM2-LM. In our historical analysis, we analysed (per GCM) one model realization of the ‘historical’ simulation, which includes anthropogenic greenhouse gas emissions, and one realization from the ‘historical-natural’ simulation, which includes only solar and volcanic climate forcing. For both the historical and historical-natural (hereafter and in the main text, ‘historical climate’ and ‘historical counterfactual’, respectively) simulations, we analysed the period 1901–2014.
To investigate the continued effect of climate change on malaria prevalence between 2015 and 2100, we analysed three CMIP6 future climate change simulations from each of the 10 GCMs. SSPs refer to the level of potential future global development (social, economic and technological) and the implication for climate change mitigation and/or adaptation actions or policy63,64. SSPs are combined with various possible future radiative forcings (RCPs) to form the climate change scenarios used in CMIP6. Of the available SSP–RCP scenarios, we selected and used three. The first two suggest enhanced human development outcomes with increased potential towards a more sustainable (SSP1)65 or a less sustainable (SSP5)66 economy. The third, SSP2 (ref. 67), is a mid-way scenario, which assumes a future that mostly follows historical trends64. We selected these scenarios in combination with a low (SSP1–RCP2.6), intermediate (SSP2–RCP4.5) and high (SSP5–RCP8.5) greenhouse gas concentration scenario.
We applied a standard quantile–quantile (Q–Q) bias-correction68,69 to the CMIP6 precipitation and temperature datasets for both of the historical simulations for the period 1901–2014, and all three future simulations for the period 2015–2100. Before the bias correction, we first remapped all simulated CMIP6 precipitation and temperature datasets to the same grid cell size (0.5° × 0.5°) as the CRU-TS observation data. We then performed for each CMIP6 model, the Q–Q bias correction at each grid point by mapping the quantile values (qi) for the empirical cumulative distribution functions for each of the 12 months over the period 1901–2014 (for each grid point) onto the corresponding quantiles in the observational dataset (CRU-TS), so that the observed precipitation or temperature values associated with qi become the bias-corrected value in the simulations. For the counterfactual (and future) simulations, we first determined, at each grid point, for each value of precipitation or temperature (for each month) over the period 1901–2014 (2015–2100), the equivalent quantile (qj) in the factual simulation and then identified the precipitation or temperature value associated with qj in the observational dataset as the bias-corrected value. We detrended both precipitation and temperature datasets before applying the bias-correction procedure, and then added the trends back after69.
Spatial data aggregation
Our statistical analysis is designed to isolate variation in the weather that is uncorrelated with other socioeconomic and/or environmental factors that influence malaria prevalence. As detailed in the next section, we build on a large body of climate econometrics research38,40,70 to do so, estimating a model that leverages variation over time in weather conditions within the same location. To estimate such a model, we required observations of malaria prevalence covering the same region in multiple time periods. By contrast, the raw prevalence data that we obtained from ref. 2 are point data observations from individual surveys conducted at different times, such that single geolocations are not observed repeatedly over time. Therefore, we aggregated the point-level data from ref. 2 by averaging PfPR2−10 observations to the first administrative level within each country (that is, state or province level, or as shorthand, ADM1), using shapefiles provided by the Database of Global Administrative Areas dataset v3.6 (www.gadm.org). This level of aggregation provides sufficient granularity to capture differences in climate impacts within countries and to control for local heterogeneity in confounders, while ensuring sufficient data coverage within these units. This aggregation scale has also been conducted in previous work that models this dataset at the same spatial resolution2. For robustness, we also show results from a statistical model that does not aggregate data, and instead uses the prevalence data at its native resolution (see below for details).
To compute average prevalence values at the scale of ADM1, we used an unweighted arithmetic mean over all prevalence surveys observed in the corresponding ADM1 month. This approach imposes minimal assumptions on the spatiotemporal process of malaria transmission and requires no additional high-resolution data (for example, population) for use as weights, which are unavailable for sub-Saharan Africa for years as early as 1901. Although previous work aiming to construct comprehensive high-resolution estimates of health outcomes using point data often uses spatiotemporal smoothing methods (for example, ref. 71), doing so here would artificially introduce spatial and temporal correlations that could bias recovered regression coefficients and threaten inference72.
We similarly aggregate monthly 0.5° grid-level weather data (from all CRU-TS and CMIP6 models) to the ADM1-month level. To do so without introducing aggregation biases, we applied methods from previous research demonstrating that it is possible to statistically recover nonlinear relationships that take place at high spatial and temporal resolution, even when the resolution of available outcome data is relatively coarse (that is, ADM1-month-level average malaria prevalence)40,73,74,75,76. In our setting, this is achieved by computing nonlinear polynomial transformations of temperature at the grid-cell-by-month level before aggregating these values across administrative units. Such an approach ensures the temperature variables used for estimation reflect the full distribution of temperatures experienced across administrative regions of varying sizes and terrains. For example, many of the 12 ADM1 regions in Ethiopia include both hot low-elevation zones and cold highlands, such that temperatures can vary substantially within an ADM1 during the same month. We computed second-order polynomials at each grid cell before aggregating across such diverse landscapes to ensure the regressor variables capture both extreme cold and extreme heat, even when they occur simultaneously within the boundaries of ADM1.
To see this method in practice, let PfPRgit denote average malaria prevalence in children 2–10 years of age in grid cell g located within administrative unit i during month t and let Tgit indicate temperature observed at the same spatiotemporal scale. Following previous studies recovering local-level quadratic responses between malaria prevalence and temperature77,78, we assumed that prevalence in grid cell g in month t is a quadratic function of the temperature experienced in that same grid cell and month (noting that we show results relaxing this assumption, such as other nonlinear functional forms and the possibility of temporal lags):
$$Pf{{\rm{P}}{\rm{R}}}_{git}={\beta }_{1}{T}_{git}+{\beta }_{2}{T}_{git}^{2},$$
(1)
where β1 and β2 are constant average coefficients. As discussed above, we cannot empirically estimate a model like equation (1) reliably because we did not observe prevalence over multiple months t for the same grid cell g. Instead, our empirical specification relies on average ADM1-month-level prevalence variables PfPRit. Thus, we must aggregate equation (1) in a manner that allows us to recover the same β1 and β2 coefficients that describe the local-level temperature response, and which we would have recovered had we been able to estimate equation (1) directly. Specifically, average ADM1-month prevalence can be written as:
$$\begin{array}{r}Pf{{\rm{PR}}}_{it}=\sum _{g\in i}Pf{{\rm{PR}}}_{git}{\omega }_{gi}=\sum _{g\in i}({\beta }_{1}{T}_{git}+{\beta }_{2}{T}_{git}^{2}){\omega }_{gi}\\ \,=\,{\beta }_{1}\sum _{g\in i}{T}_{git}{w}_{gi}+{\beta }_{2}\sum _{g\in i}{T}_{git}^{2}{w}_{gi},\end{array}$$
(2)
where ωgi denotes a grid-by-ADM1 weight. In our setting, we estimated an area-weighted average prevalence value by setting ωgi equal to the share of administrative unit i’s area that falls into grid cell g, as the lack of high-resolution population data make population weighting infeasible.
Equation (2) shows that a regression of average prevalence in ADM1 unit i and month t on variables that are ADM1-month weighted aggregates of the nonlinear temperature terms Tgit and \({T}_{git}^{2}\) will, in expectation, recover the same coefficients β1 and β2 that describe the fundamental grid-level relationship described by equation (1). This same procedure has been used to estimate the relationship between: monthly administrative-level dengue incidence and daily grid-level temperature76; annual all-cause mortality and daily grid-level temperature75,79; annual country-level crop yields and daily grid-level soil moisture80; among many other examples. As in these other cases, our approach mitigates aggregation bias up to the level of the grid resolution of the climate data. Although climate and prevalence probably vary across space within each grid cell, we cannot resolve such dynamics here, given the lack of reliable higher-resolution weather data in Africa over the extended time frame of our analysis62.
We note that one could, alternatively, construct nonlinear weather variables after aggregating to the administrative unit, thus estimating:
$$Pf{\mathrm{PR}}_{it}={\widetilde{\beta }}_{1}\sum _{g\in i}{T}_{git}{\omega }_{gi}+{\widetilde{\beta }}_{2}{(\sum _{g\in i}{T}_{git}{\omega }_{gi})}^{2},$$
(3)
where recovered parameters \({\widetilde{\beta }}_{1}\) and \({\widetilde{\beta }}_{2}\) are biased relative to the fundamental relationship in equation (1) because \({({\sum }_{g\in i}{T}_{git}{\omega }_{gi})}^{2} < {\sum }_{g\in i}{T}_{git}^{2}{\omega }_{gi}\), due to Jensen’s inequality.
Thus, we constructed a vector of ADM1-month temperature variables by computing nonlinear transformations at the grid level before aggregating across space, following equation (2). Although we could, in principle, follow the same procedure for precipitation, we instead computed drought and flood variables at the ADM1-month level, due to high rates of mismeasurement in grid-level rainfall estimates62 and due to the likelihood that prevalence–precipitation relationships occur over larger spatial scales than a single grid cell (that is, water flows through hydrological systems linking precipitation in one location to water availability and prevalence downstream). We show in a robustness exercise detailed below that using grid-level precipitation data generates similar results as our aggregated model, but increases uncertainty, consistent with evidence on rainfall mismeasurement in Africa.
All main results rely on the spatial aggregation procedure described above. However, we additionally estimated an alternative model that leverages the point-level prevalence data directly, introducing no spatial aggregation beyond the resolution of the weather data (0.5°). As we detail below and show in Supplementary Fig. 9, the recovered prevalence–temperature results are very similar using these two distinct methods.
Statistical model
The influence of climatic conditions on malaria prevalence has been heavily studied using transmission models based in vector ecophysiology and calibrated using laboratory experiments3,4. The important benefit of this approach is that the mechanistic links between a particular environmental condition (for example, temperature) and malaria prevalence in the human population, such as effects on biting rate and survival probability, can be independently isolated. However, this approach is limited in its ability to generalize to real-world contexts, in which complex socioeconomic factors interact with modelled relationships based on laboratory conditions. Clinical data, which measures malaria prevalence in human populations, have been used to validate modelled results3, but inconsistent findings arise due to challenges in statistically isolating the role of climate from the many correlated factors influencing prevalence, such as public health interventions, drug resistance, conflict and social instability, and economic shocks15,52,81,82,83.
This study seeks to provide generalizable population-scale evidence of the malaria–climate link across sub-Saharan Africa using field-collected clinical data and a statistical approach designed to isolate changing environmental conditions from spatiotemporal confounding factors. Specifically, we drew on the climate econometrics literature40, which has developed causal inference approaches to quantify and project the effects of anthropogenic climate change on a host of socioeconomic outcomes, from agricultural yields73, to civil conflict84, to all-cause mortality79. This approach is designed to approximate controlled experiments by semi-parametrically accounting for unobservable spatial and temporal confounding factors, isolating variation in the climate system that is less likely to be correlated with other socioeconomic factors85. This approach is often referred to as ‘reduced form’, as it allows for a plausibly causal interpretation of recovered relationships between socioeconomic conditions and the climate, but it does not easily enable the researcher to isolate individual mechanisms linking a changing climate to shifts in outcomes (for example, mosquito population dynamics or parasite development rates). However, causal estimates enable counterfactual simulation in which climate is changed and all other factors are held constant; this is the exercise conducted here and in many applications of climate econometric frameworks, including estimating the effects of climate change on dengue cases76, international human migration86, all-cause mortality75,79 and more. Moreover, these relationships can be used to calibrate more structured transmission models by providing empirical grounding from observational data.
We developed a statistical model using monthly survey-based malaria PfPR2−10 covering all of sub-Saharan Africa over 116 years. Our outcome variable is the average prevalence for each ADM1 i (for example, province or state) in country c during month–year t, which we denote as PfPRit. We estimated prevalence as a flexible function of monthly temperature and precipitation variables as follows:
$$\begin{array}{l}Pf{{\rm{PR}}}_{it}={\beta }_{1}\sum _{g\in i}{T}_{git}{\omega }_{gi}+{\beta }_{2}\sum _{g\in i}{T}_{git}^{2}{\omega }_{gi}\\ \,\,\,\,+\mathop{\sum }\limits_{{\ell }=0}^{L}{\rho }_{{\ell }}{\mathbb{1}}\{{{\rm{drought}}}_{i,t-{\ell }}\}+\mathop{\sum }\limits_{{\ell }=0}^{L}{\psi }_{{\ell }}{\mathbb{1}}\{{{\rm{flood}}}_{i,t-{\ell }}\}\\ \,\,\,\,+\,{\alpha }_{i}+{\gamma }_{rm}+\sum _{c\in C}[{\phi }_{1c}t+{\phi }_{2c}{t}^{2}]\\ \,\,\,\,+\,{\delta }_{1}{\mathbb{1}}{\{{\rm{intervention\; 1}}\}}_{t}+{\delta }_{2}{\mathbb{1}}{\{{\rm{intervention\; 2}}\}}_{t}+{\varepsilon }_{it},\end{array}$$
(4)
where g subscripts denote grid cells, which fall within administrative units i, and ωgi are area weights equal to the share of unit i’s area covered by grid cell g, such that ∑g∈iTgitωgi equals the area-weighted average monthly temperature across all grid cells falling within administrative unit i. Together, parameters β1 and β2 recover a quadratic response between prevalence and monthly average temperature. As described above, polynomials are computed before aggregating across grid cells to preserve local nonlinearities and avoid aggregation bias. Precipitation extremes are captured by a vector of dummy variables \({\mathbb{1}}\{{{\rm{drought}}}_{i,t-{\ell }}\}\) and \({\mathbb{1}}\{{{\rm{flood}}}_{i,t-{\ell }}\}\), which indicate whether an administrative unit’s monthly rainfall total can be categorized as drought (defined as 10% or lower of the long-run location-specific and month-specific mean) or flood (defined as 90% or higher of the long-run location-specific and month-specific mean) during month–year t − ℓ. We allowed for up to 3 months of lags (that is, L = 3) for these extreme precipitation conditions in our main specification, based on hypotheses from previous literature regarding the timescales of larvae drying and of ‘flushing’42,43. Various sensitivity analyses detailed below demonstrate that key findings are robust to: including lags for temperature as well as precipitation (Extended Data Fig. 2); the drought and flood cut-offs used for precipitation (Supplementary Figs. 2, 5 and 6); alternative functional forms of temperature (Supplementary Fig. 4); and the estimation of a grid-level regression that does not aggregate prevalence or weather across space (Supplementary Fig. 9).
Equation (4) uses a suite of semi-parametric spatiotemporal controls to isolate variation in climatological conditions that is independent from other disease transmission factors, following standard practices in the climate econometrics literature38,40. First, αi is a vector of indicator variables for each of 853 ADM1 units across our multi-country sample. These spatial ‘fixed effects’ control for all time-invariant characteristics of an administrative unit that may confound the relationship between temperature, rainfall and prevalence. For example, higher-altitude regions may exhibit cooler temperatures, but they also may be more geographically isolated communities with limited access to malaria prevention interventions. By controlling for mean conditions in each location, these spatial fixed effects avoid conflating climate conditions with other geographical correlates.
Second, γrm is a vector of region-by-month-of-year indicator variables, where regions r are defined using the Global Burden of Disease regional definitions of West, southern, central and East Africa (see figure 2 in ref. 87). Note that the subscript m indicates month of the year (for example, February), whereas the month–year index t indicates the month–year time index (for example, February, 1998). These spatiotemporal fixed effects γrm account for region-specific seasonality in prevalence that may spuriously relate to seasonally varying climatological conditions. We allowed these seasonal controls to vary by region because of large differences in climatological seasonality and in malaria cyclicality across sub-Saharan Africa88, and we show below that our main findings are robust to more stringent seasonality controls defined at the country level (Extended Data Fig. 4). Third, ϕ1c and ϕ2c are coefficients estimating, for each country c in the full set of countries C, a nonlinear, country-specific quadratic in the month–year time index, which adjusts the regression for country-specific gradual trends that may confound the malaria–climate relationship, particularly under historical conditions of anthropogenic climate change. Extended Data Fig. 4 shows that our results are robust to multiple alternative approaches to controlling for long-run trends that may vary across space.
Finally, the indicator variables \({\mathbb{1}}{\{{\rm{intervention\; 1}}\}}_{t}\) and \({\mathbb{1}}{\{{\rm{intervention\; 2}}\}}_{t}\) are equal to one when an observation falls into the 1955–1969 or 2000–2015 period, respectively. These two periods saw substantial malaria intervention programmes across the subcontinent, leading to considerable declines in malaria that were unrelated to changes in the climate2,89. These indicator variables control for shocks to prevalence during these two periods, and the coefficients δ1 and δ2 allow for differential effectiveness of the two distinct intervention periods. Although these variables are strongly correlated with average prevalence and the first is highly statistically significant (Extended Data Table 1), our main findings are robust to their exclusion (Extended Data Fig. 4).
Together, these set of flexible controls imply that the residual variation in temperature and precipitation events used to identify the coefficients β1, β2, ρℓ and ψℓ is month-to-month variation over time within the same location, after controlling for gradual country-specific trends, regional seasonality and the aggregate effects of two substantial malaria prevention intervention programmes.
We estimated equation (4) using the lfe package in R. In estimation, we clustered standard errors εimt at the country-by-5-year group level to account for spatial correlation within a country and serial correlation within a 5-year time span (see Supplementary Table 5 and the associated discussion on spatiotemporal structure in model residuals below for details on this choice). When computing probabilistic historical and future climate change simulations, we repeatedly resampled coefficients from the clustered variance–covariance matrix so that this same spatial and temporal correlation was accounted for when computing estimates of the impacts of climate change. These resampled draws of both temperature and precipitation coefficients are plotted in Fig. 2. We additionally show sensitivity to alternative methods of capturing uncertainty in Supplementary Table 3 and Supplementary Fig. 7. In Supplementary Fig. 10, we show that model residuals are close to normally distributed, although with slightly heavier tails, making the application of ordinary least squares appropriate in this context.
Statistical model robustness
In this section, we describe a set of model sensitivity analyses that probe the robustness of our empirical model. Specifically, we investigated the sensitivity of our key findings to: alternative spatiotemporal controls; inclusion of dynamic temperature effects; alternative definitions of extreme rainfall events; alternative functional forms for the prevalence–temperature relationship; and the estimation of a survey-level regression in which the point-level nature of the raw prevalence data are used directly, with minimal spatial aggregation. Finally, we have provided a set of diagnostics investigating the spatiotemporal structure of our model residuals.
Spatiotemporal controls
Our preferred empirical specification in equation (4) includes ADM1 fixed effects (that is, indicator variables), region-by-month-of-year fixed effects, country-specific quadratic time trends and two indicator variables for each of two malaria intervention periods (1955–1969 and 2000–2015). Extended Data Fig. 4 shows that our estimated prevalence–temperature relationship is highly robust to many alternative spatial and temporal controls. All panels in this figure include ADM1 fixed effects to control for time-invariant characteristics that may confound the relationship between prevalence and temperature, but each panel varies in the additional spatial and/or temporal controls included in the regression. A tabular version of these results is shown in Extended Data Table 1. Although the temperature at which prevalence peaks changes slightly across model specifications, it remains within a degree of the 24.9 °C value from our preferred specification for most models, particularly those including time trends that are spatially differentiated (note that peak temperatures indicated in Extended Data Fig. 4 are rounded to the nearest degree for display purposes). Predictably, stringent controls, such as region-by-year and country-by-month fixed effects, tend to increase statistical uncertainty. The specification without an expected inverted U shape includes country-by-year fixed effects, which absorb nearly all residual variation in temperature. However, overall, the estimated shape and magnitude of the prevalence–temperature relationship remain robust to alternative spatial and temporal controls.
Dynamic temperature effects
Our preferred empirical specification estimates contemporaneous (within 1 month) and lagged (up to 3 months) effects of extreme rainfall on malaria prevalence, but only contemporaneous effects of temperature. Although it is possible that temperature also exhibits lagged effects, we show in Extended Data Fig. 2 that the cumulative effect of temperature on PfPR2−10 is similar whether 0, 1, 2 or 3 months of lagged temperatures are accounted for. The prevalence response to temperature does become stronger with 3 months of lags, suggesting that our historical and future climate predictions shown throughout the main text may be somewhat conservative. However, overall, these findings suggest that climate change impact predictions are unlikely to change meaningfully under different assumptions of the lag structure of temperature exposure.
Definitions of extreme rainfall events
Our main empirical specification defines drought as months for which total precipitation is less than or equal to 10% of the long-run location-specific and month-specific mean. Flood is analogously defined as months for which total precipitation is greater than or equal to 90% of the long-run location-specific and month-specific mean. Here we investigated the sensitivity of our main findings to these definitions. To do so, we systematically varied both the drought and flood cut-off values, ranging from less than 1% to less than 20% for drought and from more than 85% to more than 95% for flood, respectively. Supplementary Fig. 6 shows that the relationship between malaria prevalence and temperature is insensitive to the definition of drought and flood events. Supplementary Fig. 5 shows that under most drought and flood definitions, extremely low precipitation events have a negative effect on prevalence with a lag of 1–2 months. However, this effect is rarely statistically significant. Supplementary Fig. 2 shows that extremely high rainfall events increase prevalence with a lag of 2–3 months, a result that is statistically significant and generally robust to alternative drought and flood definitions. In general, these sensitivity analyses show that our main findings are not sensitive to the specific definitions of drought and flood used in estimation of equation (4).
Temperature’s functional form
Following from theoretical and laboratory-based literature (for example, refs. 3,4), we modelled the prevalence–temperature relationship as quadratic. However, Supplementary Fig. 4 shows that this relationship is similar when more flexible functional forms are used. In particular, the temperature at which prevalence peaks changes little when higher-order polynomials are estimated. Estimating higher-order polynomials increases uncertainty, particularly in the tails of the temperature distribution, but point estimates are similar across the majority of the observed temperature range.
Estimating a grid-level regression
As described above, our main analysis relies on an ADM1-level regression in which malaria prevalence and weather data are aggregated from higher spatial resolutions to the ADM1 scale. Here we show the results from an alternative approach, in which point-level malaria prevalence survey observations are matched to corresponding 0.5° resolution CRU climate data grid cells and the regression is estimated at this grid level.
Using these disaggregated data, we estimated a regression model analogous to equation (4), but modified to fit the spatial scale of the data. Specifically, we estimated:
$$\begin{array}{l}Pf{{\rm{PR}}}_{git}={\eta }_{1}{T}_{git}+{\eta }_{2}{T}_{git}^{2}\\ \,\,\,\,+\mathop{\sum }\limits_{{\ell }=0}^{L}{\lambda }_{{\ell }}{\mathbb{1}}\{{{\rm{drought}}}_{gi,t-{\ell }}\}+\mathop{\sum }\limits_{{\ell }=0}^{L}{\xi }_{{\ell }}{\mathbb{1}}\{{{\rm{flood}}}_{gi,t-{\ell }}\}\\ \,\,\,\,+{\alpha }_{i}+{\gamma }_{rm}+\sum _{c\in C}[{\phi }_{1c}t+{\phi }_{2c}{t}^{2}]\\ \,\,\,\,+{\delta }_{1}{\mathbb{1}}{\{{\rm{intervention\; 1}}\}}_{t}+{\delta }_{2}{\mathbb{1}}{\{{\rm{intervention\; 2}}\}}_{t}+{\varepsilon }_{git},\end{array}$$
(5)
where all variables are defined as above for equation (4). In particular, PfPRgit represents average prevalence for all surveys located in grid g falling within ADM1 unit i during month t, Tgit denotes temperature in the same grid and month, and precipitation extremes are captured via dummy variables \({\mathbb{1}}\{{\mathrm{drought}}_{gi,t{\ell }}\}\) and \({\mathbb{1}}\{{\mathrm{flood}}_{gi,t{\ell }}\}\) that indicate when the monthly rainfall total of each grid cell is less than 10% of its long-run month-specific mean (drought) or more than 90% (flood). As for the main model, when estimating equation (5), we clustered standard errors at the country-by-5-year group level.
Two features render this estimating equation distinct from the main analysis. First, temperature and precipitation variables are matched exactly to the grid cell within which the malaria survey was conducted, such that no aggregation is necessary. This increases the precision of the match between weather and outcome variables. Second, the spatial ‘fixed effects’ denoted by αi—that is, indicator variables for each of the ADM1 units in our sample—are estimated at a lower spatial resolution (ADM1) than the data itself (grid). This implies that the weather variation used to identify coefficient vectors η, λ and ξ includes both variation over time within an ADM1, but also across survey locations located within the same ADM1. Thus, although this model has the benefit of leveraging higher spatial resolution in weather, it potentially suffers from omitted variables bias, as weather conditions in different survey locations may be correlated with other unobservable determinants of malaria prevalence (for example, access to healthcare, rates of poverty and proximity to water bodies).
In Supplementary Fig. 9, we show that the malaria prevalence–temperature relationship recovered from estimation of equation (5) is very similar to that estimated from our main regression model in equation (4), suggesting that aggregation of the underlying survey data does not influence the key results of the paper. By contrast, the estimated drought and flood coefficients recovered from the grid-level regression are highly imprecise, consistent with substantial measurement error in local-level precipitation datasets in Africa for much of the twentieth century62.
Correlations in model residuals
Here we evaluated the extent to which our model residuals are correlated over space and time by calculating correlations between residuals across various subsets of our data, following similar tests in ref. 90. Specifically, we calculated correlations across observations that are: (1) within the same ADM1 unit but from different time periods; and (2) from different ADM1 units but within the same time period. For temporal correlations, we computed correlations between temporally consecutive observations within windows of up to 5 years, whereas for spatial correlations, we investigated correlations within countries, across countries, within Global Burden of Disease multi-country regions and based on physical distance (using ADM1 centroids).
Supplementary Table 5 reports mean correlations, as well as the first and third quartiles of the distribution of correlations across different pairs of units. Rows labelled ‘temporal’ quantified correlations over time, whereas rows labelled ‘spatial’ quantified correlations over space. We note that these tests should be interpreted with care, as the malaria prevalence data are highly unbalanced in space and time, often leaving few observations with which to estimate correlations and/or few regional or temporal pairs over which to summarize a distribution of correlations. To ensure interpretability, we restricted analysis to correlations with at least ten observations. The last column in Supplementary Table 5, labelled N, indicates the number of pairs for which sufficient data were available to construct correlations for a given grouping. These results reveal moderate serial correlation in residuals, especially for consecutive observations (row 2; mean ρ = 0.40 for consecutive month–years within the same ADM1 unit). Spatial correlations range from negligible (mean ρ = 0.04 across ADM1s from different countries or regions) to moderate (mean ρ = 0.28 and ρ = 0.30 for ADM1s within the same country and ADM1s with centroids less than 500 km of one another, respectively). Country boundaries are critical for determining spatial correlations: ADM1s with centroids less than 500 km of one another have high mean correlations within countries (mean ρ = 0.30), but low mean correlations when crossing country borders (mean ρ = 0.08).
Given these results, we clustered standard errors at the country-by-5-year group level, accounting for correlation in model residuals across all ADM1s within the same country and across all months within the same 5-year window. However, in Supplementary Table 3 and Supplementary Fig. 7, we show sensitivity of our recovered confidence intervals to alternative approaches to standard error estimation.
Predictions
In both historical and future simulations, we applied the estimated panel regression to calculate the effect of climate change on PfPR2−10. Our predictions capture the full range of statistical uncertainty (1,000 model estimates resampled from the clustered variance–covariance matrix) and climate model uncertainty (10 climate models), producing a total of 10,000 estimates of historical or future impacts in any given scenario. Each of these 10,000 estimates was normalized to a long-run baseline (past: 1901–1930; present: 2015–2020) before estimates are averaged, creating an estimate of climate change impacts relative to that baseline. In our historical analyses, we only used these models to estimate changes in prevalence attributable to climate change: although the panel regression model accounts for other historical drivers through the fixed effects structure, these are not the focus of our analysis, and so we choose not to estimate total prevalence including these effects. Similarly, we elected not to make assumptions about non-climate drivers of malaria prevalence in the future, and thus do not apply the model to predict future trends in overall prevalence.
For overall trends (for example, reported in Figs. 2d, 3d and 4d), we generated continent-wide averages or four regional averages using the unweighted average of estimates for each ADM1 unit. This is a deliberate oversimplification, as we did not adjust averages based on either ADM1 units’ land area or the estimated population they contain; we made this decision based on the challenges of reconstructing historical population density at fine scales, as well as the need to otherwise make assumptions about how disease burden is allocated over space (for example, the distribution of transmission across rural or urban areas). For similar reasons, we chose not to estimate the effect of prevalence changes on overall malaria incidence. Although some studies have attempted this using a linear conversion with total population91, proper estimation of incidence (and the effects of treatment variables, through prevalence, on case burden) requires malaria transmission models that require substantially more demographic assumptions89. Future work could explore both of these methodologically complex directions, and potentially generate finer-scale estimates of how many cases of childhood malaria, and resulting deaths, are attributable to climate change.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

