key: cord-0681911-t6sd71b5
authors: Jaya, I. Gede Nyoman Mindra; Folmer, Henk
title: Spatiotemporal high-resolution prediction and mapping: methodology and application to dengue disease
date: 2022-02-19
journal: J Geogr Syst
DOI: 10.1007/s10109-021-00368-0
sha: 2ab332bb47653cdaae0febdbebc95488a561a3d3
doc_id: 681911
cord_uid: t6sd71b5

Dengue disease has become a major public health problem. Accurate and precise identification, prediction and mapping of high-risk areas are crucial elements of an effective and efficient early warning system in countering the spread of dengue disease. In this paper, we present the fusion area-cell spatiotemporal generalized geoadditive-Gaussian Markov random field (FGG-GMRF) framework for joint estimation of an area-cell model, involving temporally varying coefficients, spatially and temporally structured and unstructured random effects, and spatiotemporal interaction of the random effects. The spatiotemporal Gaussian field is applied to determine the unobserved relative risk at cell level. It is transformed to a Gaussian Markov random field using the finite element method and the linear stochastic partial differential equation approach to solve the “big n” problem. Sub-area relative risk estimates are obtained as block averages of the cell outcomes within each sub-area boundary. The FGG-GMRF model is estimated by applying Bayesian Integrated Nested Laplace Approximation. In the application to Bandung city, Indonesia, we combine low-resolution area level (district) spatiotemporal data on population at risk and incidence and high-resolution cell level data on weather variables to obtain predictions of relative risk at subdistrict level. The predicted dengue relative risk at subdistrict level suggests significant fine-scale heterogeneities which are not apparent when examining the area level. The relative risk varies considerably across subdistricts and time, with the latter showing an increase in the period January–July and a decrease in the period August–December. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s10109-021-00368-0.

Dengue disease is a major challenge to healthcare worldwide, potentially leading to death, especially among the poor in low-and middle-income countries (Ak et al. 2018 ). In addition, there are latent costs related to stress, productivity loss, school absence, and care taking (Wilastonegoro et al. 2020) . To reduce health impacts and treatment costs, there have been substantial efforts aimed at the prevention of dengue disease outbreaks (Kampen et al. 2014) . For these purposes, statistical models aimed at identifying the causes, transmission mechanisms, and prediction of outbreaks at a fine spatiotemporal scale are of crucial importance (Messina et al. 2019 ).

High-resolution information is needed for research on the etiology of the disease and the development of control and prevention strategies (Ak et al. 2018; Jaya et al. 2017; Folmer 2020, 2021a, b; Pokharel and Deardon 2016) . Dengue disease typically involves a spatiotemporal pattern Folmer 2020, 2021a, b; Phanitchat et al. 2019; Puggioni et al. 2020) ; therefore, high-resolution spatiotemporal models and maps of the distribution of relative risk are basic elements of an early warning system aimed at identifying where and when an outbreak will occur (Hanigan et al. 2019; Shi et al. 2013; Xu et al. 2019) .

In many areas of spatial research including epidemiology, data are often only available at different levels of aggregation (Moraga et al. 2017; Shi et al. 2013; Utazi et al. 2019) . If low-resolution information is available whereas high-resolution (cell level) information is needed but lacking, area-to-cell disaggregation through joint area-cell estimation can be applied to obtain the missing information (Moraga et al. 2017; Utazi et al. 2019; Wang et al. 2018a ). 1 Data fusion or data assimilation (Banerjee et al. 2015) and Bayesian melding (Fuentes and Raftery 2005; Liu et al. 2011) are terms used to denote integrating multiple data sources of different spatial resolutions. 2 The basic concept involves combining area and cell observations in a single statistical model. Combining data measured at different levels of aggregation can improve parameter estimation and increase prediction accuracy (Wang et al. 2018b ). However, it may lead to spatial misalignment (Moraga et al. 2017; Sahu et al. 2010; Truong et al. 2014; Utazi et al. 2019 ), which induces biased or inconsistent estimators (Liu and Bertazzon 2016; Peng and Bell 2010; Saez and López-Casasnovas 2019) .

Several kinds of correction approaches based on regression models have been developed to deal with area-to-cell misalignment problems (Banerjee and Gelfand 2002; Banerjee et al. 2015) . Moraga et al. (2017) and Utazi et al. (2019) applied Bayesian geostatistical analysis to deal with misalignment in spatial non-Gaussian 1 Below we use the notion of "cell" to denote spatial units of higher resolution. 2 Joint area-cell estimation is based on the assumption that the same factors at area and cell level drive the spatiotemporal process, although usually in various degrees (Banerjee et al. 2015; Utazi et al. 2019 ). In the present case, the same disease risk factors (weather) drive disease incidence at area and cell levels. data with linear covariates. However, methods that address area-to-cell misalignment in spatiotemporal non-Gaussian data with nonlinear covariates are less well known. This applies especially to Poisson or Negative Binomial (NB) spatiotemporal data, which are typically applied in dengue and other disease incidence modeling, but also in other kinds of spatial and regional research.

In this paper, we introduce the Fusion Area-Cell Spatiotemporal Generalized Geoadditive (GG)-Gaussian Field (GF) model, abbreviated as FGG-GF model, to generate high-resolution (cell) predictions based on observations at lower (area) resolution of the variable of interest (i.e., the number of dengue incidences) and the population at risk, and high-resolution cell data (i.e., weather variables), while controlling for misalignment. Moraga et al. (2017) and Utazi et al. (2019) have shown that the prediction performance of integrated area-cell models can outperform models that use single level data sources. Kammann and Wand (2003) introduced the Generalized Geoadditive Model (GGM), which has become popular in disease mapping (among other fields) because of its suitability for making high-resolution maps (Muleia et al. 2020; Wand et al. 2011) . The model assumes that there is a spatially continuous variable underlying all observations which can be modeled using a Gaussian process, usually denoted Gaussian Field (GF). A GF is characterized by a first-order autoregressive model with spatially correlated innovations. The GGM combines the Generalized Additive Model (GAM) and the Geostatistical Model (GM). The former was introduced by Hastie and Tibshirani (1986) to provide a flexible means of handling nonlinear and interacting covariates. GAMs are also suitable for handling complex spatial and temporal autocorrelation (French and Wand 2004; Ma et al. 2014) . GAMs are nonparametric because they do not require a priori specification of the regression function (Wang et al. 2018a) . The GM was introduced by Matheron (1963) to construct high-resolution maps over a particular geographical region based on cell data on (risk) factors associated with a (dependent) variable of interest.

The integrated area-cell observations and the combination of non-Gaussian data, a nonlinear predictor and latent model components, in particular the spatiotemporal GF, make estimation of FGG-GF model, prediction, and mapping computationally complex and time-consuming (Barber et al. 2016 ) because of the "big n" problem. This issue can be handled by transforming a GF with a dense covariance matrix to a Gaussian Markov Random Field (GMRF) with a sparse precision matrix 3 of Matérn covariances (Lindgren et al. 2011) .

The objective of this paper is to develop a high-resolution prediction and mapping procedure for spatiotemporal Poisson or Negative Binomial data applying a Fusion Area-Cell Spatiotemporal Generalized Geoadditive-Gaussian Markov Random Field (abbreviated as FGG-GMRF) model. Inference and prediction are handled in a Bayesian framework. The approach will subsequently be applied to dengue disease risk in Bandung city, Indonesia. The purpose is to predict and map the relative dengue risk at subdistrict level, given observations on dengue incidence and population at risk at district level and weather risk factors at cell level. 4 Special attention is paid to high-risk districts and subdistricts requiring public intervention (Aguayo et al. 2020) .

The structure of the remainder of this paper is as follows. Section 2 introduces the spatiotemporal GG-GF model. Section 3 presents the FGG-GF and FGG-GMRF models and the Bayesian inference framework. The link between the GF and GMRF models is summarized in Appendix 1. Section 4 applies the methodology to dengue incidence in Bandung city, Indonesia, and Sect. 5 summarizes and concludes the conducted research.

Consider region A ∈ ℝ 2 , partitioned into n A areas (e.g., districts in a city), each measured for T periods. The areas are labeled A 1 , 1 , … , A 1 , T , A 2 , 1 , … , A i , t , … , A n A , T , where A i , t denotes area i at time t, for i = 1, … , n A and t = 1, … , T . Region A is further divided into a finite set of n p cells for T periods. The set of n p cells over T periods is denoted s 1 , 1 , … , s n A 1 , T , s n A 1 +1 , 1 , … , s g , t , … , s n p , T , with n A i denoting the number of cells in area A i for g = 1, … , n p and t = 1, … , T and n p =

. Note that the notation s 1(A 1) , 1 , … , s n A 1 (A1) , T , s n A 1 +1(A 2) , 1 , … , s g(A i) , t , … , s n p A n A , T will be used to explicitly denote that cell s g(A i) belongs to area A i . If the area is not relevant, s g will be used. Moreover, the notation g for g will be incidentally used if there is no risk of misunderstanding. Finally, s g,1 and s g,2 denote the latitude and longitude coordinates of its centroid, respectively. Let y it and N it denote the number of (dengue) incidences and population at risk in area A i at time t , respectively, and y gt and N gt the number of (dengue) incidences and population at risk in cell g at time t , respectively. Note that both y gt and N gt are unobserved at the cell level. y it and y gt are assumed to follow Poisson distributions 5 with means it = E it it and gt = E gt gt , respectively, with E it and E gt denoting the expected number of (dengue) incidences and it and gt the relative (dengue) risk for area i and cell g at time t , respectively (Jaya and Folmer 2020. The expected rate E it is calculated using external standardization. It is defined based on the overall average across all areas and periods (Abente et al. 2018; Folmer 2020, 2021a, b) :

The relative risk is defined as the ratio of the local risk in a spatiotemporal unit relative to the average risk across the whole study region over the entire time period (Yin et al. 2014) . It is centered around one, meaning that the total number of incidences is equal to the expected rate. The maximum likelihood (ML) estimator of the relative risk it is (Jaya et al. 2017; Jaya and Folmer 2020): This is known as the crude risk or the standardized incidence ratio (SIR). Following Moraga et al. (2017) , Folmer (2020, 2021a, b) , and Utazi et al. (2019) , we model the relative risk as a non-separable Poisson log-linear model as follows 6 :

with it = log it and gt = log gt . In Eqs. (3a) and (3b), 0 is the overall intercept denoting the average risk across space and time, i.e., across all i = 1, … , n A , g = 1, … , n p , and t = 1, … , T . The latent 7 functions f k x k,it and f k x k,gt for k = 1, … , K , represent the (non)linear effects of the metrical area and cell risk factors, respectively. The risk factors at cell level for a given area A i and time t are fixed. However, they vary across areas and times. The latent (non)linear risk factor functions are based on observations at cell level but are predicted at area level. The risk at cell and area levels are assumed to be driven by the same factors; therefore, we adopt joint risk factor functions. For this purpose, we stack the observations on the risk factors such that risk factor k , at both (1)

for g = 1, … , n p and t = 1, … , T 6 Spatiotemporal variation (in disease outcomes) consists of the following three components: a spatial component capturing the overall spatial distribution, a temporal component capturing the overall temporal pattern, and a space-time component capturing space-time interaction. Spatiotemporal models are classified as separable or non-separable (Knorr-Held 2020). A spatiotemporal model is said to be separable if it consists of spatial and temporal components without their interaction, while a non-separable model comprises all three components (Haining and Li 2020; Knorr-Held 2000; Martinez-Beneito and Botella-Rocamora 2019) . 7 The regression functions are latent in that they do not have pre-specified functional forms. This feature makes latent regression functions suitable to accommodate complicated relationships between the risk factors and the dependent variable. area and cell level, becomes k = x k,11 , … x k,n A T , x k,11 , … , x k,n p T � for k = 1, … , K and latent function f k k . The functions f k z k are commonly centered at the mean, i.e., f k z k = 0, for identifiability reasons (Fahrmeir and Lang 2001) .

Let f (z) be the sum of the functions f k (z k ) for k = 1, … , K:

To account for spatiotemporal variation, Eq. (4) can be extended to a varying coefficients model 8 :

where the design vector = 1 , … , K ′ contains components of z or additional covariates. The vector k for k = 1, … , K , modifies the relationship between the covariate z k and the log-linear conditional expectation |z . If it is identical to the vector 1, i.e., k = (1, … , 1) � with dimension ( n A + n p )T × 1 , then f k (z k ) presents the overall (main) effect of z k . If it is different from 1, f k z k k presents the effect of k that varies along with k . In other words, f k z k k models the interaction between k and k (Fahrmeir and Lang 2001) . According to Martınez-Bello et al. (2017a; , the varying coefficients model helps refine the association between the regressors (e.g., the weather variables) and the response, thus improving predictions at a fine spatiotemporal scale. For example, if k denotes the calendar day and k is the spatiotemporal covariate temperature, then f k z k k represents the temperature effect varying by day. In this study, we apply the temporally varying coefficients model to accommodate the temporally varying nonlinear effects of risk factors on the response. The time-varying effect of, for example, the k-th covariate can be written as (Franco-Villoria et al. 2019): where k ( k,t ) for t = 1, … , T is the time-varying regression coefficient, which can be regarded as a stochastic process over k (Fahrmeir and Lang 2001) . For ease of notation, we ignore the term k and write k,t .

A time-varying coefficient can be conveniently specified as the sum of a fixed (global mean) effect and a temporal random effect of the risk factor: k,t = k + k,t for k = 1, … , K and t = 1, … , T . The fixed effect ( k ) presents the effect of the risk factor that remains constant across space or time, while the temporal random effect ( k,t ) accounts for the time-varying effect of the risk factor (Song et al. 2020) . The temporal random effect k,t can be conveniently specified as a random walk model

f k z k,t k,t = k k,t z k,t for every i and g, and for k = 1, … , K and t = 1, … , T. of order one (RW1) or two (RW2) 9 (Bernardinelli et al. 1995; Martinez-Bello et al. 2017b; Schrödle and Held 2011): with u k,t ∼ N 0, 2 k white noise, and 2 k denoting the variance of the RW process controlling the smoothness of k,t . A random walk process of order one needs an initial value of k,1 and a random walk of order two needs initial values of k,1 and k,2 .

The components i and i are the spatially structured and unstructured main effects at area level, respectively, whereas g(A i ) and g(A i ) denote the spatially structured and unstructured main effects for cell s g in the area A i to which it belongs. The components t and t are the temporally structured and unstructured main effects. For a given time t, they are equal for all areas and cells. it represents the spatiotemporal interaction effect of the unobserved risk factors at area level and g(A i )t the impact of the interaction effect it in cell s g in area A i to which it belongs (see Sect. 3.1 for details). We consider four type of space-time interactions (see Table 5 in Appendix 2).

The final component, Φ gt , in Eq. (3b) is the spatiotemporal GF in cell g for g = 1, … , n p at time t , indicating the true but unobserved relative risk (Cameletti et al. 2013; Godana et al. 2019) . Hence, Φ gt is the "own" spatiotemporal interaction effect of cell g . Because of the large number of cells, it is continuously indexed (Blangiardo and Cameletti 2015) . The component Φ it in Eq. (3a) denotes the area average of Φ gt across the cells within A i . Following Cameletti et al. (2013) and Godana et al. (2019) , we assume that Φ s g , t changes over time following a firstorder autoregressive (AR1) process with coefficient 1 , | | 1 | | < 1:

and s g , t defined as a mean square differentiable process 10 (Stein 1999) with the temporally independent but spatially correlated innovations following a zero-mean Gaussian distribution with spatiotemporal covariance function:

i.e., Var s g , t = 2 Φ for every s g and t , and R(d) the spatial autocorrelation matrix as a function of the distance d between s g and s h at time t (e.g., the Euclidean distance). Under the assumption that the covariance function only depends on d , it is a Matérn covariance

Random walk models are the Bayesian equivalents of P(enalized) Spline regression models (Wang et al. 2018a ). 10 A GF s g , t = t s g is mean square differentiable in A , if for every s g in A and for every t, function satisfying the second-order stationarity and isotropy assumptions. Consequently, the mean of the process is constant and only depends on the locations of s g and s h through the Euclidean distance d = ‖s g − s h ‖ ∈ ℝ (Song et al. 2008 ). The spatial autocorrelation function R(d) is defined as:

where Γ(.) is the gamma function, K v (.) the modified Bessel function of the second order (Abramovitz and Stegun 1965) and v > 0 the parameter controlling the smoothness of the GF (smoothness parameter). In applications, v is commonly fixed because it is usually poorly identified (Miller et al. 2019; Utazi et al. 2019) . In several software packages, including R-INLA (Integrated Nested Laplace Approximation), the default value is one(v = 1), corresponding to moderate smoothness (Lindgren et al. 2011; Utazi et al. 2019) . The scale parameter , > 0 controls the rate of decay of the correlation and is inversely related to the range parameter r of the Euclidean distance between (s g , t) and s h , t . For large r , goes to zero. Because of a lack of a simple relationship between and r , Lindgren et al. (2011) proposed the empirically derived relationship r = √ 8v∕ for spatial autocorrelation near 0.1. Substituting Eq. (10) in Eq. (9), the spatiotemporal Matérn covariance function Σ(d) for each time t is:

For the joint latent spatiotemporal GF t = Φ 1t , … , Φ n p t � at cell level, we have: with t = ( 1,t , … , n p ,t ) � and = 2 Φ R a Matérn covariance matrix. That is, the joint latent spatiotemporal GF is a second-order stationary, isotropic GF with Matérn covariance function Eq. (11) and initial value distributed as 1 ∼ N 0,

This section consists of two subsections. In the first, we present the Fusion Area-Cell Spatiotemporal Generalized Geoadditive-Gaussian Field (FGG-GF) model which integrates the sub-models (3a) and (3b) into a single statistical model. The section also presents the Bayesian statistical tools. In the second section, we discuss solving the "big n" problem resulting in the Fusion Area-Cell Spatiotemporal Generalized Geoadditive-Gaussian Markov Random Field (FGG-GMRF) model and point out that it can be estimated using the R-INLA package. Details on the link between the GF and the GMRF through the Linear Stochastic Partial Differential Equation (LSPDE) approach are discussed in Appendix 1.

(10)

(12) t = 1 t−1 + t and t ∼ N(0, ) for t = 2, … , T.

As observed above, combining low-resolution and high-resolution data to generate high-resolution predictions entails the risk of misalignment (Moraga et al. 2017; Utazi et al. 2019) . To handle misalignment, we first stack the corresponding objects of the area and cell models in Eqs. (3a) and (3b) to give the FGG-GF model 11 which is then estimated as a single model (Blangiardo and Cameletti 2015; Kifle et al. 2017; Utazi et al. 2019) . The FGG-GF model reads:

the joint k th risk factor with fixed coefficient k and temporal random coefficient

t and t defined as in Eq.

(3), and

and t = (Φ 1t , … , Φ gt , … , Φ n p t ) ′ . Note that the above vectors are ((n A + n p ) × 1) for t = 1, … , T.

The following observations apply. First, the basic components of a high-resolution spatiotemporal relative risk model are the covariates and/or the GF at cell level Φ s g , t . Either one or both are required for the estimation of the relative risk at cell level. Second, the interaction terms it and Φ gt in the non-separable models in Eqs. (3a) and (3b), respectively, have as covariance matrices the Kronecker products of the spatial and temporal covariance matrices (Blangiardo and Cameletti 2015; Fuentes et al. 2008) . See Table 5 in Appendix 2 and Sect. 3.2 for further details. For alternative approaches to handling non-separable models, see among others Bakka et al. (2020) , Gneiting (2002) , and Sherman (2011) . Third, the parameters i , i , it are estimated at area level. For the cell level, they are the corresponding area level parameters, implying that they do not vary among cells within a given area A i . Fourth, to control for misalignment, for each area A i , the model component Φ it and the area values of the risk factors x k,it are taken as the block averages of the cells within A i for a given time point t , respectively (Banerjee et al. 2015) . That is, for i = 1, … , n A and

1ds denotes the size of A i . The simplest procedure to estimate x k,it is 

, n A , and t = 1, … , T , with n A i denoting the number of cells in A i (Lawson et al. 2012; Utazi et al. 2019) .

Bayesian estimation of the FGG-GF model is initiated by defining the estimated parameter and hyperparameter vectors. Let = 0 , 1 , … , K , 1 , … , K , , , , , ,

ω , 2 , 2 , 2 , 2 δ , 2 Φ , , 1 , 2 , r � denote the parameter and hyperparameter vectors, respectively, of the FGG-GF model in Eq. (13). The joint posterior distribution of the FGG-GF model is:

where p(.) denotes the probability density function. Below, we first discuss the likelihood function p( | , ) and next the joint prior of the GF at cell level. Based on the assumption that follows a Poisson distribution at area and cell levels (see Eq. (3a) and (3b)), the likelihood function p( | , ) is given by:

The joint prior of the GF at cell level is obtained as follows. Since the GF in Eq. (8) at cell level is assumed to follow an AR1 model, the joint prior distribution of = ( 1 , … , T ) � , i.e., p | 1 , , is (Godana et al. 2019 ):

Because of the AR1 process, we have:

Thus, the joint distribution of the latent spatiotemporal Gaussian process is:

The joint distribution of the GF in Eq. (18) consists of two probability distributions: p 1 | 1 , and p t | t−1 , 1 , for t = 2, … , T . To obtain the joint distribution of , we need the joint distributions of p 1 | 1 , for g = 1, … , n p and t = 1 and p t | t−1 , 1 , for g = 1, … , n p and t = 2, … ., T . 1 is an AR1 stationary process, i.e., 1 | 1 , ∼ N 0, 1− 2 1 . It is called the initial distribution for g = 1, … , n p and reads as:

The joint distribution p t | t−1 , 1 , for g = 1, … , n p and t = 2, … , T is given by:

Finally, the joint prior distribution for the AR1 process, denoted as p | −1 , is given by multiplying Eqs. (20) and (21). It reads:

with −1 = denoting the covariance matrix of the GF in Eq. (12). The priors and joint priors of = 0 , 1 , … , K , 1 , … , K , , , , , , and the hyperpriors of Table 5 in Appendix 2. 12 Details can be found in Jaya and Folmer (2021a) and the references therein. The prior distributions are assumed to be independent, implying that:

The joint hyperparameter distribution is given by 13 :

Given the likelihood function, the joint prior distributions for the parameter vectors and the joint hyperparameter, the joint posterior distribution in Eq. (14) can be written as:

12 A joint prior may refer to multiple parameters, to the spatial units (area and cell), to multiple time periods, or to space-time. The specification of the parameter vector or the subscripts of the precision matrices indicate what kind of jointness is in order. 13 The total number of hyperparameters is J = 2K + 11 , where K is the total number of risk factors and 11 represents the total of the elements of { 2 0 , 2 ω , 2 , 2 , 2 , 2 Φ , 2 δ , , 1 , 2 , r}.

A continuously indexed GF typically has a dense covariance matrix, such as in Eq. (24), leading to complex, time-consuming numerical estimation challenges, commonly referred to as the "big n" problem. Lindgren et al. (2011) proposed to solve the big n problem by substituting a sparse, discretely indexed Gaussian Markov Random Field (GMRF) for the continuously indexed GF. 14 For a GMRF, the full conditional distribution for each component γ g,t for g = 1, … , n p , only depends on a set of neighbors N(g) as follows:

where −g,t denotes all the elements in except γ g,t , and N(g),t denotes all the elements of t in the neighborhood N(g ) of γ g,t . The vector of elements of t not in the neighborhood N(g) of γ g,t is denoted as −{g,N(g)},t . γ g,t is conditionally independent of the elements of −{g,N(g)},t . The conditional independence relationship is written as:

If Eq. (26) holds, the precision matrix = −1 of t is sparse for each t . In other words, for a pair g and h with g ≠ h and h ∉ {N(g)} , we have:

implying that the nonzero pattern in the precision matrix is given by the neighborhood structure. Conversely,

Lindgren et al. (2011) proposed the Linear Stochastic Partial Differential Equation (LSPDE) 15 approach based on a mesh of the study area, to transform a dense Matérn covariance matrix of a GF, such as in Eq. (11), into a sparse Matérn precision matrix of a GMRF (see Appendix 1). Specifically, for t = 1, … , T, the GF γ g,t in Eq. (12) with Matérn covariance function

with the initial value distributed as:

, the Kronecker product of the autoregressive temporal covariance matrix ( T ) (see Table 5 in Appendix 2) and the Matérn spatial covariance matrix ( (13)) as a GMRF ∼ , we introduce the n A + n p × L -dimensional partitioned or block matrix = 1 2 −1 that maps the GMFRs associated with the L triangulation nodes 16 to the n A areas and n p cells, respectively. The elements of 1

and zero otherwise and V i is the number of vertices in the area A i . Hence, matrix H 1 reads:

The LSPDE approach uses finite basis functions defined by a triangulation of the study region (see Online Resource 1 for application to Bandung city) to transform a GF to a GMRF. Implementation can be conveniently handled using INLA (Rue et al. 2009 ). The LSPDE approach is immediately applicable to a large class of spatiotemporal models. For instance, it applies to models with complex hierarchical structures or with non-separable covariance functions, as well as non-stationary models with time-varying parameters (Cameletti et al. 2013) . 16 The number of nodes (vertices) L is larger than the number of cells n p because we extend the domain of interest by an outer area to avoid the boundary effect (see Appendix 1 for details). 

subdistricts. While the number of dengue incidences in Bandung city is reported at district level, for efficient and effective prevention and control, figures at the subdistrict scale are needed. In Sect. 4.1, we discuss and explore the data; in Sect. 4.2, we estimate the FGG-GMRF model; and in Sect. 4.3, we use the model to predict the relative dengue risk at subdistrict level.

The data were obtained from existing databases. Annual observations at district level on the population at risk (see Online Resource 1) and monthly dengue incidence at district level (see Fig. 1 ) were obtained from the Bandung Central Statistical Bureau (2012, 2013, 2014, 2015, 2016, 2017, 2018) and the Bandung Health Department (2013 Department ( , 2014 Department ( , 2015 Department ( , 2016 Department ( , 2017 Department ( , 2018 Department ( , 2019 , respectively. From January 1, 2012, until December 31, 2018, a total of 26,095 dengue incidences (1,030 per 100,000 inhabitants) were reported. The monthly incidence pattern is highly similar from year to year and is taken as constant. Particularly, the mean Pearson correlation coefficient of monthly dengue incidence for the years 2012-2018 is approximately 0.70. In addition, there were no major shifts across the years in the annual cycle. Hence, as in Jaya and Folmer (2021a), we only consider the monthly cycle for the 30 districts. 17 Figure 1 shows a high number of incidences from January to July, followed by a sharp drop in July and a low level of incidences for the remainder of the year. The monthly incidences range from 5 to 265 (0.197-10.461 per 100,000 inhabitants, respectively). We derived the crude risk rate (i.e., the standardized incidence ratio, SIR) as the ratio of the observed to the expected number of incidences (see Eq. (2)). Figure 2 presents the monthly dengue SIR per district for the period 2012-2018. It ranges from 0.229 to 3.132. Most districts, primarily those in northern and southern Bandung, have a SIR greater than one from January to July. The districts with the highest SIR are Buah Batu ( id = 20 ), Lengkong ( id = 27 ) and Rancasari ( id = 29).

As observed by Ebi and Nealon (2016) , Jaya and Folmer (2021a) , and Zellweger et al. (2017), among others, socioeconomic and environmental conditions are the main factors influencing dengue disease risk over space and time. However, for Bandung, socioeconomic risk variables such as income, education, occupation and living conditions are unavailable for districts and cells. These factors are accounted for by the random effects, thus controlling for omitted variable bias Folmer 2020, 2021a, b) . By contrast, the monthly averages of the weather risk variables of precipitation (mm), temperature (°C), sunshine duration (kJ/m 2 day) and water vapor pressure (kPa) are available at cell level from the WorldClim2.0 database (Fick and Hijmans 2017) , obtained from 19 weather stations surrounding Bandung in West Java for 1970-2000. We selected cells of resolution 1 km 2 . Accordingly, Bandung city was divided into 179 cells. Table 1 presents the monthly average weather variables for the period 1970-2000. Figure 3 shows that precipitation and water vapor pressure are relatively high in the period November-April, temperature is relatively high in the period April-June and in October, and solar radiation is relatively high in the period August-November. Figure 4 indicates that the average temperature varies strongly over space. The minimum average temperature occurs in the northern districts, which are mountainous areas at approximately 800 m above sea level. In addition, they are densely covered with forests and have relatively high precipitation. The central districts, where the governmental facilities and businesses are located, also have high precipitation in the period November to April. Moreover, they have higher temperatures than northern Bandung because of differences in forest density and elevation. They also have high population density, high mobility and high air pollution (Jaya and Folmer 2020) .

In preparation for estimating the FGG-GMRF model in Eq. (33), we calculated the variance inflation factor (VIF) of the weather variables to check multicollinearity. Table 2 shows that the maximum VIF is 9.312 (for water vapor pressure) which is below the critical (rule of thumb) value of 10, indicating that the correlation among variables is unlikely to affect estimation (Montgomery et al. 2012 ).

The first step in estimating the FGG-GMRF model given by Eq. (33) is the construction of a triangle mesh of the study area for the application of the Finite Element Method (FEM) and LSPDE approach. As described in Appendix 1, the accuracy of the FEM calculations and the precision of the forecasts, is a function of the number of vertices in the mesh (edge length). Blangiardo and Cameletti (2015) and Utazi et al. (2019) recommended varying the edge length between the minimum distance and approximately 5-8% of the maximum distance of any two cells (18,681 m). Hence, we considered G = {966, 831, 717, 616, 548, 501} vertices, corresponding to edge lengths varying from 1000 to 1500 m, with a difference 18 of 100 m (see Online Resource 1). The data and R code are available in Online Resource 2. Before turning to the estimations of the spatiotemporal FGG-GMRF model, we make the following observations. First, as explained in Sect. 3, the covariates and/ or the state process at cell level are required for high-resolution spatiotemporal prediction using the FGG-GMRF model. Hence, either the covariates or , or both, are included in the selected model. Second, for every model, we considered Poisson and Negative Binomial model specifications for the number of incidences, a random walk of order one (RW1) and two (RW2) for the time-varying coefficients, structured and unstructured spatial and temporal random effects and their interaction, and six different edge lengths. Third, the best model was selected using the deviance information criterion (DIC), the Watanabe-Akaike information criterion (WAIC) and the marginal predictive likelihood (MPL). As a rule of thumb, the best model is the one with the smallest DIC and WAIC, and the largest MPL. Fourth, we started the estimation with the simplest models with covariates only (M1), and then, we proceeded to the model specifications with covariates and four types of interaction (see Table 5 in Appendix 2) at area and cell levels (M2) and, finally, we estimated the full models with covariates, four types of interaction at area and cell levels, and spatially and temporally structured and unstructured main effects (M3). 19 Finally, Fig. 3 Monthly variation of the mean annual weather variables a precipitation (mm), b temperature (°C), c solar radiation in 1000 (kJ/m 2 day), and d water vapor pressure (kPa) due to the large number of outcomes, we only present the estimates for interaction Type IV (spatially structured ⊗ temporally structured) which, as in Jaya and Folmer (2020) , performed best among the models with spatiotemporal interaction. 20 The estimations are presented in Table 6 (in Appendix 3 ). The table shows that the models 21 with covariates only (M1) have the worst fit and predictive performance among the three classes of models. They have the highest DICs, WAICs and smallest MPLs. Given their relatively poor fit and predictive performance, the M1 models were not considered further in the selection procedure. Introduction of the interaction effect type IV (M2) yielded substantially better predictive performance. The full models with covariates, interaction at area and cell levels and spatially and temporally structured and unstructured main effects (M3) had fit and predictive performance similar to the M2 models. Hence, the main spatially and temporally structured and unstructured effects did not improve the model fit and prediction performance, which is consistent with Jaya and Folmer (2020) . Based on these observations, and because the M2 models have a simpler structure, we selected the class of M2 models.

Next, we turned to the selection of the best model from the 24 models M2. Figure 5a shows that for model M2.1.1.6, the observed and predicted dengue relative risks are strongly correlated (Pearson correlation coefficient = 0.986) , indicating that the model fits the data well. Figure 5b shows that the PIT histogram is close to the Uniform distribution, also indicating that model M2.1.1.6 fits the data well. Table 3 summarizes various components of model M2.1.1.6 which are subsequently used to calculate the posterior means of the monthly relative risk at district and subdistrict levels. Before doing so, we discuss the components separately. To this end, we also make use of Figs. 6 and 7. Before going into detail, we make the following remarks. First, as observed in Sect. 2, p.8, the posterior means of the time-varying coefficients k,t = k + k,t for k = 1, … , K and t = 1, … , T , which are presented in Fig. 6 , consist of the fixed effect plus the temporal random effect. The minimum and maximum posterior means of the temporal random effects present the largest negative and largest positive differences of the temporal random effects relative to the global effects. Secondly, the contributions of the weather variables in explaining the spatiotemporal variation of the dengue risk are conveniently summarized by the posterior means of their hyperparameter variances and their percentage contributions as fractions of the total variance (the last column of Table 3 ). In a similar vein, the posterior means and fractions of the hyperparameter variance of the random components present their variability and strength in explaining the relative risk across space and time.

We will begin the discussion of Table 3 with the global mean effects of the weather variables. We only consider the posterior means and disregard credible intervals. Before going into the detail, it is worth noting that an indirect relationship exists between the relative risk of dengue and weather variables via the development and survival of the dengue virus and its vector (Jaya and Folmer 2021a) .

Precipitation in general has a negative impact, with a posterior global (overall) mean of − 0.0041, which is consistent with Jaya and Folmer (2021a) . For a one mm increase in the global mean, the dengue risk decreases by (exp(−0.0041) − 1)100% = −0.41% . The explanation is that heavy rainfall disrupts the Aedes-spp mosquito's reproductive cycle by washing away breeding sites (Abiodun et al. 2016; Benedum et al. 2018) .

Temperature in general has a positive impact (0.1002), which is consistent with Hurtado-Díaz et al. (2007) and Jaya and Folmer (2021a) . The relative risk increases by (exp(0.1002) − 1)100% = 10.54% for an increase of global mean temperature by 1 0 C. The explanation is that higher temperatures offer good conditions for mosquito development, particularly feeding (Hales et al. 2002; Lambrechts et al. 2012) .

Solar radiation has a negative impact (− 0.0001), which is consistent with Ekasari et al. (2018) , Jaya and Folmer (2021a) and Martínez-Bello et al. (2017b) . An increase of 1 kJ m −2 day −1 of solar radiation decreases the relative dengue risk by (exp(−0.0001) − 1)100% = −0.01% . As shown by Rasjid et al. (2019) , strong solar radiation negatively influences the breeding and spread of Aedes-spp mosquitoes. A longer spell of solar radiation implies a shortened spell of dawn and dusk, during which the Aedes-spp mosquito preys on animals and humans, particularly 20 to 30 min after sunset (Ekasari et al. 2018; Jaya and Folmer 2021a) .

Water vapor pressure in general has a negative impact, with a posterior global mean of 0.5033. An increase of the global mean by 1% increases the relative dengue risk by (exp(0.5033) − 1)100% = 65.42% due to an increase in breeding ability (Bambrick et al. 2009 ).

For the temporal random effects of the weather variables, the greatest negative temporal random effect of precipitation occurred in June (-0.0248), while the greatest positive occurred in March (0.0176). For temperature, the largest negative temporal random effect was in February (-0.0713), while the largest positive was in June (0.1104). Solar radiation has a temporal random effect close to zero, as indicated by its minimum and maximum posterior means of -0.0002 and 0.0003, respectively. The greatest negative of temporal random effect of water vapor pressure was in January (− 0.0137), while the greatest positive effect was in June (0.0161).

The variances of the weather variables together account for 31.87% of the total variance of the hyperpriors with solar radiation being the most important ( 2 3 = 0.023) , explaining 9.88%, while temperature is the smallest, accounting for 6.2%. The variance of the area level interaction effect ( 2 δ ) ) accounts for the highest fraction of the total variance (43.47%), followed by the cell level interaction effect ( 2 Φ ), with a fraction of 24.68%. These fractions indicate that the trend of dengue relative risk for each district and subdistrict is strongly affected by neighboring districts and subdistricts, respectively. The relatively low fractions of the total variance for the other effects imply that they are less important. Specifically, the low fraction of the total variance for the average temperature indicates that only a small part of the variability of the relative risk of dengue in districts and subdistricts is explained by the average temperature.

The posterior mean of hypermeter of the Leroux CAR spatial autoregressive coefficient ( ) , and the posterior means of the hyperparameters of the temporal autoregressive coefficients at area level ( 2 ) and cell level ( 1 ) are substantial (larger than 0.700), indicating strong spatial and temporal dependency. The estimated range is r equals 13.597 km . This implies that beyond the distance r = 13.597 the spatial correlation among any two cells is smaller than 0.1. Hence, the observations are spatially strongly correlated. Beyond 13.597km it is negligible.

Based on Table 3 , we now discuss the estimated parameters for the time-varying effects of the risk factors and, next, the spatiotemporal interaction effects at the area and cell levels.

The posterior means of time-varying effects of the weather variables are presented in Fig. 6 . The figure shows that the time-varying effects of precipitation, average temperature and solar radiation vary considerably over the year. The timevarying effect of water vapor pressure, in contrast, is highly constant over time.

The time-varying effect of precipitation is positive for the periods January-April, August-September and December and negative for May-July and October-November. The strongest negative effect was in June (− 0.0289), and the strongest positive effect was in March (0.0135). The negative impact for the period May-July follows after the peak of the rainy season from November-April. The negative impact for October-November is caused by the increase in precipitation after the peak of the dry season in June-July. The positive impacts for the periods January-April and August-September, and the peak in March, correspond to the relatively low rainfall one to two months before these periods.

The time-varying effect of temperature is positive for all months and ranges from 0.0289 (end of August) to 0.2106 (June-July). It increases from January-June, is at its peak in June, decreases from July until mid-September, and then starts increasing up to the global mean, where it remains for the rest of the year. Note that the monthly temperature has a delayed risk effect in that it increases from January-May, Table 3 The posterior means and 95% credible intervals of the fixed effects, the minimum and maximum posterior means and 95% credible intervals of the temporal random effects and the posterior means and 95% credible intervals of the hyperparameters of the FGG-GMRF model M2. while its impact is largest in June. The delay is due to the mosquito life cycle and incubation period (Jaya and Folmer 2021a). The time-varying coefficient of solar radiation is below zero for almost all months, except for November, and varies from − 0.0003 to 0.0002. This is due to the fact that tropical countries such as Indonesia receive a lot of solar radiation throughout the year (Handayani and Ariyanti 2012) . The time-varying effect of water vapor pressure is positive all year round and hardly varies. The effect varies from 0.4896 to 0.5194. The strongest effect is in June (0.5194). Figure 7 presents the estimated parameters of the spatiotemporal interaction effects at area and cell level which depend on their hyperparameters in Table 3 . Figure 7a presents the district level spatiotemporal interaction effects (i.e., the residual effect after accounting for the weather effects). The figure shows that the interaction effect varies across districts and time. In the northern districts, it is positive and quite high during January and February, followed by a decrease to around zero during the period March-May. From June-September it is moderately positive followed by a period of high positive interaction for the rest of the year, especially in the most north-western districts. The high interaction effect in the period September-February is related to multiple factors, in particular environmental conditions. The northern areas are ideal breeding habitats because of dense vegetation and high humidity, especially during the rainy season, with low sunshine duration and high humidity.

The central districts have high spatiotemporal interaction effects because of favorable socioeconomic conditions for the spread of the dengue virus, including high population density and high density of hotels, hostels and student apartments. For the northern central districts, there is the additional effect of spillover of mosquitoes from the northern districts. As a consequence, the interaction effect of the northern central districts follows a time pattern similar to time patterns of the northern districts, though less intense. The most central districts have high positive interaction effects all year round because they have the highest population density and density of hotels and hostels. The two most central districts have low interaction effects all year round, indicating that the main effects (risk factors) virtually fully explain the dengue risk. These districts have no special socioeconomic or environmental conditions affecting the dengue incidence rate.

In the southern districts, the spatiotemporal interaction effect is similar to that in the most central districts, although for partly different reasons. They have high population density with many residential areas that have unhygienic conditions. See Hsu et al. (2017) for details on the relationship between hygiene and dengue infection.

The situation in the eastern districts differs from that in the northern, southern and central districts. The interaction effect is negative in January and February, highly positive in March-May, slightly positive in June-August and negative for the rest of the year. The negative interaction effect in January-February is probably caused by the interaction of the weather variables and the environmental conditions. The absence of forests, heavy rainfall, and the short spells of sunshine in the period January-February keep the humidity low, which is unfavorable for the presence of dengue mosquitoes. The districts are residential areas with inadequate drainage and sanitation. The heavy rainfall until March combined with inadequate drainage and sanitation leads to large quantities of standing water, which provides favorable breeding habitats, contributing to the positive interaction effect in March.

The western districts have medium to strong negative interaction effects all year round, reducing the effects of the weather variables. The majority of the western districts have good drainage and sanitation, and the lifestyle and health behavior of the population is substantially better than in the other parts of Bandung, reducing dengue infection (Bandung Health Profile 2019). For example, the district with the highest healthy behavior index, Cicendo, is located in the western region. The negative effects also indicate that there is limited spillover of mosquitoes from the other districts. Figure 7b presents contour maps of the cell level interaction effects. In contrast with Fig. 7a , the cell level interaction effects are almost the same across the months. Hence, after accounting for the weather variables, the cell level residual varies over space but is relatively constant over time, implying that it is related more to a topographical dimension, such as elevation, than to time. Positive cell level interaction effects are found in the northern part of Bandung, which is at 800 m above sea level and has high precipitation and dense vegetation, providing an ideal breeding ground and habitat for the Aedes-spp mosquito (Arboleda et al. 2009 ). Figure 8a shows the posterior means of the relative risk estimates (based on the posterior means of the time-varying coefficients of the risk factors and the posterior means of the spatiotemporal interaction effects) at district level. The posterior means of the relative risk at subdistrict level (see Fig. 8b ) are obtained as the block average of the cell values within each subdistrict boundary and the cell values that are partly outside its boundary. 22 Accordingly, Bandung city is divided into 30 districts, 151 subdistricts and 179 cells.

Comparing Figs. 8a and 8b shows similar temporal trends, in particular, an increase in January-July and a decrease in August-November. This applies especially to the spatial units with the highest risk, notably the districts in southern Bandung. For the spatial dimension, however, we notice substantial differences. For example, according to Fig. 8a , the entire central district of Lengkong (surrounded area in Fig. 8 ) is categorized as a high-risk area in the period January-July, whereas Fig. 8b shows that this only applies to parts of the district. The explanation is that categorization based on the model given by Eq. (3a) ignores within-district heterogeneity, while this is taken into account when categorization is based on the model given by Eq. (3b). To explore this issue further, we calculated the high-risk and low-risk districts and subdistricts based on the posterior exceedance probability for the two approaches, denoted as the top-down and the bottom-up approaches, respectively. Following Sparks (2015) and Osei and Stein (2017) , we fixed the posterior exceedance probability threshold for it at 1.25 . The exceedance probability � Pr it > 1.25| over space and time is presented in Fig. 9 . Table 4 presents the classification of the subdistricts into high and low risk based on the bottom-up and top-down approaches, respectively. Table 4 shows substantial misclassification for the top-down approach for the period January-July. The overall misclassification rate is 15.6% for all periods (January-December) and 26.7% for the high-risk period (January-July). The table furthermore shows that all the misclassifications occurred in the period January-July, which partly overlaps with the rainy season in November-May, whereas there is no misclassification from August-December. The explanation for the misclassification as such is that the bottom-up approach averages out the differences in risk factors, and consequently the number of incidences, over a set of relatively small number of relatively homogenous cells within a subdistrict. The top-down approach, on the other hand, averages out the differences over a substantially larger number of relatively heterogeneous cells in a district (a district contains at least four subdistricts; see Sect. 4.1).

The misclassification is obviously concentrated in the rainy period January-July with local variation in the risk factors due variation in local characteristics such as elevation or vegetation density. Although the rainy season is from November to May, there are positive rates of misclassification in June-July and unexpected zero rates in November-December. These misclassifications and unexpected rates are due to the 22 The reason that some cells are partly outside some subdistrict boundaries is that the cells are squares. delayed responses of mating, breeding and hunting by Aedes-spp mosquitos. Mating and breeding occur mainly during the rainy season. Following this, it takes approximately two weeks to one month for the eggs to develop into adult mosquitoes and for Fig. 9 Monthly posterior exceedance probability of the relative risk at a district and b subdistrict levels (surrounded area: Lengkong district) the virus to multiply and reach the salivary glands before it is transmitted to humans. If an individual is infected, the symptoms can be observed approximately four to seven days after being bitten (Ehelepola et al. 2015) . Accordingly, there is a delayed infection response with respect to the weather conditions (Jaya and Folmer 2021a).

Effective and efficient control of a variety of spatial problems, including dengue disease abatement, requires data at a fine spatiotemporal scale. However, data availability at the same (especially fine) spatial scale is quite rare (Moraga et al. 2017; Utazi et al. 2019) . A major challenge in spatial sciences, including modeling of infectious diseases such as dengue and COVID-19, is how to align data bases of different resolutions consistently. In this study, we presented the Fusion Area-Cell Spatiotemporal Generalized Geoadditive-Gaussian Markov Random Field model as a solution to this problem. This model combines observations on the dependent variable and population at risk at the area level and covariates at the cell level to generate predictions of relative risk at the subdistrict level. Special attention was paid to the model setup to generate predictions for the target cells during model-fitting, using Bayesian Integrated Nested Laplace Approximation (INLA). The methodology was applied to monthly dengue disease data for 30 districts in the city of Bandung, Indonesia, for the period January 2012 to December 2018. The risk factors consisted of the monthly averages of precipitation, temperature, solar radiation and water vapor pressure. The analysis showed that the effects of precipitation, temperature and solar radiation varied considerably across space and time, while the effect of water vapor pressure was highly constant over time. Solar radiation was found to be the most important risk factor. The spatiotemporal interaction effect, capturing the effects of omitted variables at area level, also varied across districts and time. In contrast, the cell level interaction effect was almost constant over the months but varied substantially over space, indicating a strong spatial spillover effect. Based on the posterior means of the relative risk at cell level, we obtained the relative risk estimates at subdistrict level. We found a similar temporal pattern for district and subdistricts. Relative dengue risk was relatively high in the period January-July and relatively low during the period August-December. We further compared the risk estimates per subdistrict based on: (i) the bottom-up approach using the cell level estimates and (ii) the top-down approach assigning the district value to its subdistricts. Using the posterior exceedance probability of the relative risk, we identified high-risk and low-risk districts to find that during the high-risk period of January-July, the topdown approach misclassified 26.4% of the subdistricts as high risk, which according to the bottom-up approach was low risk. The overall misclassification rate was 15.6%.

The main conclusions of the paper are the following. First, effective and efficient policy intervention, such as the control of infectious diseases, requires data at the right level of resolution. In particular, low-resolution maps may misclassify regions. If regions are incorrectly misclassified as high-risk, unnecessary policy intervention with undue financial, social and environmental costs may result. In contrast, if regions are incorrectly misclassified as low-risk, opportunities for policy intervention may be missed which may also have costs of various kinds. Secondly, the proposed FGG-GMRF model adjusts data and maps of different resolutions consistently, and allows more data to be utilized, thus improving the statistical efficiency. Third, application of the FGG-GMRF model to the dengue disease data for Bandung from 2012 to 2018 shows that the relative infection risk is high in various cells, subdistricts and districts from January-July. The strong spatiotemporal interaction indicates that the occurrence of the dengue disease vector is highly contagious and must be detected early in order to prevent its spread. Rapid response measures such as fogging are critical in areas with high dengue incidence. Finally, based on the experiences in the present paper, it is worthwhile to investigate the suitability of the FGG-GMRF model for a variety of other spatiotemporal problems, including other infectious diseases such as COVID-19 (see Jaya and Folmer 2021b), vaccination coverage (Utazi et al. 2019) , particulate matter concentration (Cameletti et al. 2013; Lee et al. 2016) , and social issues such as unemployment and crime.

A GF with dense covariance matrix can be transformed to a GMRF with sparse covariance matrix by means of a Linear Stochastic Partial Differential Equation (LSPDE) (Lindgren et al. 2011) which reads:

where t (s) is a temporally independent GF, W t (s) a Gaussian white noise process, a positive integer related to the smoothness parameter of the Matérn covariance function in Eq. (11) by = v + 1, and ∇ 2 the Laplace operator, i.e., ∇ 2 = Using spectral decomposition, Whittle (1954 Whittle ( , 1963 showed that for > 1, the only exact stationary solution to Eq. (34) is the isotropic Matérn field, i.e., the stationary GF t (s) with the Matérn covariance function in Eq. (11).

A closed-form solution for the LSPDE in Eq. (34) is restricted to regular lattices (Lindgren et al. 2011) . For an irregular lattice, it can be approximated through a basis function representation using the Finite Element Method (FEM) defined on the domain A ∈ ℝ 2 . The basis function representation is defined on a mesh, that is, a collection of (i) vertices, (ii) the edges between the vertices, and (iii) the polygons described by the edges. A mesh consists of a minimum of three connected edges conforming to the shape of the domain. It subdivides a continuous geometric space into a finite set of discrete geometric or topological elements such as triangles or rectangles (for a two-dimensional geometric space) or tetrahedral or rectangular prisms (for three-dimensional spaces). It reduces the degrees of freedom from infinite to finite. Because the FEM calculations are based on a finite number of cells and the results are generalized through interpolation for the entire domain, the accuracy of the global solution is a function of the number of elements of the mesh (Bohn and Feischl 2021) .

Triangulation is a common FEM meshing scheme (Sloan 1993) because of its flexibility for irregular domains (Lindquist and Gilest 1989) and accuracy (due to the minimization of the discretization error) 23 (Ahmadian et al. 1998 ). Triangulation divides the domain into a set of non-intersecting triangles, where any two triangles meet in, at most, a common edge or vertex. The popular Delaunay triangulation scheme (Cheng et al. 2013) ensures that the triangulation maximizes the minimum angle of the triangles, thus avoiding sliver (i.e., long and thin) triangles and rendering the transitions between small and large triangles smooth (Lindgren et al. 2011 ). The restricted Bowyer-Watson algorithm, which is designed to conform to the domain's boundary, is a popular Delaunay triangulation algorithm (Cheng et al. 2013) .

A drawback of applying an algorithm with boundary conditions is that the variance near the boundary is inflated by a factor of two (Lindgren et al. 2011) . To avoid the boundary effect, Lindgren and Rue (2015) proposed to extend the domain of interest by an outer area at a distance of at least r , corresponding to a correlation of (34) 2 − ∇ 2 ∕2 t (s) = W t (s) for t = 1, … , T and s = (s 1 , … s g , … , s n p ) � with s g = (s g,1 , s g,2 )

approximately 0.1 between two points in the inner and outer areas. The edge length of the triangles of the outer area should be at least equal to the edge length of the triangles of the inner area (Blangiardo and Cameletti 2015) . To find the appropriate mesh for the data at hand, meshes of different sizes are usually considered, which can be evaluated using the DIC and WAIC (Righetto et al. 2018) . Given the triangular mesh, the FEM of the solution of the LSPDE in Eq. (34) (Lindgren et al. 2011 ). To show that Eq. (34) approximates the GF t (s) , we take 24 = 2 corresponding to v = 1, define the LSPDE in Eq. (34) as a variational problem, i.e., multiply it by an arbitrary test function (s) for s = (s 1 , … s g , … , s n p ) � with s g = (s g,1 , s g,2 ) , and integrate it by parts over the domain A using Green's first identity theorem 25 (Langtangen and Logg 2016) . Multiplying the LSPDE in Eq. (34) by a test function (s) gives: (Lindgren et al. 2011) . For = 2 , 2 − ∇ 2 ∕2 t (s) is a linear differential operator which is easier to solve than a fractional differential operator when is an odd number, e.g. = 1 (Miller et al. 2019) . 25 Green's first identity theorem is a multidimensional generalization of integration by parts to reduce the order of derivatives (Langtangen and Logg 2016) . Applied to Eq. (37), it reads: with γ t (s) the directional derivative:

where A is the entire spatial domain over which Eq. (37) is to be solved and ds is shorthand for ds 1 ds 2 of the two-dimensional integral. Partial integration of −∫ A (s)∇ 2 t (s)ds gives: 

The test functions are commonly taken to be equal to the basis functions, i.e., h (s) = l (s) 26 for h = 1, … , L . Hence:

The integral on the right-hand side of Eq. (42) is the Gaussian white noise distribution N 0, ∫ A 2 l (s)ds , with mean zero and covariance matrix 27 :

We can write Eq. (42) in matrix form as (Simpson et al. 2012) :

Because of the highly local nature of the piecewise linear basis functions, s and s are sparse matrices. However, −1 s is a dense matrix and will generally not be a GMRF. Lindgren et al. (2011) proposed 26 This procedure is known as the Galerkin finite element method (Langtangen and Logg 2016) . 27 By definition, a generalized GF, for example, (s) , in a domain A is a random L 2 (A) generalized function, with L 2 (A) a vector space of square-integrable functions in two-dimensional space A (i.e., L 2 (A) = { (s) ∶ A → ∫ A | (s)| 2 ds < ∞} ) such that for every finite set of test functions h ∈ L 2 (A), h = 1, … , H , the inner product ∫ A h (s)ds for h = 1, … , H is jointly Gaussian. Because Gaussian white noise W t (s) is an L 2 (A)-bounded generalized GF in the domain A , the distribution of ∫ A l W t (s)ds is Gaussian, with expectation and covariance given by (Lindgren et al. 2011) :

and with ∼ the precision matrix where κ = √ 8v∕r. Because of the serial independence assumption in Eq. (9), ∼ is constant over time. Bolin and Lindgren (2009) compared the exact FEM approach (using s with all the elements of s calculated as a Hilbert space wavalet model such as a B-spline or a Daubechies wavalet model), and the Markov approximation (replacing s with ∼ ) for spatial prediction and found that the differences between both approaches in terms prediction errors are negligible.

This appendix consists of three parts. The first part is Table 5 . which presents the priors and hyperpriors for the FGG-GMRF model given by Eq. (33). The second part consists of comments on the priors and hyperpriors and the final part deals with identification.

See Table 5 .

The following observations apply for the parameters and hyperparameters of the FGG-GMRF. Due to the lack of strong prior knowledge, we used a vague Gaussian prior distribution with a zero mean and a very large variance for the parameters 0 ∶ 0 ∼ N 0, 10 6 and k ∶ k ∼ N 0, 10 6 for k = 1, … , K(Blangiardo and Cameletti 2015; Martinez-Beneito and Botella-Rocamora 2019). A weakly informative prior 29 was assigned to the log-odds of the spatial autoregressive parameter ∶ log( ∕(1 − )) ∼ N(0, 0.45) (Utazi et al. 2019) . Note that the transformation is

for every t, 

precision matrix of the parameters k given by 

Temporally unstructured random effect Type II: combines the spatially unstructured i and the temporally structured t random effects. Implies ≠ 0 and 2 = 0 An exchangeable (Gaussian) process is a sequence of random variables that are independent and identically (normal) distributed (iid).

c Scaling is applied to the temporal structure matrices of RW1 and RW2 to make their marginal variances comparable. Scaling is implemented in R-INLA by using the option scale.model = TRUE For the scale hyperparameters (variance and standard deviation) of the temporally varying coefficients, spatially and temporally structured and unstructured random effects, spatiotemporal interaction of the random effects, and the Gaussian field (GF), the Inverse Gamma (IG), half-Cauchy (HC), Uniform and Penalized Complexity (PC) are common hyperpriors. The IG distribution is a well-known, easy to use, conjugate hyperprior for the variance (Coly et al. 2021; Hamura et al. 2021) . 30 The HC is a Cauchy distribution truncated at zero. Gelman (2006) recommended the HC with scale parameter 25 as an alternative to the IG when the variance hyperparameter is close to zero. The HC is also a frequently used hyperprior for the standard deviation (Gómez-Rubio 2020). When the limit of the scale parameter goes to infinity, the HC converges to the Uniform hyperprior (Gelman 2006; Gómez-Rubio 2020) . However, for the Uniform hyperprior, it is difficult to set ranges of values for the standard deviation using R-INLA. The reason is that R-INLA assumes that the standard deviation is unbounded (above) (Gelman et al. 2017 ; Gómez-Rubio 2020), whereas a fixed upper bound on the standard deviation is frequently needed to avoid overfitting 31 . As an alternative to the HC and Uniform hyperpriors for the standard deviation, Simpson et al. (2017) proposed the PC hyperprior. The PC allows using probability statements on values for the parameters such as the standard deviation, autocorrelation and range parameters (Franco-Villoria et al. 2019) . The parameters can be restricted using either an upper bound or a lower bound but not both. The PC is frequently used for autoregressive processes (Sørbye and Rue 2017) , varying coefficients models and high-resolution prediction models (Franco-Villoria et al. 2019) .

In the application, we assigned the IG (1, 0.01) hyperprior to the variance of the autoregressive AR1 model of the temporally structured effects ( 2 ) and to the variances of the exchangeable (iid) priors of the spatially ( 2 ) and temporally ( 2 ) unstructured random effects, respectively. Following Lawson et al. (2010) , we assigned the HC hyperprior with scale parameter to the standard deviation of the Leroux prior (Leroux et al. 2000) of the spatially structured random effect 2 .

Regarding the hyperpriors to the variances of the RW1 and RW2 priors the following applies. The marginal variances of the RW1 and RW2 priors reflect the smoothness of the k th vector of the temporal random effect k = k,1 , … , k,T � for k = 1, … , K. 32

However, before assigning hyperpriors we need to render the marginal variances as smoothness indicators comparable as they are affected by their temporal structures and the scale of the risk factors which render them inadequate as smoothness indicators (Sørbye and Rue 2014; Blangiardo and Cameletti 2015; Gómez-Rubio 2020; Kang et al. 2015) . The temporal structure matrices carry over to the hyperpriors. To overcome the incomparability problem, Sørbye and Rue (2014) proposed to scale the 30 Note that in the present setting we are dealing with the latent Gaussian model (LGM) for which the IG distribution is a conjugate hyperprior for the variance of the Gaussian prior (Hamura et al. 2021) . 31 The Uniform distribution becomes an increasingly uninformative prior as the upper bound increases (Lemoine 2019). Hence, it is unclear what upper bound is appropriate for a Uniform prior. 32 The smaller the marginal variance, the smoother the time-varying regression coefficient.

between the grid cells. For the prior belief we took r = 0.5 , i.e., Pr(r < 5) = 0.5 . Note that the two preceding decisions imply that we know the true limits ahead of the analysis.

Because the ICAR, Leroux, AR1, RW1, and RW2 priors are specified conditionally on neighboring observations and their parameters are unique up to an additive constant, their specifications implicitly include the overall intercept. Consequently, if the model would also include an additional intercept, it would not be identified due to the perfect collinearity (Goicoa et al. 2018) . To solve this identification problem, the explicit intercept can be excluded or sum-to-zero constraints on the random effects can be imposed (Goicoa et al. 2018) . When sum-to-zero constraints are imposed, the random effects become orthogonal to the explicit overall intercept. To achieve identifiability, we imposed the following sum-to-zero constraints on the temporal random of effects of the risk factors, the spatial and temporal random effects and their spatiotemporal interaction effects components:

The above-mentioned sum-to-zero constraints are imposed in INLA by setting the constr = TRUE argument to the function f (.) which defines the temporal random effects of the risk factors, the spatial, temporal, and spatiotemporal interaction effect components.

The best specification of the FGG-GMRF model in Eq. (33) was selected using the deviance information criterion (DIC), the Watanabe-Akaike information criterion (WAIC), and the marginal predictive likelihood (MPL). The results are presented in Table 6 . 35 As a rule of thumb, the best model is the one with the smallest DIC and WAIC, and the largest MPL. 

Disease mapping and spatio-temporal analysis: importance of expected-case computation criteria

Modelling the influence of temperature and rainfall on the population dynamics of Anopheles arabiensis

Identifying hotspots of cardiometabolic outcomes based on a Bayesian approach: the example of Chile

Minimization of the discretization error in mass and stiffness formulation by an inverse method

Spatiotemporal prediction of infectious diseases using structured Gaussian processes with application to Crimean-Congo hemorrhagic fever

Mapping environmental dimensions of dengue fever transmission risk in the Aburrá Valley Colombia

How to solve the stochastic partial differential equation that gives a Matérn random field using the finite element method

The diffusion-based extension of the Matérn field to space-time

Climate change could threaten blood supply by altering the distribution of vector-borne disease: An Australian case-study

Bandung Government, Bandung Bandung Health Department (2013) Health Profile of Bandung Municipality in 2012. Bandung Government, Bandung Bandung Health Department (2014) Health Profile of Bandung Municipality in 2013. Bandung Government, Bandung Bandung Health Department (2015) Health Profile of Bandung Municipality in 2014. Bandung Government, Bandung Bandung Health Department (2016) Health Profile of Bandung Municipality in 2015. Bandung Government, Bandung Bandung Health Department (2017) Health Profile of Bandung Municipality in 2016. Bandung Government

Prediction interpolation and regression for spatially misaligned data

Modelling the presence of disease under spatial misalignment using Bayesian latent Gaussian models

Statistical modeling of the effect of rainfall flushing on dengue transmission in Singapore

Overdispersion and Poisson regression

Bayesian analysis of space-time variation in disease risk

Spatial data analysis with R-INLA with some extensions

Spatial and spatio-temporal Bayesian models with R-INLA

Recurrent neural networks as optimal mesh refinement strategies

Wavelet Markov models as efficient alternatives to tapering and convolution fields

Spatio-temporal modeling of particulate matter concentration through the SPDE approach

Analysis and comparison of Bayesian methods for measurement uncertainty evaluation

Bayesian hierarchical models for disease mapping applied to contagious pathologies

Dengue in a changing climate

A study of the correlation between dengue and weather in Kandy city Sri Lanka (2003-2012) and lessons learned

Climate factors and dengue fever in Jakarta

Bayesian inference for generalized additive mixed models based on Markov random field prior

WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas

A unified view on Bayesian varying coefficient models

Generalized additive models for cancer mapping with incomplete covariates

Model evaluation and spatial interpolation by Bayesian combination of observations with outputs from numerical models

A class of nonseparable and nonstationary spatial temporal covariance functions

Intuitive joint priors for variance parameters

Prior distribution for variance parameters in hierarchical models

The prior can often only be understood in the context of the likelihood

Nonseparable, stationary covariance functions for space-time data

Dynamic spatiotemporal modeling of the infected rate of visceral leishmaniasis in humans in an endemic area of Amhara regional state Ethiopia

In spatio-temporal disease mapping models, identifiability constraints affect PQL and INLA results

Spatial models using Laplace Approximation methods

Modelling spatial and spatial-temporal data a Bayesian approach

Potential effect of population and climate changes on global distribution of dengue fever: an empirical model

Time-varying coefficient models: a proposal for selecting the coefficient drivers sets

On global-local shrinkage priors for count data

Potency of solar energy applications in Indonesia

A statistical downscaling approach for generating high spatial resolution health risk maps: A case study of road noise and ischemic heart disease mortality in Melbourne

Generalized additive models

Trend and geographic analysis of the prevalence of dengue in Taiwan 2010-2015

Short communication: impact of climate variability on the incidence of dengue in Mexico

Modeling of infectious diseases: a core research topic for the next hundred years

Bayesian spatiotemporal mapping of relative dengue disease risk in Bandung Indonesia

Identifying spatiotemporal clusters by means of agglomerative hierarchical clustering and Bayesian regression analysis with spatiotemporally varying coefficients: methodology and application to dengue disease in Bandung Indonesia

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

Geoadditive models

Disease prevention: saving lives or reducing health care costs?

A case study for modelling cancer incidence using Bayesian spatio-temporal models

Cross-covariance functions for additive and coupled joint spatiotemporal SPDE models in R-INLA

Bayesian modeling of inseparable space-time variation in disease risk

Impact of daily temperature fluctuations on dengue virus transmission by Aedes aegypti

Solving PDEs in Python. The FEniCS tutorial I. Springer Open, Cham Langtangen HP, Mardal KA (2019) Introduction to numerical methods for variational problems

Hotspot detection and clustering: ways and means

Bayesian 2-stage space-time mixture modeling with spatial misalignment of the exposure in small area health data

Spatiotemporal prediction of fine particulate matter using high resolution satellite images in the southeastern

Estimation of disease rates in small areas: a new mixed model for spatial dependence

Bayesian spatial modelling with R-INLA

An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach

A comparison of numerical schemes on triangular and quadrilateral meshes

Fine scale spatio-temporal modelling of urban air pollution

An empirical assessment of Bayesian melding for mapping ozone pollution

The use of mixed generalized additive modeling to assess the effect of temperature on the usage of emergency electrocardiography examination among the elderly in Shanghai

Relative risk estimation of dengue disease at small spatial scale

Bayesian dynamic modeling of time series of dengue disease case counts

Principles of geostatistics

The current and future global distribution and population at risk of dengue

Understanding the stochastic partial differential equation approach to smoothing

A geostatistical model for combined analysis of point level and area level data using INLA and SPDE

Spatial distribution of HIV prevalence among young people in Mozambique

Diarrhea Morbidities in small areas: accounting for non-stationarity in sociodemographic impacts using Bayesian spatially varying coefficient modelling

Approaches for dealing with various sources of overdispersion in modeling count data: scale adjustment versus modeling

Spatial misalignment in time series studies of air pollution and health data

Spatial and temporal patterns of dengue incidence in northeastern Thailand

Gaussian process emulators for spatial individual level models of infectious disease

Spatiotemporal modeling of dengue fever risk in Puerto Rico

Climate change: An overview of the prevalence of dengue hemorrhagic fever in the South Sulawesi province of Indonesia

On the choice of the mesh for the analysis of geostatistical data using R-INLA

Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations

Assessing the effects on health inequalities of di erential exposure and differential susceptibility of air pollution and environmental noise in Barcelona

Fusing point and areal level space-time data with application to wet deposition

Spatio-temporal disease mapping using INLA

The spatial and temporal scales of local dengue virus transmission in natural settings: A retrospective analysis

Mapping disease at an approximated individual level using aggregate data: a case study of mapping New Hampshire birth defects

Efficient covariance approximations for large sparse precision matrices

Think continuous: Markovian Gaussian models in spatial statistics

Penalising model component complexity: a principled, practical approach to constructing priors

A fast algorithm for generating constrained Delaunay triangulations

Spatiotemporally Varying Coefficients (STVC) model: a Bayesian local regression to detect spatial and temporal nonstationarity in variables relationships

A comparative study of Gaussian geostatistical models and Gaussian Markov random field models

Tutorial: scaling IGMRF-models in R-INLA

Penalised complexity priors for stationary autoregressive processes

Scaling intrinsic Gaussian Markov random field priors in spatial modelling

An examination of disparities in cancer incidence in Texas using Bayesian random coefficient models

Bayesian area-to-point Kriging using expert knowledge as informative priors

Computational fluid dynamics: a practical approach

A spatial regression model for the disaggregation of areal unit-based data to high-resolution grids with application to vaccination coverage mapping

Geoadditive models to assess spatial variation of HIV infections among women in Local communities of Durban South Africa

Generalized spatial fusion model framework for joint analysis of point and areal data

Stochastic processes in several dimensions

On stationary processes in the plane

Cost of dengue illness in Indonesia across hospital, ambulatory, and not medically attended settings

Global climate and weather data, version 2.1

High-resolution mapping of tuberculosis transmission: Whole genome sequencing and phylogenetic modelling of a cohort from Valencia Region

Hierarchical Bayesian modeling of spatio-temporal patterns of lung cancer incidence risk in Georgia

Socioeconomic and environmental determinants of dengue transmission in an urban setting: an ecological study in Noumea New Caledonia

temporal structure matrices by the geometric means of the diagonal elements of their inverses. However, the random walk models of order one (RW1) and two (RW2) are intrinsic Gaussian Markov random fields (IGMRFs). 33 the generalized inverses of the temporal structure matrices varying along with RW1 and RW2, respectively. Sørbye and Rue (2014) proposed to scale the temporal structure matrices of the RW1 and RW2 priors by the geometric mean of the corresponding diagonal elements of their generalized inverses − k : with 2 GV called the generalized variance of the diagonal elements of − k (see Sørbye and Rue (2014) for details). Scaling of R k gives the scaled temporal structure matrixGV and scaled generalized covariance matrix asScaling the RW1 and RW2 priors can be straightforwardly done in R-INLA by using the option scale.model = TRUE (Blangiardo and Cameletti 2015; Sørbye 2013) . After scaling, the same hyperprior IG (1, 0.01) can be assigned to the variance 2 k of the IGMRF priors RW1 and RW2. 34 We followed this procedure in the application.To avoid overfitting, we followed Fuglstad et al. (2020) and applied a weakly informative PC hyperprior to the range (r) and the standard deviation ( Φ ) of the LSPDE model. Particularly, we applied Pr Φ > Φ0 = Φ with Φ0 and Φ denoting the lower limit and lower tail probability of the PC hyperprior for , respectively. The actual value for the standard deviation was 1 and 0.001 for the tail probability, i.e., Pr Φ > 1 = 0.01 . For the range r we selected Pr r < r 0 = r with r 0 and r denoting the upper limit and upper tail probability, respectively. We took prior information r 0 = 5km which is the distance beyond which dengue disease risk is no longer spatially correlated (Sedda et al. 2018) . It corresponds to the median distanceAn IGMRF has a sparse precision matrix which is not of full rank. According to Rue and Held (2005) , the order of an IGMRF is the rank deficiency of its precision matrix, i.e., the number of its zero eigenvalues. Hence, a zero-mean IGMRF of order s is:where k is the (T × T) the semi-positive definite precision matrix of the parameters k given bywhich is equal to the product of the T − s non-zero eigenvalues of k (s = {1,2} for a RW1 and RW2, respectively). The generalized inverse of the precision matrix ( − k ) is the generalized covariance matrix Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. The second number refers to the likelihood: 1: Poisson and 2: Negative Binomial, the third to the time varying effect of the coefficients, 1: random walk of order 1 (RW1) and 2: random walk of order 2 (RW2), the fourth digit defines the edge length, 0: edge length not relevant, 1: 1000 m, 2: 1,100 m, 3: 1,200 m, 4: 1,300 m, 5: 1,400 m, and 6: 1,500 m. Note that there is no edge length for the M1 models because they do not need meshing. M1 contains 4 sub-models, and M2 and M3 contain 24 sub-models each a The FGG-GMRF models are specified as t | t b The following notation applies. M denotes Model. M1 is the sub-model with covariates only, i.e., t = β 0 1 (n A +n p ) + ∑ 4 k=1 � k + k,t � k,t , M2 the model with covariates plus area and cell level interaction, i.e., t = β 0 1 (n A +n p ) + ∑ K k=1 � k + k,t � k,t + t + t (Area interaction type IV with precision matrix = ⊗ ), .M3 the full model, i.e., t = β 0 1 (n A +n p ) + ∑ K k=1 � k + k,t � k,t +̈ +̈ + t 1 (n A +n p ) + t 1 (n A +n p ) + t + t