A combined statistical and machine learning
approach for spatial prediction of extreme
wildfire frequencies and sizes

Daniela Cisneros1, Yan Gong1, Rishikesh Yadav1,
Arnab Hazra1∗, and Raphaël Huser1

11footnotetext:  Computer, Electrical and Mathematical Sciences and
Engineering (CEMSE) Division, King Abdullah University of Science and
Technology (KAUST), Thuwal 23955-6900, Saudi Arabia.
∗E-mail<EMAIL_ADDRESS>

August 28, 2024

Abstract

Motivated by the Extreme Value Analysis 2021 (EVA 2021) data challenge we
propose a method based on statistics and machine learning for the spatial
prediction of extreme wildfire frequencies and sizes. This method is tailored
to handle large datasets, including missing observations. Our approach relies
on a four-stage high-dimensional bivariate sparse spatial model for zero-
inflated data, which is developed using stochastic partial differential
equations (SPDE). In Stage 1, the observations are categorized in zero/nonzero
categories and are modeled using a two-layered hierarchical Bayesian sparse
spatial model to estimate the probabilities of these two categories. In Stage
2, before modeling the positive observations using a spatially-varying
coefficients, smoothed parameter surfaces are obtained from empirical
estimates using fixed rank kriging. This approximate Bayesian method inference
was employed to avoid the high computational burden of large spatial data
modeling using spatially-varying coefficients. In Stage 3, the standardized
log-transformed positive observations from the second stage are further
modeled using a sparse bivariate spatial Gaussian process. The Gaussian
distribution assumption for wildfire counts developed in the third stage, is
computationally effective but erroneous. Thus in Stage 4, the predicted values
are rectified using Random Forests. Posterior inference is drawn for Stages 1
and 3 using Markov chain Monte Carlo (MCMC) sampling. A cross-validation
scheme is then created for the artificially generated gaps, and the EVA 2021
prediction scores of the proposed model are compared to those obtained using
certain natural competitors.

Keywords: Approximate Bayesian inference, Extreme wildfire frequencies and
sizes, GMRF, Random Forests, SPDE.

##  1 Introduction

Wildfires have become one of the important concerns in recent years because
they cause air pollution, extinction of floral and faunal species, significant
economic loss, irreparable damage to the environment and the atmosphere, and
threats to people’s lives. Wildfires occur because of multiple causes, such as
human intervention (e.g., agricultural activities, campfires, and smoking,
Pyne et al.,, 1996), lightning, volcanic eruption, debris burning, sparks from
rock falls, and spontaneous combustion (Scott,, 2000). Typically, wildfires
are accelerated by favorable conditions such as high flammability, extreme
weather-like conditions (e.g., a drought), and the presence of combustible
materials (e.g., forest matter). Recent wildfires in the United States (US)
have led to considerable economic losses and social stresses (Brown et al.,,
2021). Moreover, there is concern that climate change may increase the
intensity, duration, and frequency of wildfires (Abatzoglou and Williams,,
2016; Wuebbles et al.,, 2017; Brown et al.,, 2021). Wildfire prediction is an
important component of wildfire management because it impacts resource
distribution, mitigation of adverse effects, and recovery efforts, and it is
thus of crucial importance to develop resilient statistical methods that can
reliably predict extreme wildfire events over space and time.

From a probabilistic viewpoint, wildfire occurrences and sizes can be seen as
the results of random spatiotemporal processes; thus, it is important to
understand their spatiotemporal distributions and underlying risk factors
associated with them. In particular, wildfires can be linked to their spatial
coordinates, e.g., the location of the fire origin, or the center of a burnt
area, their temporal instant, and other relevant covariates. Extreme wildfires
are attracting considerable attention because they are increasingly hazardous
and often considered a more severe threat to ecosystems.

Multiple statistical approaches have been proposed in the literature to
predict the counts and sizes of wildfires using univariate probability models
(Cumming,, 2001; Preisler et al.,, 2004; Preisler and Westerling,, 2007;
Preisler and Ager,, 2013; Pereira and Turkman,, 2019). Xi et al., (2019) in
particular made multiple important contributions to modeling fire risk
components over recent decades, describing some key yet often overlooked fire
characteristics, and they highlighted various areas of recent research that
may enhance fire risk assesment models. However the spatial/temporal
dependence of wildfires in their modeling approach was disregarded. To develop
effective disaster management policies, statistical modeling of this
dependence is important. Certain papers proposed using the so-called $K$ and
$L$ functions for treating fire occurrences as spatial point-pattern datasets,
(Genton et al.,, 2006; Hering et al.,, 2009; Juan et al.,, 2012). Serra et
al., (2012) and Bivand et al., (2015) treated these point-pattern datasets as
gridded spatial datasets, in which the counts of fire occurrences (CNT) in
each grid cell were modeled as spatially dependent count data; furthermore,
the spatial dependence was modeled in terms of a Gaussian process (GP), which
is the most common tool for modeling spatial dependence because of its
attractive theoretical and computational properties (Gelfand and Schliep,,
2016). A log-Gaussian Cox process (LGCP, Møller et al.,, 1998), is a doubly
stochastic construction, which compromises a Poisson point process at the data
level, characterized with a random log-intensity modeled using a GP at the
latent level. Diggle et al., (2013) reviewed the available literature on LGCPs
in which LGCP models in different applications of spatiotemporal point pattern
analysis were reported. Many researchers used LGCPs for occurrence modeling in
a wide range of contexts; see, for example, Serra et al., (2012), Møller and
Díaz-Avalos, (2010), and Gabriel et al., (2017). Furthermore, Abdelfatah et
al., (2016) and Trucchia et al., (2018) used GPs for modeling burnt area.

For large spatial datasets on a discretized spatial domain, continuos-space
GPs are often replaced by a discrete Gaussian Markov random fields (GMRFs, Rue
and Held,, 2005). GMRFs allow using sparse precision matrices, affording
faster computations. Many researchers used joint analyses to study both fire
occurrences and sizes in large datasets. Ríos-Pena et al., (2018) used a zero-
inflated beta distribution, where a zero inflation was used to model the
absence of fires, while the fraction of the burnt area was modeled using a
beta distribution. Furthermore, the model parameters were assumed to be
spatially varying, and GMRF priors were used in a fully Bayesian analysis.
Joseph et al., (2019) compared certain probability distributions used for
modeling frequencies and sizes of large wildfires to generate a posterior
predictive distribution based on finite sample maxima for extreme events. The
best performance was achieved using a zero-inflated negative binomial model
for CNT and a lognormal model for burnt areas (BAs). Similarly, a marked LGCP
model was proposed by Pimont et al., (2021) in which the authors modeled the
occurrences using a point process and treated the fire sizes as marks.
However, their approach estimated the model components for occurrences and
fire sizes separately, which clearly limits the interaction between these two
components. More recently, Koh et al., (2021) rectified this method by
allowing joint estimation for all components using fully Bayesian inference.

Motivated by the 2021 Extreme Value Analysis conference (EVA 2021) data
challenge, in which we participated as the team named The Bedouins, we here
propose an alternative approach based on statistics and machine learning (ML)
for the spatial prediction of wildfire sizes or BAs and CNT at masked
spatiotemporal locations. The complete dataset includes monthly observations
(from March to September) at 3503 grid cells across the US between 1993 and
2015, in which each of the variables (i.e., BA and CNT) are masked at 80,000
spatiotemporal points. Because the spatial dimension is large and the dataset
includes many zeros, a model that allows a scalable inferential scheme with
high spatial dimension and zero inflation is required. Here, a four-stage
high-dimensional bivariate spatial model is proposed for zero-inflated data.
Our model is developed using stochastic partial differential equations
(SPDEs). In Stage 1, the observations are categorized in zero/nonzero
categories (zero BA indicates zero CNT and vice versa). Moreover, a two-
layered hierarchical Bayesian model is fitted, whereby the first layer is used
for defining the zero/nonzero data in two categories of a real-valued latent
process, and the second layer is used for developing an SPDE-based
construction of the latent process. This model is used to estimate the
probabilities of the categories at unobserved spatiotemporal sites. In Stage
2, before modeling the positive observations using a log-Gaussian process with
spatially-varying parameters, the model parameters are empirically estimated
(i.e., sample means and standard deviations (SDs) of the log-transformed
positive observations) at each spatial location and smoothed parameter
surfaces are obtained using fixed rank kriging (FRK, Cressie and Johannesson,,
2008). Fully Bayesian inference with spatially–varying parameters involves a
large computational burden, which can be avoided using an efficient
approximate Bayesian inference technique. In Stage 3, the standardized log-
transformed positive observations from the second stage are modeled using a
bivariate spatial GMRF. Despite the computational advantages of GP-based
modeling, the model in Stage 3 erroneously assumes that the marginal
distribution of CNT is Gaussian. Thus, the predicted values of CNT are finally
corrected in Stage 4 using Random Forests (RF) in which BAs are treated as a
covariate, and the missing BA values are imputed by the predicted BA values in
Stage 3. Recently, ML algorithms, such as RF and neural networks (NN), have
been successfully used to model fire occurrence data (Jain et al.,, 2020). In
this study, CNT data are independently calibrated at each spatial location.
Similarly, in Saha et al., (2021) was recently proposed a RF technique for
spatially dependent data. Here, we use the simulation-based Markov Chain Monte
Carlo (MCMC) method to draw posterior inference in Stages 1 and 3. This method
involves moderate computational time even for high spatial dimension because
of the sparse spatial structure implied by the SPDE. Our modeling framework
targets wildfire prediction in the US; however, it is worth noting that it can
be also adapted to other data scenarios with large spatial dimensions and zero
inflation.

The paper is structured as follows: In Section 2, an exploratory analysis of
the US wildfire dataset is discussed. Section 3 shows the development of an
approach combining statistics and RF for the joint modeling of BA and CNT. In
Section 4, a brief overview of computational details is obtained. The proposed
approach is applied to the US wildfire dataset and the results are then
discussed in Section 5. The conclusions of this study and perspectives for
future research are presented in Section 6.

##  2 The US Wildfire dataset and exploratory analysis

In this section, the primary features of the US wildfire dataset are described
and an exploratory graphical support for our modeling choices is provided.

###  2.1 Data description and the missing data pattern

The wildfire dataset for the EVA 2021 data challenge comprises monthly
observations at 3503 grid cells across the US Mainland using a spatial
resolution of $0.5^{\circ}\times 0.5^{\circ}$. The original dataset (before
masking) contains wildfire counts (CNTs) and aggregated burnt areas (BAs) in
each pixel over a period of 23 years between (1993–2015) and for seven months
per year (March to September). Information about 18 spatiotemporal land cover
covariates (e.g., proportion of urban area, shrubland, grassland), 10
spatiotemporal meteorological covariates (e.g., temperature at 2 m above the
ground, precipitation, evaporation of water), as well as certain purely
spatial and temporal covariates (mean and standard deviation of the altitude,
proportion of a pixel that is within the US Mainland, longitude and latitude
for the center of the pixel, year, and month) was available. The organizers
used the Shuttle Radar Topography Mission (SRTM) database available at a 90–m
spatial resolution to calculate altitude-related covariates. A detailed
description of the dataset is available at (Opitz,, 2022). The primary aim of
this study was to estimate the predictive distribution function of BA and CNT
at 28 severity thresholds as follows:

$\displaystyle\mathcal{U}_{\textrm{CNT}}$ $\displaystyle=$
$\displaystyle\\{0,1,2,\ldots,9,10,12,14,\cdots,30,40,50,\cdots,100\\},$
$\displaystyle\mathcal{U}_{BA}$ $\displaystyle=$
$\displaystyle\\{0,1,10,20,30,\ldots,100,150,200,250,300,400,500,1000,$ (1)
$\displaystyle 1500,2000,5000,10000,20000,30000,40000,50000,100000\\}.$

The final evaluation of all the registered teams was based on the weighted sum
of the squared error between the empirical and predictive distribution
functions, with weights given by
$\tilde{\omega}_{\textrm{BA}}(u)=\tilde{\omega}_{\textrm{CNT}}(u)=1-(1+(u+1)^{2}/1000)^{-1/4}$,
$u\in\mathcal{U}_{\text{CNT}}$ or $\mathcal{U}_{\text{BA}}$, then rescaled to
add up to 1.

The original dataset does not contain any missing values. However, the
organizers masked a total of 80,000 observations (14.18%) across all the seven
months of the alternative even years (1994, 1996, 1998, 2000, 2002, 2004,
2006, 2008, 2010, 2012, and 2014) to compare the spatial prediction
performances of models proposed by participating teams. The masked
spatiotemporal locations are not similar for BA and CNT; furthermore, the data
were removed from low to high fire-prone regions. This demonstrates the
requirement for modeling spatial dependence to borrow information from nearby
pixels. Figure 1 shows the spatial maps of BA and CNT in March 1994, where the
masked data locations are highlighted in white. The masked observations were
from small and large clusters of pixels and from regions in the southeast of
the US, where high CNT values are observed close to masked pixels. Table 1
shows the proportions of masked observations of BA and CNT. As reported,
either BA or CNT data or both are missing in 14.18% of total cases. They are
jointly missing in 8.68% cases. Thus, BA and CNT information can be borrowed
for the rest of 5.50% non-missing observations by joint modeling, which leads
to smaller standard errors of model parameters.

Figure 1: Spatial maps of BA (left) and CNT (right) in March 1994. The locations in which data are masked are highlighted in white. Table 1: The proportions of zero, nonzero, and missing values of BA and CNT in the total 563,983 observations. CNT \BA | Zero | Nonzero | Missing | Total
---|---|---|---|---
Zero | 49.60% | 0.00% | 3.34% | 52.94%
Nonzero | 0.00% | 30.70% | 2.17% | 32.87%
Missing | 3.32% | 2.18% | 8.68% | 14.18%
Total | 52.93% | 32.89% | 14.18% | 100.00%

###  2.2 Zero inflation

The high proportion of zero observations in the US wildfire dataset is one of
its important characteristics; BA being zero/nonzero at a spatiotemporal
location is equivalent to CNT being zero/nonzero. Table 1 lists the
proportions of zero and nonzero values for BA and CNT. Using this equivalence,
some of the missing data can be retrieved. CNT is zero in 3.34% (out of the
14.18%) of the cases where BA has missing value; hence, the BA values are
zeros in these instances. Similarly, BAs are zero in 3.32% (out of the 14.18%)
of the cases where CNT has missing values; hence, CNT values in these cases
are zero. Thus, by filling these masked locations with zeros, a dataset with a
smaller number of missing values can be obtained (specifically, 48,947 cases
instead of 80,000 cases). Here, out of the final available observations,
61.61% of the values are zeros. Therefore, the available information can be
divided into two parts. In the first part, a spatiotemporal dataset of binary
observations indicating whether BA/CNT is zero or not, is obtained; in the
second part, only the positive values are maintained (treating zeros as
missing data). For modeling positive BA and CNT, the logarithmic
transformation can be used, and we call the resulting transformed dataset
$\log$-BA and $\log$-CNT, respectively.

We now explore the covariate effects on the zero/nonzero indicators. The
results of a probit regression model demonstrated significance for certain
examined covariates at a significance level of 0.01. However, assuming a
single regression coefficient for the entire spatiotemporal domain is not
realistic. Despite increasing the computational burden, a cross-validation
study does not demonstrate any significant improvement in the prediction
performance after incorporating the available covariate information. Thus, to
simplify calculations we choose to ignore the covariate information.
Similarly, the spatiotemporal covariates can be ignored in the modeling of
log-BA and log-CNT, unless they do not increase significantly the
computational burden.

Figure 2: Left: empirical spatial correlation (smoothed) profiles of the
zero/nonzero indicators (black), log-BA (red), and log-CNT (blue). Right:
histogram of the empirical cross-correlation between log-BA and log-CNT across
space.

We then explore the spatial and temporal correlation profiles of the
zero/nonzero indicators, log-BA, and log-CNT, as well as the cross-correlation
between log-BA and log-CNT. The left panel of Figure 2 shows (smoothed)
empirical spatial correlation profiles as a function of distance. All three
spatial correlation profiles decrease towards zero with the increase in
geographical distance between pixels, and they demonstrate similar spatial
range and small-scale variability (nugget effect). Thus, a reasonable
separable correlation structure can be assumed for multivariate spatial
modeling of log-BA and log-CNT. For every spatial location, the Lag-1 temporal
autocorrelation was empirically calculated and was reported to be not
significant at a significance level 0.01 for a large proportion of spatial
locations (64.13% locations for the zero/nonzero indicators, 96.39% locations
for log-BA, 89.46% locations for log-CNT). A cross-validation study does not
demonstrate any significant improvement in the prediction performance by
incorporating temporal dependence. Thus, ignoring the temporal autocorrelation
for all three spatiotemporal processes is reasonable, particularly considering
the additional computational burden in a high spatial dimension. The cross-
correlation between log-BA and log-CNT is empirically calculated at every
spatial location (based on the temporal replicates), and the right panel of
Figure 2 shows the histogram of empirical cross-correlation values. The cross-
correlation is quite high for most spatial locations, which indicates the
requirement for joint statistical modeling of log-BA and log-CNT.

Finally, we explore the requirement for spatially-varying marginal
distribution parameters for log-BA and log-CNT. When log-BA and log-CNT are
modeled using a probability distribution from a location-scale family, thus
allowing a spatially-varying location profile is more common than allowing a
spatially-varying scale profile, because of the computational benefits of the
former. Figure 3 shows the empirical location-wise standard deviations of log-
BA and log-CNT. For log-BA, standard deviations are generally lower in Eastern
US than in Western US. For log-CNT, standard deviations are lower in the
middle regions of US compared to the states closer to the Atlantic or the
Pacific coasts. A similar spatial pattern is observed for the location-wise
mean values of log-BA and log-CNT. Thus, a joint analysis of log-BA and log-
CNT using a model with spatially-varying location and scale parameters is
required.

Figure 3: Pixel-wise standard deviation of log-BA (left) and log-CNT (right).
Pixels where all observations are zero or that have only one positive
observation (i.e. standard deviations are not available) are shown in white.

##  3 Joint modeling of burnt area (BA) and counts of fire occurrences (CNT)

In this section, a four-stage model based on statistics and machine learning
for the joint analysis of BA and CNT is described. As reported in Section 2,
the available information is divided into two parts. In the first part, a
spatiotemporal dataset of binary observations is obtained to determine whether
BA/CNT values are zero; in the second part, only the positive values are
maintained, and log-BA and log-CNT are modeled on the logarithmic scale. In
Stage 1, a sparse latent GP model for binary spatial data is proposed. For the
other three stages, the modeling of log-BA and log-CNT using a combination of
an approximate Bayesian inference technique and Random Forests (RFs) is
described.

###  3.1 Stage 1: A sparse latent Gaussian process model for wildfire
occurrence data

The wildfire occurrence data are assumed independent and identically
distributed (IID) across months and years; furthermore the observations across
the US mainland for each month are assumed to only be spatially dependent.
Thus, by ignoring the month–year combinations, a generic notation is used for
temporal replications.

For a spatial location $\bm{s}_{i}$ and time $t$, the BA and CNT values are
denoted $\textrm{BA}_{t}(\bm{s}_{i})$ and $\textrm{CNT}_{t}(\bm{s}_{i})$,
respectively, where $i\in\\{1,\ldots,N\\}$ and $t\in\\{1,\ldots,T\\}$, $N$ is
the total number of pixels ($N=3503$), and $T$ is the total number of months
($T=161$). The wildfire occurrence indicator $Z_{t}(\bm{s}_{i})$ at location
$\bm{s}_{i}$ and time $t$ is defined as

$\displaystyle Z_{t}(\bm{s}_{i})=\begin{cases}1,\hskip
11.38092pt\textrm{if}\begin{cases}\parbox[t]{381.5874pt}{$\textrm{BA}_{t}(\bm{s}_{i})>0,\leavevmode\nobreak\
\textrm{CNT}_{t}(\bm{s}_{i})>0$,}\\\
\text{$\textrm{BA}_{t}(\bm{s}_{i})>0,\textrm{CNT}_{t}(\bm{s}_{i})$ is missing,
$\leavevmode\nobreak\ ${or}$\leavevmode\nobreak\ \leavevmode\nobreak\ $}\\\
\text{$\textrm{CNT}_{t}(\bm{s}_{i})>0,\textrm{BA}_{t}(\bm{s}_{i})$ is
missing,}\end{cases}\\\ \vspace{-5mm}\\\ 0,\hskip
11.38092pt\textrm{if}\begin{cases}\parbox[t]{381.5874pt}{$\textrm{BA}_{t}(\bm{s}_{i})=0,\leavevmode\nobreak\
\textrm{CNT}_{t}(\bm{s}_{i})=0$,}\\\
\text{$\textrm{BA}_{t}(\bm{s}_{i})=0,\textrm{CNT}_{t}(\bm{s}_{i})$ is missing,
$\leavevmode\nobreak\ ${or}$\leavevmode\nobreak\ \leavevmode\nobreak\ $}\\\
\text{$\textrm{CNT}_{t}(\bm{s}_{i})=0,\textrm{BA}_{t}(\bm{s}_{i})$ is
missing,}\end{cases}\\\ ${NA},$\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{if \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $\textrm{BA}_{t}(\bm{s}_{i})$ and
$\textrm{CNT}_{t}(\bm{s}_{i})$ are both missing.}\end{cases}$

By assumption the replicated indicator processes, $Z_{t}(\cdot),t=1,\ldots,T$,
are IID across time $t$, and we model $Z_{t}(\cdot)$ as

$\displaystyle
Z_{t}(\bm{s}_{i})=\begin{cases}1,&\parbox[t]{95.39502pt}{if\leavevmode\nobreak\
\leavevmode\nobreak\ $X_{t}(\bm{s}_{i})>0$}\\\ 0,&\text{if\leavevmode\nobreak\
\leavevmode\nobreak\
$X_{t}(\bm{s}_{i})<0$},\end{cases}\textrm{where}\leavevmode\nobreak\
X_{t}(\bm{s}_{i})=\mu_{Z}(\bm{s}_{i})+\varepsilon_{t}(\bm{s}_{i}),$ (2)

and $\varepsilon_{t}(\cdot),t=1,\ldots,T$, are IID spatial GPs, in which
$\textrm{E}[\varepsilon_{t}(\bm{s}_{i})]=0$ and
$\textrm{Var}[\varepsilon_{t}(\bm{s}_{i})]=1$ for all $i\in\\{1,\ldots,N\\}$.
The process $\varepsilon_{t}(\cdot)$ is assumed to follow an isotropic Matérn
spatial correlation (with nugget effect) given by

$\displaystyle\rho_{\varepsilon}(\bm{s}_{i},\bm{s}_{j})={\textstyle{r_{\varepsilon}\over\Gamma(\nu)2^{\nu-1}}}\left({\textstyle{d(\bm{s}_{i},\bm{s}_{j})\over\phi_{\varepsilon}}}\right)^{\nu}K_{\nu}\left({\textstyle{d(\bm{s}_{i},\bm{s}_{j})\over\phi_{\varepsilon}}}\right)+(1-r_{\varepsilon})\mathbb{I}(\bm{s}_{i}=\bm{s}_{j}),$
(3)

where $d(\bm{s}_{i},\bm{s}_{j})$ is the Euclidean distance between
$\bm{s}_{i}$ and $\bm{s}_{j}$, $\phi_{\varepsilon}>0$, $\nu>0$ and
$r_{\varepsilon}\in[0,1]$ are the range, smoothness, and ratio of the spatial
to total variation, respectively. In (3), $K_{\nu}$ is the modified Bessel
function of the degree $\nu$, and $\mathbb{I}(\bm{s}_{i}=\bm{s}_{j})=1$ if
$\bm{s}_{i}=\bm{s}_{j}$, and 0 otherwise. When $r_{\varepsilon}=1$,
$\epsilon_{t}(\cdot)$ is mean-square differentiable if $\nu$ is an integer.
For practical applications, identifying $\nu$ is difficult, and thus, it is
generally fixed a priori. Here, the process $\varepsilon_{t}(\cdot)$ is not
observable and hence estimating $\nu$ is more challenging. Therefore, we here
set $\nu$ to one. To overcome the high computational burden due to the large
spatial dimension, $\varepsilon_{t}(\cdot)$ is defined as a Gaussian Markov
random field (GMRF) that has an (approximately) equivalent spatial covariance
structure to the dense GP with a spatial Matérn
$\rho_{\varepsilon}(\cdot,\cdot)$. As described in Lindgren et al., (2011),
the equivalence is derived from the one-to-one link between dense isotropic
Matérn GPs and GMRFs. We next briefly summarize this link.

Suppose that $\varepsilon(\cdot)$ is a dense GP with correlation structure (3)
and $r_{\varepsilon}=1$. Then, $\varepsilon(\cdot)$ is the solution to the
SPDE
$(8\phi_{\varepsilon}^{-2}-\Delta)\varepsilon(\bm{s})=\mathcal{W}(\bm{s})$,
where
$\Delta={\textstyle{\delta^{2}\over\delta^{2}x}}+{\textstyle{\delta^{2}\over\delta^{2}y}}$
is the Laplacian operator, $(8\phi_{\varepsilon}^{-2}-\Delta)$ is a pseudo-
differential operator, and
$\mathcal{W}(\bm{s})\overset{\textrm{IID}}{\sim}\textrm{Normal}(0,1)$. We can
solve this SPDE using finite element methods (Ciarlet,, 2002) over a
triangular mesh in $\mathbb{R}^{2}$, where the triangles are formed following
a Delaunay triangulation. Let the set of mesh nodes be denoted by
$\mathcal{S}^{*}=\\{\bm{s}^{*}_{1},\ldots,\bm{s}^{*}_{N^{*}}\\}$. We construct
a finite element representation of the solution to
$\varepsilon(\bm{s})=\sum_{j=1}^{N^{*}}\psi_{j}(\bm{s})\varepsilon^{*}_{j}$
for some chosen basis functions $\psi_{j}(\cdot)$ and normally distributed
weights $\varepsilon^{*}_{j}$ defined at the mesh nodes $\mathcal{S}^{*}$. We
calculate the inner products $\langle\,\psi_{j}(\cdot),1\rangle$ and
$\langle\,\nabla\psi_{j_{1}}(\cdot),\nabla\psi_{j_{2}}(\cdot)\rangle$, where
$\langle\,f,g\rangle=\int f(\bm{s})g(\bm{s})d\bm{s}$, and obtain three
$(N^{*}\times N^{*})$-dimensional finite element matrices $\bm{C}$,
$\bm{G}_{1}$, and $\bm{G}_{2}$. Here, $\bm{C}$ is a diagonal matrix, of which
the $(j,j)^{th}$ entry $C_{j,j}=\langle\,\psi_{j}(\cdot),1\rangle$,
$\bm{G}_{1}$ is a sparse matrix, of which the $(j_{1},j_{2})^{th}$ entry
$G_{1_{j_{1},j_{2}}}=\langle\,\nabla\psi_{j_{1}}(\cdot),\nabla\psi_{j_{2}}(\cdot)\rangle$,
and $\bm{G}_{2}=\bm{G}_{1}\bm{C}^{-1}\bm{G}_{1}$. Further theoretical details
are discussed in Bakka et al., (2018). The vector
$\bm{\varepsilon}^{*}=[\varepsilon^{*}_{1},\ldots,\varepsilon^{*}_{N^{*}}]^{\prime}\sim\textrm{Normal}_{N^{*}}(\bm{0},\bm{Q}_{\phi_{\varepsilon}}^{-1})$,
where the precision matrix is
$\bm{Q}_{\phi_{\varepsilon}}=(4\pi)^{-1}\phi_{\varepsilon}^{2}\left[\phi_{\varepsilon}^{-4}\bm{C}+2\phi_{\varepsilon}^{-2}\bm{G}_{1}+\bm{G}_{2}\right]$.
To project $\bm{\varepsilon}^{*}$ back to the data locations $\mathcal{S}$, we
evaluate $a_{ij}=\psi_{j}(\bm{s}_{i})$ for each location $\bm{s}_{i}$ and mesh
node $\bm{s}_{j}^{*}$. The $(N\times N^{*})$-dimensional matrix $\bm{A}$, of
which $(i,j)^{th}$ entry is $a_{ij}$, is called the SPDE projection matrix,
and
$\bm{A}\bm{\varepsilon}^{*}\sim\textrm{Normal}_{N}(\bm{0},\bm{A}\bm{Q}_{\phi_{\varepsilon}}^{-1}\bm{A}^{\prime})$.
The covariance matrix $\bm{A}\bm{Q}_{\phi_{\varepsilon}}^{-1}\bm{A}^{\prime}$
approximates the Matérn correlation matrix obtained by evaluating (3) at
$\mathcal{S}$ (for $r_{\varepsilon}=1$).

Figure 4: Left: triangulated mesh over the United States, which is used to
develope the spatial process $\varepsilon_{t}(\cdot)$, using SPDEs. Right:
comparison of a possible Matérn correlation structure (red) and the pairwise
covariances (as a function of distance) between two spatial locations obtained
from the SPDE approximation (black dots).

Suppose that for each $t\in\\{1,\ldots,T\\}$, $\bm{\varepsilon}^{*}_{t}$ is an
independent copy of $\bm{\varepsilon}^{*}$. We construct a vector
$\bm{\varepsilon}_{t}=[\varepsilon_{t}(\bm{s}_{1}),\ldots,\varepsilon_{t}(\bm{s}_{N})]^{\prime}$
from $\bm{\varepsilon}^{*}_{t}$ as

$\displaystyle\bm{\varepsilon}_{t}=\sqrt{r_{\varepsilon}}\bm{A}\bm{\varepsilon}^{*}_{t}+\sqrt{1-r_{\varepsilon}}\tilde{\bm{\varepsilon}}_{t},$
(4)

where
$\tilde{\bm{\varepsilon}}_{t}=[\tilde{\varepsilon}_{t}(\bm{s}_{1}),\ldots,\tilde{\varepsilon}_{t}(\bm{s}_{N})]^{\prime}$
with
$\tilde{\varepsilon}_{t}(\bm{s}_{i})\overset{\textrm{IID}}{\sim}\textrm{Normal}(0,1)$.
The final covariance matrix of $\bm{\varepsilon}_{t}$ is
$\bm{\Sigma}_{\varepsilon}=r_{\varepsilon}\bm{A}\bm{Q}_{\phi_{\varepsilon}}^{-1}\bm{A}^{\prime}+(1-r_{\varepsilon})\bm{I}_{N}$,
and it approximates the Matérn correlation matrix obtained by evaluating (3)
at $\mathcal{S}$, for any $r_{\varepsilon}\in[0,1]$. We then discuss the GMRF
approximation accuracy in our data application, for the values of
$\phi_{\varepsilon}$ and $r_{\varepsilon}$ set to 3 and 0.8 respectively.
These values are similar to estimates obtained in our data application. The
left panel of Figure 4 shows the mesh that we used for the US wildfire data
analysis; here, there are $N^{*}=1027$ mesh nodes. In this SPDE mesh, we
calculate the covariance between every pair of spatial locations (the elements
of the matrix $\bm{\Sigma}_{\varepsilon}$) and the corresponding true Matérn
correlation based on (3). The results of these calculations are presented as a
function of distance in the right panel of Figure 4. The true correlation
structure of dense GP is well approximated by the corresponding GMRF.
Moreover, the sparsity of $\bm{Q}_{\phi_{\varepsilon}}$ can be exploited to
allow quick computations. Moreover, while categorizing the latent process
given in (2), the conditional spatial independence structure
$\bm{\varepsilon}_{t}|\bm{\varepsilon}^{*}_{t}\sim\textrm{Normal}_{N}(\sqrt{r_{\varepsilon}}\bm{A}\bm{\varepsilon}^{*}_{t},(1-r_{\varepsilon})\bm{I}_{N})$
allows the univariate imputation of latent variables, which is exploited for
spatial prediction.

After marginalization with respect to $X_{t}(\bm{s}_{i})$, we obtain that
$\textrm{Pr}(Z_{t}(\bm{s}_{i})=1)\approx\Phi(\mu_{Z}(\bm{s}_{i}))$, where
$\Phi(\cdot)$ is the standard normal distribution function. The conditional
distribution of $Z_{t}(\bm{s}_{i})$ given $\bm{\varepsilon}^{*}_{t}$ is
$\textrm{Pr}(Z_{t}(\bm{s}_{i})=1|\bm{\varepsilon}^{*}_{t})=\Phi((\mu_{Z}(\bm{s}_{i})+\sqrt{r_{\varepsilon}}\bm{a}^{\prime}_{i}\bm{\varepsilon}^{*}_{t})/\sqrt{1-r_{\varepsilon}})$,
where $\bm{a}_{i}$ is the $i^{th}$ row of $\bm{A}$, and
$\textrm{Pr}(Z_{t}(\bm{s}_{i})=0|\bm{\varepsilon}^{*}_{t})=1-\textrm{Pr}(Z_{t}(\bm{s}_{i})=1|\bm{\varepsilon}^{*}_{t})$.
The variables $Z_{t}(\bm{s}_{i})$ and $Z_{t}(\bm{s}_{j})$ are conditionally
independent given $\bm{\varepsilon}^{*}_{t}$. There is no closed form
expression for the joint distribution of $Z_{t}(\bm{s}_{i})$ and
$Z_{t}(\bm{s}_{j})$, after marginalizing with respect to $X_{t}(\bm{s}_{i})$.

###  3.2 Stage 2: Approximate Bayesian inference to smooth spatially-varying
parameters

Despite the discreteness of CNTs, we model log-BA and log-CNT using a GP with
spatially-varying location and scale terms, because of the computational
attractiveness of GPs, as follows

$\displaystyle\log\text{-BA}_{t}(\bm{s}_{i})$ $\displaystyle=$
$\displaystyle\mu_{1}(\bm{s}_{i})+\sigma_{1}(\bm{s}_{i})W_{t1}(\bm{s}_{i}),$
$\displaystyle\log\text{-CNT}_{t}(\bm{s}_{i})$ $\displaystyle=$
$\displaystyle\mu_{2}(\bm{s}_{i})+\sigma_{2}(\bm{s}_{i})W_{t2}(\bm{s}_{i}),$
(5)

where
$\bm{W}_{t}(\bm{s}_{i})=[W_{t1}(\bm{s}_{i}),W_{t2}(\bm{s}_{i})]^{\prime}$ is a
bivariate standard GP (zero mean and unit variance for all marginal
distributions). A full Bayesian inference is computationally challenging, and
thus, in Stage 2, we focus on estimating the parameter surfaces
$\mu_{1}(\cdot)$, $\sigma_{1}(\cdot)$, $\mu_{2}(\cdot)$, and
$\sigma_{2}(\cdot)$ only, while the parameters within the bivariate spatial GP
$\bm{W}_{t}(\cdot)$ are treated as nuisance parameters. Here, the procedure
for $\mu_{1}(\cdot)$, which is the same procedure used to estimate the
surfaces $\log[\sigma_{1}(\cdot)]$, $\mu_{2}(\cdot)$, and
$\log[\sigma_{2}(\cdot)]$ is detailed; all surfaces evaluated at any
$\bm{s}_{i}$ are defined over the whole real line.

An approximate Bayesian inference scheme, similar to Max-and-Smooth
(Hrafnkelsson et al.,, 2021; Johannesson et al.,, 2021), is applied to obtain
parameter surfaces in two steps. In the first step, we estimate
$\mu_{1}(\bm{s}_{i})$, $\sigma_{1}(\bm{s}_{i})$, $\mu_{2}(\bm{s}_{i})$, and
$\sigma_{2}(\bm{s}_{i})$ at each $\bm{s}_{i}$ separately, using location-wise
empirical means and standard deviations, ignoring any spatial/temporal trend
or dependence. We denote the estimates by ${\widehat{\mu}}_{1}(\bm{s}_{i})$,
$\log[{\widehat{\sigma}}_{1}(\bm{s}_{i})]$, ${\widehat{\mu}}_{2}(\bm{s}_{i})$,
and $\log[{\widehat{\sigma}}_{2}(\bm{s}_{i})]$.

In the second step we smooth parameter surfaces by treating the preliminary
estimates as noisy measurements of the true underlying parameters.
Specifically, for the parameter surface $\mu_{1}(\cdot)$, we assume that
${\widehat{\mu}}_{1}(\bm{s}_{i})=\mu_{1}(\bm{s}_{i})+e(\bm{s}_{i})$, where
$\mu_{1}(\bm{s}_{i})$ is the true parameter value; furthermore the estimate
${\widehat{\mu}}_{1}(\bm{s}_{i})$ is perturbed from the true parameter by a
pure nugget term $e(\bm{s}_{i})\sim\textrm{Normal}(0,\sigma^{2}_{e})$. Because
of the large spatial dimension, the prior for $\mu_{1}(\bm{s}_{i})$ is assumed
to follow a low-rank structure as

$\mu_{1}(\bm{s}_{i})=\beta_{0}+\beta_{1}\textrm{lon}(\bm{s}_{i})+\beta_{2}\textrm{lat}(\bm{s}_{i})+\sum_{r=1}^{3}\sum_{k=1}^{K_{r}}h_{rk}(\bm{s}_{i})\omega^{*}_{rk}+\xi(\bm{s}_{i}),\leavevmode\nobreak\
\leavevmode\nobreak\ \bm{s}_{i}\in\mathcal{S},$ (6)

where $h_{rk}(\cdot)$ is the $k^{th}$ spatial Gaussian kernel at the $r^{th}$
resolution, the spatial random effects are
$\bm{\omega}^{*}_{r}=[\omega^{*}_{rk_{1}},\ldots,\omega^{*}_{rk_{r}}]^{\prime}\overset{\textrm{Indep}}{\sim}\textrm{Normal}_{K_{r}}(\bm{0},\bm{\Sigma}(\bm{\theta}_{r}))$,
and
$\xi(\bm{s})\overset{\textrm{IID}}{\sim}\textrm{Normal}(0,\sigma^{2}_{\xi})$.
By selecting flat priors for hyperparameters, the estimation of
$\mu_{1}(\bm{s}_{i}),\bm{s}_{i}\in\mathcal{S}$, lies in the setting of fixed
rank kriging (FRK) using a frequentist approach, which is readily implemented
using the R package FRK (Zammit-Mangion and Cressie,, 2021). For certain
pixels, there was no available positive observation. In these instances
${\widehat{\mu}}_{1}(\bm{s}_{i})$ is treated as missing and the corresponding
values of $\mu_{1}(\bm{s}_{i})$ are predicted based on the available first-
step estimates. We repeat the same procedure for other parameter surfaces and
obtain smoothed estimates, say, $\tilde{\mu}_{1}(\cdot)$,
$\tilde{\sigma}_{1}(\cdot)$, $\tilde{\mu}_{2}(\cdot)$, and
$\tilde{\sigma}_{2}(\cdot)$. Finally, we also obtain
$\bm{W}_{t}(\bm{s}_{i}),\bm{s}_{i}\in\mathcal{S},t=1,\ldots,T$, by plugging
the smoothed estimates in (3.2).

###  3.3 Stage 3: Bivariate spatial modeling of standardized log-BA and log-
CNT

In this stage, we model the standardized variables
$\widehat{W}_{t1}(\bm{s}_{i})=\\{\log\text{-BA}_{t}(\bm{s}_{i})-\mu_{1}(\bm{s}_{i})\\}/\sigma_{1}(\bm{s}_{i})$
and
$\widehat{W}_{t2}(\bm{s}_{i})=\\{\log\text{-CNT}_{t}(\bm{s}_{i})-\mu_{2}(\bm{s}_{i})\\}/\sigma_{2}(\bm{s}_{i})$,
obtained in Stage 2 using (3.2). Suppose that
$\widehat{\bm{W}}_{tp}=[\widehat{W}_{tp}(\bm{s}_{1}),\ldots,\widehat{W}_{tp}(\bm{s}_{N})]^{\prime}$
for $p=1,2$, and
$\widehat{\bm{W}}_{t}=[\widehat{\bm{W}}^{\prime}_{t1},\widehat{\bm{W}}^{\prime}_{t2}]^{\prime}$.
We model $\widehat{\bm{W}}_{t}$ as
$\widehat{\bm{W}}_{t}=\left[\bm{I}_{2}\otimes\bm{A}\right]\bm{\eta}^{*}_{t}+\tilde{\bm{\eta}}_{t}$,
where

$\displaystyle\bm{\eta}^{*}_{t}\sim\textrm{Normal}_{2N}\left(\bm{0}_{2N},r_{\eta}\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\
1\end{array}\right)\otimes\bm{Q}^{-1}_{\phi_{\eta}}\right),\leavevmode\nobreak\
$ (9)

$\tilde{\bm{\eta}}_{t}\sim\textrm{Normal}_{2N}\left(\bm{0}_{2N},(1-r_{\eta})\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\ 1\end{array}\right)\otimes\bm{I}_{N}\right).$

Similarly to Stage 1, we model the processes $\widehat{W}_{tp}(\cdot)$ using
GMRFs. The same SPDE mesh and SPDE projection matrix as those used in Stage 1
are used. The spatial correlation of each component $\widehat{W}_{tp}(\cdot)$
is approximately equal to (3), with $\phi_{\varepsilon}$ and $r_{\varepsilon}$
replaced by $\phi_{\eta}$ and $r_{\eta}$, respectively. This is is confirmed
by similar empirical spatial correlation profiles to those of log-BA and log-
CNT shown in Figure 2. The marginal standard deviation of each
$\widehat{W}_{tp}(\bm{s}_{i})$ is approximately equal to one (as described in
Stage 1); thus, the cross-covariance between $\widehat{W}_{t1}(\bm{s}_{i})$
and $\widehat{W}_{t2}(\bm{s}_{i})$ is approximately equal to the cross-
correlation $\rho_{\eta}\in[-1,1]$ for each $\bm{s}_{i}$ and $t$. To
summarize, an approximately separable bivariate GMRF, where the marginal
distributions have zero mean and unit variance, is defined.

At a spatiotemporal prediction location $(\bm{s}_{i},t)$, we need to simulate
the missing process $\widehat{W}_{tp}(\bm{s}_{i})$, possibly for both $p=1,2$.
Let the elements of $2N$-length vector
$\left[\bm{I}_{2}\otimes\bm{A}\right]\bm{\eta}^{*}_{t}$ be denoted by
$[\bm{\tilde{\eta}}^{*^{\prime}}_{t1},\bm{\tilde{\eta}}^{*^{\prime}}_{t2}]^{\prime}$,
where
$\bm{\tilde{\eta}}^{*}_{t1}=[\tilde{\eta}^{*}_{t1}(\bm{s}_{1}),\ldots,\tilde{\eta}^{*}_{t1}(\bm{s}_{N})]^{\prime}$
and
$\bm{\tilde{\eta}}^{*}_{t2}=[\tilde{\eta}^{*}_{t2}(\bm{s}_{1}),\ldots,\tilde{\eta}^{*}_{t2}(\bm{s}_{N})]^{\prime}$.
Then,

$\displaystyle\widehat{\bm{W}}_{t}(\bm{s}_{i})\sim\textrm{Normal}_{2}\left(\left(\begin{array}[]{c}\tilde{\eta}^{*}_{t1}(\bm{s}_{i})\\\
\tilde{\eta}^{*}_{t2}(\bm{s}_{i})\end{array}\right),(1-r_{\eta})\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\ 1\end{array}\right)\right).$ (13)

The predicted values are denoted by $\tilde{W}_{tp}(\bm{s}_{i})$, possibly for
both $p=1$ and $p=2$. Then, the predicted values of log-BA and log-CNT are
obtained by plugging $\tilde{\mu}_{1}(\cdot)$, $\tilde{\sigma}_{1}(\cdot)$,
$\tilde{\mu}_{2}(\cdot)$, and $\tilde{\sigma}_{2}(\cdot)$ (obtained in Stage
2) and $\tilde{W}_{tp}(\bm{s}_{i})$ (obtained in Stage 3) into (3.2). The
distribution functions of BA and CNT at $(\bm{s}_{i},t)$ are approximately
$\textrm{Lognormal}(\mu_{1}(\bm{s}_{i}),\sigma^{2}_{1}(\bm{s}_{i}))$ and
$\textrm{Lognormal}(\mu_{2}(\bm{s}_{i}),\sigma^{2}_{2}(\bm{s}_{i}))$
distribution functions, and they are denoted $F_{LN}(\cdot\leavevmode\nobreak\
;\leavevmode\nobreak\ \mu_{1}(\bm{s}_{i}),\sigma^{2}_{1}(\bm{s}_{i}))$ and
$F_{LN}(\cdot\leavevmode\nobreak\ ;\leavevmode\nobreak\
\mu_{2}(\bm{s}_{i}),\sigma^{2}_{2}(\bm{s}_{i}))$, respectively. Furthermore,
by incorporating zero inflation, $F_{\textrm{BA}}$ and $F_{\textrm{CNT}}$, the
distribution functions of BA and CNT at $(\bm{s}_{i},t)$, respectively, are
given by (up to the SPDE approximation),

$\displaystyle F_{\textrm{BA}}(x;\leavevmode\nobreak\
\mu_{Z}(\bm{s}_{i}),\mu_{1}(\bm{s}_{i}),\sigma^{2}_{1}(\bm{s}_{i}))$
$\displaystyle\approx$ $\displaystyle 1-\Phi(\mu_{Z}(\bm{s}_{i}))$
$\displaystyle+$
$\displaystyle\Phi(\mu_{Z}(\bm{s}_{i}))F_{LN}(x;\mu_{1}(\bm{s}_{i}),\sigma^{2}_{1}(\bm{s}_{i})),$
$\displaystyle F_{\textrm{CNT}}(x;\leavevmode\nobreak\
\mu_{Z}(\bm{s}_{i}),\mu_{2}(\bm{s}_{i}),\sigma^{2}_{2}(\bm{s}_{i}))$
$\displaystyle\approx$ $\displaystyle 1-\Phi(\mu_{Z}(\bm{s}_{i}))$
$\displaystyle+$
$\displaystyle\Phi(\mu_{Z}(\bm{s}_{i}))F_{LN}(x;\mu_{2}(\bm{s}_{i}),\sigma^{2}_{2}(\bm{s}_{i})),x\geq
0.$

###  3.4 Stage 4: Random Forests (RFs)

Until Stage 3, only the statistical approaches for joint modeling of BA and
CNT, were considered. In the joint modeling, the stage-specific models did not
incorporate any covariate information (except in certain priors that are
discussed later). As reported in Section 2, incorporating covariates in the
statistical modeling framework and allowing for spatially-varying regression
coefficients is computationally challenging and do not significantly improve
the prediction performance. Moreover, assuming that the marginal distributions
of log-CNT are Gaussian is computationally beneficial; however, it is
inaccurate because of the discrete nature of CNT. Thus, in Stage 4, a
rectification of the predicted values of CNT obtained from Stages 1 through 3
is proposed. Typically, machine learning (ML) approaches (such as RFs) do not
have any distributional assumptions and involve low computational burden,
while allowing a straightforward incorporation of covariate information.

As reported, the prediction performance is not improved using a simple
regression setting. Nevertheless, for interpretation and prediction purposes,
it is important to incorporate meaningful covariates, such as some landcover
types and climate/weather conditions, because they play a key role in the
occurrence of wildfires, as shown in several studies (Fusco et al.,, 2019;
Nadeem et al.,, 2020). These covariates are interrelated and show seasonal
dependence, but classical statistical methods, including generalized linear
models, cannot account for the sophistication in this process. Jain et al.,
(2020) reviewed extensive literature (until 2019) on ML application in
wildfire science and management, in which they reported RFs to be a more
common method for predicting fire occurrences since 2012, before models based
on artificial neural networks (NNs) and support vector machines (SVMs) were
widely reported in the literature. It is reasonable to use RFs in the
rectification of the US wildfire data analysis, because tree-based models are
often suitable for classification problems, hence they are suitable for the
discrete data CNT, and RF generates nonlinear regression models, which
facilitates correct formulations of covariates.

The RF algorithm (Breiman,, 2001) involves an ensemble of many decision trees,
in which individual trees are trained based on a random subset of the data,
and are drawn with replacement. A random subset (resampled with replacement)
of covariates is selected at each node of every decision tree. For
classification problems, each individual tree forms a class; furthermore, the
predicted class is selected by the highest votes. The success of RFs are
attributed to the low correlation among trees (each tree is trained
independently from others); thus, it results in low prediction variance.
Importantly, RFs minimize the correlation between trees; hence it provides
higher accuracy than compared to individual trees, which explains the success
of classification and regression trees (CART).

For the EVA 2021 data challenge, the prediction performance is evaluated based
on score functions with high weights to the data categorized into higher
severity levels as described in Section 2. Considering the size of the dataset
and underlying computational cost, this specific task (i.e., that of the
challenge) was converted in classification problem, in which CNT was
classified into 29 categories, (based on the 28 severity thresholds provided
for model evaluation in 2.1) and labeled at each spatiotemporal location. Then
CNT is converted into categorical data, which inevitably decreased in
resolution and essentially did not affect the ability to achieve high
prediction accuracy. Moreover, BA is here included as a covariate, in which
the missing values of BA are imputed by the joint statistical modeling of BA
and CNT.

###  3.5 Modeling based on log-Gaussian Cox processes (LGCPs): another
competing approach for CNT

Log-Gaussian Cox processes (LGCPs), or Cox process, are a commonly flexible
approach for analyzing point pattern data and are obtained by assuming a
hierarchical Bayesian structure. At the first level of this structure (data
level), the response at every discretized location is assumed to follow a
Poisson distribution conditioned on its random intensity measure; at the
second level (process level), the intensity measure is assumed to follow a
log-Gaussian process defined over the spatial domain. In particular, a spatial
point process $Y(\cdot)$ defined over $\mathcal{S}\subset\mathbb{R}^{2}$ is
called a Cox process if
$Y(\bm{s})|\Lambda(\bm{s})\overset{\textrm{Indep}}{\sim}\textrm{Poisson}(\exp[\Lambda(\bm{s})])$
and the log-intensity process $\\{\Lambda(\bm{s});\bm{s}\in\mathcal{S}\\}$ is
a GP.

To fit the LGCP model in a tractable manner, a common approach is to divide
the whole spatial region in a uniform grid and assume a constant log-intensity
in each grid cell. Here, the US wildfire dataset is already gridded, and
$\\{\bm{s}_{i},i=1,\ldots,N\\}$ denote the spatial locations (the centroids of
the grid cells) as defined in Section 3.1. The log-intensity vector of the
$t$-th time point is denoted
$\bm{\Lambda}_{t}=[\Lambda_{t}(\bm{s}_{1}),\ldots,\Lambda_{t}(\bm{s}_{N})]^{\prime}$.
We have
$\bm{\Lambda}_{t}\sim\textrm{Normal}_{N}(\bm{X}_{\Lambda}\bm{\beta}_{\Lambda},\sigma_{\Lambda}^{2}\bm{C}_{\Lambda})$,
where $\bm{X}_{\Lambda}$ is the design matrix, $\bm{\beta}_{\Lambda}$ is the
vector of regression coefficients, $\sigma_{\Lambda}^{2}$ is the common
variance term, and $\bm{C}_{\Lambda}$ is the corresponding correlation matrix.
For a dense correlation matrix $\bm{C}_{\Lambda}$, the estimation procedure
involves a high computational burden and thus, similar to Section 3.1, we use
an SPDE approximation-based construction of $\bm{C}_{\Lambda}$. We denote the
CNT data at spatial location $\bm{s}_{i}$ and time $t$ by
$\textrm{CNT}_{t}(\bm{s}_{i})$ and propose the following LGCP model:

$\displaystyle\textrm{CNT}_{t}(\bm{s}_{i})\mid\Lambda_{t}(\bm{s}_{i})$
$\displaystyle\overset{\textrm{Indep}}{\sim}$
$\displaystyle\text{Poisson}\left(\exp[\Lambda_{t}(\bm{s}_{i})]\right),\quad
i=1,\ldots,N,\,t=1,\ldots T,$ $\displaystyle\Lambda_{t}(\bm{s}_{i})$
$\displaystyle=$ $\displaystyle\mu_{t}(\bm{s}_{i})+\zeta_{t}(\bm{s}_{i}),$
(14)

where $\zeta_{t}(\cdot),t=1,\ldots,T$, are IID copies from a GMRF that has an
(approximately) equivalent isotropic Matérn spatial correlation given by (3)
(see Section 3.1 for more details), with $\phi_{\varepsilon}$ and
$r_{\varepsilon}$ replaced by $\phi_{\zeta}$ and $r_{\zeta}$, respectively,
and $\mu_{t}(\bm{s}_{i})$ denotes the spatiotemporal mean process at spatial
location $\bm{s}_{i}$ and time $t$ and is defined in terms of fixed
covariates. Multiple combinations of covariates (i.e., all covariates, only
spatial covariates, including/excluding BA) are compared and discussed in
Section 5.

##  4 Computation

###  4.1 Computational details for Stages 1 and 3

Inferential statistical analysis is conducted on the model parameters
mentioned in Section 3.1 and 3.3 based on MCMC sampling. Conjugate priors are
selected whenever possible. The full posterior distributions of the model
parameters and hyperparameters are provided in the Supplementary Materials. We
here briefly outline the MCMC steps. In Stage 2, a fixed rank kriging model
was directly fitted using the R package FRK; hence, this was skipped here.

In Stage 1, the parameters and hyperparameters are
$\bm{\Theta}_{1}=\left\\{\\{\bm{X}_{t}\\}_{t=1}^{T},\bm{\mu}_{Z},\bm{\theta}_{\mu},\tau_{\mu},\\{\bm{\varepsilon}^{*}_{t}\\}_{t=1}^{T},\phi_{\varepsilon},r_{\varepsilon}\right\\}$.
By an abuse of notation, we reparametrize
$\sqrt{r_{\varepsilon}}\bm{\varepsilon}^{*}_{t}$ by
$\bm{\varepsilon}^{*}_{t}$. The full posterior distribution of the latent
variables $X_{t}(\bm{s}_{i})$ depends on $Z_{t}(\bm{s}_{i})$. If
$Z_{t}(\bm{s}_{i})$ is missing, the posterior distribution of
$X_{t}(\bm{s}_{i})$ is normal; otherwise, if $Z_{t}(\bm{s}_{i})$ is zero or
one, the posterior of $X_{t}(\bm{s}_{i})$ is a truncated normal distribution,
supported on the negative or positive side of the real line, respectively. The
prior distribution we choose for $\bm{\mu}_{Z}$ is
$\bm{\mu}_{Z}\sim\textrm{Normal}_{N}(\bm{D}\bm{\theta}_{\mu},\tau^{-1}_{\mu}\bm{I}_{N})$.
Here, $\bm{D}$ is a $(N\times 6)$-dimensional design matrix with its columns
representing an intercept term, longitude, latitude, mean altitude, standard
deviation of altitude, and the proportion of a pixel that is within the US
Mainland. Because the aim of this study was to predict the underlying spatial
process at a new set of locations, it is reasonable to assume an unstructured
covariance for the prior of $\bm{\mu}_{Z}$. For the hyperparameters
$\bm{\theta}_{\mu}$ and $\tau_{\mu}$, we choose weakly-informative conjugate
priors $\bm{\theta}_{\mu}\sim\textrm{Normal}_{6}(\bm{0},10^{2}\bm{I}_{6})$ and
$\tau_{\mu}\sim\textrm{Gamma}(0.1,0.1)$. The unconditional distribution of
$\bm{\varepsilon}^{*}_{t}$ is
$\bm{\varepsilon}^{*}_{t}\sim\textrm{Normal}_{N^{*}}(\bm{0},r_{\varepsilon}\bm{Q}_{\phi_{\varepsilon}}^{-1})$
and the conditional distribution of $\bm{X}_{t}$ given
$\bm{\varepsilon}^{*}_{t}$ is
$\bm{X}_{t}|\bm{\varepsilon}^{*}_{t}\sim\textrm{Normal}_{N}(\bm{\mu}_{Z}+\bm{A}\bm{\varepsilon}^{*}_{t},(1-r_{\varepsilon})\bm{I}_{N})$.
Thus, the full conditional posterior distribution of
$\bm{\varepsilon}^{*}_{t}$ is an $N^{*}$-variate normal distribution;
furthermore, the calculation of the mean vector and covariance matrix is
straightforward. Based on the remaining parameters and hyperparameters,
$\bm{\varepsilon}^{*}_{t};t=1,\ldots,T$, are conditionally independent and
thus are updated in parallel. In case of the parameters $\phi_{\varepsilon}$
and $r_{\varepsilon}$, any existence of conjugate priors is not known and
hence independent priors are selected as
$\phi_{\varepsilon}\sim\textrm{Uniform}(0,2\Delta_{\mathcal{S}})$ and
$r_{\varepsilon}\sim\textrm{Uniform}(0,1)$, where $\Delta_{\mathcal{S}}$ is
the largest Euclidean distance between two data locations. The posterior
samples from $\phi_{\varepsilon}$ and $r_{\varepsilon}$ are drawn using the
well-established Metropolis-Hastings (M-H) algorithm.

In Stage 3, the parameters and latent variables are
$\bm{\Theta}_{2}=\left\\{\\{\bm{\eta}^{*}_{t}\\}_{t=1}^{T},\phi_{\eta},r_{\eta},\rho_{\eta}\right\\}$.
The unconditional distribution of $\bm{\eta}^{*}_{t}$ is a $2N^{*}$-variate
normal distribution; furthermore, the conditional distribution of
$\widehat{\bm{W}}_{t}$ given $\bm{\eta}^{*}_{t}$ is a $2N$-variate normal
distribution, (the detailed expressions are presented in Section 3.3). The
calculation of the full conditional posterior distribution of
$\bm{\eta}^{*}_{t}$ is straightforward and it is again a $2N^{*}$-variate
normal distribution. When the remaining parameters and hyperparameters,
$\bm{\eta}^{*}_{t};t=1,\ldots,T$, are conditionally independent and thus are
updated in parallel. For the parameters $\phi_{\eta}$, $r_{\eta}$, and
$\rho_{\eta}$, any existence of conjugate priors is not known and we choose
independent priors $\phi_{\eta}\sim\textrm{Uniform}(0,2\Delta_{\mathcal{S}})$,
$r_{\eta}\sim\textrm{Uniform}(0,1)$, and
$\rho_{\eta}\sim\textrm{Uniform}(0,1)$. Posterior samples from $\phi_{\eta}$,
$r_{\eta}$, and $\rho_{\eta}$ are drawn using an M-H algorithm.

Each MCMC chain was run for 60,000 iterations with the first 10,000 iterations
discarded as burn-in. The post-burn-in samples were then thinned maintining
one in each five samples. Convergence of the chains was monitored via trace
plots. The computation of Stage 1 and Stage 3 were undertaken on a desktop
with an Intel Xeon CPU E5-2680, a 2.40GHz processor and 128GB RAM , and the
corresponding computational times for Stage 1 and Stage 3 were 447 minutes and
744 minutes respectively. Note that these two stages can be run in parallel.

###  4.2 Computational details for Stage 4

We used the function randomForest from the R package randomForest, which
implements Breiman’s random forest algorithm. The constructed RF models are
trained and their performances tested using a cross-validation study detailed
in Section 5.1. For classification problems, the optimal number of covariates
used at each splitting node is $\sqrt{P}$, where $P$ is the total number of
covariates. However, multiple models were fitted using different number of
covariates, and the model with the best prediction performance was selected
based on a cross-validation scheme. The running time on a computer with the
same configuration described in Section 4.1, was approximately 20 minutes.
After some trial-and-error, the tuning parameter configurations were finally
chosen as follows: $mtry=36$ (all available covariates), where $mtry$ is the
number of covariates randomly sampled as candidates for each split; and
$ntree=200$, where $ntree$ denotes the number of trees to grow. Note that
$ntree$ should not be extremely small to ensure that every input row is
predicted a few times. These two tuning-parameters are selected based on some
exploratory experiments and a cross-validation study.

###  4.3 Computational details for fitting the LGCP model

Posterior inference is obtained from the LGCP model (3.5) based on a
stochastic gradient-based MCMC method (For more details refer to Welling and
Teh, (2011) and Algorithm 1 in Yadav et al., (2021)). In LGCP models of the
form (3.5), the set of parameters, hyperparameters, and latent variables is
given by
$\bm{\Theta}_{\textrm{LGCP}}=\left\\{\\{\bm{\Lambda}_{t}\\}_{t=1}^{T},\bm{\beta},\\{\bm{\zeta}^{*}_{t}\\}_{t=1}^{T},\phi_{\zeta},r_{\zeta}\right\\}$,
where $\bm{\beta}$ is the vector of regression coefficients associated with
the known covariates, and
$\left\\{\\{\bm{\zeta}^{*}_{t}\\}_{t=1}^{T},\phi_{\zeta},r_{\zeta}\right\\}$
has a similar interpretation as
$\left\\{\\{\bm{\varepsilon}^{*}_{t}\\}_{t=1}^{T},\phi_{\varepsilon},r_{\varepsilon}\right\\}$
in Section 3.1. We use conjugate priors whenever possible, and these
parameters are updated using Gibbs sampling. For $\bm{\beta}$, we choose a
weakly-informative conjugate Gaussian prior with mean zero and variance 100
when a small number of covariates are in the model, and for the case of a
large number of covariates, we use an informative Gaussian prior with mean 0
and variance equal to 0.1 that allows appropriate penalization. For the Matérn
correlation parameters $\phi_{\zeta}$ and $r_{\zeta}$, we use the same priors
as for $\phi_{\varepsilon}$ and $r_{\varepsilon}$, respectively. We update
them within MCMC using a M-H algorithm, similar to updating
$\phi_{\varepsilon}$ and $r_{\varepsilon}$. For the latent vectors
$\\{\bm{\Lambda}_{t}\\}_{t=1}^{T}$, we do not have closed-from posteriors, and
thus we update them jointly using the stochastic gradient Langevin dynamics,
which is similar to Algorithm 1 in Yadav et al., (2021). For the latent
vectors $\\{\bm{\zeta}^{*}_{t}\\}_{t=1}^{T}$, we have closed-form full
posteriors, and thus they are updated using Gibbs sampling. For the stochastic
gradient MCMC algorithm, we set the batch size to be $b=10$ (i.e., we update
10 out of the $T$ vectors $\bm{\Lambda}_{t}$ at a time). Because the remaining
parameters and hyperparameters, $\bm{\zeta}^{*}_{t};t=1,\ldots,T$, are
conditionally independent, we update them in parallel.

The MCMC chains were run for a total of 250,000 iterations, and the first
200,000 samples were discarded as burn-in samples. The chains were thinned by
keeping one for each 25 samples. Thus, all the summary statistics were
calculated based on the final 2000 samples. The computation is approximately
37 hours, when only spatial covariates are included in the models, and
approximately 110 hours when we use all the covariates (both spatial and
spatio-temporal).

##  5 Data Application

###  5.1 Cross-validation schemes and model comparison

The prediction performance of the four-stage model proposed in Sections 3.1 to
3.4 is compared to a few sub-models, in addition to the alternative LGCP model
described in Section 3.5. For comparison, a cross-validation study was
performed by dividing the available data after masking into training and test
sets. To ensure that the used cross-validation scheme is compatible to that
used by the EVA 2021 data challenge organizers, we attempt to replicate the
original missingness pattern in the test set, and we choose the same metric of
model comparison as that used for the data challenge described in Section 2.1.

As reported in Section 2.1, the observations are available for 23 years
(1993–2015), from March to September. Thus, data for a total of 161 months
were available. Of the 23 years, complete observations for 12 years (84
months) are available, and a significant proportion of the data are missing
(NA) for the remaining of 11 years (77 months). Here, two types of test sets
were developed by replicating the spatial patterns of the missing observations
for 77 months out of 84 months (chosen randomly) with complete data. The first
type is called ‘Fixed Month’ scheme, where missingness is created for a month
using the pattern from the same month index; e.g., for March 1993 (has no
missing data), the month March is randomly selected for one out of 11 years
where certain data are missing, e.g. 1994, and then the missingness pattern
for March 1993 is set to be the same as that for March 1994. The second type
is called ‘Random Month’ scheme where missingness for a month is created using
the pattern from any randomly selected month with certain observations
missing. Figure 5 shows the process of generating the new test set for CNT
using the principle of the ‘Fixed Month’ scheme. The original test set created
by the organizers contains 80,000 spatiotemporal locations, and we masked data
at additional 80,000 observations. These two types of cross-validation schemes
are primarily selected to confirm the validity of the prediction performance
of the proposed models while accounting for seasonality and clustering of
masked data in space.

Figure 5: A graphical schematic of the cross-validation schemes. The top
panel, we illustrates binary time series (across the X-axis) for the first 50
pixels (across the Y-axis), where data availability and missingness are
presented in black and white, respectively. The pixels within two consecutive
white dashed lines correspond to seven months within a year. The 14th column
(September 1994), which highlighted within the red vertical lines, is used in
the bottom panel for illustration. The bottom panel, shows the binary time
series for data availability and missingness, after masking additional 80,000
spatiotemporal cases under the ‘Fixed month’ scheme. Here, column 133
(September 2011) has the same missingness pattern as column 14.

For BA prediction, five competing models, closely linked to the final proposed
model, are compared. These include the benchmark, which is a linear regression
model for BA on a logarithmic scale, with all covariates including all filled
CNT (achieved thanks to a Poisson regression). The spatial prediction
performances, of the five models are compared using the same metric employed
by the data challenge organizers. For the first competitor, a scenario where
there is no data available for CNT, is assumed; thus, the zero/nonzero
indicator created in Stage 1 is only based on the BA information.
Consequently, the original 80,000 missing observations (rather than 48,947
cases, as mentioned in Section 2.1) in Stage 1 are obtained. Moreover, rather
than a bivariate modeling of BA and CNT, a similar univariate SPDE-based
spatial model is used for BA only. For other competing models, the CNT data is
assumed to be available. For the second competitor, CNT data are used only for
the zero/nonzero indicator part in Stage 1, for which data of only 48,947
spatiotemporal locations are missing. However, for the positive part of BA,
the same univariate SPDE-based spatial model is fitted as in the case of the
first competitor.

For the third competitor, all the 35 spatiotemporal, purely spatial, or purely
temporal covariates mentioned in Section 2.1 are used in Stage 1. However, the
same univariate SPDE-based model in Stage 3 as for the first two competitors,
is fitted. In Stage 1, we replace (2) as follows

$\displaystyle
Z_{t}(\bm{s}_{i})=\begin{cases}1,&\parbox[t]{95.39502pt}{if\leavevmode\nobreak\
\leavevmode\nobreak\ $X_{t}(\bm{s}_{i})>0$}\\\ 0,&\text{if\leavevmode\nobreak\
\leavevmode\nobreak\
$X_{t}(\bm{s}_{i})<0$},\end{cases}\textrm{where}\leavevmode\nobreak\
X_{t}(\bm{s}_{i})=\alpha(\bm{s}_{i})+\sum_{l=1}^{35}\gamma_{l}D_{l,t}(\bm{s}_{i})+\varepsilon_{t}(\bm{s}_{i}),$

where $\alpha(\cdot)$ is a spatially-varying intercept term,
$D_{l,t}(\bm{s}_{i}),l=1,\ldots,35$, are the 35 covariates,
$\gamma_{l},l=1,\ldots,35$, are the corresponding spatially/temporally
invariant regression coefficients, and $\varepsilon_{t}(\bm{s}_{i})$ are the
same as in (2). Here,
$\textrm{Pr}(Z_{t}(\bm{s}_{i})=1)\approx\Phi(\alpha(\bm{s}_{i})+\sum_{l=1}^{35}\gamma_{l}D_{l,t}(\bm{s}_{i}))$.
Table 2 lists the evaluation scores for both cross-validation schemes. Under
both schemes, the final model described in Stages 1 to 3 exhibited better
performance than the alternative univariate spatial models and the benchmark.

Table 2: EVA 2021 data challenge evaluation scores for different models and
settings under the two proposed cross-validation schemes, i.e., Fixed Month
and Random Month. (A smaller value indicates better performance). Cross-
Validation for BA
---
Model | Fixed Month | Random Month
| Benchmark model (Log-Gaussian
---
regression, with filled CNT)
3468.85 | 3574.88
| Univariate spatial modeling,
---
with BA zero/nonzero indicator
2842.88 | 2972.47
| Univariate spatial modeling,
---
with BA/CNT zero/nonzero indicator
2842.64 | 2972.25
| Univariate spatial modeling,
---
with BA/CNT zero/nonzero
indicator and all covariates
3841.73 | 3971.32
| Bivariate spatial modeling,
---
with BA/CNT zero/nonzero indicator
2796.93 | 2923.08
Extra Space
Cross-Validation for CNT
Model | Fixed Month | Random Month
| Benchmark model (Poisson regression,
---
with log-link, excluding BA)
5663.02 | 5218.30
| LGCP with SPDE basis with
---
no covariates (only intercept)
4369.05 | 4822.24
| LGCP with SPDE basis and
---
purely spatial covariates
4358.78 | 4798.53
| LGCP with SPDE basis and
---
purely spatial covariates and BA
4469.38 | 4835.74
| LGCP with SPDE basis and
---
all covariates and BA
4423.26 | 4795.33
| Bivariate spatial modeling, with BA/CNT
---
zero/nonzero indicator, without Stage 4
3821.24 | 3570.91
| RF with mtry = 3 and ntree = 200
---
2784.01 | 2990.90
| RF with mtry = 6 and ntree = 200
---
2638.03 | 2827.02
| RF with mtry = 36 and ntree = 200
---
2612.01 | 2785.05

In CNT prediction, three types of competing models were selected with
different settings for each type, including the benchmark, which is a Poisson
regression model with all covariates except BA. The first type is a univariate
LGCP model constructed using a latent SPDE, as in (3.5). In the first setting,
it is assumed that the term $\mu_{t}(\bm{s}_{i})$ in (3.5) is constant across
space and time. In the second setting, $\mu_{t}(\bm{s}_{i})$ is assumed to be
only spatially-varying and it is written as a linear combination of purely
spatial covariates, including an intercept term (the columns of the design
matrix $\bm{D}$ in Section 4.1). In the third setting, $\mu_{t}(\bm{s}_{i})$
is assumed to be both spatially and temporally varying; it is written as a
linear combination of the purely spatial covariates, as well as BA, where the
missing values are filled using Stages 1 to 3. Finally, in the fourth setting,
$\mu_{t}(\bm{s}_{i})$ is written as a linear combination of all covariates and
filled BA data. The second model type is the purely statistical model
described in Stages 1 to Stage 3, without incorporating any rectification of
the erroneous Gaussian assumption for CNT. The third model type is based on
RFs, in which all the spatiotemporal covarites, as well as filled BA, are
used. Under this setting, we fix $ntree=200$ as mentioned in Section 4.2, and
set the number of covariates randomly sampled as candidates for each split
(mtry) at 3, 6, and 36, under three different settings. The final evaluation
scores for the two cross-validation schemes are presented in Table 2. Under
both cross-validation schemes, the final model described in Stages 1 to 4
performs better than the alternative univariate LGCP and the bivariate model
without the rectification using RF, as well as the benchmark.

###  5.2 Results

In this section, the posterior means and posterior standard deviations of the
model parameters and hyperparameters in Stages 1 and 3 are discussed. The
fixed rank kriging estimates of parameter surfaces $\mu_{1}(\cdot)$,
$\mu_{2}(\cdot)$, $\sigma_{1}(\cdot)$, and $\sigma_{2}(\cdot)$ in Stage 2, and
the prediction performance of the final model described in Stage 1 through
Stage 4 are also reported and discussed.

The mixing and convergence of MCMC chains in Stages 1 and 3 were evaluated
using trace plots. For scalar parameters, Figure 6 shows the trace plots of
the thinned MCMC chains. All plots exhibit good mixing and convergence
diagnostics. For other parameter vectors and latent variables, the convergence
and mixing are also confirmed though, not shown here. The computations in
Stages 2 and 4 were conducted using the R packages FRK and randomForest,
respectively.

Figure 6: Trace plots for some parameters and hyperparameters in Stages 1
(top) and 3 (bottom).

Figure 7 shows the spatial maps of the posterior mean and posterior standard
deviation of $\mu_{Z}(\cdot)$ in Stage 1. The values of the posterior mean of
$\mu_{Z}(\cdot)$ are higher near the Southeastern (the state of Georgia) and
Southwestern (the state of California) parts of the US and are generally lower
in the middle parts of the US. Among available temporal replications, the
posterior standard deviation is extremely high (more than 0.4) for 54 spatial
locations, mostly in the northeastern parts of the US, where the zero/nonzero
indicator is either always zero (for 53 locations) or always one (for 1
location).

Figure 7: Posterior mean and standard deviation profiles of $\mu_{Z}(\cdot)$.

Table 3 lists the posterior mean and standard deviation of non-spatial
parameters and hyperparameters in Stages 1 and 3. In a frequentist sense, all
parameters $\theta_{\mu,1}$ to $\theta_{\mu,6}$ are significant (the absolute
value of the ratio of posterior mean and posterior standard deviation is
larger than 2 for each of the six cases). The posterior means of
$\phi_{\varepsilon}$ and $r_{\varepsilon}$ are 3.0491 and 0.5319,
respectively, for the latent Gaussian process, which indicates a correlation
of 0.050 at a spatial distance of 10 degrees. In Stage 3, the posterior means
of $\phi_{\eta}$ and $r_{\eta}$ are 3.6408 and 0.3442, respectively, for each
component of the bivariate spatial Gaussian process, which indicates a
correlation of 0.052 at a spatial distance of 10 degrees. These results show
the requirement for modeling the spatial dependence in fire occurrences and
sizes. The posterior mean of $\rho_{\eta}$ is 0.4575, which shows a strong
positive correlation between BA and CNT.

Table 3: Posterior mean and standard deviation (SD) of the non-spatial parameters and hyperparameters in Stage 1 and Stage 3. Here, $\theta_{\mu,1},\ldots,\theta_{\mu,6}$ denote the six components of $\bm{\theta}_{\mu}$ that represent an intercept term and the regression coefficients of longitude, latitude, mean altitude, standard deviation of altitude, and the proportion of a pixel that is within the mainland US, respectively. | Stage 1 |
---|---|---
Parameter | Posterior mean | Posterior SD
$\theta_{\mu,1}$ | -0.3805 | 0.0262
$\theta_{\mu,2}$ | 0.0894 | 0.0239
$\theta_{\mu,3}$ | -0.2596 | 0.0192
$\theta_{\mu,4}$ | -0.2583 | 0.0280
$\theta_{\mu,5}$ | 0.4540 | 0.0289
$\theta_{\mu,6}$ | 0.1826 | 0.0154
$\tau_{\mu}$ | 1.5398 | 0.1276
$\phi_{\varepsilon}$ | 3.0491 | 0.9754
$r_{\varepsilon}$ | 0.5319 | 0.0194
| Stage 3 |
Parameter | Posterior mean | Posterior SD
$\phi_{\eta}$ | 3.6408 | 1.2179
$r_{\eta}$ | 0.3442 | 0.0042
$\rho_{\eta}$ | 0.4575 | 0.0022

Figure 8 shows the fixed rank kriging estimates of $\mu_{1}(\cdot)$,
$\mu_{2}(\cdot)$, $\sigma_{1}(\cdot)$, and $\sigma_{2}(\cdot)$, in Stage 2.
For $\mu_{1}(\cdot)$, the FRK estimates are generally higher in a large
portion of the southeastern US, while the values are lower in the mid-west,
north-west and the north-east. For $\mu_{2}(\cdot)$, the values are higher in
two small regions in the southeast and southwest of the US. In some regions
(mid-north), the estimated profile is not highly smooth indicating high local
nonstationarity. For $\sigma_{1}(\cdot)$, large estimates are visible in a
small portion of the mid-US, while moderately large values are observed in a
large region in the western US. For $\sigma_{2}(\cdot)$, large values are
observed in two large regions of the eastern and western US. In the mid-US,
the estimated profile is not smooth in some regions, particularly, near the
mid-north of the US, similar to the $\mu_{2}(\cdot)$ profile. The spatial maps
of $\sigma_{1}(\cdot)$ and $\sigma_{2}(\cdot)$ illustrate the underlying
spatial heteroscedasticity and sharp local variability.

Figure 8: Approximate Bayesian estimates of $\mu_{1}(\cdot)$,
$\mu_{2}(\cdot)$, $\sigma_{1}(\cdot)$, and $\sigma_{2}(\cdot)$ using the fixed
rank kriging procedure described in Section 3.2.

After incorporating the Stage 2 estimates, we obtain the residuals
$\widehat{W}_{tp}(\bm{s}_{i})$ (Section 3.3). Despite fitting a Gaussian model
to log-BA and log-CNT in Stage 3 for computational suitability, we study the
histograms of the standardized log-BA and log-CNT values to check the validity
of our Gaussian assumption, and we present them in Figure 9. For the
standardized log-BA, the histogram is bell-shaped and symmetric around zero
and the Gaussian assumption appears to be reasonable. However, for the
standardized log-CNT, the histogram is bimodal and the shapes of the two modes
are considerably different. Furthermore, the histogram appears to be right-
skewed; thus, a Gaussian assumption for $\widehat{W}_{tp}(\bm{s}_{i})$ in
Stage 3 for CNT is questionable and this justifies the rectification using RFs
in Stage 4.

Figure 9: Histograms of the standardizing $\log$-BA and $\log$-CNT based on
the Stage 2 estimates.

In Stage 4, the RF algorithm, is used. In this algorithm, the final model
based on cross-validation includes all the covariates and BA (missing values
are imputed using Stages 1–3), as described in Section 5.1.

Figure 10 shows the variable importance plot (or mean decrease accuracy plot).
The plot expresses how much accuracy the model looses by excluding each
variable at a time during classification. The mean decrease in Gini
coefficient is a measure of how each variable contributes to the homogeneity
of nodes and leaves in the resulting RF. The increase in the value mean
reduction in accuracy or mean reduction Gini score, increases the importance
of the variable in the model. Here, BA is reported to be the most important
covariate, followed by clim5 (Evaporation) and clim1 (wind speed in Eastern
direction), whereas lc6 (tree broadleaved deciduous closed) is the least
important covariate. The covariates were included into the models in the order
of their variable importance.

Figure 10: Variable importance (or mean decrease accuracy plot) for the RF
algorithm in Stage 4, in which all the available covariates and filled BA were
used for classification.

Finally, we compare the estimated predictive distribution functions for the
80,000 test observations, for each BA and CNT, evaluated at 28 levels, with
the empirical distribution functions calculated based on the test
observations. Figure 11 shows boxplots of the absolute differences between the
two cumulative distribution functions (CDFs) at 28 evaluation levels,
specified in (2.1). At lower levels certain boxes and upper endpoints of error
bars are significantly different from zero; however, the median of the
absolute differences is close to zero in most cases, particularly for BA. For
higher levels, which are assigned larger weights for model evaluation, the
boxes are close to zero indicating that the proposed model performs well for
predicting the distribution of BA and CNT at unobserved spatiotemporal
locations.

Figure 11: Boxplots of absolute differences between the estimated predictive
distribution functions and empirical distribution functions evaluated at 28
levels given in (2.1).

##  6 Discussions and Conclusions

###  6.1 Summary

Motivated by the EVA 2021 data challenge, in which the first four authors of
this paper participated as the The Bedouins team, a four-stage high-
dimensional zero-inflated bivariate spatial model based on statistics and
machine learning was proposed for the prediction of BA and CNT at masked
spatiotemporal locations. Here, a spatial dependence structure was developed
using SPDEs, which reduces the computational burden by allowing sparsity in
the precision matrices. In Stage 1, the data were categorized into
zero/nonzero categories and a two-layered Bayesian hierarchical model was
fitted for estimating the probabilities of the two categories at unobserved
sites. In Stage 2, the parameter surfaces were estimated using two-step
approximate Bayesian inference technique that bypasses high computational
burden.

In Stage 3, the standardized log-transformed positive observations from Stage
2 were modeled using a bivariate spatial GMRF. The log-Gaussian assumption for
modeling positive wildfire frequencies was computationally helpful; however,
it was erroneous because the observations are discrete-valued. Thus, in Stage
4, the predicted values of wildfire frequencies were rectified using RFs, in
which BA was treated as a covariate, and the missing BA values were imputed by
the predicted BA values in Stage 3. MCMC sampling was used to draw posterior
inference in Stages 1 and 3. The computation in Stages 2 and 4 were done
directly using the R packages FRK and randomForest, respectively. Our final
model was shown to outperform some alternatives in a well-designed cross-
validation study and to effectively predict low to high quantiles of BA and
CNT at unobserved sites. Table 4 lists the final scores for the three best
participating teams. While the teams BlackBox and Kohrrelation focused on
purely ML-based approaches, our method only depends on ML algorithms for
rectifying predictions given by statistical models.

A cross-validation scheme was created to effectively compare the models. For
modeling discrete spatial data, an alternative approach using LGCPs was
discussed. The fitted LGCP was constructed using a latent SPDE spatial effect;
however, in case of spatial prediction of the US wildfire data, our proposed
model outperforms the LGCP approach. While fully machine learning algorithms
do not have any distributional assumption and are thus more robust, our
proposed statistical method can better quantify the underlying uncertainty,
specifically in Stages 1 and 3.

Table 4: Final evaluation scores for the three best performing teams in the EVA 2021 data challenge and the benchmark. (A smaller score indicates better performance). | BlackBox | Kohrrelation | The Bedouins | Benchmark
---|---|---|---|---
BA | 3315.65 | 3446.02 | 3408.31 | 4244.36
CNT | 2804.95 | 2989.85 | 3145.81 | 5565.15
Total | 6120.60 | 6435.87 | 6554.12 | 9809.51

###  6.2 Drawbacks and possible solutions

Despite its elegant performance in predicting the distribution functions at
masked sites, our approach also has some limitations. First, we selected a
multistage approach in which information was not borrowed from one stage to
another. For example, in Stage 1, the data are treated as zero/nonzero
indicators; however, in other stages, zeros are treated as missing data. Joint
modeling of zero and nonzero data is a possible solution. Second, in Stages 1
and 3, it was assumed that the underlying dependence structures are of the
isotropic Matérn form given in (3). Because the spatial domain is large, it
would be reasonable to assume a nonstationary dependence structure, e.g.,
using empirical orthogonal functions (Wikle,, 2010; Hazra and Huser,, 2021),
that are more realistic for modeling the wildfire process over the entire US
Mainland. In Stage 2, we obtained smoothed parameter estimates in two steps;
in the first step, we calculated empirical estimates, and in the second step,
we treated such estimates as spatial data. Moving from Step 1 to Step 2, we
ignored the uncertainty of the estimates in the first step. While the
posterior coverage probabilities are only mildly affected by such a choice as
shown in Hazra et al., (2019), a superior method for borrowing the uncertainty
information from the first step to the second step has been recently discussed
in Hrafnkelsson et al., (2021). In Stage 4, we used the univariate random
forest algorithm of Breiman, (2001) independently at each site. However,
recently, a random forest algorithm for spatially-dependent data has been
proposed by Saha et al., (2021). Because we ignored covariate information in
Stage 1 to Stage 3, the proposed model is not suitable to draw inference about
the significance of a specific predictor in the context of fire modeling. A
possible solution would be to fit a regression model for the mean surfaces.

###  6.3 Other applications

The model proposed here is a general tool for all bivariate zero-inflated
spatiotemporal datasets, although our methodology has been motivated by the
joint analysis of fire occurrences and sizes. The model supplements spatial
marked point processes; thus, it is suitable for multiple scenarios. For
example, Penttinen et al., (1992) discussed a marked point process approach
for forest statistics, where ‘points’ are the tree positions, whereas the
‘marks’ are tree characteristics such as stem diameters and tree species.
After summarizing the data over a grid, we can obtain a bivariate spatial
dataset with two components representing the number of trees and average stem
diameter within each grid cell. In such a scenario, our proposed inference
approach would be suitable. Another example is the joint modeling of the
number of rainy days within a month and the total monthly precipitation within
a pixel/region. In addition to the marked point process scenario, our model
can also be used to fit temperature and rainfall data over a spatiotemporal
domain (Gelfand et al.,, 2005), for example. Temperature is a real-valued
continuous variable for which a Gaussian assumption is common (and thus a log-
transformation is not necessary), and rainfall data are non-negative and zeros
for the dry periods; here, we can ignore Stage 4 if the log-transformed
nonzero precipitation amounts follow a normal distribution.

###  6.4 Possible extensions

The method proposed in Sections 3.1 and 3.4 can be extended in multiple
directions. The distribution function estimation problem for the EVA 2021 data
challenge can be treated as a classification problem. Thus, by categorizing
the data based on the 28 levels, a set of spatial indicators can be obtained,
as described in Agarwal et al., (2021). Furthermore, a multivariate version of
Stage 1 would be a reasonable model for drawing inference. The proposed model
assumes temporal independence; thus, a natural extension would be to assume
that the latent processes have temporal autocorrelation. By replacing the mean
profiles as functions of the covariates, we can extend the proposed model to
study the significance of a specific predictor in the context of fire
modeling. The underlying dependence structures are assumed to be stationary
for Stage 1 and Stage 3; however, some approaches available in the literature
allow nonstationary spatial modeling for large datasets (Katzfuss,, 2013;
Banerjee,, 2020), and we can extend the proposed model in that direction.
Although the main focus in the data challenge is the accurate estimation of
the upper tail of CNT and BA observations, we build our model using Gaussian
processes that have been criticized for modeling spatial extremes (Davison et
al.,, 2013). While the extreme-value theory-justified models typically entail
a huge computational burden even in low dimensions, a simpler alternative
would be to consider scale mixture models (Huser et al.,, 2017; Huser and
Wadsworth,, 2020; Hazra et al.,, 2021). Although we used an adaptive
Metropolis-Hastings algorithm for the computations in Stages 1 and 3, we can
extend it using some faster and more recently developed algorithms, like the
stochastic gradient-based (Welling and Teh,, 2011) or deterministic
transformation-based (Dutta and Bhattacharya,, 2014) MCMC algorithms.

## Acknowledgments

The first three authors (Cisneros, Gong, Yadav) contributed equally to this
work by implementing some of the methods and writing parts of the paper. The
last two authors (Hazra, Huser) oversaw the whole project, with Hazra having a
leading role throughout all practical aspects of the data competition
(supervision of the Bedouins Team, methods’ implementation, results’
interpretation, writing).

We would like to thank Thomas Opitz for organizing this very interesting data
competition for the EVA 2021 Conference, as well as Thomas Mikosch for
welcoming a Special Issue in Extremes about this topic. This publication is
based upon work supported by the King Abdullah University of Science and
Technology (KAUST) Office of Sponsored Research (OSR) under Award No. OSR-
CRG2020-4394.

## References

  * Abatzoglou and Williams, (2016) Abatzoglou, J. T. and Williams, A. P. (2016).  Impact of anthropogenic climate change on wildfire across western US forests.  Proceedings of the National Academy of Sciences, 113(42):11770–11775.
  * Abdelfatah et al., (2016) Abdelfatah, K., Bao, J., and Terejanu, G. (2016).  Environmental modeling framework using stacked Gaussian processes.  arXiv preprint:1612.02897.
  * Agarwal et al., (2021) Agarwal, G., Sun, Y., and Wang, H. J. (2021).  Copula-based multiple indicator kriging for non-Gaussian random fields.  Spatial Statistics, page 100524.
  * Bakka et al., (2018) Bakka, H., Rue, H., Fuglstad, G.-A., Riebler, A., Bolin, D., Illian, J., Krainski, E., Simpson, D., and Lindgren, F. (2018).  Spatial modeling with R-INLA: A review.  Wiley Interdisciplinary Reviews: Computational Statistics, 10(6):e1443.
  * Banerjee, (2020) Banerjee, S. (2020).  Modeling massive spatial datasets using a conjugate bayesian linear modeling framework.  Spatial Statistics, 37:100417.
  * Bivand et al., (2015) Bivand, R., Gómez-Rubio, V., and Rue, H. (2015).  Spatial data analysis with R-INLA with some extensions.  Journal of Statistical Software, 63(20):1–31.
  * Breiman, (2001) Breiman, L. (2001).  Random Forests.  Machine Learning, 45(1):5–32.
  * Brown et al., (2021) Brown, E. K., Wang, J., and Feng, Y. (2021).  US wildfire potential: a historical view and future projection using high-resolution climate data.  Environmental Research Letters, 16(3):034060.
  * Ciarlet, (2002) Ciarlet, P. (2002).  The Finite Element Method for Elliptic Problems.  Society for Industrial and Applied Mathematics, SIAM.
  * Cressie and Johannesson, (2008) Cressie, N. and Johannesson, G. (2008).  Fixed rank kriging for very large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(1):209–226.
  * Cumming, (2001) Cumming, S. (2001).  A parametric model of the fire-size distribution.  Canadian Journal of Forest Research, 31(8):1297–1303.
  * Davison et al., (2013) Davison, A. C., Huser, R., and Thibaud, E. (2013).  Geostatistics of dependent and asymptotically independent extremes.  Mathematical Geosciences, 45(5):511–529.
  * Diggle et al., (2013) Diggle, P. J., Moraga, P., Rowlingson, B., and Taylor, B. M. (2013).  Spatial and spatio-temporal log-Gaussian Cox processes: extending the geostatistical paradigm.  Statistical Science, 28(4):542–563.
  * Dutta and Bhattacharya, (2014) Dutta, S. and Bhattacharya, S. (2014).  Markov chain Monte Carlo based on deterministic transformations.  Statistical Methodology, 16:100–116.
  * Fusco et al., (2019) Fusco, E. J., Finn, J. T., Balch, J. K., Nagy, R. C., and Bradley, B. A. (2019).  Invasive grasses increase fire occurrence and frequency across US ecoregions.  Proceedings of the National Academy of Sciences, 116(47):23594–23599.
  * Gabriel et al., (2017) Gabriel, E., Opitz, T., and Bonneu, F. (2017).  Detecting and modeling multi-scale space-time structures: the case of wildfire occurrences.  Journal de la Société Française de Statistique, 158(3):86–105.
  * Gelfand et al., (2005) Gelfand, A. E., Banerjee, S., and Gamerman, D. (2005).  Spatial process modelling for univariate and multivariate dynamic spatial data.  Environmetrics, 16(5):465–479.
  * Gelfand and Schliep, (2016) Gelfand, A. E. and Schliep, E. M. (2016).  Spatial statistics and Gaussian processes: A beautiful marriage.  Spatial Statistics, 18:86–104.
  * Genton et al., (2006) Genton, M. G., Butry, D. T., Gumpertz, M. L., and Prestemon, J. P. (2006).  Spatio-temporal analysis of wildfire ignitions in the St. Johns River water management district, Florida.  International Journal of Wildland Fire, 15(1):87–97.
  * Hazra and Huser, (2021) Hazra, A. and Huser, R. (2021).  Estimating high-resolution Red Sea surface temperature hotspots, using a low-rank semiparametric spatial model.  Annals of Applied Statistics, 15(2):572–596.
  * Hazra et al., (2021) Hazra, A., Huser, R., and Bolin, D. (2021+).  A sparse gaussian scale mixture process for short-range extremal dependence and long-range independence.  arXiv preprint 2112.10248.
  * Hazra et al., (2019) Hazra, A., Reich, B. J., Reich, D. S., Shinohara, R. T., and Staicu, A.-M. (2019).  A spatio-temporal model for longitudinal image-on-image regression.  Statistics in Biosciences, 11(1):22–46.
  * Hering et al., (2009) Hering, A. S., Bell, C. L., and Genton, M. G. (2009).  Modeling spatio-temporal wildfire ignition point patterns.  Environmental and Ecological Statistics, 16(2):225–250.
  * Hrafnkelsson et al., (2021) Hrafnkelsson, B., Siegert, S., Huser, R., Bakka, H., and Jóhannesson, Á. V. (2021).  Max-and-smooth: a two-step approach for approximate Bayesian inference in latent Gaussian models.  Bayesian Analysis, 16(2):611–638.
  * Huser et al., (2017) Huser, R., Opitz, T., and Thibaud, E. (2017).  Bridging asymptotic independence and dependence in spatial extremes using Gaussian scale mixtures.  Spatial Statistics, 21:166–186.
  * Huser and Wadsworth, (2020) Huser, R. and Wadsworth, J. L. (2020).  Advances in statistical modeling of spatial extremes.  Wiley Interdisciplinary Reviews: Computational Statistics, page e1537.
  * Jain et al., (2020) Jain, P., Coogan, S. C., Subramanian, S. G., Crowley, M., Taylor, S., and Flannigan, M. D. (2020).  A review of machine learning applications in wildfire science and management.  Environmental Reviews, 28(4):478–505.
  * Johannesson et al., (2021) Johannesson, Á. V., Siegert, S., Huser, R., Bakka, H., and Hrafnkelsson, B. (2021).  Approximate bayesian inference for analysis of spatio-temporal flood frequency data.  Annals of Applied Statistics.
  * Joseph et al., (2019) Joseph, M. B., Rossi, M. W., Mietkiewicz, N. P., Mahood, A. L., Cattau, M. E., St. Denis, L. A., Nagy, R. C., Iglesias, V., Abatzoglou, J. T., and Balch, J. K. (2019).  Spatiotemporal prediction of wildfire size extremes with Bayesian finite sample maxima.  Ecological Applications, 29(6):e01898.
  * Juan et al., (2012) Juan, P., Mateu, J., and Saez, M. (2012).  Pinpointing spatio-temporal interactions in wildfire patterns.  Stochastic Environmental Research and Risk Assessment, 26(8):1131–1150.
  * Katzfuss, (2013) Katzfuss, M. (2013).  Bayesian nonstationary spatial modeling for very large datasets.  Environmetrics, 24(3):189–200.
  * Koh et al., (2021) Koh, J., Pimont, F., Dupuy, J.-L., and Opitz, T. (2021).  Spatiotemporal wildfire modeling through point processes with moderate and extreme marks.  arXiv preprint:2105.08004.
  * Lindgren et al., (2011) Lindgren, F., Rue, H., and Lindström, J. (2011).  An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(4):423–498.
  * Møller and Díaz-Avalos, (2010) Møller, J. and Díaz-Avalos, C. (2010).  Structured spatio-temporal shot-noise Cox point process models, with a view to modelling forest fires.  Scandinavian Journal of Statistics, 37(1):2–25.
  * Møller et al., (1998) Møller, J., Syversveen, A. R., and Waagepetersen, R. P. (1998).  Log-Gaussian Cox processes.  Scandinavian Journal of Statistics, 25(3):451–482.
  * Nadeem et al., (2020) Nadeem, K., Taylor, S., Woolford, D. G., and Dean, C. (2020).  Mesoscale spatiotemporal predictive models of daily human-and lightning-caused wildland fire occurrence in British Columbia.  International Journal of Wildland Fire, 29(1):11–27.
  * Opitz, (2022) Opitz, T. (2022).  Editorial: Eva 2021 data competition on spatio-temporal prediction of wildfire activity in the united states. extremes.  Extremes (to appear).
  * Penttinen et al., (1992) Penttinen, A., Stoyan, D., and Henttonen, H. M. (1992).  Marked point processes in forest statistics.  Forest Science, 38(4):806–824.
  * Pereira and Turkman, (2019) Pereira, J. and Turkman, K. (2019).  Statistical models of vegetation fires: Spatial and temporal patterns.  In Handbook of Environmental and Ecological Statistics, pages 401–420. Chapman and Hall/CRC, Taylor & Francis.
  * Pimont et al., (2021) Pimont, F., Fargeon, H., Opitz, T., Ruffault, J., Barbero, R., Martin-StPaul, N., Rigolot, E., Rivière, M., and Dupuy, J.-L. (2021).  Prediction of regional wildfire activity in the probabilistic Bayesian framework of Firelihood.  Ecological applications, page e02316.
  * Preisler and Ager, (2013) Preisler, H. and Ager, A. (2013).  Forest-fire models.  Encyclopedia of Environmetrics, 3:2181–2185.
  * Preisler et al., (2004) Preisler, H. K., Brillinger, D. R., Burgan, R. E., and Benoit, J. (2004).  Probability based models for estimation of wildfire risk.  International Journal of Wildland Fire, 13(2):133–142.
  * Preisler and Westerling, (2007) Preisler, H. K. and Westerling, A. L. (2007).  Statistical model for forecasting monthly large wildfire events in western United States.  Journal of Applied Meteorology and Climatology, 46(7):1020–1030.
  * Pyne et al., (1996) Pyne, S., Andrew, P., and Laven, R. (1996).  Introduction to wildland and rural fire.  Princeton University Press, Princeton, NJ (654p).
  * Ríos-Pena et al., (2018) Ríos-Pena, L., Kneib, T., Cadarso-Suárez, C., Klein, N., and Marey-Pérez, M. (2018).  Studying the occurrence and burnt area of wildfires using zero-one-inflated structured additive beta regression.  Environmental Modelling & Software, 110:107–118.
  * Rue and Held, (2005) Rue, H. and Held, L. (2005).  Gaussian Markov random fields: theory and applications.  Chapman and Hall/CRC, Taylor & Francis.
  * Saha et al., (2021) Saha, A., Basu, S., and Datta, A. (2021).  Random forests for spatially dependent data.  Journal of the American Statistical Association, to appear, pages 1–46.
  * Scott, (2000) Scott, A. C. (2000).  The pre-quaternary history of fire.  Palaeogeography, palaeoclimatology, palaeoecology, 164(1-4):281–329.
  * Serra et al., (2012) Serra, L., Saez, M., Varga, D., Tobías, A., Juan, P., and Mateu, J. (2012).  Spatio-temporal modelling of wildfires in Catalonia, Spain, 1994–2008, through log-Gaussian Cox processes.  Modelling, Monitoring and Management of Forest Fires III, page 11139\.
  * Trucchia et al., (2018) Trucchia, A., Egorova, V., Pagnini, G., and Rochoux, M. C. (2018).  Surrogate-based global sensitivity analysis for turbulence and fire-spotting effects in regional-scale wildland fire modeling.  arXiv preprint:1809.05430.
  * Welling and Teh, (2011) Welling, M. and Teh, Y. W. (2011).  Bayesian learning via stochastic gradient Langevin dynamics.  In Proceedings of the 28th International Conference on Machine Learning (ICML-11), pages 681–688.
  * Wikle, (2010) Wikle, C. K. (2010).  Low-rank representations for spatial processes.  In Handbook of spatial statistics, pages 114–125. CRC Press, Taylor & Francis.
  * Wuebbles et al., (2017) Wuebbles, D. J., Fahey, D. W., Hibbard, K. A., Arnold, J. R., DeAngelo, B., Doherty, S., Easterling, D. R., Edmonds, J., Edmonds, T., Hall, T., et al. (2017).  Climate science special report: Fourth national climate assessment (NCA4), Volume I.
  * Xi et al., (2019) Xi, D. D., Taylor, S. W., Woolford, D. G., and Dean, C. (2019).  Statistical models of key components of wildfire risk.  Annual Review of Statistics and its Application, 6:197–222.
  * Yadav et al., (2021) Yadav, R., Huser, R., and Opitz, T. (2021+).  A flexible Bayesian hierarchical modeling framework for spatially dependent peaks over-threshold-data.  arXiv preprint:2112.09530.
  * Zammit-Mangion and Cressie, (2021) Zammit-Mangion, A. and Cressie, N. (2021).  Frk: An R package for spatial and spatio-temporal prediction with large datasets.  Journal of Statistical Software, 98(4):1–48.

Supplementary Material

##  7 MCMC details

###  7.1 Stage 1: Spatial indicator kriging for wildfire occurrence data

We define the indicator variables $Z_{t}(\bm{s}_{i})$ as

$\displaystyle Z_{t}(\bm{s}_{i})=\begin{cases}1,\hskip
11.38092pt\textrm{if}\begin{cases}\parbox[t]{381.5874pt}{$\textrm{BA}_{t}(\bm{s}_{i})>0,\leavevmode\nobreak\
\textrm{CNT}_{t}(\bm{s}_{i})>0$,}\\\
\text{$\textrm{BA}_{t}(\bm{s}_{i})>0,\textrm{CNT}_{t}(\bm{s}_{i})$ is missing,
$\leavevmode\nobreak\ ${or}$\leavevmode\nobreak\ \leavevmode\nobreak\
$}\text{$\textrm{CNT}_{t}(\bm{s}_{i})>0,\textrm{BA}_{t}(\bm{s}_{i})$ is
missing,}\end{cases}\\\ \vspace{-5mm}\\\ 0,\hskip
11.38092pt\textrm{if}\begin{cases}\parbox[t]{381.5874pt}{$\textrm{BA}_{t}(\bm{s}_{i})=0,\leavevmode\nobreak\
\textrm{CNT}_{t}(\bm{s}_{i})=0$,}\\\
\text{$\textrm{BA}_{t}(\bm{s}_{i})=0,\textrm{CNT}_{t}(\bm{s}_{i})$ is missing,
$\leavevmode\nobreak\ ${or}$\leavevmode\nobreak\ \leavevmode\nobreak\
$}\text{$\textrm{CNT}_{t}(\bm{s}_{i})=0,\textrm{BA}_{t}(\bm{s}_{i})$ is
missing,}\end{cases}\\\ ${NA},$\leavevmode\nobreak\ \leavevmode\nobreak\
\text{if \leavevmode\nobreak\ \leavevmode\nobreak\
$\textrm{BA}_{t}(\bm{s}_{i})$ and $\textrm{CNT}_{t}(\bm{s}_{i})$ both are
missing.}\end{cases}$

Further, we assume that $Z_{t}(\bm{s}_{i})$ are independent and identically
distributed (iid) across $t$, and model it using a latent Gaussian process
(GP) $X_{t}(\cdot)$ as

$\displaystyle
Z_{t}(\bm{s}_{i})=\begin{cases}1&\parbox[t]{95.39502pt}{if\leavevmode\nobreak\
\leavevmode\nobreak\ $X_{t}(\bm{s}_{i})>0$}\\\ 0&\text{if\leavevmode\nobreak\
\leavevmode\nobreak\
$X_{t}(\bm{s}_{i})<0$},\end{cases}\textrm{where}\leavevmode\nobreak\
X_{t}(\bm{s}_{i})=\mu_{Z}(\bm{s}_{i})+\varepsilon_{t}(\bm{s}_{i}),$

and $\varepsilon_{t}(\cdot)$ are spatial GPs that are iid across $t$, with
$\textrm{E}[\varepsilon_{t}(\bm{s})]=0$ and
$\textrm{Var}[\varepsilon_{t}(\bm{s})]=1$ (approximately). We assume that the
process $\varepsilon_{t}(\cdot)$ follows an isotropic Matérn spatial
correlation with range parameter $\phi_{\varepsilon}$ and smoothness parameter
fixed at one. Further, we model
$\bm{\varepsilon}_{t}=[\varepsilon_{t}(\bm{s}_{1}),\ldots,\varepsilon_{t}(\bm{s}_{N})]^{\prime}$
as
$\bm{\varepsilon}_{t}\sim\textrm{Normal}_{N}(\bm{A}\bm{\varepsilon}^{*}_{t},(1-r_{\varepsilon})\bm{I}_{N})$
where
$\bm{\varepsilon}^{*}_{t}\sim\textrm{Normal}_{N}(\bm{0},r_{\varepsilon}\bm{Q}_{\phi_{\varepsilon}}^{-1})$.
We denote
$\bm{\mu}_{Z}=[\mu_{Z}(\bm{s}_{1}),\ldots,\mu_{Z}(\bm{s}_{N})]^{\prime}$ and
choose the prior
$\bm{\mu}_{Z}\sim\textrm{Normal}_{N}(\bm{D}\bm{\theta}_{\mu},\tau^{-1}_{\mu}\bm{I}_{N})$.
Also, we denote
$\bm{X}_{t}=[X_{t}(\bm{s}_{1}),\ldots,X_{t}(\bm{s}_{N})]^{\prime}$. The set of
latent processes, parameters and hyper-parameters in the model are

$\displaystyle\bm{\Theta}_{1}=\left\\{\\{\bm{X}_{t}\\}_{t=1}^{T},\bm{\mu}_{Z},\bm{\theta}_{\mu},\tau_{\mu},\\{\bm{\varepsilon}^{*}_{t}\\}_{t=1}^{T},\phi_{\varepsilon},r_{\varepsilon}\right\\}.$

The MCMC steps for updating the parameters in $\bm{\Theta}_{1}$ are as
follows. For a specific a parameter (or a set of parameters), we call “$rest$”
the observed data, all the parameters and hyper-parameters in
$\bm{\Theta}_{1}$ except that parameter (or that set of parameters).

Throughout the algorithm, we update the parameters and hyper-parameters
consecutively following the full conditional posterior distributions at each
iteration as follows:

$\bm{\mu}_{Z}|rest$
The prior distribution of $\bm{\mu}_{Z}$ is
$\bm{\mu}_{Z}\sim\textrm{Normal}_{N}(\bm{D}\bm{\theta}_{\mu},\tau^{-1}_{\mu}\bm{I}_{N})$.
The full conditional posterior distribution of $\bm{\mu}_{Z}$ is

$\displaystyle\bm{\mu}_{Z}|rest\sim\textrm{Normal}_{N}\left(\left[{\textstyle{\displaystyle
T\over\displaystyle
1-r_{\varepsilon}}}+\tau_{\mu}\right]^{-1}\left[{\textstyle{\displaystyle
1\over\displaystyle
1-r_{\varepsilon}}}\sum_{t=1}^{T}(\bm{X}_{t}-\bm{A}\bm{\varepsilon}^{*}_{t})+\tau_{\mu}\bm{D}\bm{\theta}_{\mu}\right]^{-1},\left[{\textstyle{\displaystyle
T\over\displaystyle
1-r_{\varepsilon}}}+\tau_{\mu}\right]^{-1}\bm{I}_{N}\right).$

$\bm{\theta}_{\mu}|rest$
We consider the conjugate noninformative prior
$\bm{\theta}_{\mu}\sim\textrm{Normal}_{P}(\bm{0},10^{2}\bm{I}_{P})$. The full
conditional posterior distribution is

$\displaystyle\bm{\theta}_{\mu}|rest\sim\textrm{Normal}_{P}\left(\left[\tau_{\mu}\bm{D}^{\prime}\bm{D}+10^{-2}\bm{I}_{P}\right]^{-1}\tau_{\mu}\bm{D}^{\prime}\bm{\mu}_{Z},\leavevmode\nobreak\
\left[\tau_{\mu}\bm{D}^{\prime}\bm{D}+10^{-2}\bm{I}_{P}\right]^{-1}\right).$

$\tau_{\mu}|rest$
We consider the prior $\tau_{\mu}\sim\textrm{Gamma}(0.1,0.1)$. The full
conditional posterior is

$\displaystyle\tau_{\mu}|rest\sim\textrm{Gamma}\left(0.1+N/2,\leavevmode\nobreak\
0.1+(\bm{\mu}_{Z}-\bm{D}\bm{\theta}_{\mu})^{\prime}(\bm{\mu}_{Z}-\bm{D}\bm{\theta}_{\mu})/2\right).$

$\\{\bm{\varepsilon}^{*}_{t}\\}_{t=1}^{T}|rest$
The unconditional distribution of $\bm{\varepsilon}^{*}_{t}$ is
$\bm{\varepsilon}^{*}_{t}\sim\textrm{Normal}_{N^{*}}(\bm{0},r_{\varepsilon}\bm{Q}_{\phi_{\varepsilon}}^{-1})$.
The full conditional posterior distribution of $\bm{\varepsilon}^{*}_{t}$ is
$\bm{\varepsilon}^{*}_{t}|rest\sim\textrm{Normal}_{N^{*}}(\bm{\mu}^{\ast}_{\varepsilon^{*}},\bm{\Sigma}^{\ast}_{\varepsilon^{*}})$,
where

$\displaystyle\bm{\Sigma}^{\ast}_{\varepsilon^{*}}=\left[(1-r_{\varepsilon})^{-1}\bm{A}^{\prime}\bm{A}+r_{\varepsilon}^{-1}\bm{Q}_{\phi_{\varepsilon}}\right]^{-1},\leavevmode\nobreak\
\leavevmode\nobreak\
\bm{\mu}^{\ast}_{\varepsilon^{*}}=\bm{\Sigma}^{\ast}_{\varepsilon^{*}}[(1-r_{\varepsilon})^{-1}\bm{A}^{\prime}(\bm{X}_{t}-\bm{\mu}_{Z})].$

The observations are assumed to be independent across $t=1,\ldots,T$, and
hence, the full conditional posterior distributions of
$\bm{\varepsilon}^{*}_{t};t=1,\ldots,T$ are independent. Thus,
$\bm{\varepsilon}^{*}_{t};t=1,\ldots,T$ are updated within MCMC in parallel.

$\phi_{\varepsilon}|rest$
We consider the prior
$\phi_{\varepsilon}\sim\textrm{Uniform}(0,2\Delta_{\mathcal{S}})$, where
$\Delta_{\mathcal{S}}$ is the largest Euclidean distance between two data
locations. Let the $m$-th MCMC sample from $\phi_{\varepsilon}$ be denoted by
$\phi_{\varepsilon}^{(m)}$. Considering a logit transformation, we obtain
$\phi_{\varepsilon}^{*(m)}\in\mathbb{R}$ from $\phi_{\varepsilon}^{(m)}$, and
simulate $\phi_{\varepsilon}^{*(c)}$ from
$\textrm{Normal}(\phi_{\varepsilon}^{*(m)},s_{\phi_{\varepsilon}}^{2})$
distribution, where $s_{\phi_{\varepsilon}}$ is the standard deviation of the
candidate normal distribution. Subsequently, using an inverse-logit
transformation, we obtain $\phi_{\varepsilon}^{(c)}$ from
$\phi_{\varepsilon}^{*(c)}$, and either accept $\phi_{\varepsilon}^{(c)}$ as a
candidate from the posterior distribution of $\phi_{\varepsilon}$ or reject
it. The acceptance ratio is

$\displaystyle\mathcal{R}$ $\displaystyle=$
$\displaystyle{\textstyle{\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{N^{*}}}\left(\bm{\varepsilon}^{*}_{t};\bm{0},r_{\varepsilon}\bm{Q}_{\phi_{\varepsilon}^{(c)}}^{-1}\right)\over\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{N^{*}}}\left(\bm{\varepsilon}^{*}_{t};\bm{0},r_{\varepsilon}\bm{Q}_{\phi_{\varepsilon}^{(m)}}^{-1}\right)}}\times{\textstyle{\displaystyle\phi_{\varepsilon}^{(c)}\left(2\Delta_{\mathcal{S}}-\phi_{\varepsilon}^{(c)}\right)\over\displaystyle\phi_{\varepsilon}^{(m)}\left(2\Delta_{\mathcal{S}}-\phi_{\varepsilon}^{(m)}\right)}},$

where $f_{\textrm{Normal}_{n}}(\cdot;\bm{\mu},\bm{\Sigma})$ denotes the
$n$-variate normal density with mean $\bm{\mu}$ and covariance matrix
$\bm{\Sigma}$. The candidate is accepted with probability
$min\\{\mathcal{R},1\\}$.

$r_{\varepsilon}|rest$
We consider the prior $r_{\varepsilon}\sim\textrm{Uniform}(0,1)$. Suppose
$r_{\varepsilon}^{(m)}$ denotes the $m$-th MCMC sample from $r_{\varepsilon}$.
We simulate a candidate sample $r_{\varepsilon}^{(c)}$ from
$r_{\varepsilon}^{(m)}$ following a procedure similar to simulating
$\phi_{\varepsilon}^{(c)}$ from $\phi_{\varepsilon}^{(m)}$. The Metropolis-
Hastings acceptance ratio is

$\displaystyle\mathcal{R}={\textstyle{\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{N}}\left(\bm{X}_{t};\bm{\mu}+\bm{A}\bm{\varepsilon}^{*}_{t},(1-r^{(c)})\bm{I}_{N}\right)\over\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{N}}\left(\bm{X}_{t};\bm{\mu}+\bm{A}\bm{\varepsilon}^{*}_{t},(1-r^{(m)})\bm{I}_{N}\right)}}\times{\textstyle{\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{N^{*}}}\left(\bm{\varepsilon}^{*}_{t};\bm{0},r_{\varepsilon}^{(c)}\bm{Q}_{\phi_{\varepsilon}}^{-1}\right)\over\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{N^{*}}}\left(\bm{\varepsilon}^{*}_{t};\bm{0},r_{\varepsilon}^{(m)}\bm{Q}_{\phi_{\varepsilon}}^{-1}\right)}}\times{\textstyle{\displaystyle
r_{\varepsilon}^{(c)}\left(1-r_{\varepsilon}^{(c)}\right)\over\displaystyle
r_{\varepsilon}^{(m)}\left(1-r_{\varepsilon}^{(m)}\right)}}.$

$X_{t}(\bm{s}_{i}),i=1,\ldots,N,t=1,\ldots,T|rest$
Let the $i$-th entry of $\bm{A}\bm{\varepsilon}^{*}_{t}$ be
$\tilde{\varepsilon}^{*}_{t}(\bm{s}_{i})$. The full conditional distribution
of $X_{t}(\bm{s}_{i})$ is

$\displaystyle X_{t}(\bm{s}_{i})|Z_{t}(\bm{s}_{i})=0,rest$ $\displaystyle\sim$
$\displaystyle\textrm{Truncated-
Normal}_{(-\infty,0)}\left(\mu_{Z}(\bm{s}_{i})+\tilde{\varepsilon}_{t}(\bm{s}_{i}),(1-r_{\varepsilon})\right),$
$\displaystyle X_{t}(\bm{s}_{i})|Z_{t}(\bm{s}_{i})=1,rest$ $\displaystyle\sim$
$\displaystyle\textrm{Truncated-
Normal}_{(0,\infty)}\left(\mu_{Z}(\bm{s}_{i})+\tilde{\varepsilon}_{t}(\bm{s}_{i}),(1-r_{\varepsilon})\right),$
$\displaystyle X_{t}(\bm{s}_{i})|Z_{t}(\bm{s}_{i})=\texttt{NA},rest$
$\displaystyle\sim$
$\displaystyle\textrm{Normal}\left(\mu_{Z}(\bm{s}_{i})+\tilde{\varepsilon}_{t}(\bm{s}_{i}),(1-r_{\varepsilon})\right).$

###  7.2 Stage 2: Bivariate GMRF modeling for positive CNT and BA data

We model the positive CNT and BA data jointly using a bivariate model, where
we model log-CNT and log-BA as

$\displaystyle\log\\{\text{CNT}_{t}(\bm{s}_{i})\\}=$
$\displaystyle\mu_{1}(\bm{s}_{i})+\sigma_{1}(\bm{s}_{i})W_{t1}(\bm{s}_{i})$
$\displaystyle\log\\{\text{BA}_{t}(\bm{s}_{i})\\}=$
$\displaystyle\mu_{2}(\bm{s}_{i})+\sigma_{2}(\bm{s}_{i})W_{t2}(\bm{s}_{i}),$
(15)

where
$\bm{W}_{t}(\bm{s}_{i})=[W_{t1}(\bm{s}_{i}),W_{t2}(\bm{s}_{i})]^{\prime}$ is a
bivariate standard Gaussian process with separable spatial covariance
structure and independent across time. We approximate the underlying Gaussian
process by a GMRF, detailed in the main paper. We draw approximate Bayesian
inference using a two-stage method, where first we obtain smoothed estimates
of $\mu_{1}(\bm{s}_{i})$, $\sigma_{1}(\bm{s}_{i})$, $\mu_{2}(\bm{s}_{i})$, and
$\sigma_{2}(\bm{s}_{i})$ at each $\bm{s}_{i}$, using fixed rank kriging.

In the second step, we model the standardized variables $W_{t1}(\bm{s}_{i})$
and $W_{t2}(\bm{s}_{i})$. We use the same SPDE mesh and the same projection
matrix $\bm{A}$ as in Stage 1. Let
$\bm{W}_{tp}=[W_{tp}(\bm{s}_{i}),\ldots,W_{tp}(\bm{s}_{N})]^{\prime}$ for
$p=1,2$, and
$\bm{W}_{t}=[\bm{W}^{\prime}_{t1},\bm{W}^{\prime}_{t2}]^{\prime}$. We model
$\bm{W}_{t}$ as
$\bm{W}_{t}=\left[\bm{I}_{2}\otimes\bm{A}\right]\bm{\eta}^{*}_{t}+\tilde{\bm{\eta}}_{t}$,
where

$\displaystyle\bm{\eta}^{*}_{t}\sim\textrm{Normal}_{2N^{*}}\left(\bm{0},r_{\eta}\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\
1\end{array}\right)\otimes\bm{Q}^{-1}_{\phi_{\eta}}\right),\leavevmode\nobreak\
\textrm{and}\leavevmode\nobreak\ $ (18)
$\displaystyle\tilde{\bm{\eta}}_{t}\sim\textrm{Normal}_{2N}\left(\bm{0},(1-r_{\eta})\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\ 1\end{array}\right)\otimes\bm{I}_{N}\right).$ (21)

The set of latent process, parameters, and hyper-parameters in the model are

$\displaystyle\bm{\Theta}_{2}=\left\\{\\{\bm{\eta}^{*}_{t}\\}_{t=1}^{T},\phi_{\eta},r_{\eta},\rho_{\eta}\right\\}.$

The MCMC steps for updating the parameters in $\bm{\Theta}_{2}$ and simulating
posterior samples from missing $W_{tp}(\bm{s}_{i})$ are as follows. For a
specific a parameter (or a set of parameters), we call “$rest$” the observed
data, all the parameters and hyper-parameters in $\bm{\Theta}_{2}$ except that
parameter (or that set of parameters).

Throughout the algorithm, we update the parameters and hyper-parameters
consecutively following the full conditional posterior distributions at each
iteration as follows:

$\\{\bm{\eta}^{*}_{t}\\}_{t=1}^{T}|rest$
The full conditional posterior distribution of $\bm{\eta}^{*}_{t}$ is
$\bm{\eta}^{*}_{t}|rest\sim\textrm{Normal}_{2N^{*}}(\bm{\mu}^{\ast}_{\eta^{*}},\bm{\Sigma}^{\ast}_{\eta^{*}})$,
where

$\displaystyle\bm{\mu}^{\ast}_{\eta^{*}}=\left(\bm{I}_{2}\otimes\left[(1-r_{\eta})^{-1}\bm{A}^{\prime}\bm{A}+r_{\eta}^{-1}\bm{Q}_{\phi_{\eta}}\right]^{-1}\right)(1-r_{\varepsilon})^{-1}[\bm{I}_{2}\otimes\bm{A}^{\prime}]\bm{W}_{t},$
(22)
$\displaystyle\bm{\Sigma}^{\ast}_{\eta^{*}}=\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\
1\end{array}\right)\otimes\left[(1-r_{\eta})^{-1}\bm{A}^{\prime}\bm{A}+r_{\eta}^{-1}\bm{Q}_{\phi_{\eta}}\right]^{-1}.$
(25)

The vectors $\bm{W}_{t},t=1,\ldots,T$ are assumed to be independent across
$t=1,\ldots,T$, and hence, the full conditional posterior distributions of
$\bm{\eta}^{*}_{t};t=1,\ldots,T$ are independent. Thus,
$\bm{\eta}^{*}_{t};t=1,\ldots,T$ are updated within MCMC in parallel.

$\phi_{\eta}|rest$
Similar to $\phi_{\varepsilon}$, we consider the prior
$\phi_{\eta}\sim\textrm{Uniform}(0,2\Delta_{\mathcal{S}})$. Let the $m$-th
MCMC sample from $\phi_{\eta}$ be denoted by $\phi_{\eta}^{(m)}$. We simulate
a candidate sample $\phi_{\eta}^{(c)}$ from $\phi_{\eta}^{(m)}$ following a
procedure similar to simulating $\phi_{\varepsilon}^{(c)}$ from
$\phi_{\varepsilon}^{(m)}$. The Metropolis-Hastings acceptance ratio is

$\displaystyle\mathcal{R}$ $\displaystyle=$
$\displaystyle{\textstyle{\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N^{*}}}\left(\bm{\eta}^{*}_{t};\bm{0},r_{\eta}\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\
1\end{array}\right)\otimes\bm{Q}^{-1}_{\phi_{\eta}^{(c)}}\right)\over\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N^{*}}}\left(\bm{\eta}^{*}_{t};\bm{0},r_{\eta}\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\
1\end{array}\right)\otimes\bm{Q}^{-1}_{\phi_{\eta}^{(m)}}\right)}}\times{\textstyle{\displaystyle\phi_{\eta}^{(c)}\left(2\Delta_{\mathcal{S}}-\phi_{\eta}^{(c)}\right)\over\displaystyle\phi_{\eta}^{(m)}\left(2\Delta_{\mathcal{S}}-\phi_{\eta}^{(m)}\right)}}.$
(30)

$r_{\eta}|rest$
We consider the prior $r_{\eta}\sim\textrm{Uniform}(0,1)$. Suppose
$r_{\eta}^{(m)}$ denotes the $m$-th MCMC sample from $r_{\eta}$. We simulate a
candidate sample $r_{\eta}^{(c)}$ from $r_{\eta}^{(m)}$ following a procedure
similar to simulating $\phi_{\varepsilon}^{(c)}$ from
$\phi_{\varepsilon}^{(m)}$. The Metropolis-Hastings acceptance ratio is

$\displaystyle\mathcal{R}={\textstyle{\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N}}\left(\bm{W}_{t};\left[\bm{I}_{2}\otimes\bm{A}\right]\bm{\eta}^{*}_{t},(1-r_{\eta}^{(c)})\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\
1\end{array}\right)\otimes\bm{I}_{N}\right)\over\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N}}\left(\bm{W}_{t};\left[\bm{I}_{2}\otimes\bm{A}\right]\bm{\eta}^{*}_{t},(1-r_{\eta}^{(m)})\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\ 1\end{array}\right)\otimes\bm{I}_{N}\right)}}\times$ (34)
$\displaystyle{\textstyle{\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N^{*}}}\left(\bm{\eta}^{*}_{t};\bm{0},r_{\eta}^{(c)}\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\
1\end{array}\right)\otimes\bm{Q}^{-1}_{\phi_{\eta}}\right)\over\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N^{*}}}\left(\bm{\eta}^{*}_{t};\bm{0},r_{\eta}^{(m)}\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\
1\end{array}\right)\otimes\bm{Q}^{-1}_{\phi_{\eta}}\right)}}\times{\textstyle{\displaystyle
r_{\eta}^{(c)}\left(1-r_{\eta}^{(c)}\right)\over\displaystyle
r_{\eta}^{(m)}\left(1-r_{\eta}^{(m)}\right)}}.$ (37)

$\rho_{\eta}|rest$
We consider the prior $\rho_{\eta}\sim\textrm{Uniform}(0,1)$. Suppose
$\rho_{\eta}^{(m)}$ denotes the $m$-th MCMC sample from $\rho_{\eta}$. We
simulate a candidate sample $\rho_{\eta}^{(c)}$ from $\rho_{\eta}^{(m)}$
following a procedure similar to simulating $\phi_{\varepsilon}^{(c)}$ from
$\phi_{\varepsilon}^{(m)}$. The Metropolis-Hastings acceptance ratio is

$\displaystyle\mathcal{R}={\textstyle{\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N}}\left(\bm{W}_{t};\left[\bm{I}_{2}\otimes\bm{A}\right]\bm{\eta}^{*}_{t},(1-r_{\eta})\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}^{(c)}\\\
\rho_{\eta}^{(c)}\leavevmode\nobreak\ \leavevmode\nobreak\
1\end{array}\right)\otimes\bm{I}_{N}\right)\over\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N}}\left(\bm{W}_{t};\left[\bm{I}_{2}\otimes\bm{A}\right]\bm{\eta}^{*}_{t},(1-r_{\eta})\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}^{(m)}\\\
\rho_{\eta}^{(m)}\leavevmode\nobreak\ \leavevmode\nobreak\
1\end{array}\right)\otimes\bm{I}_{N}\right)}}\times$ (42)
$\displaystyle{\textstyle{\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N^{*}}}\left(\bm{\eta}^{*}_{t};\bm{0},r_{\eta}\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}^{(c)}\\\
\rho_{\eta}^{(c)}\leavevmode\nobreak\ \leavevmode\nobreak\
1\end{array}\right)\otimes\bm{Q}^{-1}_{\phi_{\eta}}\right)\over\displaystyle\prod_{t=1}^{T}f_{\textrm{Normal}_{2N^{*}}}\left(\bm{\eta}^{*}_{t};\bm{0},r_{\eta}\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}^{(m)}\\\
\rho_{\eta}^{(m)}\leavevmode\nobreak\ \leavevmode\nobreak\
1\end{array}\right)\otimes\bm{Q}^{-1}_{\phi_{\eta}}\right)}}\times{\textstyle{\displaystyle\rho_{\eta}^{(c)}\left(1-\rho_{\eta}^{(c)}\right)\over\displaystyle\rho_{\eta}^{(m)}\left(1-\rho_{\eta}^{(m)}\right)}}.$
(46)

Further, we also need to simulate missing $W_{tp}(\bm{s}_{i})$ from their full
conditional posterior distributions. Let the elements of $2N$-length vector
$\left[\bm{I}_{2}\otimes\bm{A}\right]\bm{\eta}^{*}_{t}$ be denoted by
$[\bm{\tilde{\eta}}^{*^{\prime}}_{t1},\bm{\tilde{\eta}}^{*^{\prime}}_{t2}]^{\prime}$,
where
$\bm{\tilde{\eta}}^{*}_{t1}=[\tilde{\eta}^{*}_{t1}(\bm{s}_{1}),\ldots,\tilde{\eta}^{*}_{t1}(\bm{s}_{N})]^{\prime}$
and
$\bm{\tilde{\eta}}^{*}_{t2}=[\tilde{\eta}^{*}_{t2}(\bm{s}_{1}),\ldots,\tilde{\eta}^{*}_{t2}(\bm{s}_{N})]^{\prime}$.
The full conditional distribution of
$\bm{W}_{t}(\bm{s}_{i})=[W_{t1}(\bm{s}_{i}),W_{t2}(\bm{s}_{i})]^{\prime}$ is

$\displaystyle\bm{W}_{t}(\bm{s}_{i})\sim\textrm{Normal}_{2}\left(\left(\begin{array}[]{c}\tilde{\eta}^{*}_{t1}(\bm{s}_{i})\\\
\tilde{\eta}^{*}_{t2}(\bm{s}_{i})\end{array}\right),(1-r_{\eta})\left(\begin{array}[]{c}1\leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{\eta}\\\ \rho_{\eta}\leavevmode\nobreak\
\leavevmode\nobreak\ 1\end{array}\right)\right).$ (50)