COVID-19 transmission risk in Surabaya and Sidoarjo: an inhomogeneous marked Poisson point process approach

Marked Poisson point processes

We denote a point pattern data as a set of coordinate locations $\mathbf{x}=\{x_1,\cdots ,x_n\}$ representing all n recorded home addresses of COVID-19 new confirmed cases observed within a spatial domain $D\subset {\mathbb{R}}^2$ with an area |D|. In this study, D represents the Surabaya and Sidoarjo regions. Each coordinate $x_i,i=1,\cdots ,n$ is accompanied by a mark $m\in \mathcal{M}=\{1,\dots ,M\}$ denoting the time period when the new cases occurred, where we consider $M=6$ (see Table 1). The point pattern can also be written as $\mathbf{x}=\{{\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(M)}\}$, where ${\mathbf{x}}^{(m)}=\{x_1,\cdots ,x_{n_m}\}$ is the pattern of home coordinates belonging to period m, and $n_m$ is the total new COVID-19 cases at that period $m=1,\dots ,M$. The corresponding stochastic marked point process is then denoted by $\mathbf{X}=\{{\mathbf{X}}^{(1)},\dots ,{\mathbf{X}}^{(M)}\}$. More details on marked point processes can be found, for example, in Baddeley et al. (2015); Illian et al. (2008), although literature in not prominent when marks are taken into the modeling framework.

If $\mathbf{X}$ is a marked Poisson point process, the model is fully determined by its intensity function. The intensity represents the probability of having a new COVID-19 case at a given coordinate u and marked time m. To study the COVID-19 transmission in Surabaya and Sidoarjo, we focus in assessing the intensity in relation with the spatial covariates detailed in Sect. 2. In this study, the relation is captured by a log-linear intensity form (e.g.Waagepetersen 2007; Choiruddin et al. 2018; Husain and Choiruddin 2021)

where $m=1,2,\cdots ,6$ is the mark representing the different times of COVID-19 phase, $z_j(u)$ is the jth-spatial covariate measured at u, $\mathit{\beta }={({\beta }_0,{\beta }_1,\cdots ,{\beta }_4)}^{\top }$, ${\beta }_{0,m}$ stands for the intercept parameter, and ${\beta }_{j,m}$ represents the corresponding regression parameter. Note that both ${\beta }_0$ and ${\beta }_j$ are indexed by mark m to allow the effect of spatial covariates to differ from time to time. This is essential for the goals of the paper, as authorities need to evaluate the information of covariates depending on the period in which restrictions are taking place. To fit $\mathtt{ppm}$ function of the $\mathtt{spatstat}$ package, we rearrange the intensity model (1) by

where $i=2,\cdots ,6$ and $\mathtt{I}(i=m)$ is an indicator function. Equation (2) is equivalent to that of (1) by equating ${\beta }_{0,m}={\alpha }_{0,1}+{\alpha }_{0,i}\mathtt{I}(i=m)$ and ${\beta }_{j,m}={\alpha }_{j,1}+{\alpha }_{j,i}\mathtt{I}(i=m)$. Note that parameter estimation of $\mathit{\beta }$ and model interpretation detailed in Sects. 3.2–4 are conducted through $\mathit{\alpha }$.

Parameter estimation and model assessment

To estimate the parameter $\mathit{\beta }$, we maximize the corresponding log-likelihood function (Baddeley et al. 2015)

where in particular the Berman-Turner scheme is employed to discretize the integral term in (3), as

where $u_i,i=1,\dots ,n_m+d_m$ are points consisting of the $n_m$ data points and $d_m$ dummy points, and where the $w(u_i)>0$ are quadrature weights such that ${\sum }_{i}w(u_i)=|D|$. Using this technique, (3) is then approximated by

where $w_i=w(u_i),{\lambda }_{i,m}(\mathit{\beta })=\lambda (u_i,m;\mathit{\beta }),$ and $y_i=w_i^{-1}\mathtt{I}(u_i\in {\mathbf{x}}^{(m)},m\in \mathcal{M})$. The approximation in (4) is equivalent to the weighted likelihood function of Poisson generalized linear models (GLMs), and it is efficiently implemented in e.g. $\mathtt{ppm}$ function of the $\mathtt{spatstat}$ package (Baddeley et al. 2015). We denote the resulting estimator obtained from (4) by $\widehat{\mathit{\beta }}$ (and denote the estimator of $\mathit{\alpha }$ by $\widehat{\mathit{\alpha }}$). To handle the multicollinearity problems (See Table 2 and Sect. 2), we perform model selection based on the Bayesian Information Criterion (Choiruddin et al. 2021).

To perform model validation, we consider simulation envelopes. Envelopes define the critical bounds of the statistical test of a summary function, such as K-function or cross-K-function, which validates suitability point pattern data to point process model (Baddeley et al. 2015; Choiruddin et al. 2021).

Model assessment

Due to high bivariate correlations among covariates, we perform model selection (Choiruddin et al. 2021) to find the best intensity model (1) or (2) and to tackle the multicollinearity problem. For $i=2,\cdots ,6$, the resulting model is

where among the four covariates, only enterprise ($z_1(u)$), place of worship ($z_3(u)$), and population ($z_4(u)$) densities are retained in the model. Note that the regression coefficient ${\widehat{\beta }}_{k,m},k=0,1,3,4,m=1,\cdots ,6$ is represented by ${\widehat{\alpha }}_{k,m}$ through ${\widehat{\beta }}_{k,m}={\widehat{\alpha }}_{k,1}+{\widehat{\alpha }}_{k,i}\mathtt{I}(i=m),i=2,\cdots ,6$. The resulting estimates $\widehat{\mathit{\alpha }}$ are depicted in Table 3. Following the intensity model (5), we are able to interpret the effect of covariates for each period by assessing the magnitude of ${\widehat{\beta }}_{j,m}$ as well as to evaluate the evolution of the covariates effect for each period compared to the one at pre-PSBB by ${\widehat{\alpha }}_{j,m}$. Table 3 additionally shows the standard error associated to $\widehat{\mathit{\alpha }}$ as an indication of the uncertainty magnitude of such estimates. See Sects. 4.2–4.3 for further details.

We then evaluate the performance of our inhomogeneous marked Poisson point process model with intensity (5) using envelopes of inhomogeneous K-functions (Fig. 3). In general, the point process designed in this study fits the COVID-19 pattern well since for each time the empirical K-function (black solid line) lies within the envelope interval (grey area). Furthermore, the interval is tight (especially in Fig. 3c–f where more cases are observed) and the average of the simulated K-functions (red dashed line) is pretty close to the empirical K-function. We observe a wider envelope in Fig. 3a, b due to a larger variability prompted by a lower number of cases. In addition, the empirical K-functions of Fig. 3a, b tend to be slightly lower than what mark Poisson point process expects after a radius of 8 km. One reason may come from the lack of information as the number of cases in pre-PSBB and PSBB 1 is much lower (See also Sect. 4.5). We note that the K-function is a second-order characteristic describing the model structure in terms of the possible existence of clustering due to latent interactions. The small variability of the envelopes based on simulations of the null model is an indication of the low uncertainty of our fitted models.

Effect of sources of crowd

Table 3 presents the parameter estimates $\widehat{\mathit{\alpha }}$ of the mark Poisson point process with intensity (5). Although we denote the model in terms of $\mathit{\alpha }$, the interpretation can also be conducted through $\mathit{\beta }$. We prefer to consider $\mathit{\alpha }$ in order to better understand how the effect of each ones of the covariates evolves with time, and evaluate the large-scale social distancing (PSBB) regulated by the government. For example, the ${\widehat{\alpha }}_{3,i}$ explains the effect of increasing/decreasing (1 unit per km${}^2$) of the worship place density ($z_3$) to the relative risk of COVID-19 at period $i,i=2,\cdots ,6$ compared to the pre-PSBB (period 1). We additionally provide in Table 4 the effect of each of the spatial covariates (in terms of $exp(\widehat{\beta })$) to understand the different effects over each period of time.

The intercept provides an insight of the predicted probability of observing a new COVID-19 case in an area with no enterprise, place of worship, and population within area of 1 km${}^2$; such a typical area might represent the southeastern part of the Surabaya-Sidoarjo region (conservation area, see Sect. 2). The probability of observing a new case is very small at the pre-PSBB (0.067). The relative risk does not significantly change during the PSBB 1, but then largely increases during the PSBB 2 to the New Normal, by 2.8 to 11 times compared to the one in the pre-PSBB. By observing the number of cases and its spatial distribution over time (Table 1 and Fig. 5a), it is clearly seen that the disease risk inflates over time, except for PSBB 1.

According to (5), two sources of crowd are involved in the model: enterprise and place of worship. Enterprise density shows a positive effect. For example, at the pre-PSBB, $exp({\widehat{\beta }}_{1,1})=exp({\widehat{\alpha }}_{1,1})=1.353$ implies that for every increase in one unit of enterprise within an area of 1 km${}^2$ during the pre-PSBB the risk of COVID-19 cases increases up to 36%. The relative risk remains increasing by 14% for every increase of an enterprise at the PSBB 1 relative to the one at pre-PSBB. To measure the enterprise density effect in PSBB 1 itself, one could quantify $exp({\widehat{\beta }}_{1,2})=exp({\widehat{\alpha }}_{1,1}+{\widehat{\alpha }}_{1,2})=exp(0.302+0.133)=1.545$ (also presented in Table 4), meaning that for each addition of one enterprise within a radius of 1 km during the PSBB 1, the risk increases by 54%. This also means that every reduction of one enterprise within 1 km${}^2$ (interpreted as an enterprise closing during the PSBB 1) will reduce the transmission risk by 54%. Next, during the PSBB 2 to New Normal period, the effect of enterprise density is not significantly different from the effect of that in the pre-PSBB, i.e., $exp({\widehat{\beta }}_{1,1})\approx exp({\widehat{\beta }}_{1,3})\approx \cdots \approx exp({\widehat{\beta }}_{1,6})$, see Table 4. In this context, enterprise density might be viewed as a measure of the urbanization of an area. The denser the enterprise in an area, the more human interaction would take place and therefore higher risk of virus spread. This perspective was also previously used by Kwok et al. (2021), associating building geometry (height) with urbanization level and higher risk of COVID-19. To our knowledge, no previous study has investigated the direct effect of enterprise density on the risk of COVID-19. The positive effect of the enterprise suggests that human interaction in the workplace plays a significant role in spreading COVID-19 in Surabaya-Sidoarjo region.

On the other hand, place of worship shows a stronger effect with a decreased risk of about 84% (at the pre-PSBB) when reducing/closing 1 unit of place of worship in an area of 1 km${}^2$. The effect of place of worship during the PSBB 1 and PSBB 2 remains as high as the one at the pre-PSBB, and even grows up to 41% at the PSBB 3. A similar result was reported by Scarpone et al. (2020), arguing that it is not the place of worship per se that is the source of transmission, but it rather serves as an environment where more social interpersonal interactions may take place. Nevertheless, higher effects of place of worship relative to enterprise in our data is somewhat expected. In a region with a Muslim majority population, the public activity in Surabaya and Sidoarjo was hugely affected by Islamic practices. Ramadan, a month based on the Lunar calendar, is a special month in Islamic faith and therefore is usually profited by Muslims to go to the mosques more frequently. It started on 23 April and ended on 23 May 2020, a day in which mosques usually hold a public prayer attended by a relatively larger number of people than the daily routine prayer. In addition, this religious holiday is generally enriched until a week after Ramadan with cultural habits of gathering with family and friends. This period coincides with PSBB 1-PSBB 3. This shows that the action taken by government to close place of worship during these periods was effective to slow down the transmission risk.

It is important to note that for both enterprise and place of worship, positive effects are observed during PSBB (at any stage) when compared to Pre-PSBB, so that closing one unit of enterprise or place of worship at any PSBB stage provides higher effects on the transmission risk reduction in PSBB than in Pre-PSBB. This result reflects the effectiveness of social distancing in reducing transmission risk of COVID-19 related to sources of crowd since enterprises and places of worship were strictly closed during the PSBB 1, and then such a regulation was relaxed gradually in the PSBB 2 and 3. Table 4 also highlights this finding.

Effect of population density

Although the population density is significant, the effect is low (increasing 5% of risk at the pre-PSBB for each addition of 1000 population in a surrounding 1 km${}^2$ area). Interestingly, the effect of population density to relative risk keeps decreasing over time, which might also reflect the effectiveness of social distancing (PSBB) to reduce human interaction during PSBB stage, leading to reducing the COVID-19 transmission risk. The arguments about whether or not the population density affects COVID-19 risk are heterogeneous in the literature. After human-to-human transmission of the virus was declared, it is tempting to assume that higher density tends to aggravate COVID-19 transmission. This view has been empirically evidenced by previous studies (e.g. Rocklöv and Sjödin 2020; Kadi and Khelfaoui 2020). However, most of these studies did not account for source of crowds or other possible confounding factors (Kadi and Khelfaoui 2020) or used definitions of population density other than total population over total area (Rocklöv and Sjödin 2020). Carozzi (2020) found a similar correlation, but further argued that population density affects the COVID-19 spread more through the timing of the outbreaks, in that having a denser population is associated with higher risk of earlier outbreaks (therefore higher cases over time), and not higher rate of the spread. Furthermore, Hamidi et al. (2020) reported that connectivity of the population was more important than the density, with areas highly linked through socioeconomic relationships to each other being the most vulnerable to the outbreaks. Our data is in agreement with this view, showing that population density is a significant factor, but the sources of crowd provide stronger effects. However, it has to be noted that our data correspond to the first wave of the pandemic and that as the diseases started to spread and more cases are confirmed, it is possible that population density can be more influential at the later stages, as suggested by Wong and Li (2020).

Predicted risk map

Based on the resulting model, maps of the predicted intensity at each time period are built to provide an overview of the predicted COVID-19 transmission risk in each location all over the Surabaya-Sidoarjo region, accounting for the estimated effect of the sources of crowd (enterprise and place of worship) and population densities at New Normal period. The resulting map can serve as a basis for the COVID-19 risk assessment at the early stage.

For the sake of interpretation and recommendation, we provide the district names in the predicted risk map at New Normal period depicted in Fig. 4. The resulting map takes the risk variation within district into account and does not assume a constant risk. The districts associated with highest risk are Sawahan, Simokerto, Tegalsari, Tambaksari, Mulyorejo, Sukolilo, and Gunung Anyar. This perspective could serve as an additional information to the currently used demographic-based risk analysis. This map is expected to be useful for the government to prioritize regions with higher risk of COVID-19 transmission for testing, tracing, and vaccination programs.

The full predicted risk maps for each period depicted in Fig. 5b reflect the fitted intensity in Eq. (5) with the estimated parameter values as in Table 3. The uncertainty associated to these risk maps is inherited from the standard errors of the $\widehat{\mathit{\alpha }}$ estimates (column 5). This means that the uncertainty of our model for risk of infection is kept to low values indicating the accuracy of our model. With our fitted intensity (Fig. 5b), we can also provide in Fig. 5c simulated patterns coming from the corresponding inhomogeneous mark Poisson process. For comparative purposes, we show in Fig. 5 the real data, the risk maps and simulated data from our model. We note that real and simulated data show close patterns of behaviour.

Methodological consideration

Considering point locations data, our study allows to analyze the distribution of COVID-19 at the home address level using a point process modeling framework. As the dynamic movement of COVID-19 infected people is almost impossible to obtain, both for practical and ethical reasons, our analysis covers one of the best proxies to get the spatial information of COVID-19. Within the point process framework, our study is one of few which includes the covariates effect to link with the transmission risk. In addition, by taking the mark into account, we are able to explore how the effect of sources of crowds and population densities alters with different phases of COVID-19 in Surabaya and Sidoarjo. This is indeed of great importance as we can evaluate potential sources of spread and provide evidence for authorities to conduct more focused actions.

Nonetheless, these results must be interpreted with caution, and a number of limitations should be borne in mind. First, some potential confounding factors, such as the location of healthcare facilities, and the number of diagnostic tests per unit time, were not included due to limited data availability. As major cities have more access to healthcare facilities, it is possible that some cases in rural areas were not detected and included in the analysis. For the same reason, more tests are available from time to time, and thus some portion of the increasing confirmed cases over time could be explained by this factor.

Second, the time gap between a real-time infection, case confirmation, data recording and reporting could drive to a potential temporal bias. To address this issue, we separate the point patterns into six times (each period lasts approximately two weeks, see Table 1), classified according to different phases of PSBB, and assume that all cases confirmed in one period happen in that particular period. This strategy might lead to independence among cases within and between periods, depicted by the envelopes of K and cross-K functions (see Figs. 3 and 7), so the degree of clustering as natural pattern of COVID-19 transmission cannot be exploited. However, this setting benefits us to evaluate the spatial trend of COVID-19 cases at different periods as a function of health policy, which is one main focus of this paper. If one has more concern on investigating the clustering pattern, the point pattern should not be divided into different phases, and then a cluster-type point process could be an alternative for modeling. For example, modeling the point point pattern depicted in Fig. 1 using a Thomas cluster point process (e.g. Baddeley et al. 2015) demonstrates a good performance in terms of the envelope K-function (Fig. 6) and presents the clustering pattern described by $\widehat{\omega }$ and $\widehat{\kappa }$ in Table 5 (see e.g. Baddeley et al. 2015; Choiruddin et al. 2021 for the interpretation of the clustering estimates.) Again, this overall model can not offer two-weeks snapshots of the spatial distribution of infected cases in relation to some spatial covariates, aspect that was mainly needed by authorities to drive their restrictions to society.

Competing Interests

The authors declare that they have no conflict of interest.