Multivariate “Bayesian” regression via a shared component model has gained popularity in recent years, particularly in modeling and mapping the risks associated with multiple diseases. This method integrates joint outcomes, fixed effects of covariates, and random effects involving spatial and temporal components and their interactions. A shared spatial–temporal component considers correlations between the joint outcomes. Notably, due to spatial–temporal variations, certain covariates may exhibit nonlinear effects, necessitating the use of semiparametric regression models. Sometimes, choropleth maps based on regional data that is aggregated by administrative regions do not adequately depict infectious disease transmission. To counteract this, we combine the area-to-point geostatistical model with inverse distance weighted (IDW) interpolation for high-resolution mapping based on areal data. Additionally, to develop an effective and efficient early warning system for controlling disease transmission, it is crucial to forecast disease risk for a future time. Our study focuses on developing a novel multivariate Bayesian semiparametric regression model for forecasting and mapping HIV and TB risk in West Java, Indonesia, at fine-scale resolution. This novel approach combines multivariate Bayesian semiparametric regression with geostatistical interpolation, utilizing population density and the Human Development Index (HDI) as risk factors. According to an examination of annual data from 2017 to 2021, HIV and TB consistently exhibit recognizable spatial patterns, validating the suitability of multivariate modeling. The multivariate Bayesian semiparametric model indicates significant linear effects of higher population density on elevating HIV and TB risks, whereas the impact of the HDI varies over time and space. Mapping of HIV and TB risks in 2022 using isopleth maps shows a clear HIV and TB transmission pattern in West Java, Indonesia.