Slow slip events (SSEs) at subduction zone plate boundaries sometimes trigger earthquake swarms and megathrust earthquakes. The causal relationship between SSEs and seismicity has been studied worldwide, but the epidemic-type aftershock-sequence (ETAS) model, which is a standard statistical model of seismicity, does not explicitly consider the seismicity-triggering effect of SSEs. Therefore, if an SSE occurs at a plate boundary, probabilistic earthquake forecasts based on the ETAS model fail to predict observed seismicity. Here, we constructed a statistical model named the SSE-modulated ETAS model by incorporating SSE moment rates estimated from observation data from the global navigation satellite system into the original ETAS model. Our model assumes a linear or power-law relationship between the SSE moment rates and seismicity rates and estimates its proportionality constant as a new ETAS parameter. We applied this new model to three SSEs and M 2.5 or greater earthquakes in the shallow part of the Hikurangi Trench, New Zealand. The results show that it is better than the original ETAS model, giving a significant reduction in the Akaike information criterion. In addition, we examined the functional forms (e.g., lag time and power exponent) of the equation relating the moment rate of the SSEs to the seismicity rate. The results imply that, in addition to SSE-induced stress changes, crustal fluid migration may be related to SSE-induced seismicity. We also examine the influence of SSEs on aftershock productivity. Our model can improve short-term forecasts of seismicity associated with SSEs and is useful for quantifying characteristics of the seismicity.