Urban heat vulnerability (UHV) caused by anthropogenic activities and climate changes has given rise to heat health issues in urban areas worldwide. Previous studies have extensively revealed a simple linear relationship between heat vulnerability indices (HVIs) and morbidity or mortality of heat-related illnesses, but the nonlinear relationship and interactions between main HVIs have not yet been fully explored. Based on vulnerability assessment framework, this paper selected fifteen indicators from built environment, sociodemographic and socioeconomic attributes, resource accessibility and residential thermal comfort, obtained from multisource data. Through the evaluation and analysis of composite HVI and its dimensions, we found that Qingshan district and East Lake scenic area contain more high to very high heat vulnerability communities. The performances of the ordinary least squares (OLS) and gradient boosting decision trees (GBDT) were compared, and results indicate GBDT outperforms the OLS model and captures the nonlinear relationship more efficiently in study areas with higher accuracy. When analyzing HVIs’ contributions and interactions with the GBDT model and the SHAP algorithm, nighttime light (NTL), building year (BY), PM2.5, floor area ratio (FAR), number of elderly (≥65 years) (NE) and urban surface roughness (USR) are six key indicators of morbidity of heat-related diseases (mean SHAP value>2.5), and they have an evident nonlinear relationship with the threshold effect and spatially heterogeneous contributions for the morbidity variation of heat-related diseases. Our study provides insights into machine learning (ML) model for the effect of heat vulnerability on city residential health and mitigation and adaptation strategies for governments and urban planners to develop heat resilience cities.