Empirical dynamic modeling (EDM) 164 EDM infers a system's mechanistic underpinnings and predicts its dynamics using 165 time series data of one or more variables to construct an attractor in state space 166 ( Figure S1 ). This procedure is called univariate (using lagged versions of a single To test for nonlinear state dependence of a variable-the motivation behind EDM-185 we used the S-map test for nonlinearity (Sugihara 1994) (Figures S2b, c and S5;  186 Supporting Information). 187
We used an EDM approach called convergent cross-mapping (CCM) (Sugihara et al. 189 2012) to identify drivers of dengue incidence. If two variables are causally related, 190 then a multivariate attractor-where each variable in the system represents a 191 dimension that traces the dynamics of the system-can be reconstructed (up to a 192 practical limit) using lagged versions of just one of the variables ( Figure S1 ). Based 193 on Takens' Theorem, this univariate "shadow attractor" preserves the structural and 194 dynamic properties of the original multivariate attractor (Takens 1981; Sugihara et 195 al. 2012 ). The concept behind CCM is that if temperature causes dengue incidence, 196 then information about past temperature will be embedded in the dynamics of 197 dengue, such that the shadow attractor produced using only incidence data allows 198 us to accurately reconstruct temperature in the past. However, the converse 199 scenario would not be true: since dengue does not cause temperature, the shadow 200 attractor constructed using temperature data should not contain information to 201 accurately reconstruct past dengue incidence (Supporting Information). 202
The critical criterion for estimating causal (directional) associations between 203 two variables using CCM is checking that the cross-mapping skill (i.e., Pearson's 204 correlation coefficient, ρ, between predicted driver values using the univariate SSR 205 of the response variable, and the observed driver values) monotonically increases 206 12 anomalies (Deyle et al. 2016a). In the second, more conservative "Ebisuzaki" null 229 model, we conserved any periodicity (beyond seasonal) and randomized the phases 230 of Fourier-transformed time series (Ebisuzaki 1997). We tested for statistically 231 significant differences in cross-mapping skill between the model that used the data 232 versus the null models by performing Kolmogorov-Smirnov (K-S) tests after 233 convergence. 234
We also repeated CCM in the nonsensical, reverse-causal direction (e.g., to 235 test whether incidence drives climate) as a control for potential spurious 236 relationships generated by non-causal covariation (e.g., due to seasonality). 
We examined the predictive power of the drivers on dengue incidence by assessing 242 how well we can predict dengue dynamics using temperature, rainfall, susceptibles 243 index, and their combined effects. We used a combination of univariate SSR (i.e., 244 with incidence data) and multivariate SSR to build forecasting models and to 245 determine the improvement of forecasting using simplex projection when including 246 (Supporting Information). We built the SSR forecasting models/attractors using the We assessed model forecasting performance using leave-one-out cross-validation. 250
Next, we evaluated out-of-sample forecasting performance of these models 251 using testing data from four additional seasons (2009/2010-2012/2013 we followed the procedure as directed in the challenge (Supporting Information).
In nonlinear systems, drivers generally have an effect that is state-dependent: the 268 strength and direction of the effect depends on the current state of the system. 269
Scenario exploration with multivariate EDM allowed us to assess the effect of a 270 small change in temperature or rainfall on dengue incidence, across different states 271 of the system. The outcome of these small changes allowed us to deduce the 272 relationship between each climate driver and dengue incidence and how they 273 depend on the system state. For each time step t we used S-maps (Sugihara 1994; 274 Deyle et al. 2016a) to predict dengue incidence using a small increase (+ΔX) and a 275 small decrease (-ΔX) of the observed value of driver ( ) (temperature or rainfall). 276
For each putative climate driver, the difference in dengue predictions between these 277 small changes is Δ = ( + 1) J ( ) + for both temperature and rainfall to recover their approximate relationships with 281 dengue incidence at different states of the system. Scenario exploration analyses 282
were repeated across several model parameterizations to address potential
EDM showed that temperature, rainfall, and the susceptibles index drive dengue 287 incidence since the convergence criterion was met (Kendall's > 0, P < 0.01) in all 288 three cases (Figure 3a We cannot rule out the possibility that the apparent forcing of temperature on 295 dengue is due to a seasonal confounder. However, if no such confounder exists, then 296 the seasonal trend in temperature, which accounts for most temperature variation 297 in San Juan, drives the seasonal trend observed in dengue incidence. Compared to 298 the other drivers, the converging cross-mapping skill of the temperature null 299 models were relatively high (Figures 3 and S8 As expected, EDM tests for putative causality in the nonsensical directions-304 incidence driving temperature or rainfall-were not significant (i.e., no convergence; Figure S7 , black lines). This result further supports the finding that temperature and rainfall drive dengue incidence, because their causal relationships 307
were not confounded by spurious bidirectionality. The null models for the 308 nonsensical directions of causality ( Figure S7 , grey lines) also displayed no 309 convergence (completely flat), as expected (i.e., seasonality of dengue incidence 310 does not drive seasonality of temperature or rainfall). However, seasonality (or any 311 periodicity) of temperature, rainfall and susceptibles index drive dengue dynamics, 312
shown by convergence of the seasonal and Ebisuzaki null models (grey lines in 313
The multivariate SSR model using only temperature and rainfall data did not predict 316 predicting peak incidence, peak week, and seasonal incidence for all seasons on 344 average (Tables S1-S2, Figures S9-S12 ). This demonstrates the benefit of the EDM 345 approach for capturing the mechanistic, nonlinear, interdependent relationships 346 among drivers over both equation-based mechanistic models and phenomenological 347 models. 348
As expected, we find state-dependent effects of temperature and rainfall with non-350 zero median effects. We found that temperature had a small positive median effect 351
(2.88 cases/°C, Wilcox P < 0.001) on dengue incidence (Figure 5a) finding, since evidence of climate functional responses for disease dynamics is rare 386 due to the difficulty of obtaining appropriately informative field data. It is possible 387 that if we had temperature data ranging across a larger spectrum-possibly by 388 assembling data across multiple climates-that the empirical functional response 389 derived from EDM would also look unimodal. Further, when the susceptibles index 390 was high, the slope of the relationship between rainfall and dengue incidence 391 became more negative as rainfall increased, suggesting a concave-down effect of rainfall on incidence (Figure 6g, h) We found that rainfall, susceptible availability, and plausibly temperature 408 (via its seasonality) interact to drive dengue incidence. Combined, these three 409 drivers predicted dengue incidence with high accuracy (Figure 4c Even when accounting for susceptible availability, the effects of temperature 461 and rainfall on dengue were strongly state-dependent (Figure 6d, g) . This result is 462 potentially due to nonlinear effects of each climate driver (Figure 6e 