Abstract: As an important member of the East Asian summer monsoon system, the western North Pacific subtropical high (WPSH) significantly affects precipitation in China and tropical cyclone activity in the western North Pacific. However, predicting the WPSH’s summer behavior remains challenging. Using Liang’s (2014) causality of information flow method, we have selected predictors for the dominant modes of the summer WPSH, utilizing monthly air/sea datasets and indices from 1979 to 2022 to develop a seasonal statistical model. We first applied empirical orthogonal function (EOF) analysis to identify the two EOF modes and their corresponding principal components (PCs) of the summer WPSH from 1979 to 2020. The first two dominant modes are closely related to the fading of El Niño events, while the latter correlate with the emergence of La Niña. To identify predictors for these patterns, we used information flow analysis and multiple linear stepwise regression during the model training period of 1980–2015. We identified seven predictors for each PC. The predictors for PC1 mainly include the Niño 4 index and sea surface temperatures (SSTs) in the northern equatorial Central Pacific. For PC2, predictors mainly include SSTs in the northern tropical Central Pacific and the Atlantic multidecadal oscillation (AMO) index, suggesting that La Niña development from the Pacific meridional mode (PMM) in spring and North Atlantic SSTs significantly affect EOF2. During the training period (1980–2015), the correlation between predicted and actual PC time series was 0.91 for PC1 and 0.88 for PC2. For the forecasting period (2016–2020), these correlations were 0.74 and 0.85, respectively. Over the 1980–2022 period, the regional average time correlation coefficient and the multiyear average pattern correlation coefficient of the reconstructed fields compared to observations were 0.63 and 0.48, both exceeding the 99% significance level. The forecast success largely depends on its limitations. We also conduct similar year forecasts based on PC values to address weak anomalies in the reconstructed fields.