BACKGROUND : Pediatric hand, foot, and mouth disease (HFMD) has generally been found to be associated with climate. However, knowledge about how this association varies spatiotemporally is very limited, especially when considering the influence of local socioeconomic conditions. This study aims to identify multi-sourced HFMD environmental factors and further quantify the spatiotemporal nonstationary effects of various climate factors on HFMD occurrence. METHODS : We propose an innovative method, named spatiotemporally varying coefficients (STVC) model, under the Bayesian hierarchical modeling framework, for exploring both spatial and temporal nonstationary effects in climate covariates, after controlling for socioeconomic effects. We use data of monthly county-level HFMD occurrence and data of related climate and socioeconomic variables in Sichuan, China from 2009 to 2011 for our experiments. RESULTS : Cross-validation experiments showed that the STVC model achieved the best average prediction accuracy (81.98%), compared with ordinary (68.27%), temporal (72.34%), spatial (75.99%) and spatiotemporal (77.60%) ecological models. The STVC model also outperformed these models in the Bayesian model evaluation. In this study, the STVC model was able to spatialize the risk indicator odds ratio (OR) into local ORs to represent spatial and temporal varying disease-climate relationships. We detected local temporal nonlinear seasonal trends and spatial hot spots for both disease occurrence and disease-climate associations over 36 months in Sichuan, China. Among the six representative climate variables, temperature (OR = 2.59), relative humidity (OR = 1.35), and wind speed (OR = 0.65) were not only overall related to the increase of HFMD occurrence, but also demonstrated spatiotemporal variations in their local associations with HFMD. CONCLUSIONS : Our findings show that county-level HFMD interventions may need to consider varying local-scale spatial and temporal disease-climate relationships. Our proposed Bayesian STVC model can capture spatiotemporal nonstationary exposure-response relationships for detailed exposure assessments and advanced risk mapping, and offers new insights to broader environmental science and spatial statistics.