Bayesian spatiotemporal forecasting and mapping of COVID‐19 risk with application to West Java Province, Indonesia

Abstract The coronavirus disease (COVID‐19) has spread rapidly to multiple countries including Indonesia. Mapping its spatiotemporal pattern and forecasting (small area) outbreaks are crucial for containment and mitigation strategies. Hence, we introduce a parsimonious space–time model of new infections that yields accurate forecasts but only requires information regarding the number of incidences and population size per geographical unit and time period. Model parsimony is important because of limited knowledge regarding the causes of COVID‐19 and the need for rapid action to control outbreaks. We outline the basics of Bayesian estimation, forecasting, and mapping, in particular for the identification of hotspots. The methodology is applied to county‐level data of West Java Province, Indonesia.

. Relevant risk factors are still emerging with great uncertainties, and many are difficult to measure (Subramanian et al., 2020). Hence, modeling COVID-19 has proven to be a complex and difficult endeavor.
There is growing evidence, however, that COVID-19 affects different demographic and socioeconomic groups in distinct ways. Healthcare workers, older age groups, and people with poor health are more susceptible to infection and vulnerable to negative health outcomes than younger people and those with good health (Saadat et al., 2020;Williamson et al., 2020). Risk factors in particular include being over the age of 70 or having medical problems such as chronic respiratory diseases, diabetes, cardiovascular diseases, obesity, and cancer (Fang et al., 2020;Giannisa et al., 2020;Zheng et al., 2020). Socioeconomic factors such as occupation, education, lifestyle, and socioeconomic status have also been found to affect the COVID-19 infection risk (Khalatbari-Soltani et al., 2020;Stojkoski et al., 2020). Weather variables, such as temperature and humidity, and environmental conditions, notably air pollution, have been designated as risk factors that affect transmission . However, other studies reported contradictory results and claimed that the relevance of weather and environmental conditions as risk factors is questionable (Tosepu et al., 2020;Xie & Zhu, 2020). Hence, a comprehensive overview of the determinants of the spatiotemporal distribution of COVID-19 is still lacking. However, even if knowledge of the risk factors is available, it is debatable if updated data is readily available for estimation and prediction using complex models. According to M. Wang and Flessa (2020), the formulation and computation of even the best COVID-19 models are under severe uncertainty. This applies particularly to developing countries such as Indonesia, where data availability and accuracy is a common issue (Djalante et al., 2020).
The impacts of the unobserved factors affecting COVID-19 risk transmission can be inferred from changes in the number of incidences in space and time captured by hierarchical random effects statistical models, denoted as pure models (Lopez-Quılez & Munoz, 2009). As in time series analysis, when there is lack of solid information on the covariates or uncertainty in the model specification, univariate, that is, pure models are commonly applied as alternatives to multivariate models to generate short-term forecasts (Chayamaa & Hirata, 2016). Additionally, Martinez-Beneito and Botella-Rocamora (2019) observed that random effects models are a suitable tool for modeling correlated data in general. A spatiotemporal pure model consists of structured and unstructured spatial and temporal random effects and their interactions, thus capturing the unobserved factors driving the spatiotemporal distribution of the relative risk (Jaya & Folmer, 2020a;Kazembe, 2007;Wakefield, 2007). A hierarchical pure model provides robust estimation of the spatiotemporal relative risk parameters. However, care is required in interpreting the model, because it may not reflect any causal mechanism due to the fact that it does not contain explicit information about the risk factors (Wakefield, 2007). It is an exercise in mapping and forecasting in the first place. The spatiotemporal pure model introduced in this paper is a generalization of the pure time series model, although in a Bayesian setting. As its time series analog, it generates unbiased and precise forecasts with minimal data requirements. It is generally accepted that when the goals are mapping and forecasting, a simple model without exogenous variables suffices as a first approximation (Leone, 1979;Naylor et al., 1972;Preez & Witt, 2003). Even more so, the more parsimonious the model, the wider its applicability.
The pure model does not contain explicit information on the risk factors. It can, however, be used to produce maps to generate new hypotheses on the disease etiology, particularly to identify relevant risk factors (Huque et al., 2016;Jaya & Folmer, 2020a;Wakefield, 2007). Note that great care is required in modeling covariates at an early stage of research on an unexplored topic because of uncertainty about nonlinear or interacting effects of the covariates (Wakefield, 2007).
Rather than just observing relative risk estimates, it is important to define "unusual" counties where public intervention is needed (Aguayo et al., 2020). The unusual counties are denoted "hotspots," that is, isolated areas in which the relative risk is unusually high. Bayesian methods are convenient to define hotspots via a classification decision rule. The rule defines a reference threshold on the posterior probability distribution of the relative risk, and a cutoff value on the posterior probability distribution that the relative risk is above the reference threshold. Then, a county is classified as a hotspot if the posterior probability that the relative risk is above the reference threshold is greater than the cut-off value (Richardson et al., 2004). Hotspots can provide clues to disease etiology, suggesting, for example, local environmental or social characteristics that promote increased risk (Mclafferty, 2015).
Herein, we present a parsimonious county-level space-time model and corresponding maps of the distribution of COVID-19, outline the basics of Bayesian mapping and forecasting, and apply it to West Java Province, Indonesia, as an illustration. As a first step, we estimated several Bayesian spatiotemporal models and selected the most appropriate one using Bayesian evaluation criteria. Next, we interpreted the posterior means (intercept, random effects) and standard errors which were used to plot maps for the relative risk by county and week. The maps were subsequently used to identify high-risk counties (hotspots) using the selected exceedance probability.
The remainder of this paper is structured as follows. In Section 2 we present the spatiotemporal pure model to describe and forecast the relative risk distribution. Section 3 presents the estimations, forecasts, and maps of COVID-19 relative risk for West Java, and Section 4 contains the discussion and conclusions. Supporting information on methods and material, including supplements on data and detailed relative risk estimates for each county, is presented in Appendices A and B.

| THE RELATIVE RISK AS A PURE MODEL
The spatiotemporal pure model presented below only requires the number of incidences and population size per region and period. It assumes that the number of COVID-19 incidences in county i in period t, y it , follows a Poisson distribution (or, in the case of overdispersion when the data contains a large number of zeros, a negative binomial (NB) distribution (Berk & MacDonald, 2008;Payne et al., 2017) or a zero-inflated Poisson (ZIP) distribution (Lewsey & Thomson, 2004) or a zero-inflated negative binomial (ZINB) distribution (Agarwal et al., 2002) as follows: where λ it denotes the mean and variance of y it . To capture differences in the population at risk, λ it is divided by the expected number of incidences, E it defined as follows (Jaya & Folmer, 2020a): 1, , and 1, , ,  (Abente et al., 2018;Jaya & Folmer, 2020a). Accordingly, y it follows a Poisson distribution with the mean and variance equal to λ θ = E it it it , as follows: The ratio of the observed number of incidences (y it ) to the expected number of incidences (E it ) is the SIR. It is the unbiased maximum likelihood estimator of the relative risk (Jaya et al., 2017) defined as 1, , and 1, , .
Although the SIR is an unbiased estimator of the relative risk, it can be unreliable for sparse count data.
Particularly, a spatiotemporal unit with a small number of observed incidences with a small population size, and hence a low expected number of cases, may be incorrectly classified as a high-risk unit (Jaya & Folmer, 2020a;Yin et al., 2014). To reduce the amount of variation due to population size and heterogeneity, smoothing can be applied (Jaya et al., 2017). Spatiotemporal smoothing effectively borrows strength across spatiotemporal units by introducing spatiotemporal dependence and heterogeneity into the Poisson or NB regression model (Yin et al., 2014). Taking the natural logarithm of λ it yields The offset in Equation (5) is assumed to have a regression coefficient fixed at 1.
To analyze the spatiotemporal distribution of COVID-19, we consider the loglinear relative risk, and model it as a generalized linear mixed model accounting for the spatiotemporal dependence and heterogeneity as follows (For inclusion of age-structure effects, see Appendix A1): Here, η θ α = log( ), it it denotes the intercept representing the overall relative risk; ω υ ϕ , , , i i t and φ t are the spatially structured (spatial autocorrelation), spatially unstructured (spatial heterogeneity), temporally structured (serial autocorrelation), and temporally unstructured (temporal heterogeneity) random effects, respectively, and δ it is the interaction of a pair of the four random effects above. The random effects implicitly capture the determinants of the relative risk (Jaya & Folmer, 2020a;Kazembe, 2007).
We propose Bayesian inference to analyze model (6) which takes all the parameters as random variables with prior probability distributions. We propose a vague Gaussian prior distribution with a zero mean and a large variance for α α :~(0, 10 ) 6 N (Martinez-Beneito & Botella-Rocamora, 2019). In Bayesian (disease) mapping, spatial autocorrelation in the data is commonly modeled at the second level of a hierarchical model by a set of random effects with CAR prior distributions belonging to the class of Markov random fields. Besag et al. (1991) proposed the intrinsic and convolution CAR priors. However, these CAR priors do not allow clear distinguishing between the structured and unstructured random effects, which are thus not identifiable (MacNab, 2011). Leroux et al. (2000) proposed the Leroux before make the difference between the unstructured and structured spatial variation explicit. This prior is also more flexible in representing a range of spatial correlation scenarios (Lee, 2011). Through a simulation study, Aswi et al. (2020) found that spatial estimation based on the Leroux prior performed best on the Watanabe-Akaike information (WAIC) criterion. Krisztin et al. (2020) proposed to use spatial econometrics to quantify spatial spillover based on the SAR model (in spatial econometrics known as the spatial error model).
Rather than using latent random effects as in the Bayesian approach, spatial econometrics models spatial dependence as an explicit spatial autoregressive error structure.
In addition to the LCAR prior proposed by Aswi et al. (2020), we also considered the CAR and SAR priors for the application in Section 3. The results presented in Table B4 show that the LCAR prior performed best. For the spatially structured random effect (ω) of county i across time t the LCAR prior reads as follows (Leroux et al., 2000): where = w W ( ) ij is the spatial weights matrix, ρ is the spatial autoregressive parameter, and σ ω 2 the variance share a vertex or border, and = w 0 ij otherwise) are generally an optimal choice for modeling spatial spillover and spatial smoothing. Nevertheless, we compared it with various types of weights matrices based on the geographical distance d ( ) ij between areas i and j measured as the Euclidean distance between their respective centroids. The results are presented in Table B5. Table B5 shows that there is no uniformly best model. However, the first-order queen contiguity spatial weights matrix with ρ = 0.347 is the best, or close to the best, weights matrix. Hence, we selected it for further analysis.
A random walk prior of order one (RW1) or order two (RW2) is appropriate for the temporal trend we propose the following exchangeable prior: where σ φ 2 is the variance hyperparameter of φ. For the four interaction components δ , the priors are the products of the priors of the corresponding components (see Boulos & Geraghty, 2020;Jaya & Folmer, 2020a;Knorr-Held, 2000 for details).
In addition to the priors above, the variance hyperparameters  (Gelman, 2006).

Hence, it follows that the number of COVID-19 incidences
defined as random variables with priors Ω ψ p ( | ) and hyperparameters The Bayesian pure spatiotemporal Model (6)  In a Bayesian setting, forecasts of the relative risk are based on the posterior predictive distribution ) p y|y (( X. Wang et al., 2018). In INLA, prediction can be conveniently accomplished by fitting a model with missing observations. Specifically, one combines past observations from the previous periods and missing or not available (NA) observations for the periods one wants to forecast. For details on forecasting using INLA see Appendix A2.
Choropleth maps are commonly used to visualize the spatiotemporal distribution of the relative risk of a disease (Indrayan & Kumar, 1996).
To identify spatiotemporal hotspots, we employed exceedance probability criteria (Lawson, 2010). The spatiotemporal exceedance probability can be computed from the posterior spatiotemporal distribution of the relative risk. It is defined as the probability that the estimated posterior mean of the relative risk of area i at time t is higher than a threshold value c, that is, θ > c|y Pr ( ).
it It is estimated as is the cumulative probability of θ it with threshold value c. It can be estimated using the Laplace approximation (Blangiardo & Cameletti, 2015). The use of the exceedance probability to identify hotspots requires two parameters that must be fixed a priori. The first one is the threshold value c for θ it . The value 1 indicates that a county has average relative risk, whereas values such as 2 or 3 indicate extreme risk. The second parameter is the cut-off value γ of the exceedance probability. Common values for γ are 0.90, 0.95 and 0.99 (Lawson & Rotejanaprasert, 2014). See Jaya and Folmer (2020b) for a detailed discussion of hotspots.

| SPATIOTEMPORAL MAPPING AND FORECASTING OF COVID-19 RELATIVE RISK IN WEST JAVA
The data were obtained from Statistics West Java (2019) and COVID-19 Information Center West Java (2020), respectively. The following observations apply. First, the number of tests may not be constant overtime. Indeed, it increased in many, especially industrialized countries as the duration of the pandemic progressed (Papastefanopoulos et al., 2020;Subramanian et al., 2020). The number of confirmed COVID-19 cases tends to increase with the number of tests (Bertozzi et al., 2020;Roda et al., 2020). However, the number of tests in West Java over the study period (March 06-July 09, 2020) was approximately constant at around 25,000 tests per week due to limited testing capacity (Nabila, 2020). Second, we did not take into account the age effect, although there is growing evidence that the relative risk and impacts including mortality vary by age (Goldstein et al., 2021;Mallapaty, 2020). Particularly, older people have a substantially larger chance of getting infected and developing symptoms (Mallapaty, 2020). We did not include the age effect because our data refers to the beginning of the pandemic in Indonesia, at which point incidences were not registered by age (Djalante et al., 2020).  Figure 1b shows a rapid increase in the number of confirmed incidences from Week 1 till Weeks 9-11 (with a dip in Week 10), followed by a substantial decline and leveling off from Week 12 till Week 18 with a short peak in Week 17. Figure 1c shows that the outbreak started in the populated counties in the northeast and center, and started rapidly spreading to their surrounding counties. From Week 12 onward, the disease started leveling off in the periphery whereas the populated counties and, to a lesser extent, their neighboring counties remained affected.
To identify the most appropriate spatiotemporal model of COVID-19, we estimated 48 different submodels of Model (4) with the RW1 and RW2 trends for the Poisson and NB distributions and their zero-inflated variants. We compared their fits using the DIC and WAIC, and evaluated their predictive performances using the MPL, MAE, RMSE, and r based on cross-validation with the first 14 weeks as the training set and Weeks 15-18 as the testing set. The results are presented in Table B2.
Apart from a few minor exceptions, the evaluation statistics of the RW1 submodels outperformed their RW2 counterparts. Hence, the latter were not considered any further. Similarly, the ZIP and the ZINB distributions were rejected in favor of the Poisson and NB distribution, respectively, thereby reducing the selection of the optimal model to the six variants of the Poisson and NB models. As a next step, we assessed overdispersion by means of the deviance index of the Poisson models and the overdispersion parameter of the NB models (Mohebbi et al., 2014). A deviance index greater than one is evidence of overdispersion. Table B3 shows that the deviance index of the six Poisson models ranged from 1.49 to 2.47, and that the overdispersion parameters of the NB models were larger than zero. Therefore, the Poisson model was rejected in favor of the NB model. counties. (We also considered the CAR and SAR prior and several distance-based contiguity matrices as alternatives to the first-order queen contiguity matrix. The results are presented in Table B4).
The contribution of each random effect and the interaction term to the total variance, FV h , is calculated as ,ˆ,ˆ,ˆ,2 2 2 2 2 }. Table 1 shows that the structured spatial (ω) and unstructured spatial random effects (v) account for (0.139 + 0.069) × 100% = 20.8% of the spatiotemporal COVID-19 variation, the structured (ϕ) and unstructured temporal effects (φ) for (0.710 + 0.057) × 100% = 76.6%, and the interaction effect (δ ) for 2.6%. Accordingly, the temporal effects dominate the spatiotemporal variation of COVID-19.  Table 1. A relative risk larger than one means that the corresponding posterior mean is larger than average across space and time.  intensity (relative risk >1, in red) in one county in the northwest in Week 3, then began spreading with light (relative risk 0.5-0.8, light green) to high intensity in the northwest, and emerged in the center in Week 4. During Weeks 5-12, it spread further and intensified in the latter regions while pockets emerged in the southwest. During Weeks 13-18, a substantial fluctuation occurred. Initially, a rapid decline was observed in the northwest and center, followed by a gradual increase in these areas for the remainder of the period, ranging from light to high intensities. Figure 2b shows that for Weeks 19 and 20, a significant decrease was predicted in these counties, except in the northwest counties and in the center where the outbreak started.
In terms of policy, the identification of hotspots is crucial because, as epicenters, close monitoring and effective managing of them are required to prevent proliferation. Figure  Note the difference between Figures 2 and 3 with respect to counties with relative risk probability larger than one. Taking Week 20 as an example, Figure 2 indicates that, for example, the county in the center has relative risk probability larger than one (colored red) but Figure 3 shows that it is not a hotspot. The explanation is that the exceedance probability is smaller than 0.95 and thus does not qualify as a hotspot.  Tables B6 and B7 for detailed estimates and predictions). In the remainder of the province, the relative risk was low (less than one) and no hotspots were discovered. Note that the intensive commuting between Jakarta and the northwest and central West Java counties was hypothesized but not estimated because Jakarta does not belong to West Java Province.
The concentration of COVID-19 in the northwest and center necessitates local measures, that is, inter-county lockdowns, in particular restrictions for commuting to Jakarta as well as limitations for interactions of the hotspots with the other, not (yet) affected West Java counties. The persistence of the hotspots necessitates intra-county lockdowns, such as closing public facilities, offices, and public transportation, limiting religious gatherings, social, and cultural activities, mandating wearing of face masks, and imposing social distancing measures. However, a province-wide intra-county lockdown is not required.
Based on Figures 2 and 3, we noticed that the high-risk areas are dominated by urban areas (cities) with high population density. To further examine the relationship we calculated the Pearson correlation between the population by land area in square kilometers and the average predicted relative risk over the study period. The result is presented in Figure 4. The Pearson correlation coefficient of 0.78 supports the hypothesis of a close relationship.
This finding is in line with the growing evidence elsewhere that COVID-19 infection is more likely in urban areas, which are densely populated, than in rural areas, which are less populated (Bhadra et al., 2020;Kadi & Khelfao, 2020;Saadat et al., 2020).
Regarding mobility and relative risk we refer to Table 1, where we observed a strong structured temporal effect indicating that the relative risk is an intra-county issue in the first place. Regarding the structured spatial effect, we found a Leroux conditional autoregressive distribution with first-order queen contiguity matrix and spatial autoregressive parameter ρ = 0.347, indicating that there is low spatial spillover among neighboring regions with a maximum distance of 42 km. The low spatial spillover effect is likely to be related to the fact that the city of JAYA AND FOLMER | 859 Jakarta-the epicenter of the COVID-19 disease in Indonesia-was not included in the analysis whereas there is intensive commuting between Jakarta and the provincial hotspots in the north-east and the provincial capital Bandung in the center. Other reasons for the low spatial spillover effect are the large-scale social restrictions (PSBB) that the provincial government of West Java implemented from April 22 to July 2, 2020 (Pangestika, 2020a(Pangestika, , 2020b. PSBB allows the county governments to restrict the mobility of people and goods to and from their locations (Andriani, 2020).

| DISCUSSION AND CONCLUSIONS
An early warning system for the control and management of the COVID-19 pandemic requires mapping of the spatiotemporal dynamics of the disease and forecasting of the locations and times of future outbreaks. Herein, we propose the relative risk indicator as the key concept of such a warning system as it accounts for differences in the population at risk, such as the population size, which varies across space and time. Relative risk is based on the assumption that the number of incidences follows a Poisson distribution (or the Negative Binomial (NB) distribution or the zero-inflated variants of the Poisson or NB distribution for the case of overdispersion). We observed that registration and analysis by age is a prerequisite to ensure a good understanding of the pandemic. The methodology presented in this paper can be straightforwardly adapted for characteristics of the population at risk, such as the age or sex structure.
We modeled the relative risk as a spatiotemporal pure model of the structured and unstructured spatial and temporal random effects and their interactions to capture unknown or unobserved risk factors. The model is parsimonious because it only requires the number of cases and the population at risk across space and time. Additionally, the pure model allows adequate mapping and forecasting of future outbreaks by location (county) and time (week).
A frequently occurring problem at fine spatiotemporal scales is that the estimators of the parameters have large standard errors because of data sparseness or spatiotemporal heterogeneity. The problem can be solved by smoothing, that is, by borrowing strength across spatiotemporal units such that the observations for a specific spatiotemporal unit are weighted toward the observations of neighboring units.
We demonstrated that the pure model can be conveniently estimated using Bayesian approaches, thereby allowing Leroux conditional spatial autoregressive based on first-order queen contiguity, and random walk temporal priors to capture spatiotemporal dependence and heterogeneity. In addition, the Bayesian approach enables the relative risk uncertainty to be quantified via the posterior distributions, particularly the posterior exceedance probability (Adin et al., 2019) which can be used to identify hotspots, that is, isolated areas where the relative risk is unusually high. Hotspots are candidates for policy intervention.
The application to the West Java counties for the period March 6-July 9, 2020 demonstrated that the spatiotemporal distribution of COVID-19 was primarily temporal, while there was relatively little spatial interaction between the West Java counties. The persistence of the hotspots over time indicated the necessity for intracounty lockdowns. The pure model can be straightforwardly used to identify hotspots within a hotspot, such as districts or residential blocks. Local lockdowns or the imposition of mobility restrictions are effective policy instruments to combat the spread of the disease within a county.
The relatively low impact of the spatial component was probably due to (i) the fact that the city of Jakarta-the epicenter of the COVID-19 disease in Indonesia-does not belong to the province of West Java, whereas there is intensive commuting between Jakarta and the provincial hotspots in the bordering counties in north-east West Java and the provincial capital Bandung in the center and (ii) large-scale social restrictions (PSBB) that West Java implemented from April 22 to July 2, 2020.
Taking into consideration the distance that the disease spread from Jakarta to the provincial hotspots, restrictions on ingoing and outgoing commute to Jakarta were needed. In a similar vein, restrictions on mobility between the West Java hotspots and the other counties in the province including low-risk counties were necessary. Further analysis of the relative risk by county revealed that there is a close relationship with population density in line with the hypothesis that COVID-19 infection is transmitted via interaction of noncontaminated people with infected people or via contact with contaminated objects and surfaces.
There are several limitations to the present study related to the data quality. First, this study is based on COVID-19 cases reported by county. Due to different monitoring standards, this may have led to differences in incidence recording. Uniform data collection using international standards, multidisciplinary research, and combining statistical estimates with the outcomes of laboratory and field experiments are required for future research. Second, although the number of tests in West Java over the study period (March 06-July 09, 2020) was approximately constant at around 25,000 tests per week, it is important for future research to take changes in testing intensity into account.
Third, improved monitoring could be achieved using longer time series data thus accommodating the different stages of the evolution of the pandemic (Dobricic et al., 2020). This kind of information presents clues to the effects of behavioral changes such as social distancing, working from home, closing schools, closing public facilities, offices, and public transportation, limiting religious gatherings, social events, and cultural activities, and wearing face masks.
Finally, the analysis was restricted to the Province of West Java and excluded the heavy regular travel to Jakarta, the national epicenter of COVID-19. An adequate early warning system requires data collection in functional regions and, if necessary, policy coordination among administrative regions (Yang, 2017).
In spite of the above-mentioned shortcomings, the results of this study can still contribute to pandemic-related policymaking at the county level because the inaccurate reporting of confirmed cases applies to each county.
Additionally, predictions were based on the notion of relative risk, which is a more accurate measure to identify highrisk counties than the number of confirmed cases, especially in the case of rare or sparse data (Tango, 2010;Tenny & Hoffman, 2020). In addition, supplementary analysis such as on the theoretically plausible and elsewhere confirmed correlation between population density and the location of the hotspots gives further support to the findings.
In summary, the hierarchical Bayesian pure spatiotemporal model is an easy to apply, widely applicable model for the accurate forecasting of short-term spatiotemporal COVID-19 relative risk. Additionally, it can be used to generate and test hypotheses regarding control, such as environmental or socioeconomic conditions.

ACKNOWLEDGMENTS
We thank the West Java province COVID-19 task force for providing the data. Financial support was received from Universitas Padjadjaran, Bandung, Indonesia.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

A1. Age-structured effects
If age-structured count data is available, the methodology presented in this paper can be straightforwardly adapted to include the age-structured effects as follows. Assume that the region under study is divided into n counties and T periods with K age groups. Let N itk denote the population at risk in county i at time t, and age group k and y itk the number of confirmed COVID-19 cases. Accordingly, y itk follows a Poisson distribution with the mean and variance equal to λ θ = E itk itk itk , as follows: with E itk and θ itk the expected number of cases and the relative risk in area i at time t for age group k, respectively.
The expected number of cases E itk is defined as (Wakefield, 2007): with p k the probability of COVID-19 for age group k defined as  Here, η θ = log( ) itk itk . The additional components β k , ∂ , ik tk g , and ζ ijk denote the age-structured effect for age group k, age-space interaction, age-time interaction, and age-space-time interaction, respectively. The last components ζ ijk is rarely used because it greatly increases the complexity of the model (Lopez-Quılez & Munoz, 2009). The age-structured effects can be modeled through a fixed effects or random effects model. Sun et al. (2000) and Xia and Carlin (1998) incorporated the age-structured effects as fixed-effects covariates through dummy variables. On the other hand, Knorr-Held and Besag (1998) and Goicoa et al. (2016) treated the agestructured effects as random effects, and applied the pure spatiotemporal model (Schmid & Held, 2004). According to Adin et al. (2019) and Jaya and Folmer (2020b), the random effects model is more appropriate to model agespecific effects due to the identifiability problems of the fixed effects model caused by the large number of agegroup variables. This also applies to other group-specific effects. In a Bayesian setting, the random effects component of the age-structured effects β k can be modeled using random walk or exchangeable priors (Goicoa et al., 2016;Knorr-Held & Besag, 1998). (Note that the age-groups are indexed with an integer so that they can be directly modeled using random walk or exchangeable priors). A random walk prior of order one (RW1) β k reads: )is the likelihood function of the number of cases, = … ′ = … i n y (y , , y ) , 1, , nT 11 , and = … t T 1, , which is based on the Poisson or Negative Binomial distribution, Ω ψ | p ( ) is the prior distribution of Ω given hyperparameter vector ψ, ψ p ( ) the joint hyperprior distribution of the hyperparameters ψ and ψ p y ( | )is the marginal likelihood of the model which is a normalizing constant and usually ignored because of computational costs. The number of cases y is assumed independent given the parameters Ω and hyperparameter ψ. The likelihood function is then: where each data point y it is connected to only one element )is assumed to be the normal distribution. Hence: Because Ω kit is the posterior mode of̃Ω    owing to the property that the integral of a density function over the entire space yields 1. Hence, Equation (A2.16) becomes: Because the expectation of a Gaussian distribution is equal to the mode we have: | )ˆ, 1 , 6, 1 , a n d 1 , , . Forecasting in a Bayesian framework is out-of-bounds prediction based on the posterior predictive distribution of the unobserved quantity, + ŷ T h ( ) , given the observed data y, with = … h 1, 2, steps ahead from T . In particular (Jaya & Folmer, 2020a): where the marginal posterior distribution for a h-step-ahead forecast ψ  residual spatiotemporal autocorrelation (see Jaya and Folmer (2020b) for details on the MoranST). We ran the 15 models in INLA, and report the results in Table B5. The deviance index is defined as: