Statistical Analysis of the Mechanisms of COVID-19 Spread Among Japan's Prefectures

1. Introduction

During the COVID-19 pandemic, Japan rapidly and extensively collected data on new infections and related information, providing insights into infection trends across its prefectures. Analyzing this data is crucial not only for understanding the dynamics of COVID-19 but also for devising strategies to prevent the spread of future infectious diseases. This paper focuses on analyzing daily new positive case data of COVID-19 in Japan’s prefectures.

The daily data on new positive cases comprises a trend component reflecting long-term fluctuations, a weekly variation component with recurring patterns, and short-term fluctuations influenced by temporary factors such as weather, public holidays, seasons, local events, and regional prevention measures. To address these complexities, this study first decomposes the daily time series data of new positive cases into a trend component, a weekly variation component, and a short-term component. Based on this decomposition, we aim to analyze the mechanisms driving the spread of COVID-19 across Japan and identify factors influencing these trends.

Time series decomposition methods, such as those proposed by [3], are commonly employed in such analyses. However, the uneven fluctuations in new positive case data present challenges for traditional methods like Kitagawa’s. To overcome these limitations, we adopt the extended moving linear model (EMLM) approach proposed by [4], an extension of the moving linear model (MLM) introduced by [5], and confirm its effectiveness.

Among the components of new positive case data, the weekly variation component is primarily driven by the periodicity of socio-economic activities and is unrelated to regional infection dynamics. Thus, our analysis emphasizes the trend and short-term fluctuation components. The trend component reveals regional connections in the spread of infections, while short-term fluctuations, influenced by complex and random factors, are harder to interpret. Nevertheless, public holidays, being nationwide, provide a unique opportunity to examine their impact on infection dynamics.

Recent studies on Japan’s COVID-19 pandemic, such as [9,10,11,12,13], have reported that human mobility accelerates infection spread and that both trend and periodic fluctuations strongly characterize case variations. For instance, [10] highlight the impact of government campaigns on travel demand during the pandemic, while [11] examines infection patterns before major events like the Tokyo Olympics. Similarly, studies by [7] focuses on global trends, with [14] demonstrating the effectiveness of the Prophet model for predicting cases in the United States. Notably, [12,13] forecasts the spread of COVID-19 in Japan using epidemiological models. However, our study differs in its focus on uncovering the mechanisms behind regional transmission patterns within Japan.

In addition, [2] reviewed deep learning methods for COVID-19 forecasting, identifying 53 quality studies from over 150 reviewed publications. Their work provides valuable insights into machine learning’s potential for outbreak prediction. [8] proposed the SVIHRD model for COVID-19, employing a physics-informed neural network (PINN) to estimate parameters using real data from Japan. While machine learning methods such as these are effective for predicting outbreak trajectories, their utility for post-event analysis remains limited. This limitation arises from their opaque nature, as they often fail to explain the underlying mechanisms, posing challenges for policymaking aimed at future epidemic prevention.

Recently, [15] analyzed the spread of COVID-19 across prefectures in Japan’s Kanto region by decomposing daily time series data and performing time-lag correlation analysis on trend components. Their findings suggest that fluctuations in Tokyo’s cases lead those in surrounding prefectures. Similarly, [6] extended this method to the Kansai region, showing that cases in Osaka and Kyoto precede those in nearby prefectures. However, time-lag correlation methods have limitations, such as inconsistency when applied broadly across all prefectures. For instance, if prefecture A leads prefecture B by n days and prefecture B leads prefecture C by m days, it does not necessarily follow that prefecture A leads prefecture C by $n+m$ days.

To address these limitations, this study first decomposes the daily time series data of new positive cases in all prefectures into trend, weekly variation, and short-term fluctuation components. Using the trend component, we analyze inter-prefecture transmission dynamics by defining Euclidean distances based on time-lagged trend data and estimating lags that minimize total distances. This approach ensures consistent results across all prefectures. We further analyze factors affecting these lag values and apply cluster analysis to the lag-adjusted trend data to examine transmission mechanisms. Additionally, we investigate the holiday effect on infection spread using the short-term fluctuation component. By studying all prefectures comprehensively, this research aims to provide a deeper understanding of the mechanisms driving COVID-19 spread in Japan.

The ultimate goal of this study is to analyze the factors influencing the spread of COVID-19 from diverse perspectives through post-event analysis, even if some precision is compromised. Using straightforward methods, we aim to clarify infection dynamics and provide valuable insights for future epidemic prevention policies.

The remainder of this paper is organized as follows: Section 2 presents the hypotheses and describes the preprocessed data. Section 3 introduces the analytical methods and procedures. Section 4 presents the analysis results, along with related discussions. Finally, Section 5 offers a discussion of the findings, conclusions, and suggestions for future research.

2. Hypotheses and Data Preprocessing

2.1. Hypotheses

Regarding the spread of COVID-19 in Japan, we establish the following assumptions. First, the daily time series data of new COVID-19 cases in each prefecture in Japan, which is the subject of analysis, can be decomposed into trend components, weekly variation components, and short-term variation components. As for the trend component, we assume that its fluctuations lead in prefectures with major airports and tourist spots, such as Tokyo and Osaka, and then spread to surrounding prefectures. Next, for the short-term variation component, although it is influenced by many complex factors, we assume that there are some regularities in fluctuations due to factors common to all prefectures, such as holidays and seasons. These assumptions require statistical analysis to provide supporting evidence.

2.2. Data

This study analyzes data on new COVID-19 cases in Japan, with the aim of elucidating the mechanisms underlying the spread of the virus across the country’s prefectures. To achieve this, we utilize daily time series data of new cases from all 47 prefectures. For reference, Table 1 lists the serial numbers and names of each prefecture, while Figure 1 presents a map of Japan labeled with prefecture numbers.

The data used in this analysis were obtained from the Ministry of Health, Labour and Welfare (MHLW) website of Japan (see https://www.mhlw.go.jp/stf/english/ index.html). The analysis covers all prefectures of Japan. The daily changes in new COVID-19 cases for each prefecture span from January 20, 2020, to May 7, 2023, covering a total of 1204 days, equivalent to 172 weeks. The original time series is denoted as $x_t$, and for the analysis, we transform it into $y_t=log(x_t+1)$. The purpose of applying the logarithmic transformation was to make the amplitude of cyclical fluctuations in the time series approximately uniform. However, since the original time series $x_t$ contains zeros, we added 1 before applying the logarithmic transformation. Therefore, the minimum value of the transformed time series is 0.

The dataset is extensive, and the time series data of new cases exhibit similar fluctuation patterns across all prefectures. Therefore, Tokyo is selected as an example, with its data presented in Figure 2. Panel (a) displays the original time series data, $x_t$, while panel (b) shows the transformed time series data, $y_t$. In the subsequent analysis, the transformed time series, $y_t$, for each prefecture will be analyzed.

2.3. Decomposition of Time Series Components

Here, the time series $y_t$ for each prefecture is decomposed into trend, weekly fluctuation, and short-term components using the extended moving linear model (EMLM) approach1. Due to the subtle differences in the decomposition results across prefectures, which are difficult to recognize at a glance, and the limitations of space, only the decomposition results for Tokyo are presented as a typical example in Figure 3.

The following observations can be made from Figure 3. First, the trend component is smooth, reflecting not only the stable fluctuation trends in the number of new COVID-19 cases but also long-term waves. The stable fluctuations are primarily attributed to the mechanisms of infection spread within the prefecture, while the long-term waves are likely influenced by external factors, such as preventive measures and policies. Furthermore, the fluctuation patterns of the trend component are highly similar across prefectures, making it possible to analyze the progression relationships of COVID-19 spread among prefectures based on the similarity of their trend components.

Next, the weekly fluctuation component exhibits a nearly consistent pattern of variation on a weekly basis in the short term, but its fluctuation pattern changes over time in the long term, warranting further analysis. Additionally, the short-term component fluctuates irregularly around zero. These fluctuations are caused by temporary shocks such as weather, holidays, seasons, regional events, and local preventive measures. Due to the complexity of these influencing factors, a detailed analysis is required to obtain meaningful results. However, since the effects of holidays and seasons are nearly uniform nationwide, it is possible to analyze their impact.

It should be noted that the weekly fluctuation component not only changes its pattern over time but also shows variations in amplitude. However, this paper does not delve into a detailed analysis of these aspects.

3. Methods and Procedures

3.1. Estimating Time Lags of Infection Spread Among Prefectures

Based on the consideration that the waves of COVID-19 infection spread among prefectures may exhibit time lags (hereafter referred to as lags), we analyze the transition relationships of infection spread among prefectures. However, the mutual influence relationships among prefectures form a complex network, making it difficult to obtain consistent results when estimating lags on a pairwise basis. Additionally, prior studies have indicated that the trend components provide more information regarding the transitions of infection spread among prefectures. Consequently, the time series of the trend components are utilized in this analysis. Therefore, in this paper, an analysis is conducted for each prefecture based on the total distance calculated from the trend component time series in relation to other prefectures.

For convenience, each prefecture is identified by the number shown in Table 1. Let the time series of trend components for prefecture i and prefecture j be denoted as ${\tilde{z}}_t^{\left(i\right)}$ and ${\tilde{z}}_t^{\left(j\right)}$, respectively. To eliminate dependency on the scale of each prefecture, the time series data for all prefectures are standardized such that their mean is zero and variance is one. The standardized time series corresponding to ${\tilde{z}}_t^{\left(i\right)}$ and ${\tilde{z}}_t^{\left(j\right)}$ are denoted as $z_t^{\left(i\right)}$ and $z_t^{\left(j\right)}$, respectively. The data incorporating lags into their indices are expressed as

$$\{z_{t+L_i}^{\left(i\right)};t=1,2,\dots ,N\},\{z_{t+L_j}^{\left(j\right)};t=1,2,\dots ,N\}(i,j=1,2,\dots ,47).$$

Here, N denotes the length of the time series, and $L_i$ and $L_j$ are the lags applied to the respective time series, both of which are integers. The Euclidean distance between prefectures i and j is defined as follows:

$$d_{ij}(L_i,L_j)=\sqrt{\sum _{t=1}^N{(z_{t+L_i}^{\left(i\right)}-z_{t+L_j}^{\left(j\right)})}^2}.$$

(1)

Since the data for $z_{t+L_i}^{\left(i\right)}$ and $z_{t+L_j}^{\left(j\right)}$ are standardized, the greater the similarity between their fluctuation patterns, the smaller $d_{ij}(L_i,L_j)$ becomes, and the higher the correlation coefficient between them. In other words, for prefectures i and j, the distance $d_{ij}(L_i,L_j)$ can be used as a measure of dissimilarity, while the correlation coefficient can be used as a measure of similarity, both reflecting the similarity in trend component fluctuations between the two. However, when used in conjunction with cluster analysis, the distance measure is more convenient.

Here, for prefecture i, the total distance with other prefectures is defined as follows:

$${\tilde{d}}_i\left(L_i\right)=\sum _{j\ne i}{d}_{ij}(L_i,L_j).$$

(2)

As is evident, the total distance ${\tilde{d}}_i\left(L_i\right)$ defined in Eq. (2) is a function of the lag $L_i$. Therefore, by minimizing ${\tilde{d}}_i\left(L_i\right)$ for $i=1,2,\dots ,47$, we can obtain an estimate of the lag $L_i$ denoted as ${\widehat{L}}_i$. Specifically, in general,

$${\tilde{d}}_i\left({\widehat{L}}_i\right)=min\{{\tilde{d}}_i\left(L_i\right);-M\le L_i\le M\}$$

where ${\widehat{L}}_i$ is obtained discretely. Here, M is a suitably chosen positive integer.

However, there is a dilemma hidden in the above steps. As can be seen from Eq. (2), the total distance ${\tilde{d}}_i\left(L_i\right)$ depends on the lags of other prefectures, meaning that unless all the lags are estimated simultaneously for every prefecture, the estimation results will not be unique. However, estimating all lags simultaneously is a difficult task.

In this paper, the estimation of the lags is performed recursively, and the procedure is repeated to address the above issue. In the first step, the initial value of the lags is set as ${\widehat{L}}_i^{\left(0\right)}=0(i=1,2,\dots ,47)$, and instead of Eq. (2), the total distance is newly defined as follows:

$${\tilde{d}}_i^{\left(1\right)}\left(L_i\right)=\sum _{j\ne i}{d}_{ij}(L_i,{\widehat{L}}_j^{\left(0\right)}).$$

Then, for $i=1,2,\dots ,47$, by minimizing ${\tilde{d}}_i^{\left(1\right)}\left(L_i\right)$, the value of the lag ${\widehat{L}}_i^{\left(0\right)}$ is updated to ${\widehat{L}}_i^{\left(1\right)}$. Furthermore,

$${\tilde{d}}_i^{\left(I\right)}\left(L_i\right)=\sum _{j\ne i}{d}_{ij}(L_i,{\widehat{L}}_j^{(I-1)})$$

and in the next step, for $i=1,2,\dots ,47$, the updated lag estimates ${\widehat{L}}_i^{\left(1\right)}$ are used in place of the initial values ${\widehat{L}}_i^{\left(0\right)}$, and the total distance ${\tilde{d}}_i^{\left(1\right)}\left(L_i\right)$ is updated to ${\tilde{d}}_i^{\left(2\right)}\left(L_i\right)$. The lag is then re-estimated by minimizing ${\tilde{d}}_i^{\left(2\right)}\left(L_i\right)$. This process is repeated for $I=2,3,\dots$, and the lag estimates will eventually converge. The final lag estimates obtained are denoted as ${\widehat{L}}_i(i=1,2,\dots ,47)$. Empirically, for Japan, the results can be obtained within five iterations.

The estimated lag results allow for the analysis of the transition relationships of infection spread between prefectures. For example, if ${\widehat{L}}_i<{\widehat{L}}_j$ for $i\ne j$, this indicates that the wave of COVID-19 infection in prefecture i is preceding that in prefecture j. However, this does not necessarily imply that the COVID-19 infection spread in prefecture i is a causal factor for that in prefecture j. Such an interpretation cannot be made conclusively without considering the similarity in the trend component fluctuations between the prefectures. Therefore, a comprehensive analysis that integrates the results of the trend component cluster analysis is needed for a more complete understanding.

3.2. Cluster Analysis of Prefectures

Using the estimated lag values ${\widehat{L}}_i(i=1,2,\dots ,47)$ obtained through the procedure described in the previous section, the inter-prefectural relationships in the spread of COVID-19 infections can be analyzed. However, the lag values can also be interpreted relatively. In such cases, the lags can be converted into relative lags ${\tilde{L}}_i$ using the following formula:

$${\tilde{L}}_i={\widehat{L}}_i-L_{\mathtt{max}}(i=1,2,\dots ,47),$$

(3)

where $L_{\mathtt{max}}=max\{{\widehat{L}}_i;i=1,2,\dots ,47\}$ represents the maximum value of ${\widehat{L}}_i(i=1,2,\dots ,47)$. Clearly, the relative lags have a maximum value of zero.

Here, the time series data of the trend components are lag-adjusted as follows:

$$\{z_{t+{\tilde{L}}_i}^{\left(i\right)};t=1,2,\dots ,N\}(i=1,2,\dots ,47).$$

(4)

For the first few values of the time series $z_{t+{\tilde{L}}_i}^{\left(i\right)}$, if their indices satisfy $t+{\tilde{L}}_i<1$, these values become missing and are supplemented with zeros. Since the initial values in the time series of the trend components for each prefecture are mostly zero or close to zero, the advantage of adjusting lags using relative lags in this way is that it preserves the length of the time series at N without compromising its accuracy.

By applying the lag-adjusted time series data described in Eq. (4) to Eq. (1), the distance $d_{ij}({\tilde{L}}_i,{\tilde{L}}_j)$ for each pair of prefectures, where $i=2,3,\dots ,47$ and $j=1,\dots ,i-1$, is calculated, thereby defining an inter-prefectural distance matrix. This resulting distance matrix corresponds to the minimum total distance between prefectures, reflecting their similarity based on trend components. The matrix is then utilized for cluster analysis among prefectures2.

The following describes the procedure for grouping the prefectures nationwide based on the results of the cluster analysis. The number of groups to be formed is equal to the number of cluster branches that are cut horizontally, and by determining this, the group membership of each prefecture is decided. For example, if the desired number of groups is $K>1$, the set of prefecture numbers belonging to each group, denoted as $G_k$ for $k=1,2,\dots ,K$, is defined, and the number of prefectures in each group is denoted as $n\left(G_k\right)$. The following time series are then defined for the average and total variance of the trend components for all prefectures:

$${\overline{z}}_t={\displaystyle \frac{1}{47}}\sum _{i=1}^{47}{z}_{t+{\tilde{L}}_i}^{\left(i\right)},V_t={\displaystyle \frac{1}{46}}\sum _{i=1}^{47}{(z_{t+{\tilde{L}}_i}^{\left(i\right)}-{\overline{z}}_t)}^2(t=1,2,\dots ,N).$$

Next, the average time series within group $G_k$ is defined as:

$${\overline{z}}_t^{\left(k\right)}={\displaystyle \frac{1}{n\left(G_k\right)}}\sum _{i\in G_k}{z}_{t+{\tilde{L}}_i}^{\left(i\right)}(k=1,2,\dots ,K;t=1,2,\dots ,N).$$

Then, the time series of the variance between groups is defined as:

$$W_t={\displaystyle \frac{1}{K-1}}\sum _{k=1}^K{({\overline{z}}_t^{\left(k\right)}-{\overline{z}}_t)}^2(t=1,2,\dots ,N).$$

The ratio of the total variance to the variance between groups is calculated as:

$$r={\displaystyle \frac{{\sum }_{t=1}^N{W}_t}{{\sum }_{t=1}^N{V}_t}}$$

Finally, by maximizing the ratio of the variance between groups, r, the number of groups, K, is determined, and the grouping of the prefectures based on the results of the cluster analysis is naturally decided.

Using the results of the cluster analysis and prefectural grouping, the similarity in the COVID-19 infection spread patterns between prefectures is analyzed, and grouping of the prefectures is performed accordingly. By associating this with the analysis of the transition relationships in the spread of COVID-19 between prefectures, the mechanisms of the spread of COVID-19 across prefectures are brought to light.

3.3. Further Analysis Based on Trend Components

The framework for analyzing the trend components described above is primarily based on the similarity of fluctuation patterns in the trend components of each prefecture. In addition to this, the trend components can also be used for the following analyses.

3.3.1. Mutual Influences Between Prefectures

The time series data of lag-adjusted trend components are defined as

$${\tilde{z}}_{t+{\tilde{L}}_i}^{\left(i\right)};t=1,2,\dots ,N(i=1,2,\dots ,47).$$

Here, it should be noted that ${\tilde{z}}_t$ represents the original trend time series, that is, the trend time series before standardization. The first-order difference series for each of these is defined as $q_{t+{\tilde{L}}_i}^{\left(i\right)}={\tilde{z}}_{t+{\tilde{L}}_i}^{\left(i\right)}-{\tilde{z}}_{t-1+{\tilde{L}}_i}^{\left(i\right)}(i=1,2,\dots ,47)$. These difference series are considered to be approximately mean-stationary time series.

Let $C_{ij}$ denote the covariance between $q_{t+{\tilde{L}}_i}^{\left(i\right)}$ and $q_{t+{\tilde{L}}_j}^{\left(j\right)}$. If the value of $C_{ij}$ is large, it can be inferred that there is a strong interconnection in the spread of COVID-19 infections between prefectures i and j. This indicates strong mutual influence in the spread of infections, providing crucial information for analyzing the propagation relationships of infection spread among prefectures.

3.3.2. Factors Influencing the Lag Values

The lag values can be interpreted as information regarding the propagation relationships of COVID-19 between prefectures. Identifying the factors influencing these lags is crucial for understanding inter-prefecture propagation dynamics. Factors that may affect the lag include population, land area, prefectural boundary length, and the number of airline passengers. Although railway usage could also influence the lag, consistent data were unavailable and, therefore, excluded from the analysis.

To investigate this, we employ a regression model where the estimated lag value for each prefecture serves as the objective variable. The explanatory variables include population, land area, perimeter length, the number of domestic airline passengers, and the number of international airline passengers (or their transformed values, where applicable) for each prefecture.

Through the estimation results for the regression coefficients, we can analyze the factors influencing the lag.

3.4. Analysis of Factors Influencing Short-Term Components

As mentioned earlier, short-term components are caused by temporary shocks such as weather, holidays, seasons, local events and festivals, as well as regional preventive measures. Due to the complexity of these influencing factors, detailed analysis is required. However, since the effects of holidays and seasons are almost uniform nationwide, it is possible to analyze their impacts.

In this study, we analyze the impact of holidays and seasonal factors on short-term components using regression analysis. Denoting the time series of short-term components for prefecture i as $u_t^{\left(i\right)}$, we compile the time series data for short-term components across all prefectures to obtain the following 47-dimensional time series data of length N.

$$\{u_t^{\left(i\right)};t=1,2,\dots ,N\}(i=1,2,\dots ,47).$$

Furthermore, for this 47-dimensional data, we normalize the data for each component individually and apply principal component analysis to obtain the first principal component score data $\{v_t;t=1,2,\dots ,N\}$. Assuming that the effects of holidays and seasons are common to all prefectures, it is plausible that the time series $v_t$ contains variations attributable to these effects.

Holidays may differ significantly from regular days in terms of people’s daily routines and social activity patterns, potentially exerting a widespread influence that increases the likelihood of COVID-19 transmission. However, given that there is likely a lag between transmission and the confirmation of positive COVID-19 cases, it is necessary to estimate this lag using the available data. First, for holidays that do not overlap with Saturdays or Sundays, we define $J+1$ holiday dummy variables as follows, covering the holiday itself and up to J days after the holiday.

$$\begin{array}{ccc}& & g_t^{\left(1\right)}=\left\{\begin{array}{cc}1,\hfill & \mathrm{if\; day} t \mathrm{is\; a\; holiday}\hfill \\ 0,\hfill & \mathrm{if\; day} t \mathrm{is\; not\; a\; holiday}\hfill \end{array}\right.,\hfill \\ & & g_t^{(j+1)}=\left\{\begin{array}{cc}1,\hfill & \mathrm{if} j \mathrm{days\; before\; day} t \mathrm{was\; a\; holiday}\hfill \\ 0,\hfill & \mathrm{if} j \mathrm{days\; before\; day} t \mathrm{was\; not\; a\; holiday}\hfill \end{array}\right.(j=1,2,\dots ,J).\hfill \end{array}$$

The seasonal effects are represented using 11 monthly seasonal dummy variables as follows:

$$h_t^{\left(m\right)}=\left\{\begin{array}{cc}1,\hfill & \mathrm{if\; day} t \mathrm{is\; in\; month} m\hfill \\ 0,\hfill & \mathrm{if\; day} t \mathrm{is\; not\; in\; month} m\hfill \end{array}\right.(m=1,2,\dots ,11).$$

Based on these, we establish the following regression model:

$$v_t=a+\sum _{j=0}^J{\alpha }_j{g}_t^{\left(j\right)}+\sum _{m=1}^{11}{\beta }_m{h}_t^{\left(m\right)}+{\epsilon }_t(t=1,2,\dots ,N).$$

(5)

In this model, ${\alpha }_j(j=0,1,\dots ,J)$ represents the holiday effect for each lag of holidays, while ${\beta }_m(m=1,2,\dots ,11)$ represents the seasonal effect, and ${\epsilon }_t$ denotes the error term. Note that this setup assumes the seasonal effect for December to be zero, and thus the values of ${\beta }_m(m=1,2,\dots ,11)$ should be interpreted relative to December.

Note that not all explanatory variables in the model given by Eq. (5) need to be included. In other words, variable selection must be performed before using this model. In this study, the minimum AIC method is primarily employed3.

4. Results

4.1. Estimation of Lags

In this study, the interval limit for lag estimation was set to $M=30$, and the minimum total distance method proposed in Section 3.1 was used to estimate the trend component variation lags among prefectures. Furthermore, to normalize the maximum lag to zero, these lags were converted into relative lags ${\tilde{L}}_i$ based on Eq. (3). The results are shown in Table 2.

The following observations can be drawn from Table 2. Firstly, Okinawa (No. 47) and Tokyo (No. 13) have estimated relative lags of $-10$ and $-8$, respectively, indicating small values. These regions have international airports with numerous flights connecting Japan to various countries worldwide, especially in Asia, which results in a high volume of foreign tourists. Consequently, it is believed that COVID-19 initially emerged in these areas and subsequently spread across Japan.

Next, Hokkaido (No. 1) and Gunma (No. 10) have relative lags estimated at $-7$, while the surrounding prefectures of Tokyo, such as Miyagi (No. 4), Saitama (No. 11), Chiba (No. 12), and Kanagawa (No. 14), as well as the central regions of the Kansai area, including Kyoto (No. 26), Osaka (No. 27), and Nara (No. 29), all have relative lags estimated at $-6$. Most of these prefectures are located near international airports, have active business environments, or are home to famous tourist destinations, which leads to high levels of movement and interaction.

In contrast, Akita (No. 5), Gifu (No. 21), Mie (No. 24), and Tokushima (No. 36) have the maximum relative lag value of 0, while Toyama (No. 16), Okayama (No. 33), Ehime (No. 38), Kochi (No. 39), and Oita (No. 44) have relative lags estimated at $-1$. These prefectures likely have a relatively large rural population, which results in less frequent interaction with external areas. Surprisingly, Iwate (No. 3), which is well known for having the latest first positive case of COVID-19 in Japan, has a relative lag of $-4$, indicating that its relative lag is not particularly high.

4.2. Cluster Analysis

Using the estimated relative lags mentioned above, a cluster analysis of the prefectures was performed using the method introduced in Section 3.2. Note that the Ward method in the free software R was used for the cluster analysis. The resulting dendrogram is shown in Figure 4. The dendrogram resembles an inverted tree. From this figure, hierarchical similarities in the fluctuation patterns of the trend components between prefectures can be observed. The numbers at the bottom of the leaves represent the prefecture numbers, and as you move upward, the branches bind hierarchically. Prefectures (or groups of prefectures, called clusters) that bind earlier have a higher degree of similarity. Additionally, the height of the branches represents the distance between prefectures (or clusters).

Based on the results in Figure 4, the grouping of prefectures was performed. By testing the number of groups from 2 to 17, it was found that the variance ratio between groups was maximized at $K=13$, with a maximum value of 1.56. Therefore, all prefectures can be divided into 13 groups.

One notable feature is the presence of two single-prefecture groups: Hokkaido (No. 1) and Okinawa (No. 47). These prefectures are estimated to have small relative lags. Both Hokkaido and Okinawa have international airports, attract many foreign tourists, and are located on remote islands. Therefore, it is understood that the spread of COVID-19 in these prefectures followed a unique pattern. In Hokkaido, the spread of COVID-19 progressed relatively early, but it was particularly rapid in Okinawa, where the spread was the fastest in Japan. However, geographically, these prefectures did not have a strong influence on other regions.

Next, there are three groups consisting of two prefectures. The first group is made up of Miyagi (No. 4) and Gunma (No. 10). In these two prefectures, the spread of COVID-19 occurred relatively early. The second group includes the two neighboring prefectures of Toyama (No. 16) and Ishikawa (No. 17), both of which experienced a delay in the spread of COVID-19. The third group is composed of Gifu (No. 21) and Oita (No. 44). While these two prefectures are similar in that the spread of COVID-19 was delayed, it is interesting that they are geographically distant from each other.

Additionally, there is one group containing three prefectures: Fukushima (No. 7), Ibaraki (No. 8), and Mie (No. 24). These prefectures are characterized by a relatively delayed spread of COVID-19.

For the remaining groups, they will be explained in relation to the progression of COVID-19 spread among prefectures, as observed in the previous section, following the order from left to right in the dendrogram shown in Figure 3.

First, the group on the far left contains four prefectures from the Kanto region, structured as follows:

$$Tokyo(No.13)\Rightarrow \{Kanagawa(No.14)\iff \{Saitama(No.11)\iff Chiba(No.12)\left\}\right\}$$

Here, ⇒ represents a unidirectional progression relationship, ⇔ indicates a simultaneous progression relationship, and {*} denotes the same cluster. This result suggests that the spread of COVID-19 progressed most rapidly in Tokyo, reflecting the active human interactions with its surrounding prefectures.

Next, the group consisting of four prefectures in the Kansai region, along with Aichi and Fukuoka, has the following structure:

$$\begin{array}{ccc}& & Fukuoka(No.40)\Leftarrow \left\{\right\{Kyoto(No.26)\iff Nara(No.29)\}\hfill \\ & \leftrightarrow & \{Osaka(No.27)\Rightarrow \{Aichi(No.23)\iff Hyogo(No.28)\left\}\right\}\}\hfill \end{array}$$

Here, ⇐ indicates a progression relationship from the latter to the former, and ↔ signifies a complex bidirectional progression relationship. This represents the mutual relationships in the spread of COVID-19 among the prefectures in the Kansai region, Aichi, and Fukuoka. Additionally, Figure 3 shows that these two groups together with Okinawa form a larger cluster.

Additionally, the next group includes only four prefectures, but it spans across the Tohoku, Chubu, and Chugoku regions. Its structure is as follows:

$$\begin{array}{ccc}& & Yamagata(No.6)\leftrightarrow \{Okayama(No.33)\hfill \\ & \Leftarrow & \{Nagano(No.20)\iff Yamaguchi(No.45)\}\}\hfill \end{array}$$

Here, ↔ represents a complex bidirectional progression relationship, ⇐ indicates a progression relationship from the latter to the former, and ⇔ denotes a simultaneous progression relationship. This result highlights the spread of COVID-19 among these prefectures, revealing unique interregional interactions despite their geographical distances.

Other groups have complex structures, and detailed analyses of these groups are omitted here.

Synthesizing the results above, the following insights can be inferred: The initial source of COVID-19 infections likely entered Japan from abroad through Tokyo. After the spread of infections in Tokyo, they dispersed to surrounding prefectures. This impact extended to Osaka, Kyoto, and Nara. Subsequently, the spread reached Aichi Prefecture, centered on Nagoya, and Fukuoka Prefecture, centered on Fukuoka City. This progression formed the primary pattern of COVID-19’s spread in Japan, closely related to the fact that these metropolitan areas are well-connected by Shinkansen and flights, ensuring easy access via both land and air routes.

4.3. Further Analysis Based on Trend Components

4.3.1. Mutual Influences Between Prefectures

As shown in Section 3.3.1, we calculated the covariance of differenced time series for the trend components between prefectures to analyze the strength of interconnections in the spread of COVID-19 infections. For clarity, we multiplied the covariance values by 1000 and included only those exceeding 2.30. The results are presented in Table 3. The first column and the first row of Table 3 indicate the prefecture numbers i and j, while the values in the table represent the covariance $C_{ij}$ between the prefectures i and j.

The following observations can be made from Table 3. The table includes 15 prefectures, among which are five prefectures from the Kanto region (Nos. 10–14) and four from the Kansai region (Nos. 26–29). These two regions are well-known as the political and economic centers of Japan, and they are also famous for their commercial areas and tourist destinations, making them hubs of human interaction. Additionally, Toyama (No. 16), Gifu (No. 21), and Aichi (No. 23) are included. These prefectures are adjacent to the Kansai region and are known for famous tourist spots. Furthermore, Fukuoka (No. 40) and Oita (No. 44) from the Kyushu region are also listed. These prefectures are hubs in Kyushu, well-connected by domestic transportation such as bullet trains and airplanes, and are close to China and South Korea, which makes them closely tied with international connections. They are also famous tourist destinations, contributing to vibrant exchanges with foreign regions.

In terms of the number of prefectures with strong connections, Fukuoka (No. 40) ranks first with 14 prefectures, followed by Saitama (No. 11) in second place with 10. This suggests a wide range of human interactions involving these two prefectures. Additionally, regarding the covariance values, the covariance between Fukuoka (No. 40), Osaka (No. 27), and Aichi (No. 23) exceeds 2.98, marking the highest values observed. These results provide further insight into the clustering analysis, which revealed that Fukuoka (No. 40) is grouped with the Kansai region cluster.

4.3.2. Factors Influencing the Lag Values

As mentioned in Section 3.3.2, to investigate the factors influencing the lag values, we consider the following regression model:

$${\tilde{L}}_i=b+{\gamma }_{1}log\left(s_1^{\left(i\right)}\right)+\sum _{j=2}^5{\gamma }_j{s}_j^{\left(i\right)}+e_i,$$

(6)

where i denotes the prefecture number, ${\tilde{L}}_i$ represents the estimated lag value, and $s_1^{\left(i\right)}$ represents the population (in millions of people). To enhance the model’s performance, logarithmic values were used as explanatory variables based on trial and error. Additionally, $s_2^{\left(i\right)}$ to $s_5^{\left(i\right)}$ represent the land area (in 10,000 square kilometers), the perimeter length length (in 10,000 kilometers), the number of domestic airline passengers (in millions of people), and the number of international airline passengers (in millions of people), respectively4. Here, b is a constant, ${\gamma }_1$ to ${\gamma }_5$ are the corresponding coefficients, and $e_i$ is the error term.

For the selection of explanatory variables, we used a stepwise regression procedure. As a result, the AIC value of the model in Eq. (6), which included all explanatory variables, was 68.87. The stepwise regression yielded a model that included the logarithm of the population, $log\left(s_1^{\left(i\right)}\right)$, and the number of domestic airline passengers, $s_5^{\left(i\right)}$, achieving a minimum AIC value of 63.26. Table 4 presents the estimated regression coefficients and the corresponding p-values. From Table 4, the following observations can be made: there is a significant negative effect from both the population and the number of domestic airline passengers. It implies that in prefectures with larger populations, the lag tends to be more negative, indicating that the spillover effect to other prefectures is greater. A similar observation is made for the number of domestic airline passengers.

4.4. Analysis of Holiday and Seasonal Effects

In this section, we analyze the effects of holidays and seasons included in the short-term component using regression analysis. First, principal component analysis (PCA) was performed based on the correlation matrix of the short-term components for all combinations of prefectures. The contribution ratio of the first principal component (PC) was 42.78%, and its time series, $v_t$, is shown in Figure 5.

Next, we set the number of holiday dummy variables to $J+1=20$ and incorporated 20 holiday dummy variables and 11 seasonal dummy variables to construct a regression model as shown in Eq. (5). The AIC value for the model using all explanatory variables was -596.13. Based on this, stepwise regression was performed for variable selection, resulting in a regression model that included 2 holiday dummy variables and 8 seasonal dummy variables, with a minimum AIC value of -628.18. Furthermore, using a 5% significance level as the criterion, t-tests confirmed that all regression coefficients for these dummy variables were significant. Table 5 presents the estimated regression coefficients and the upper bounds of their corresponding p-values. From Table 5, the following observations can be made: there is a significant negative effect on the day after a holiday and a positive effect eight days later. This suggests that despite active interpersonal interactions during holidays, the number of positive cases increases about a week later.

Regarding seasonal effects, assuming the effect for December is zero, the largest positive effect is observed in November, with significant positive effects also appearing in April, May, and August. The effect in October, while relatively weaker, remains positive. On the other hand, the effects for March and June are negative, with the effect in June being weaker than that in March. The adjusted coefficient of determination for the regression model is 0.4119, which is not particularly high. However, it is meaningful to have identified such significant holiday and seasonal effects amidst complex disturbance factors.

5. Discussion

In this study, we examined the mechanism of COVID-19 spread across prefectures in Japan by analyzing daily time series data of new COVID-19 positive cases. Through a comprehensive approach, we first performed time series decomposition, which allowed us to break down the complex fluctuation patterns in the data into trend components, weekly fluctuation components, and short-term components. We then estimated the relative lags in the progression of infection spread among prefectures using the shortest distance method, based on the decomposition results of the trend components. Finally, we classified the prefectures through cluster analysis, providing deeper insights into the regional dynamics of infection spread.

The key findings of the study are as follows. First, the spread of COVID-19 in Okinawa and Tokyo exhibited earlier fluctuations compared to other prefectures, suggesting that these regions were among the first to experience significant increases in new infections. Kyoto and Osaka followed Tokyo with a slight delay, but they also demonstrated earlier infection spread relative to other prefectures. The cluster analysis further classified the prefectures into 13 groups, highlighting regional differences in the progression of the pandemic. Notably, Hokkaido and Okinawa demonstrated unique fluctuation patterns, each forming its own distinct group. Okinawa, despite experiencing the earliest progression of the pandemic, had its influence geographically limited, as expected due to its isolation. The cluster analysis results also revealed clear progression structures in the groups centered around the Kanto region, led by Tokyo, and the Kansai region, led by Osaka and Kyoto. These findings suggest that the spread of COVID-19 in Japan likely originated in Tokyo, subsequently influenced Osaka and Kyoto, and then spread nationwide, with clear regional patterns emerging over time.

In further examining the trend components and the interconnections between prefectures, we analyzed the strength of regional ties by calculating the covariance between differenced time series of the trend components for different prefectures. The results revealed significant mutual influences, particularly between the Kanto and Kansai regions, which are central to Japan’s economy and human interactions. For example, Fukuoka stood out with the highest number of strong connections to other prefectures, reinforcing its pivotal role in the spread of infections. These findings, when combined with the cluster analysis results, highlight the crucial role of regional dynamics in the pandemic’s progression, with Tokyo acting as the primary epicenter influencing surrounding areas like Osaka, Kyoto, and Fukuoka. The clustering analysis further corroborates the idea that the spread followed distinct regional patterns, emphasizing the significance of inter-prefecture interactions in understanding the geographical spread of COVID-19 across Japan.

A significant contribution of this study was the regression analysis performed to understand the effects of holidays and seasonal factors on short-term fluctuations in infection rates. The analysis revealed a significant negative effect the day after a holiday and a positive effect eight days later, implying that the social interactions during holidays played a key role in the transmission of the virus. This aligns with the incubation period of COVID-19, where the effects of increased interactions become more evident approximately a week later. The analysis of seasonal effects further revealed that the largest positive effect occurred in November, followed by positive effects in April, May, and August. Conversely, March exhibited a significant negative effect. These findings underscore the importance of accounting for both holidays and seasonal patterns in infection control strategies, as they contribute to spikes in cases following periods of heightened social activity.

Further analysis of factors influencing the lag in infection spread, as described in Section 4.3.2, provided additional insights. Regression analysis showed that population size and the number of domestic airline passengers had a significant negative effect on the lag value, indicating a faster spread in more populated prefectures and those with higher domestic travel traffic. This suggests that social interaction, especially in regions with dense populations and high connectivity, is a key factor driving the spread. The AIC value from the stepwise regression suggested that incorporating these variables into a predictive model could refine our understanding of lag patterns in infection dynamics. These findings provide a foundation for developing more targeted interventions, especially in regions with high population density or significant transportation hubs.

Despite the valuable insights gained from analyzing the holiday and seasonal effects, this study has several limitations. One key limitation is the focus on regional trends and the lack of in-depth analysis of infection dynamics within individual prefectures. Future research could extend the current study by examining intra-prefecture transmission dynamics and identifying local factors that may influence the spread of infections. Additionally, while the regression analysis identified significant holiday and seasonal effects, the model’s adjusted coefficient of determination ($R^2$) was relatively low (0.4119), indicating that the model only explained a portion of the variability in the data. Future work could aim to refine the model by incorporating other potential explanatory variables, such as demographic factors, healthcare infrastructure, and public policy measures, to improve the model’s explanatory power.

Another limitation of this study is the exclusion of other possible factors influencing the spread of COVID-19, such as mobility patterns, public behavior, and the role of interventions such as lockdowns and vaccination campaigns. While we focused on the short-term components related to holidays and seasonal variations, a more comprehensive model would include these additional factors to provide a more holistic understanding of the spread of the virus. Moreover, the study primarily relied on statistical techniques to identify trends and relationships in the data. While these techniques are powerful, they do not account for all possible causal relationships and may overlook complex interactions between different variables. Future research could benefit from incorporating more advanced machine learning or causal inference methods to uncover deeper insights into the mechanisms of infection spread.

In conclusion, this study provides valuable insights into the regional dynamics of COVID-19 spread in Japan, highlighting the importance of holidays, seasonal variations, and geographic factors in shaping the pandemic’s progression. The findings from the cluster analysis and regression models offer practical implications for infection control strategies, emphasizing the need for targeted interventions during holidays and seasonal peaks. Furthermore, the covariance analysis between trend components, along with the investigation of lag values, reveals important regional interdependencies that can help guide future policy decisions. However, the study also identifies several areas for future research, particularly the need for a more detailed analysis of local transmission dynamics and the inclusion of additional factors in predictive models. By addressing these challenges, future studies can contribute to more effective public health strategies in managing infectious diseases.