This is a preprint of the paper. This is a preprint of the paper. This is a
preprint of the paper.

Villejo et al.

[$\ast$]Corresponding author<EMAIL_ADDRESS>
0Year 0Year 0Year

# A Data Fusion Model for Meteorological Data using the INLA-SPDE method

Stephen Jun Villejo 0000-0002-0510-3143  Sara Martino  Finn Lindgren  Janine B
Illian  School of Mathematics and Statistics, University of Glasgow, Glasgow,
United Kingdom  Department of Mathematical Sciences, Norwegian University of
Science and Technology, Trondheim, Norway  School of Mathematics, College of
Science and Engineering, University of Edinburgh, Edinburgh, United Kingdom
School of Statistics, University of the Philippines Diliman, Quezon City,
Philippines

(2023; 2023; Date; Date; Date)

###### Abstract

This work aims to combine two primary meteorological data sources in the
Philippines: data from a sparse network of weather stations and outcomes of a
numerical weather prediction model. To this end, we propose a data fusion
model which is primarily motivated by the problem of sparsity in the
observational data and the use of a numerical prediction model as an
additional data source in order to obtain better predictions for the variables
of interest. The proposed data fusion model assumes that the different data
sources are error-prone realizations of a common latent process. The outcomes
from the weather stations follow the classical error model while the outcomes
of the numerical weather prediction model involves a constant multiplicative
bias parameter and an additive bias which is spatially-structured and time-
varying. We use a Bayesian model averaging approach with the integrated nested
Laplace approximation (INLA) for doing inference. The proposed data fusion
model outperforms the stations-only model and the regression calibration
approach, when assessed using leave-group-out cross-validation (LGOCV). We
assess the benefits of data fusion and evaluate the accuracy of predictions
and parameter estimation through a simulation study. The results show that the
proposed data fusion model generally gives better predictions compared to the
stations-only approach especially with sparse observational data.

###### keywords:

data fusion, integrated nested Laplace approximation (INLA), Bayesian model
averaging

##  1 Introduction

A common aim in spatial statistical modelling is to predict values of a
spatial variable at unsampled locations, based on measurements taken at a
finite – and often relatively small – number of locations. In the context of
environmental sciences and meteorology in particular, data on e.g. air-
quality, environmental pollution, surface temperature, or rainfall are
collected through a network of monitoring stations (lawson2016handbook, ;
lee2017rigorous, ; blangiardo2016two, ; greven2011approach, ; chien2014impact,
; arab2014modelling, ; jaya2022spatiotemporal, ; lee2015impact, ). These data
are used for prediction, and for understanding the spatio-temporal dynamics of
the underlying processes and the impact of these variables on potential
outcomes of interest, such as health outcomes. However, due to the high
maintenance costs, the monitoring networks are typically spatially sparse
lawson2016handbook . Increasingly, data from additional sources derived from
satellite images or outcomes of numerical models with a high spatial
resolution are available. These can be used jointly with the monitoring data
to improve the accuracy of predictions in a process that combines information
from different data sources and is often referred to as data fusion or data
assimilation (bauer2015quiet, ; gettelman2022future, ; lawson2016handbook, ).
The goal is to exploit the better spatial resolution of the additional data to
fill gaps in areas only sparsely covered by the monitoring stations in order
to predict and map the variables into space with more accuracy at smaller
scales than based on the monitoring stations data alone. However, a general
issue that is common in attempts to combine data from more than one source is
that the various data streams differ in their quality. Numerical or satellite
data, specifically, are often biased due to calibration issues, and these
biases have to be accounted for in the modelling process. Motivated by a data
application of meteorological data from the Philippines, we develop a joint
modeling methodology that not only flexibly accounts for differences in data
quality, but also allows us to gauge the quality of the different data
sources.

In particular, this work combines two primary meteorological data sources for
the Philippines: observational data from a sparse network of weather synoptic
stations, and simulated outcomes from a numerical weather prediction model
called the Global Spectral Model (PHClimate_PAGASA, ), which provides broad
spatial coverage but typically biased due to sensitivity to model
initialization and parameterization. The observational data from the weather
stations are likely to be less biased, but they provide limited spatial
coverage since key areas in space are under-sampled, as is evident in Figure
1. It is therefore potentially advantageous to simultaneously use both data
sources through data fusion. Unlike approaches discussed in the literature so
far, which are discussed in Section 1.2, this work accounts for calibration
biases in the data fusion model, specifies a flexible bias formulation,
defines a single and interpretable latent process for the different data
sources, and specifies a measurement error for the monitoring stations data. A
simulation study compares our proposed data fusion approach to the stations-
only model. In the data application, we compare the predictions derived from
the proposed data fusion approach, a model based on the data from the stations
only, and another approach for data fusion called regression calibration using
leave-group-out cross-validation (liu2022leave, ; adin2023automatic, ).
Although the data application of the proposed data fusion model in this work
only considers two data sources, the model and the framework extend to more
than two data sources and are relevant beyond this specific scenario.

###  1.1 Meteorological data from the Philippines

The Philippines is an archipelagic country, covering an area of ca. 300
thousand km2 (see Figure 1), situated in tropical Southeast Asia. The eastern
part and some southern parts of the country are mostly classified as tropical
rainforest and is characterized by no distinct wet or dry season and a
relatively high rainfall all year round (hess2011mcknight, ). On the other
hand, most of the country’s western section is classified as tropical monsoon
or tropical Savannah and is characterized by pronounced dry and wet seasons
(beck2018present, ). Such spatio-temporal interaction is the basis for the
climate types of the Philippines (PHClimate_PAGASA, ).

Figure 1: Meteorological data sources for PH: a sparse network of weather
synoptic stations and an outcome of a numerical weather prediction model
called Global Spectral Model. The measurements are monthly aggregated values
of temperature for August 2019.

The Philippine Atmospheric, Geophysical and Astronomical Services
Administration (PAGASA) maintains a network of 56 weather synoptic stations
which regularly record several meteorological variables such as temperature,
relative humidity, and precipitation. The spatial locations of the weather
stations are shown in Figure 1. It is a sparse network relative to the total
surface area of the country. Some regions, especially in the northern part,
are heavily undersampled. Reconstructing meteorological variable surfaces
based only on the weather stations would therefore result in high uncertainty.
To augment this sparsity problem in the stations data, PAGASA utilizes the
outcomes of the Global Spectrum Model (GSM), which is a numerical weather
prediction model maintained by the Japanese Meteorological Agency. The GSM
provides forecast outputs of up to 132 hours four times a day (with initial
times 0000, 0600, 1200, and 1800 UTC) within 4 hours of the initial time, and
up to 264 hours twice a day (with initial time 0000 and 1200 UTC) within 7
hours of the initial time. As an example, Figure 1 shows the map of mean
temperatures from the GSM for August 2019 at a spatial resolution of 0.5
degree, which corresponds to approximately 55km $\times$ 55 km. Although the
figure shows the outcomes are gridded, these outcomes are interpreted by
PAGASA as point-referenced at the centroids PHClimate_PAGASA .

Figure 2: Scatterplot of the observed values at the weather stations versus
interpolated outcomes of the GSM for three meteorological variables:
temperature, relative humidity, and log-transformed rainfall. The plot shows
the discrepancies in the values between the two data sources.

The GSM outcomes have a much larger spatial coverage, but it is known that
outputs from numerical weather models are biased (lawson2016handbook, ;
bauer2015quiet, ). To assess this, we initially do spatial interpolation using
the GSM values on three meteorological variables of interest – temperature,
relative humidity, and log-transformed rainfall amounts – and then predict the
values at the locations of the weather stations. This is done by fitting the
following geostatistical model:

$\displaystyle w(\bm{s},t)=\beta_{0}+\beta_{1}z(\bm{s},t)+\xi(\bm{s},t),$

where $w(\bm{s},t)$ is the observed value of the meteorological variable from
the GSM at the spatial grid with centroid $\bm{s}$, $z(\bm{s},t)$ is a known
covariate, and $\xi(\bm{s},t)$ is a spatio-temporal random effect. For
temperature and relative humidity, we used the elevation variable as
covariate. For log-transformed rainfall, we used the relative humidity as the
covariate. We assume that $\xi(\bm{s},t)$ evolves in time as an autoregressive
process of order 1, i.e.,

$\xi(\bm{s},t)=\phi\xi(\bm{s},t-1)+\omega(\bm{s},t),$

where $|\phi|<1$ and $\omega(\bm{s},t)$ is a time-independent Gaussian process
with a Matérn covariance function specified in Equation (22) and which follows
the stationary distribution of the process at time 1. We fit this model using
INLA and the SPDE method rue2009approximate ; lindgren2011explicit . We then
predict the values of $w(\bm{s},t)$ at the locations of the weather stations
and then compare the predicted and the observed values. The resulting scatter
plots in Figure 2 show a general agreement of the two sets of values, but a
clear bias is visible. Figure 2 highlights specific weather stations and it is
apparent, especially for temperature, that there is a spatially-varying
additive bias and a constant multiplicative bias parameter since the GSM
outcomes from a specific weather station seem parallel to the regression line.
For the other two meteorological variables, it is also clear that there is a
spatially varying additive bias, but accounting for spatially-varying or time-
varying multiplicative bias might also be necessary.

Figure 2 also shows that the quality of the GSM outcomes for the three
variables varies. The discrepancy in the outcomes between the GSM and weather
stations for the rainfall data is bigger than for the other two meteorological
variables. The proposed data fusion model, which is discussed in Section 2, is
able to gauge the relative quality of the GSM outcomes for the three
meteorological variables.

PAGASA provided the aggregated data at the monthly level in terms of the mean
temperature, mean relative humidity, and the total precipitation from January
2019 to December 2020. In particular, the GSM outcomes are simulated from an
initial time of 0000 UTC. The reason for using monthly temporal scale is that
we will use the model predictions as an input to an epidemiological model,
particularly with the Dengue fever cases as the response variable, which are
usually available as counts at the monthly level (naish2014climate, ;
abdullah2022association, ). The specific goal in this work, therefore, is to
reconstruct monthly surfaces of meteorological variables, particularly
temperature, relative humidity, and precipitation.

###  1.2 Current data fusion approaches

One statistical model for data fusion called Bayesian melding was proposed by
fuentes2005model  and discussed in lawson2016handbook . The primary motivation
in their work is in modelling the concentration of air pollutants and its
public health implication. The Bayesian melding model is based on the
assumption that all data sources are error-prone realizations of a common
latent spatial process. Suppose $w_{1}(\bm{s})$ is the observed value of the
quantity of interest from a station at a spatial location $\bm{s}$, and
$w_{2}(B)$ is the realized outcome at a grid cell B for the numerical
prediction model. The melding model assumes the following structure:

$\displaystyle x(\bm{s})$ $\displaystyle=\mu(\bm{s})+\xi(\bm{s})$ (1)
$\displaystyle w_{1}(\bm{s})$ $\displaystyle=x(\bm{s})+e(\bm{s})$ (2)
$\displaystyle w_{2}(\bm{s})$
$\displaystyle=\alpha_{0}(\bm{s})+\alpha_{1}(\bm{s})x(\bm{s})+\delta(\bm{s})$
(3) $\displaystyle w_{2}(B)$
$\displaystyle=\dfrac{1}{|B|}\int_{B}w_{2}(\bm{s})d\bm{s}$ (4)

In Equation (1), $x(\bm{s})$ is the latent process of interest, $\mu(\bm{s})$
is the mean of the process which is typically a function of fixed covariates,
and $\xi(\bm{s})$ is a residual component which can be spatially correlated.
Equation (2) links the observed outcomes at the stations $w_{1}(\bm{s})$ and
the latent process $x(\bm{s})$. It follows the classical error model, i.e.,
$e(\bm{s})\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{e})$ where
$e(\bm{s})$ is the error component. Data from the numerical forecast model is
treated as a gridded data where $w_{2}(B)$ is the spatial average over grid
cell B of a conceptual point-referenced model $w_{2}(\bm{s})$ as shown in
Equation (4). Equation (3) specifies the biases in the numerical prediction
model via the model components $\alpha_{0}(\bm{s})$ and $\alpha_{1}(\bm{s})$
which are interpreted as additive and multiplicative biases, respectively. It
also involves another random noise component
$\delta(\bm{s})\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{\delta})$.
The main advantage of this specification is that it allows the two spatially
misaligned data to jointly inform about the latent process $x(\bm{s})$. Due to
model identifiability problems, $\alpha_{0}(\bm{s})$ and $\alpha_{1}(\bm{s})$
are typically parameterized as fixed effects instead of spatial random fields.
In particular, $\alpha_{0}(\bm{s})$ is a polynomial function of $\bm{s}$ while
$\alpha_{1}(\bm{s})$ is an unknown constant.

A similar idea was proposed in (moraga2017geostatistical, ). Given a zero-mean
process $\xi(\bm{s})$ with a stationary covariance function, the model for the
data outcomes $w_{1}(\bm{s})$ and $w_{2}(B)$ is given by

$\displaystyle w_{1}(\bm{s})|\xi(\bm{s})$
$\displaystyle\sim\mathcal{N}\Big{(}\mu(\bm{s})+\xi(\bm{s}),\sigma^{2}_{e}\Big{)}$
(5) $\displaystyle w_{2}(B)$
$\displaystyle=\dfrac{1}{|B|}\int_{B}\Big{(}\mu(\bm{s})+\xi(\bm{s})\Big{)}d\bm{s}.$
(6)

However, the model in Equation (6) does not account for the measurement error
in $w_{2}(B)$. This is also the same model specification in
(zhong2023bayesian, ). Although it uses the same idea as the melding model in
Equations $\eqref{eq:melding_latent}$ to $\eqref{eq:melding_w2s}$ where both
equations have a common latent process, it does not incorporate bias
parameters such as $\alpha_{0}(\bm{s})$ and $\alpha_{1}(\bm{s})$. In a joint
modelling framework, the data coming from $w_{2}(B)$ can dominate the
parameter estimation since there are considerably more outcomes from this data
source compared to $w_{1}(\bm{s})$ (lawson2016handbook, ). The calibration
parameters in Equation (3) impose a restriction on this regard by accounting
for higher measurement error from this data source. These calibration
parameters were accounted for in villejo2023data  and they also proposed the
use of the same estimation strategy as (moraga2017geostatistical, ). A
limitation of their model is that they assumed that both $\alpha_{0}(\bm{s})$
and $\alpha_{1}(\bm{s})$ are constant in space and time. This specification
may not be very flexible in calibrating the numerical forecast model.

forlani2020joint  proposed another data fusion model, but instead of assuming
a common single latent process for the observed outcomes, they assumed several
latent processes which are shared across all the observation data. But these
different latent processes do not have a clear interpretation; in fact, these
processes were simply referred to as spatial random effects. For instance,
suppose $w_{1}(\bm{s})$, $w_{2}(\bm{s})$, and $w_{3}(\bm{s})$ are three data
sources with mean $\mu_{1}(\bm{s})$, $\mu_{2}(\bm{s})$, and $\mu_{3}(\bm{s})$,
respectively, then the proposed model is as follows:

$\displaystyle\mu_{1}(\bm{s})$ $\displaystyle=\beta_{1}+\xi_{1}(\bm{s})$ (7)
$\displaystyle\mu_{2}(\bm{s})$
$\displaystyle=\beta_{2}+\lambda_{2}\xi_{1}(\bm{s})+\xi_{2}(\bm{s})$ (8)
$\displaystyle\mu_{3}(\bm{s})$
$\displaystyle=\beta_{3}+\lambda_{3}\xi_{1}(\bm{s})+\lambda_{4}\xi_{2}(\bm{s})+\xi_{3}(\bm{s})$
(9)

The parameters $\beta_{1}$, $\beta_{2}$, and $\beta_{3}$ are fixed effects
while $\xi_{i}(\bm{s}),i=1,2,3$ are spatial random effects which are shared
among the observational data, and with $\lambda_{j},j=2,3,4$ as unknown
scaling parameters. This approach is also closely related to the so-called
coregionalization model (schmidt2003bayesian, ). The original melding model in
Equations (1) to (3) assumes a single common latent process and has a clear
interpretation as the true process and that the different observation data are
error-prone realizations of the process with varying levels and degree of
measurement error. The model in Equations (7) to (9) deviates from this
general principle.

One bottleneck with the melding model in Equations (1) to (3) is that it
requires considerable computational effort because of the change-of-support
integral in Equation (4). A model which tries to overcome this difficulty was
proposed in (sahu2010fusing, ), given by:

$\displaystyle x(\bm{s})$ $\displaystyle=\tilde{x}(B)+\nu(\bm{s})$ (10)
$\displaystyle w_{1}(\bm{s})$ $\displaystyle=x(\bm{s})+e(\bm{s})$ (11)
$\displaystyle w_{2}(B)$
$\displaystyle=\alpha_{0}+\alpha_{1}\tilde{x}(B)+\psi(B).$ (12)

$\tilde{x}(B)$ is considered as the true areal process and is defined on the
same grid resolution as $w_{2}(B)$. The term $\nu(\bm{s})$ is a Gaussian error
while $\psi(B)$ is a discrete spatial effect which is commonly estimated using
a conditionally autoregressive model. A limitation of this model is that it
considers the underlying true process as discrete and that it can only provide
gridded predictions of the true process.

Another statistical approach for doing data fusion, called the regression
calibration approach, fits a regression model using the outcomes of the
numerical prediction model as a predictor and the observational data as the
response variable. The model is given by

$w_{1}(\bm{s})=\alpha_{0}(\bm{s})+\alpha_{1}(\bm{s})w_{2}(B_{\bm{s}})+e(\bm{s}),\;\;e(\bm{s})\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{e}).$
(13)

Here, $\alpha_{0}(\bm{s})$ and $\alpha_{1}(\bm{s})$ are spatially varying
additive and multiplicative biases, while $w_{2}(B_{\bm{s}})$ is the value of
the numerical forecast model at the grid $B$ which contains the point location
$\bm{s}$. chen2021novel  used this approach to estimate chlorophyll-A
concentration over eutrophic lakes, lee2017rigorous  to model air quality, and
in berrocal2010spatio  to model ozone concentration. Obviously, there is a
computational advantage with this approach since it only uses the values of
$w_{2}(B)$ linked to a stations data $w_{1}(\bm{s})$ to fit the model, and
then uses the full data $w_{2}(B)$ for predictions. A limitation of the
calibration framework is that it assumes that the observed measurements
$w_{1}(\bm{s})$ are the gold standard, even though it is very likely to have
instrumental errors (lawson2016handbook, ). Also, since $w_{2}(B)$ is used as
a predictor, then it cannot contain missing values. However, remote-sensed
data can be missing due to cloud cover and highly reflective surfaces
(lawson2016handbook, ). In addition, the resolution of the predicted latent
surface for the quantity of interest depends on the resolution of $w_{2}(B)$.
(forlani2020joint, ) performed a comparison of the melding model and the
regression calibration approach in a specific data application. The results
show that the Bayesian melding model gave better model predictions compared to
the regression calibration approach.

Unlike the approaches discussed so far, this work accounts for calibration
biases in the data fusion model, specifies a flexible bias formulation,
defines a single and interpretable latent processes for the different data
sources, and specifies a measurement error for the stations data. In
particular, we specify a random field for the additive bias
$\alpha_{0}(\bm{s})$ which we call an error field. Furthermore, we use a
flexible specification by assuming that the error field is a time-varying
random field and that $\alpha_{1}$ is constant in space and time. Finally, our
proposed model treats the outcomes of numerical prediction models as point-
referenced at the centroids of the grids $B$, which is the same strategy as in
(forlani2020joint, ; villejo2023data, ; lee2017rigorous, ). A comparison
between the predictions from using the proposed data fusion model, the
stations-only model, and the regression calibration approach in the real data
application was done using the leave-group-out crossvalidation. We also
conducted a simulation study in order to assess the benefits of doing data
fusion under specific scenarios.

Monte Carlo methods have been the standard approach to doing inference with
Bayesian melding model (falk2010estimating, ; vsevvcikova2011uncertain, ;
liu2011empirical, ; poole2000inference, ). In recent years, the use of the
integrated nested Laplace approximation (INLA) paired with the stochastic
partial differential equations (SPDE) approach has been widely used
(rue2009approximate, ; lindgren2011explicit, ). The INLA method is a
deterministic approach for doing Bayesian inference for latent Gaussian
models. Moreover, the SPDE method is another efficient computational method to
estimate Gaussian fields of the Matérn type (lindgren2011explicit, ). One
challenge in the proposed model, which is discussed in more detail in Section
2.2, is that it involves non-linear components primarily due to the
$\alpha_{1}$ multiplicative bias parameter as shown in Equation 3. To overcome
this issue, we primarily use Bayesian model averaging with INLA for model
inference.

Section 2 discusses the framework and the data fusion model, while section 3
discusses the model estimation approaches. Sections 4 and 5 discusses the
results from a simulation study and the application to real data,
respectively.

##  2 Data Fusion Framework and Model

###  2.1 Framework

We assume that the process of interest is denoted by $x(\bm{s},t)$ which is
the true outcome at a spatial location $\bm{s}$ and time $t$,
$\bm{s}\in\mathcal{S}$, $t=1,\dots,T$. This is observed via two different data
sets. The first one is given by

$\bm{w}_{1t}^{\intercal}=\begin{pmatrix}w_{1}(\bm{s}_{1},t)&w_{1}(\bm{s}_{2},t)&\ldots&w_{1}(\bm{s}_{n_{M}},t)\end{pmatrix},\;\;\;t=1,\ldots,T;\;\;s_{i}\in\mathcal{S},$

which consists of observations from a set of $n_{M}$ stations in locations
$\bm{s}_{i}$, $i=1,\dots,n_{M}$, and provides observations at time points
$t=1,\dots,T$. The second one is given by

$\bm{w}_{2t}^{\intercal}=\begin{pmatrix}w_{2}(\bm{g}_{1},t)&w_{2}(\bm{g}_{2},t)&\ldots&w_{2}(\bm{g}_{n_{G}},t)\end{pmatrix},\;\;\;t=1,\ldots,T;\;\;g_{j}\in\mathcal{S}$

which consists of the gridded output of a numerical model, for example the GSM
model. We assume that $w_{2}(\bm{g}_{j},t)$ is the value at the grid with
centroid $\bm{g}_{j}$, $j=1,\dots,n_{G}$, at time $t=1,\dots,T$. Here we
assume that the two data sources are aligned in time. If this is not the case,
it is always possible to aggregate the higher resolution data. We further
assume that $\bm{w}_{2t}$ has a much wider spatial coverage than $\bm{w}_{1t}$
($n_{M}<<n_{G}$) but with a larger bias, and that both $\bm{w}_{1t}$ and
$\bm{w}_{2t}$ are error-prone realizations of the same process of interest
$\bm{x}_{t}$. This implies that:

$\displaystyle\bm{w}_{1t}$
$\displaystyle=f_{1}(\bm{x}_{t},\bm{\theta}_{1})+\bm{e}_{1t}$ (14)
$\displaystyle\bm{w}_{2t}$
$\displaystyle=f_{2}(\bm{x}_{t},\bm{\theta}_{2})+\bm{e}_{2t},$ (15)

where $f_{1}(\cdot)$ and $f_{2}(\cdot)$ are some deterministic functions of
the process $\bm{x}_{t}$ with bias parameters $\bm{\theta}_{1}$ and
$\bm{\theta}_{2}$. The terms $\bm{e}_{1t}$ and $\bm{e}_{2t}$ are assumed
independent error components. Notice that expanding the above framework to
more than two data sources is straightforward as a new data set would be
treated as yet another error-prone realization of the process of interest.

###  2.2 Proposed Model

Using Equations (14) and (15) to represent the two data sources, we propose
the following data fusion model for the meteorological data:

$\displaystyle x(\bm{s},t)$
$\displaystyle=\bm{\beta}^{\intercal}\bm{z}(\bm{s},t)+\xi(\bm{s},t)$ (16)
$\displaystyle w_{1}(\bm{s}_{i},t)$
$\displaystyle=x(\bm{s}_{i},t)+e_{1}(\bm{s}_{i},t),\;\;\;i=1,\ldots,n_{M},$
(17) $\displaystyle w_{2}(\bm{g}_{j},t)$
$\displaystyle=\alpha_{0}(\bm{g}_{j},t)+\alpha_{1}x(\bm{g}_{j},t)+e_{2}(\bm{g}_{j},t),\;\;\;j=1,\ldots,n_{G}.$
(18)

The process of interest, $x(\bm{s},t)$, is modelled as a linear function of
some known covariates $\bm{z}(\bm{s},t)$ (including an intercept) and a random
field $\xi(\bm{s},t)$. The observed data $w_{1}(\bm{s}_{i},t)$ from the
weather stations are assumed to be unbiased observation of the latent values
$x(\bm{s}_{i},t)$ with an additive error term
$e_{1}(\bm{s}_{i},t)\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{e_{1}}),i=1,\ldots,n_{M}$.
The iid assumption is justified by the fact that the weather stations are
sparsely located in the spatial domain and operate independently of each
other. The data $w_{2}(\bm{g}_{j},t)$ from the numerical weather prediction
model are assumed to be biased observation of the latent values
$x(\bm{g}_{j},t)$, $j=1,\ldots,n_{G}$. We specify both an additive and a
multiplicative bias for $w_{2}(\bm{g}_{j},t)$. While the multiplicative bias
$\alpha_{1}$ is assumed to be constant in time and space, the additive bias or
error field $\alpha_{0}(\bm{g}_{j},t)$ is a random field and also varies in
time. A constant multiplicative bias is justified based on Figure 2. This is
very apparent for the mean temperature since the points for the highlighted
stations are strongly parallel to the line. For the other two meteorological
variables, this is less apparent, and so the assumption of a constant
$\alpha_{1}$ provides an approximation of a possibly spatially-varying or
time-varying $\alpha_{1}$. Finally, the model for $w_{2}(\bm{g}_{j},t)$
contains another unstructured error
$e_{2}(\bm{g}_{j},t)\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{e_{2}})$.

Both the spatio-temporal field $\xi(\cdot,t)$ in the latent process in
Equation (16) and the error field $\alpha_{0}(\cdot,t)$ in Equation (18) are
modelled using a Matérn Gaussian space-time field. In practice, we assume that

$\displaystyle\xi(\bm{s},t)=\phi_{1}\ \xi(\bm{s},t-1)+\omega_{1}(\bm{s},t)$
(19) $\displaystyle\alpha_{0}(\bm{s},t)=\phi_{2}\
\alpha_{0}(\bm{s},t-1)+\omega_{2}(\bm{s},t),$ (20)

where $|\phi_{1}|<1$ and $|\phi_{2}|<1$ model the temporal dependence in an
autoregressive (AR) fashion of order 1 while $\omega_{1}(\bm{s},t)$ and
$\omega_{2}(\bm{s},t)$ are two independent Gaussian innovation processes with
Matérn covariance structure, i.e.,

$\displaystyle\text{Cov}\Big{(}\omega_{h}(\bm{s}_{i},t),\omega_{h}(\bm{s}_{j},u)\Big{)}=\begin{cases}0&t\neq
u\\\ \Sigma^{(h)}_{i,j}&t=u\end{cases},$ (21)
$\displaystyle\Sigma^{(h)}_{i,j}=\frac{\sigma^{2}_{h}}{2^{\nu_{h}-1}\Gamma(\nu_{h})}(\kappa_{h}\left\lVert\bm{s}_{i}-\bm{s}_{j}\right\rVert)^{\nu_{h}}K_{\nu_{h}}(\kappa_{h}(\left\lVert\bm{s}_{i}-\bm{s}_{j}\right\rVert)),$
(22)

for $h=1,2$. Here, $\left\lVert\cdot\right\rVert$ is the Euclidean distance in
$\mathbb{R}^{2}$ between two locations $\bm{s}_{i}$ and $\bm{s}_{j}$, and
$K_{\nu_{h}}(\cdot)$ is the modified Bessel function of the second kind and
order $\nu_{h}>0$. The Matérn field is parameterized by the marginal variance
$\sigma^{2}_{h}$, a scaling parameter $\kappa_{h}$, and a smoothness parameter
$\nu_{h}$ which is related to the mean-square differentiability of the
process. The smoothness parameter is typically fixed at some value because it
is poorly identified in many applications (lindgren2011explicit, ). Moreover,
the scaling parameter is related to the range parameter via the empirically
derived relationship $\rho_{h}\approx\dfrac{\sqrt{8\nu_{h}}}{\kappa_{h}}$.
Given the model structure in Equations (19) to (22), we have
$\bm{\xi}_{t}\sim\mathcal{N}(\phi_{1}\bm{\xi}_{t-1},\bm{\Sigma}^{(1)})$ and
$\bm{\alpha}_{0_{t}}\sim\mathcal{N}(\phi_{2}\bm{\alpha}_{0_{t-1}},\bm{\Sigma}^{(2)})$
for $t=1,\ldots,T$, where
$\bm{\xi}_{t}=\begin{pmatrix}\xi(\bm{s}_{1},t)&\ldots&\xi(\bm{s}_{n_{M}},t)\end{pmatrix}^{\intercal}$,
$\bm{\alpha}_{0_{t}}=\begin{pmatrix}\alpha_{0}(\bm{g}_{1},t)&\ldots&\alpha_{0}(\bm{g}_{n_{G}},t)\end{pmatrix}^{\intercal}$,
and that $\bm{\Sigma}^{(1)}$ and $\bm{\Sigma}^{(2)}$ are dense covariance
matrices whose elements are given in Equation (22). Both $\xi(\bm{s},t)$ and
$\alpha_{0}(\bm{s},t)$ follows the stationary distribution at time 1, i.e.,
$\xi(\bm{s},1)\sim\mathcal{N}\big{(}0,\sigma^{2}_{1}/(1-\phi_{1}^{2})\big{)}$
and
$\alpha_{0}(\bm{g},1)\sim\mathcal{N}\big{(}0,\sigma^{2}_{2}/(1-\phi_{2}^{2})\big{)}$.

##  3 Model Estimation

The estimation approach uses the integrated nested Laplace approximation
(INLA) which is a method for doing deterministic Bayesian inference for
hierarchical models which are latent Gaussian (rue2009approximate, ). The
joint model in Equations (16) to (18) is a latent Gaussian model which becomes
obvious after rewriting it, in vector form, as

$\displaystyle\bm{w}_{1t}$
$\displaystyle=\bm{Z}_{t}\bm{\beta}+\bm{\xi}_{t}+\bm{e}_{1t},\;\;\;\bm{e}_{1t}\sim\mathcal{N}(\bm{0},\sigma^{2}_{e_{1}}\mathbb{I})\;\;\;,t=1,\ldots,T,$
(23) $\displaystyle\;\;\;\bm{\xi}_{t}=\phi_{1}\bm{\xi}_{t-1}+\bm{\omega}_{1t}$
(24) $\displaystyle\bm{w}_{2t}$
$\displaystyle={\bm{\alpha}}_{0_{t}}+\alpha_{1}(\bm{Z}_{t}\bm{\beta}+\bm{\xi}_{t})+\bm{e}_{2t},$
(25)
$\displaystyle\;\;\;\bm{\alpha}_{0_{t}}=\phi_{2}\bm{\alpha}_{0_{t-1}}+\bm{\omega}_{2t},\;\;\;\bm{e}_{2t}\sim\mathcal{N}(\bm{0},\sigma^{2}_{e_{2}}\mathbb{I}).$
(26)

The model specification involves two likelihood components: $\bm{w}_{1t}$ and
$\bm{w}_{2t}$. The latent part of the model includes the fixed effects
$\bm{\beta}$, the space-time effects $\bm{\xi}_{t}$, and the error field
${\bm{\alpha}}_{0_{t}}$, $t=1,\ldots,T$. The fixed effects $\bm{\beta}$ are
given a non-informative prior which is Gaussian with mean zero and a large
variance, while the random fields $\bm{\xi}_{t}$ and ${\bm{\alpha}}_{0_{t}}$
follow a Gaussian autoregressive structure as described in Equations (19) to
(22). The hyperparameters include the parameters of the Matérn field
$\bm{\omega}_{1t}$, particularly $\sigma_{1},\rho_{1},$ and $\phi_{1}$, and
the parameters of the error field field, namely $\sigma_{2},\rho_{2}$, and
$\phi_{2}$. The measurement error variance parameters $\sigma^{2}_{e_{1}}$ and
$\sigma^{2}_{e_{2}}$ are also considered hyperparameters. The aforementioned
hyperparameters are given penalized complexity (PC) priors
fuglstad2019constructing ; simpson2017penalising . PC priors are weakly
informative priors which penalize the complexity of Gaussian random fields by
shrinking the range towards infinity and the marginal variance towards zero
fuglstad2019constructing . These are defined and expressed through probability
statements of the type $\mathbb{P}(\sigma>\sigma_{\text{o}})=\nu_{1}$ and
$\mathbb{P}(\rho<\rho_{\text{o}})=\nu_{2}$, where $\nu_{1},\nu_{2}\in(0,1)$
are probability values chosen by the user, while $\sigma_{\text{o}}$ and
$\rho_{\text{o}}$ are also user-defined values of the parameters. The variance
parameters $\sigma^{2}_{e_{1}}$ and $\sigma^{2}_{e_{2}}$ are also given PC
priors.

A potential issue is that $\alpha_{1}$ is a scaling parameter of an unknown
Gaussian field, $\bm{Z}_{t}\bm{\beta}+\bm{\xi}_{t}$. Although this violates
the latent Gaussian assumption of INLA, this model specification could still
be estimated in the INLA framework by doing data augmentation, which in
consequence considers $\alpha_{1}$ as a model hyperparameter
(muff2015bayesian, ; martins2013bayesian, ; ruiz2012direct, ; villejo2023data,
). However, scaling parameters are typically difficult to estimate and the
model could run into numerical issues in such a scenario. Hence, we explore in
this work the use of a Bayesian model averaging (BMA) approach with INLA
(gomez2020bayesian, ). This approach allows us to fit the data fusion model
conditional on $\alpha_{1}$. When $\alpha_{1}$ is fixed, the proposed data
fusion model in Equations (23) to (26) is latent Gaussian.

Fitting the data fusion model conditional on $\alpha_{1}$ is a viable approach
since we have a sensible set of values for $\alpha_{1}$ which is on a grid
centered at 1. When $\alpha_{1}=1$, this implies that there is no
multiplicative bias in $\bm{w}_{2t}$; and the further the value of
$\alpha_{1}$ is away from 1, the more serious the multiplicative bias in
$\bm{w}_{2t}$.

Suppose all observed data is denoted by $\mathbf{Y}$, a latent model component
is denoted by $x_{j}$, and the model hyperparameters are denoted by
$\bm{\theta}$. INLA computes the posterior marginals based on the following
nested integrals:

$\displaystyle\pi(\theta_{i}|\mathbf{Y})$
$\displaystyle=\int\pi(\bm{\theta}|\mathbf{Y})d\bm{\theta}_{-i}$ (27)
$\displaystyle\pi(x_{j}|\mathbf{Y})$
$\displaystyle=\int\pi(x_{j}|\bm{\theta},\mathbf{Y})\pi(\bm{\theta}|\mathbf{Y})d\bm{\theta}.$
(28)

We express
$\bm{\theta}=\begin{pmatrix}\alpha_{1}&\bm{\theta}_{-\alpha_{1}}\end{pmatrix}^{\intercal}$,
where $\bm{\theta}_{-\alpha_{1}}$ includes all hyperparameters excluding
$\alpha_{1}$. In this regard, the posterior marginals in Equations (27) and
(28) can be expressed as

$\displaystyle\pi(\cdot|\mathbf{Y})$
$\displaystyle=\int\pi(\cdot|\alpha_{1},\mathbf{Y})\pi(\alpha_{1}|\mathbf{Y})d\alpha_{1}.$
(29)

The density $\pi(\cdot|\alpha_{1},\bm{Y})$ is the marginal posterior
conditional on $\alpha_{1}$. This conditional marginal posterior can be
computed using INLA for a fixed $\alpha_{1}$. The density
$\pi(\alpha_{1}|\mathbf{Y})$ in Equation (29) can be expressed as

$\displaystyle\pi(\alpha_{1}|\mathbf{Y})$
$\displaystyle\propto\pi(\mathbf{Y}|\alpha_{1})\pi(\alpha_{1}),$

where $\pi(\mathbf{Y}|\alpha_{1})$ is the conditional marginal likelihood and
$\pi(\alpha_{1})$ is the prior for $\alpha_{1}$. The computation is done by
specifying a grid of values for $\alpha_{1}$, say
$\alpha_{1}^{(k)},k=1,\ldots,K$, and then fitting the INLA conditional on each
$\alpha_{1}^{(k)}$. Given this ensemble of INLA models, the weights for model
averaging are computed as:

$w_{k}=\dfrac{\pi(\mathbf{Y}|\alpha_{1}^{(k)})\pi(\alpha_{1}^{(k)})}{\sum_{k=1}^{K}\pi(\mathbf{Y}|\alpha_{1}^{(k)})\pi(\alpha_{1}^{(k)})}.$
(30)

The marginal posteriors given in Equation (29) are then computed as

$\pi(\cdot|\mathbf{Y})\approx\sum_{k=1}^{K}\pi(\cdot|\alpha_{1}^{(k)},\mathbf{Y})w_{k}.$
(31)

All other posterior quantities of interest and all predictions are computed
using model averaging. A disadvantage of the model averaging approach is that
it requires us to fit the models conditional on each $\alpha_{1}$ value which
can be inefficient especially for large spatio-temporal datasets.

Moreover, we use the stochastic partial differential equations (SPDE) approach
to fit the Matérn fields in the model, particularly, $\bm{\omega}_{1t}$ and
$\bm{\omega}_{2t}$. The SPDE method approximates a Gaussian field with Matérn
covariance function using a Gaussian Markov random field which has a sparse
precision matrix and therefore allows fast model inference
(lindgren2011explicit, ).

##  4 Simulation Study

In this section, we perform a simulation study with the aim to show that data
fusion allows us to recover more accurate predictions of the latent process of
interest and better estimates of model parameters. To simulate our data, we
use the shapefile of the Belo horizonte region in Brazil, available in the R
package spdep (bivand2015comparing, ). This region has a total surface area of
ca 330 $\text{km}^{2}$. We perform the simulation study on a purely spatial
context and simulate the process of interest as:

$x(\bm{s})=\beta_{0}+\beta_{1}z(\bm{s})+\xi(\bm{s})$ (32)

where $\xi(\bm{s})$ is a Matérn random field with the following parameters:
effective range $\rho_{\xi}=2$ degrees, marginal standard deviation
$\sigma_{\xi}=3.16$, and smoothness parameter equal to 1. The value of the
range is such that the spatial correlation becomes negligible at a distance of
ca 222 km which corresponds to half of the maximum distance in the study
region. Moreover, $z(\bm{s})$ is a known covariate and is simulated from a
Matérn process with effective range of 3 degrees, a marginal variance of 1,
and a smoothness parameter equal to 1. The fixed effects are $\beta_{0}=10$
and $\beta_{1}=3$.

The two observed datasets are simulated as follows:

$\displaystyle w_{1}(\bm{s}_{i})$
$\displaystyle=x(\bm{s}_{i})+e_{1}(\bm{s}_{i}),$ $\displaystyle
i=1,\ldots,n_{M}$ (33) $\displaystyle w_{2}(\bm{g}_{j})$
$\displaystyle=\alpha_{0}(\bm{g}_{j})+\alpha_{1}x(\bm{g}_{j})+e_{2}(\bm{g}_{j}),$
$\displaystyle j=1,\ldots,n_{G},$ (34)

wherein the locations of the $n_{M}$ weather stations are randomly distributed
over the study area while the $n_{G}$ numerical model values are observed on a
regular grid. In the data simulation model, we set
$e_{1}(\bm{s}_{i})\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{e_{1}}=.25)$
and
$e_{2}(\bm{g}_{j})\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2}_{e_{2}}=.01)$.
The values of $\sigma^{2}_{e_{1}}$ and $\sigma^{2}_{e_{2}}$ are chosen based
on the empirical results from the temperature model in Section 5.1. This
implies that the noise in $w_{2}(\bm{g}_{j})$ is mainly attributed to the
error field $\alpha_{0}(\bm{g})$. The error field $\alpha_{0}(\bm{g})$ is a
Matérn random field with range $\rho_{\alpha_{0}}=1$, marginal standard
deviation $\sigma_{\alpha_{0}}=1$, and smoothness parameter of 1. The range
parameter and marginal variance of $\xi(\bm{s})$ is chosen to be higher
compared to $\alpha_{0}(\bm{g})$ based on the empirical results from the real
data application discussed in Section 5. Finally, the constant multiplicative
bias parameter is $\alpha_{1}=1.1$.

Figure 3a shows the dense simulation grid. Figure 3b shows a simulated “true”
surface of interest while Figure 3c shows the simulated data from the
numerical forecast model at a coarser resolution than the simulation grid. The
centroids of the grids in Figure 3c is a coarse subset of the points in Figure
3a. The corresponding simulated bias or error field $\alpha_{0}(\bm{g})$ is
shown in Figure 3d. To illustrate the difference in the severity of the biases
between $w(\bm{s}_{i})$ and $w(\bm{g}_{j})$, we generate a specific data
example from Equations (33) and (34), with $n_{M}=10$ stations and with the
$\bm{w}_{2}$ centroids shown in Figure 3c. The point locations of the 10
stations are shown in Figure 5a. Figure 4 shows a comparison of the biases
between the outcomes of $w_{1}(\bm{s}_{i})$ and $w_{2}(\bm{g}_{j})$. The
observed data at the 10 stations closely corresponds to the truth as shown in
Figure 4a. On the other hand, the outcomes of $w_{2}(\bm{g}_{j})$ are more
biased and also overestimates the truth because of the multiplicative bias
parameter $\alpha_{1}>1$ as shown in Figure 4b. Figure 4c shows that there is
a big discrepancy between the outcomes at the stations and the outcomes from
the numerical prediction model at the locations of the stations. This plot is
similar to Figure 2 which shows the discrepancies in the values between the
two data sources for the real data application.

(a)

(b)

(c)

(d)

Figure 3: (a) dense simulation grid, (b) a simulated true field $x(\bm{s})$,
(c) a simulated outcome of the numerical model $w_{2}(\bm{g}_{j})$, (d) a
simulated error field $\alpha_{0}(\bm{g})$

(a)

(b)

(c)

Figure 4: (a) simulated observed values at 10 stations versus true values (b)
simulated outcomes from numerical prediction model versus true values (c)
simulated observed values at 10 stations versus outcomes from numerical
prediction model. The simulated values for $w_{2}(\bm{g}_{j})$ are more biased
than $w_{1}(\bm{s}_{i})$.

Our main interest here is to understand if jointly modeling the two data
sources offers advantages and how these change with the density or sparsity of
the weather stations. We therefore consider three different cases. The first
scenario consists of only $n_{M}=10$ stations and with severe undersampling
(Figure 5a). In the second case, we have $n_{M}=25$ stations, and a large area
that is undersampled (Figure 5b). The third case has $n_{M}=40$ stations which
are more uniformly distributed over the study area (Figure 5c). The location
of the stations are held constant for all the data replicates so that the
configuration of the spatial locations does not influence the results.

(a)

(b)

(c)

Figure 5: Spatial location of stations over the study area: (a) a sparse
network with undersampled regions, (b) a denser network but with an
undersampled region, (c) a dense uniformly distributed network.

###  4.1 Model definition and estimation

We compare two approaches for estimating the latent process of interest. The
first approach uses only the data from the stations and fits the following
model:

$\displaystyle x(\bm{s})$
$\displaystyle=\beta_{0}+\beta_{1}z(\bm{s})+\xi(\bm{s}),$ $\displaystyle
w_{1}(\bm{s}_{i})$ $\displaystyle=x(\bm{s}_{i})+e_{1}(\bm{s}_{i}),\qquad
e_{1}(\bm{s}_{i})\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma_{e_{1}}^{2}),\qquad
i=1,\dots,n_{M}$

The second approach considers the two data sources and fits the following
joint model:

$\displaystyle x(\bm{s})$
$\displaystyle=\beta_{0}+\beta_{1}z(\bm{s})+\xi(\bm{s}),$ $\displaystyle
w_{1}(\bm{s}_{i})$ $\displaystyle=x(\bm{s}_{i})+e_{1}(\bm{s}_{i}),$
$\displaystyle
e_{1}(\bm{s}_{i})\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma_{e_{1}}^{2}),\qquad
i=1,\dots,n_{M}$ $\displaystyle w_{2}(\bm{g}_{j})$
$\displaystyle=\alpha_{0}(\bm{g}_{j})+\alpha_{1}x(\bm{g}_{j})+e_{2}(\bm{g}_{j}),$
$\displaystyle
e_{2}(\bm{g}_{j})\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma_{e_{2}}^{2}),\qquad
j=1,\dots,n_{G}$

where $\xi(\bm{s})$ and $\alpha_{0}(\bm{g}_{j})$ are Matérn random fields with
marginal standard deviations $\sigma_{\xi}$ and $\sigma_{\alpha_{0}}$ and
practical ranges $\rho_{\xi}$ and $\rho_{\alpha_{0}}$, respectively.

To finalize the Bayesian model, we need to specify the priors for all model
parameters. We assign $\beta_{0}$ and $\beta_{1}$ a vague Gaussian prior with
zero mean. We use penalized complexity (PC) priors for the variance parameters
of $e_{1}(\bm{s}_{i})$ and $e_{2}(\bm{g}_{j})$ and the parameters of the two
Matérn fields: $\rho_{\xi}$, $\sigma_{\xi}$, $\rho_{\alpha_{0}}$, and
$\sigma_{\alpha_{0}}$. In this simulation experiment, we set
$\nu_{1}=\nu_{2}=0.5$. For the values assigned to $\sigma_{\text{o}}$ and
$\rho_{\text{o}}$, we have two scenarios. The first scenario, which we
indicate as “matching priors”, uses the actual values used to simulate the
datasets. The second scenario, which we call “non-matching priors”, uses
$\sigma_{{e_{1}}_{\text{o}}}=1.5,\sigma_{{e_{2}}_{\text{o}}}=0.5,{\sigma}_{\xi_{\text{o}}}=1$,
${\rho}_{\xi_{\text{o}}}=0.5$, ${\sigma}_{\alpha_{0\text{o}}}=.5$,
${\rho}_{\alpha_{0\text{o}}}=.5$, which are arbitrarily chosen values. Note
that the matching priors are not necessarily more informative than the non-
matching priors and that both cases are a weakly informative specification of
the priors.

We use the inlabru library to fit the models. The first approach which only
uses the data from the stations is straightforward to fit. For the data fusion
model, we use the Bayesian model averaging approach with INLA as discussed in
Section 3. We define a regular grid of $\alpha_{1}$ values centered on 1.
Defining the grid around 1 is reasonable since $\alpha_{1}=1$ implies that
there is no multiplicative bias in the numerical forecast model. The weights
of the different models for model averaging are computed using Equation (30).
We use a uniform prior for $\alpha_{1}$ in computing the weights.

###  4.2 Model assessment

We use the fitted models to estimate the values of the true process
$x(\bm{s})$ over a set of locations $\bm{s}$ which defines a prediction grid
in the study domain $\mathcal{S}$. We used the simulation grid shown in Figure
3a as the prediction grid. The estimated value for $x(\bm{s})$ is the
posterior mean $\mathbb{E}[x(\bm{s})|\mathbf{Y}]$ with uncertainty given by
the posterior standard deviation $\sqrt{\mathbb{V}[x(\bm{s})|\mathbf{Y}]}$.
Here the vector $\mathbf{Y}$ indicates the observed data used to fit the
model, so either $\mathbf{Y}=\bm{w}_{1}$ or
$\mathbf{Y}=\\{\bm{w}_{1},\bm{w}_{2}\\}$ for the stations-only model and the
data fusion model, respectively.

We use the following metrics for model assessment:

  1. 1.

Average squared error of the estimated field over $\mathcal{S}$:
$\dfrac{1}{|\mathcal{S}|}\mathop{{\int}_{\mathcal{S}}}\Big{(}x(\bm{s})-\mathbb{E}[x(\bm{s})|\bm{Y}]\Big{)}^{2}d\bm{s}$

  2. 2.

Average posterior uncertainty of the estimated field over $\mathcal{S}$:
$\dfrac{1}{|\mathcal{S}|}\mathop{{\int}_{\mathcal{S}}}\sqrt{\mathbb{V}[x(\bm{s})|\mathbf{Y}]}d\bm{s}$

  3. 3.

Average Dawid-Sebastiani (DS) score which is a measure of the closeness
between an observed quantity of interest and the prediction distribution, say
$\mathcal{F}$. The DS score is based on a coherent design criteria and is
appropriate for predictive decision problems, i.e., to predict some future
observation of the process dawid1999coherent . The DS score for a prediction
on $y$ is given by

$\dfrac{\Big{(}y-\mathbb{E}_{\mathcal{F}}(y)\Big{)}^{2}}{\mathbb{V}_{\mathcal{F}}(y)}+\log[\mathbb{V}_{\mathcal{F}}(y)],$
(35)

where $\mathbb{E}_{\mathcal{F}}(y)$ and $\mathbb{V}_{\mathcal{F}}(y)$ is the
mean and variance of the predictive distribution $\mathcal{F}(y)$.

The first metric is a measure of the average discrepancy between the estimated
field and the true values of the latent field, $\bm{x}(\bm{s})$, while the
second metric assesses the average uncertainty in the estimated field. We
approximate both integrals above using the estimated values on the prediction
grid. The third metric is another proper scoring rule which depends on the
predictive mean and variance of the observed data (gneiting2007strictly, ). We
are particularly interested in the predictions for the observational data
$\bm{w}_{1}$. Similar to the first two metrics, a lower value for the average
DS score is preferred. We also look at the accuracy of the parameter estimates
by computing the relative error and the posterior uncertainty. All simulation
results are computed based on 600 independent data replicates.

###  4.3 Simulation study results

Figure 6: Comparison of squared-errors between the stations-only model and
data fusion model for the simulated data in Figure 3 and 4. The errors from
the stations-only model are much higher.
Figure 7: Comparison of the posterior uncertainty between the stations-only
model and data fusion model for the simulated data in Figure 3 and 4. The
posterior uncertainty from the stations-only model are much higher.

(a)

(b)

(c)

Figure 8: Plots of the (a) average squared-errors, (b) average posterior
uncertainty, and (c) average DS score from the 600 simulated datasets with
respect the sparsity of the number of stations, the priors used, and the use
of either stations-only model or data fusion model. The errors and the
posterior uncertainty for the data fusion model is smaller especially for the
case of sparse data from the stations. The DS score is also smaller on the
average for the data fusion model.

Figures 6 and 7 show a comparison of the squared-errors and the posterior
uncertainty in an estimated field, respectively, between the stations-only
model and the data fusion model on the data example in Figure 3 and 4. The
spatial location of the stations are shown as white points in Figures 6 and 7.
The squared errors from the stations-only model, shown in Figure 6, are
generally much larger than the squared errors from the data fusion model. The
posterior uncertainty of the latent field is also much larger for the
stations-only model as shown in Figure 7. Also, as expected, the posterior
uncertainty at the location of the stations is lower than elsewhere, which is
very apparent for the stations-only model. The average squared error for the
stations-only model and data fusion model for this specific case is 0.53 and
0.27, respectively, while the average posterior uncertainty is 0.80 and 0.48,
respectively. The average DS scores are 1.02 and -0.65, respectively.

Figure 8a shows a plot of the average squared errors of the estimated latent
field over $\mathcal{S}$ based on the 600 data replicates for the different
simulation scenarios. The data fusion approach clearly gives smaller average
squared errors, especially when the data on the stations is very sparse. The
same plot also shows a comparison of the average squared errors between the
use of matching and non-matching priors. The results show that there is no
substantial difference with respect to the priors. Figure 8b shows the results
for the average posterior uncertainty. It shows that the data fusion approach
also gives lower uncertainty estimates and that higher number of stations is
associated with smaller uncertainty. In addition, the same figure shows that
there is no difference in the average posterior uncertainty between the use of
matching and non-matching priors. Finally, Figure 8c shows the average DS
scores. The average DS scores are also smaller with the data fusion approach.
Furthermore, it also shows that the average DS scores tend to decrease with
the number of stations, and that the use of either matching or non-matching
priors gives roughly the same values on the average.

Figure 9 shows a comparison of the average relative error and average
posterior uncertainty for two model hyperparameters: the marginal standard
deviation of the spatial field $\xi(\bm{s})$ (in Figure 9a), and the standard
deviation of the measurement error in the stations $e_{1}(\bm{s})$ (in Figure
9b). There is a clear benefit from doing data fusion for all scenarios,
especially when the data from the stations is very sparse. The average
relative errors are generally smaller when matching priors are used. The
average relative error tends to decrease with more stations, especially for
the stations-only model. Figures 19 and 20 in Appendix section 10.1 shows the
results for the fixed effects, $\beta_{0}$ and $\beta_{1}$, and the the range
parameter of the spatial field $\xi(\bm{s})$, respectively. The results show
that the average relative error and average posterior uncertainty are also
smaller for the data fusion approach. Also, the errors and the uncertainty in
the estimates from the stations-only model is higher when the number of
stations is very sparse. For $\rho_{\xi}$, the use of matching priors seem to
give higher relative errors compared to the use of non-matching priors.

Finally, Figure 21 in Appendix section 10.1 shows a summary of the model
averaging weights of the INLA models for each $\alpha_{1}$ value. The weight
for $\alpha_{1}=1.1$ is the highest for all the scenarios. Moreover, the
weights decrease as the value of $\alpha_{1}$ goes further away from 1.1. The
use of either matching or non-matching PC priors does not seem to
substantially affect the weights.

(a)

(b)

Figure 9: Plot of average relative errors and average posterior uncertainty
from 600 simulated datasets for two hyperparameters: (a) marginal standard
deviation of the spatial field, $\sigma_{\xi}$ and (b) standard deviation of
measurement error at the stations, $\sigma_{e_{1}}$.

##  5 Results for Meteorological Data from the Philippines

This section presents in detail the results of the proposed data fusion model
applied to meteorological data in the Philippines. Section 5.1, 5.2, 5.3
discuss the results for temperature, relative humidity, and rainfall,
respectively, while section 5.4 presents the results of the leave-group-out
cross-validation.

###  5.1 Temperature

In modelling temperature, we have the following form for the fixed effects in
the latent process:

$\text{Temperature}(\bm{s},t)=\beta_{0}+\beta_{1}\log\Big{(}\text{Elevation}(\bm{s},t)\Big{)}+\beta_{2}\text{Cool}+\beta_{3}\text{ClimateType}$
(36)

There are two prevailing seasons in the Philippines: rainy or dry season
(PHClimate_PAGASA, ). Rainy season lasts from June to November, while the rest
of the year is generally considered dry. The rainy season is also considered a
hot episode in the country. On the other hand, the dry season is further
categorized into either a cool dry or a hot dry period, where the former lasts
from December to February and the latter from March to May. Hence, we define
the ‘Cool’ variable in Equation (36) as a binary variable which takes a value
of ‘1’ for months December to February and a value of ‘0’ for the other
months. The climate type variable is also a binary variable which takes ‘1’
for the eastern section of the country and ‘0’ for the western section. The
eastern part of the country is characterized as having no distinct wet or dry
season and with relatively high rainfall all year round. On the other hand,
most of the western section of the country are characterized as having
pronounced dry and wet seasons (beck2018present, ).

Penalised complexity (PC) priors are used for the Matérn field parameters
fuglstad2019constructing ; simpson2017penalising . The parameter values for
the Matérn PC priors are as follows:
$\rho_{1\text{o}}=\rho_{2\text{o}}=300\;\text{km},\sigma_{1\text{o}}=0.07,\sigma_{2\text{o}}=0.01$.
The value for the range parameters is one-third the maximum distance of the
spatial domain. The value for $\sigma_{1\text{o}}$ is the standard deviation
of the temperature values, while the value of $\sigma_{2\text{o}}$ is chosen
to be some value smaller than $\sigma_{1\text{o}}$ based on preliminary model
results. The variance parameters of $e_{1}(\bm{s}_{i},t)$ and
$e_{2}(\bm{g}_{j},t)$ are also given PC priors, with
$\sigma_{{e_{1}}_{\text{o}}}=0.2$ and $\sigma_{{e_{2}}_{\text{o}}}=0.01$. The
probability value of all PC priors is set to be equal to 0.50. The rest of the
model parameters are given the default non-informative priors.

We define a grid of values from 0.5 to 1.5 with a length step of 0.1 for the
multiplicative bias parameter, $\alpha_{1}$. A uniform prior was used for
$\alpha_{1}$. An ensemble of INLA models were fitted for a fixed $\alpha_{1}$,
and the BMA weights are computed using Equation 30. The marginal log-
likelihoods $\log\pi(\mathbf{Y}|\alpha_{1}^{(k)})$ and the corresponding BMA
weights $w_{k}$ are shown in Appendix 10.2 Table 4. The results show that the
weight of the model with $\alpha_{1}=1$ is approximately equal to 1, while the
rest of the models do not have a significant contribution in the model
averaging.

Table 1: Posterior estimates of fixed effects for the temperature model – stations-only model versus proposed data fusion model | Stations only | Data fusion
---|---|---
Parameter | Mean | SD | P2.5% | P97.5% | Mean | SD | P2.5% | P97.5%
$\beta_{0}$ | 28.6644 | 2.6296 | 23.5104 | 33.8184 | 28.9186 | 4.6029 | 19.8972 | 37.9401
$\beta_{1}$, log(Elevation) | -0.6308 | 0.0942 | -0.8154 | -0.4462 | -0.7085 | 0.0506 | -0.8076 | -0.6093
$\beta_{2}$, Cool | -0.6830 | 0.1982 | -1.0715 | -0.2946 | -0.6178 | 0.1769 | -0.9645 | -0.2712
$\beta_{3}$, Climate Type | 2.1827 | 0.6991 | 0.8126 | 3.5529 | 0.6062 | 0.3367 | -0.0536 | 1.2661

Shown in Table 1 are the posterior estimates of the fixed effects for the
stations-only model and the proposed data fusion model. The estimated fixed
effects for both models are quite similar, except for the climate type
variable which is significant in the stations-only model but not in the
proposed data fusion model. The uncertainty in the estimates for
$\beta_{1},\beta_{2}$, and $\beta_{3}$ are smaller for the data fusion model
as expected. Elevation is negatively associated with temperature, i.e., it is
colder at higher elevation. The ‘Cool’ binary variable is negatively
associated with temperature, i.e., colder months (December, January, and
February) have lower mean temperature compared to the rest of the months. The
climate type variable is positively associated with temperature, which means
that areas which are located in the eastern section of the country are warmer
on the average than areas in the western section of the country.

The posterior estimates of the hyperparameters from the stations-only model
and the proposed data fusion model are also quite similar (see Table 5 in
Appendix section 10.2). In the data fusion model, the range parameter
$\rho_{1}$ of the spatial field is higher than the range parameter $\rho_{2}$
of the error field. This means that the spatial field $\xi(\bm{s},t)$ is
smoother than the error field $\alpha_{0}(\bm{g},t)$. The spatial correlation
in the temperature spatial field and the error field becomes negligible at a
distance of around 765 km and 113 km, respectively. Also, the marginal
standard deviation $\sigma_{1}$ of the spatial field is much larger than that
of the error field $\sigma_{2}$. The autocorrelation parameters $\phi_{1}$ and
$\phi_{2}$ in both the spatial field and error field are close to 1 which
suggests a high degree of dependence in time.

Figure 10a shows a comparison of the estimated temperature fields for August
2019 between the stations-only model and the proposed data fusion model. The
estimated fields look very similar, although the estimated field from the data
fusion model is smoother. The dark spot in the northern part of the country is
primarily mountainous which makes it cooler than the other parts of the
country. Figure 10b shows the corresponding uncertainty estimates or standard
deviation (SD) in the predicted fields. As expected, the standard deviations
are smaller at the locations of the weather stations (in green points) which
is more apparent for the stations-only model. Also, there is more uncertainty
in the predictions from the stations-only model especially in the islands in
the lower left portion. The average of the posterior standard deviation of the
estimated fields from using the stations-only model and data fusion model is
1.2388 and 0.7711, respectively. Figures 22 in Appendix section 10.2 shows the
estimated spatial fields, $\xi(\bm{s},t)$, for the stations-only model and
data fusion model for the same month. The spatial structure looks quite
similar as the estimated temperature fields in Figures 10a.

(a)

(b)

Figure 10: Comparison of the estimated temperature fields and corresponding
posterior uncertainties for August 2019 between the stations-only model and
the proposed data fusion model. The posterior uncertainties from the proposed
data fusion model are substantially smaller.

Figure 11: Estimated error fields for the temperature model for August 2019.
The estimated error fields at the specific stations correspond to the additive
bias shown in Figure 2.

Figure 11 shows the estimated error field $\hat{\alpha}_{0}(\bm{g}_{j},t)$ for
the GSM data for the same month, August 2019. We can relate this plot with
Figure 2 which shows the discrepancies in the values between the weather
stations and the GSM outcomes. From Figure 11, the estimated error outcomes
around the neighborhood of Baguio station (in blue) is the highest, i.e., the
GSM has overestimated the temperature value in this area. This is consistent
with Figure 2, which shows that the GSM outcomes are larger than the observed
data at the station. Similarly, the estimated error values in the neighborhood
of Malaybalay station (in red) is negative, and is consistent with Figure 2,
which shows that there is a negative bias in the outcomes from the GSM in this
area. For Tanay station (in cyan), the estimated error is around zero which is
also expected since Figure 2 shows that there seems to be no additive bias in
the GSM outcomes.

We also explored the regression calibration approach for doing data fusion
which uses the GSM outcomes as a predictor in the model, as explained in
Equation (13). The model estimates are shown in Appendix 10.2 Table 6. The
corresponding predicted temperature field for August 2019, the posterior
uncertainty, and the estimated spatial field for the same month are shown in
Appendix 10.2 Figure 23. The results from the regression calibration approach
are quite similar to the stations-only model. The proposed data fusion model
has the smallest posterior uncertainty among the three approaches. A better
comparison of the three approaches is discussed in section 5.4 which shows the
cross-validation results.

Figure 12a shows a scatterplot between the observed values at the stations
$w_{1}(\bm{s}_{i},t)$ and the corresponding predicted latent values
$\hat{x}(\bm{s}_{i},t)$ using the proposed data fusion model. The plot shows a
close correspondence between the two sets of values. Figure 12b, which shows a
scatter plot between the observed GSM values $w_{2}(\bm{g}_{j},t)$ and the
corresponding predicted latent values, indicates a strong bias in the GSM
values, with several points that are either overestimated or underestimated.
Finally, Figure 12c shows a scatter plot between the observed GSM data
$w_{2}(\bm{g}_{j},t)$ and the predicted values
$\hat{w}_{2}(\bm{g}_{j},t)=\hat{\alpha}_{0}(\bm{g}_{j},t)+\hat{\alpha}_{1}x(\bm{g}_{j},t)$,
indicating a very close correspondence between the two values. This shows that
the error field can be effectively used to calibrate the GSM outcomes for
temperature.

(a)

(b)

(c)

Figure 12: Plot of observed temperature values versus predicted values using
the proposed data fusion model for (a) weather stations, (b) GSM data, and (c)
calibrated GSM data

###  5.2 Relative Humidity

In modelling relative humidity (RH), we have the following form for the fixed
effects in the latent process:

$\log\Big{(}\text{RH}(\bm{s},t)\Big{)}=\beta_{0}+\beta_{1}\log\text{Temp}(\bm{s},t)+\beta_{2}\Big{(}\log\text{Temp}(\bm{s},t)\Big{)}^{2}+\beta_{3}\log\Big{(}\text{Elevation}(\bm{s},t)\Big{)}+\beta_{4}\text{ClimateType}$
(37)

Elevation and climate type are also used as predictors for the relative
humidity model. In addition, we use log temperature and its quadratic term as
a predictor in the model as recommended by (PHClimate_PAGASA, ). This non-
linear relationship is also established in the atmospheric science literature
(goody1995principles, ). A log transformation on temperature is possible in
the data since the range of its values for the Philippines which is a tropical
country is from 16∘C to 32∘C. This is also done in the literature, such as in
sanders2007temperature . In particular, we use the predictions from the
temperature model in Equation (36), in both the locations of the weather
stations and the point locations of the GSM outcomes, as the input in the
linear predictor in Equation (37).

PC priors are used for the Matérn field parameters. The values used to define
the Matérn PC priors for the two range parameters are the same as the
temperature model ($\rho_{1\text{o}}=300,\rho_{2\text{o}}=50$). For the
marginal standard deviations of the fields, we have $\sigma_{1\text{o}}=0.08$
and $\sigma_{2\text{o}}=0.01$. The former is the standard deviation of the
log-transformed relative humidity values, while the value of
$\sigma_{2\text{o}}$ is chosen to be some value smaller than
$\sigma_{1\text{o}}$. The variance parameters of $e_{1}(\bm{s}_{i},t)$ and
$e_{2}(\bm{g}_{j},t)$ are also given PC priors, with
$\sigma_{{e_{1}}_{\text{o}}}=0.01$ and $\sigma_{{e_{2}}_{\text{o}}}=0.004$.
The probability value in the PC priors are also set equal to 0.5. The rest of
the model parameters are given the default non-informative priors.

From the model averaging results, the model with $\alpha_{1}=1$ also gave the
highest marginal log-likelihood value and with a weight close to 1. The
marginal log-likelihoods $\log\pi(\mathbf{Y}|\alpha_{1}^{(k)})$ and the
corresponding BMA weights $w_{k}$ are shown in Appendix 10.3 Table 7. Table 2
shows the estimates of the fixed effects for the stations-only model and the
proposed data fusion model. The posterior means of the fixed effects from the
two approaches are comparable, although the climate type variable is not
significant in the stations-only model. The estimates of the hyperparameters
are shown in Appendix 10.3 Table 8.

Table 2: Posterior estimates of fixed effects for the relative humidity model - stations-only model versus proposed data fusion model | Stations only | Data fusion
---|---|---
Parameter | Mean | SD | P2.5% | P97.5% | Mean | SD | P2.5% | P97.5%
$\beta_{0}$ | 4.4506 | 0.0299 | 4.3920 | 4.5091 | 4.4505 | 0.0499 | 4.3526 | 4.5484
$\beta_{1}$, log(Temp) | 0.5665 | 0.0480 | 0.4725 | 0.6606 | 0.7900 | 0.0386 | 0.7143 | 0.8656
$\beta_{2}$, log(Temp)2 | -0.1732 | 0.0144 | -0.2014 | -0.1450 | -0.2374 | 0.0115 | -0.2599 | -0.2149
$\beta_{3}$, log(Elevation) | -0.0079 | 0.0029 | -0.0135 | -0.0022 | -0.0131 | 0.0017 | -0.0165 | -0.0098
$\beta_{4}$, Climate Type | 0.0282 | 0.0179 | -0.0070 | 0.0634 | 0.0259 | 0.0085 | 0.0092 | 0.0426

Table 2 shows that there is a significant non-linear relationship between
temperature and relative humidity. There is a positive correlation between the
two variables for lower temperature values, but a negative correlation for
high temperature values. The elevation variable is negatively associated with
relative humidity. Moreover, the climate type variable is positively related
with relative humidity which means that the areas in the eastern section of
the country have higher relative humidity on the average than the western
section. Similar to the results of the temperature data fusion model, the
estimate of the range parameter $\rho_{1}$ of the spatial field in the latent
process $\xi(\bm{s},t)$ is larger than the range parameter $\rho_{2}$ of the
error field $\alpha_{0}(\bm{s},t)$ (see Appendix 10.3 Table 8). The spatial
correlation in the relative humidity values becomes negligible at a distance
of around 589 km, while the spatial correlation in the error field decays to a
negligible value at a distance of 117 km. Moreover, the spatial variance of
the spatial field is larger than that of the error field. The estimates of the
AR parameter are both close to 1, although the estimated value of the
parameter for the spatial field is higher than that of the error field.

Figures 13a and 13b show the estimated relative humidity fields for two
different months, August 2019 and January 2020, using the proposed data fusion
model and the stations-only model, respectively. These two specific months
were chosen since August is characterized as a rainy month while January is
considered a dry month PHClimate_PAGASA . The spatial structure in the
predicted fields between the proposed data fusion model and the stations-only
model is similar, although it is apparent that there is more smoothing in the
estimates from the stations-only model. Figures 13a and 13b show that in the
eastern section of the country, the level of relative humidity is similar for
the two months. In the western section, particularly in the northwestern
section, relative humidity is very high during August, which is a rainy month,
and is very low during January, which is a dry month.

Figure 14 shows the uncertainty estimates of the estimated relative humidity
fields in Figure 13. As expected, there is more uncertainty in the predictions
from the stations-only model. The average of the posterior standard deviation
in the predicted fields from using the stations-only model and the proposed
data fusion model is 0.0394 and 0.0220, respectively. Also, the uncertainty
estimates are smaller at the location of the weather stations shown in green
points.

Figures 24 and 25a in Appendix 10.3 show the estimated spatial field of the
latent process, $\xi(\bm{s},t)$, for the same two months using the stations-
only model and the proposed data fusion model, respectively. The same spatial
patterns observed in the estimated relative humidity field as shown in Figure
13 are also evident in the estimated spatial fields. The spatial effects for
the eastern part of the country are fairly constant for the two different
months, although there is a slight increase in the values for the northeastern
parts during August which is a rainy month. On the other hand, the spatial
effects for the western part are higher during August and lower during
January. Figure 25b in Appendix 10.3 shows the estimated error fields for the
GSM outcomes on relative humidity for the same two months.

A regression calibration modelling approach was also explored for relative
humidity. The model estimates are shown in Appendix 10.3 Table 9. The
estimated relative humidity fields for August 2019 and January 2020 are shown
in Appendix 10.3 Figure 27. The same spatial structure observed from Figure 13
are also observed from the estimated fields using the regression calibration
approach. The corresponding posterior uncertainty and the estimated spatial
fields for the same two months are shown in Appendix 10.3 Figure 28. Similar
to the results for temperature, the estimated relative humidity fields and the
corresponding posterior uncertainty using the regression calibration approach
are quite similar to the stations-only model. The proposed data fusion model
gave the smallest posterior uncertainty among the three approaches. The
results from the cross-validation, which is a better way of doing model
comparison, are discussed in section 5.4.

(a)

(b)

Figure 13: Comparison of estimated relative humidity fields for August 2019
and January 2020 between (a) proposed data fusion model and (b) stations-only
model. The spatial structure in the fields from the two approaches are
similar. There is more smoothing in the estimated fields using the stations-
only model.

(a)

(b)

Figure 14: Posterior uncertainty of the estimated relative humidity fields in
Figure 13. The posterior uncertainty in the estimated field from using the
stations-only model is much higher.

Figure 26a in the Appendix section 10.3 shows a scatterplot between the
observed values at the stations, $w_{1}(\bm{s}_{i},t)$, and the corresponding
predicted values, $\hat{x}(\bm{s}_{i},t)$. There is a strong correspondence
between the two values, although there are slight over- and underestimation.
Figure 26b (Appendix section 10.3) shows the scatter plot between the outcomes
from the GSM, $w_{2}(\bm{g}_{j},t)$, and the predicted latent values,
$\hat{x}(\bm{g}_{j},t)$. The plot shows a clear bias in the GSM outcomes.
Lastly, Figure 26c (Appendix section 10.3) shows a plot of the observed GSM
outcomes, $w_{2}(\bm{g}_{j},t)$, and the corresponding predicted values
$\hat{w}_{2}(\bm{g}_{j},t)=\hat{\alpha}_{0}(\bm{g}_{j},t)+\hat{\alpha}_{1}\hat{x}(\bm{g}_{j},t)$.
The plot shows a very close correspondence between the two values. Similar to
the temperature model, the scatter plot shows that the estimated error fields
can be effectively used to calibrate the GSM outcomes.

###  5.3 Rainfall

In modelling rainfall, the fixed effects in the latent process are as follows:

$\log\Big{(}\text{Rainfall}\Big{)}=\beta_{0}+\beta_{1}\log\text{Temp}(\bm{s},t)+\beta_{2}\Big{(}\log\text{Temp}(\bm{s},t)\Big{)}^{2}+\beta_{3}\text{Season}+\beta_{4}\text{ClimateType}+\beta_{5}\text{ClimateType}\times\text{Season}$
(38)

The ‘Season’ variable is binary which takes a value of ‘1’ for months June to
November (characterized as a wet period), and a value of ‘0’ for the rest of
the months which is the dry season in the country. Just like the relative
humidity model, the log temperature and its squared term are also predictors
in the rainfall model as recommended by PHClimate_PAGASA . We also used the
predictions from the temperature model in Equation (36) as the input in the
linear predictor in Equation (38). An interaction effect between climate type
and season was considered to capture the climate dynamics of the country, also
as recommended by PHClimate_PAGASA .

Similar to the previous models, PC priors are used for the Matérn field
parameters. The same values from previous models are used to define the PC
priors for the range parameters $\rho_{1}$ and $\rho_{2}$. For the marginal
standard deviations, we have $\sigma_{1\text{o}}=1.35$, which is the standard
deviation of the log-transformed rainfall values; and
$\sigma_{2\text{o}}=0.01$, which is the same for the previous models. The
variance parameters of $e_{1}(\bm{s}_{i},t)$ and $e_{2}(\bm{g}_{j},t)$ are
also given PC priors, with $\sigma_{{e_{1}}_{\text{o}}}=0.5$ and
$\sigma_{{e_{2}}_{\text{o}}}=0.26$. The probabily value in the PC priors is
equal to 0.5. The rest of the model parameters are given the default non-
informative priors.

Figure 15: Estimated marginal posterior of $\alpha_{1}$,
$\pi(\alpha_{1}|\mathbf{Y})$, for the rainfall data fusion model. The
posterior mean is 0.6733, while the 95% credible interval estimate is
$(0.5607,0.8353)$.

Shown in Figure 15 is the estimated marginal posterior distribution of
$\alpha_{1}$, $\pi(\alpha_{1}|\mathbf{Y})$. The posterior mean is 0.6733 while
the 95% credible interval estimate is (0.5607, 0.8353). Unlike the previous
two meteorological variables, the multiplicative bias parameter for the GSM
outcomes for rainfall is significantly different from 1, which implies a more
severe bias for rainfall outcomes. This is expected since Figure 2 shows a big
discrepancy between the interpolated GSM outcomes at the locations of the
weather stations and the observed values at the weather stations. This also
agrees with the insights from the LGOCV results discussed in Section 5.4.

Table 3 shows the estimates of the fixed effects for the stations-only model
and the proposed data fusion model. The results show that log temperature has
a non-linear association with log rainfall amounts. The correlation is
positive for small values of log temperature and is negative for large log
temperature values. Moreover, there is a significant interaction between
season and climate type. In particular, the results show that the western
section of the country has a pronounced dry and wet season. On the other hand,
the eastern section of the country has a less pronounced dry and wet season
and with more or less evenly distributed rainfall the whole year. This can be
confirmed in Figures 16a and 16b which shows the predicted log rainfall fields
using the proposed data fusion model and stations-only model, respectively,
for two months - August 2019 (rainy month) and January 2020 (dry month). This
climatic pattern is consistent with theory (beck2018present, ). Also, these
are the same seasonal dynamics observed for relative humidity. The uncertainty
estimates of the fields are shown in Figure 29 in Appendix 10.4. The log
rainfall fields from the two approaches look quite similar but the uncertainty
in the predictions from the proposed data fusion model is substantially
smaller as expected. The estimated spatial fields and error fields for the
same two months using the proposed data fusion model are shown in Figures 31a
and 31b in Appendix 10.4, respectively, while Figure 30 in Appendix 10.4 shows
the estimated spatial fields $\xi(\bm{s},t)$ for the stations-only model.

Table 3: Posterior estimates of fixed effects for the rainfall model - stations-only model versus proposed data fusion model | Stations only | Data fusion
---|---|---
Parameter | Mean | SD | P2.5% | P97.5% | Mean | SD | P2.5% | P97.5%
$\beta_{0}$ | 4.4272 | 0.3598 | 3.7221 | 5.1323 | 4.7589 | 0.3771 | 3.9310 | 5.4842
$\beta_{1}$,  log(Temp) | 2.1863 | 0.4541 | 1.2963 | 3.0763 | 1.6722 | 0.4297 | 0.9110 | 2.5408
$\beta_{2}$,  log(Temp)2 | -0.6989 | 0.1339 | -0.9613 | -0.4366 | -0.5703 | 0.1217 | -0.8090 | -0.3536
$\beta_{3}$,  Season | 0.7948 | 0.3061 | 0.1949 | 1.3948 | 0.4438 | 0.2610 | -0.0196 | 0.9732
$\beta_{4}$,  Climate Type | 1.1831 | 0.1456 | 0.8977 | 1.4686 | 0.6572 | 0.1119 | 0.4577 | 0.8730
$\beta_{5}$,  Climate Type $\times$ Season | -0.8441 | 0.1617 | -1.1611 | -0.5271 | -0.2868 | 0.0994 | -0.4607 | -0.1055

Table 10 in Appendix 10.4 shows the estimates of the hyperparameters for both
the stations-only model and the proposed data fusion model. Similar to the
previous meteorological variables, the estimated range $\rho_{1}$ of the
spatial field in the latent process is larger than the estimated range
$\rho_{2}$ of the error field. This is also true for the estimated marginal
variances of the two fields. The estimated autocorrelation parameters of the
two fields are very different, with the AR parameter for the error field being
much larger.

(a)

(b)

Figure 16: Comparison of estimated log rainfall fields for August 2019 (wet
season) and January 2020 (dry season) between (a) proposed data fusion model
and (b) stations-only model. The figures show that the western section of the
country has a pronounced dry and wet season.

A regression calibration modelling approach was also explored for rainfall.
The model estimates are shown in Appendix 10.4 Table 11. The predicted log
rainfall fields for August 2019 and January 2020 are shown in Appendix 10.4
Figure 33. The same spatial structure observed from Figures 16a and 16b are
also observed from Figure 33. The estimated log rainfall fields and the
corresponding posterior uncertainty using the regression calibration approach
(see Figure 33 and 34 in Appendix 10.4) are similar to the stations-only
model. Although the proposed data fusion model gave the smallest posterior
uncertainty among the three approaches, a better comparison are shown in the
LGOCV results which are discussed in section 5.4.

Figure 32a in Appendix 10.4 shows a scatter plot of the observed log rainfall
outcomes and the predicted values at the weather stations, i.e.,
$w_{1}(\bm{s},t)$ versus $\hat{x}(\bm{s},t)$. The deviation between the
values, especially for small log rainfall values, is large. The scatter plot
between $w_{2}(\bm{g},t)$ and $\hat{x}(\bm{g},t)$, as shown in Figure 32b in
Appendix 10.4, shows the severe bias in the GSM outcomes. The bias is also
higher for small log rainfall values. Finally, Figure 32c in Appendix 10.4
shows a scatter plot of the observed GSM log rainfall data versus the
predicted values calibrated by the error field, i.e., $w_{2}(\bm{g},t)$ versus
$\hat{\alpha}_{0}(\bm{g},t)+\hat{\alpha}_{1}\hat{x}(\bm{g},t)$. The plot shows
that there is a close correspondence between the two.

###  5.4 Leave-group-out cross-validation

This section presents the results of the leave-group-out cross-validation
(LGOCV). The LGOCV is a better alternative than the leave-one-out cross-
validation for structured models which includes multi-level models, time
series models, and spatial models (liu2022leave, ; adin2023automatic, ). We
usually use the spatial models to do an extrapolation task, i.e., to predict
on a spatial location far away from the observations’ locations. The LGOCV
approach is appropriate to evaluate such models since it makes the unobserved
data as represented by the validation set or hold-out set less dependent on
the observed data which is represented by the training set.

liu2022leave  proposed an approximation of the LGOCV for INLA models. They
also proposed an automatic group construction for doing LGOCV. In this work,
we manually define the groups and use a strategy tailored to the application.
When making a prediction at a particular location, which is called the testing
point, we leave-out that point and also all the observations within its
proximity which altogether are referred to as the leave-out set. This is
illustrated in Figure 17. The leave-out set consists of all the points within
the circle, while the training set consists of all the points excluding the
leave-out set. The testing point is technically considered part of the leave-
out set, and then we predict the value at the testing point.

Figure 17: Illustration of the LGOCV approach. We fit the model on the
training set (red) after excluding the leave-out set (blue and green), and
then make predictions on the testing point (green).

We use each weather synoptic station $\bm{s}_{i}$ as a testing point. Suppose
$I_{\bm{s}_{i}}$ is the leave-out set, which consists of the testing point
$\bm{s}_{i}$ and all weather stations located within a given radius from the
testing point. Also, suppose $-I_{\bm{s}_{i}}$ is the training set. The
predicted value at a testing point $\bm{s}_{i}$ is the mean of the posterior
predictive distribution
$\pi\big{(}x(\bm{s}_{i},t)|\mathbf{Y}_{-I_{\bm{s}_{i}}}\big{)}$. To compare
the competing models, we look at different posterior prediction scores. We
primarily look at the Dawid-Sebastiani score whose formula is given in
Equation (35). Since there are several leave-out sets, we compute the mean
Dawid-Sebastiani scores (MDS) for all testing points. We also compute the root
mean squared error (RMSE), mean absolute error (MAE), and the mean absolute
percentage error (MAPE) of the predictions. In addition, we look at the
standard deviation of
$\pi\big{(}x(\bm{s}_{i},t)|\mathbf{Y}_{-I_{\bm{s}_{i}}}\big{)}$ and the
Kullback-Leibler divergence (KLD) between
$\pi\big{(}x(\bm{s}_{i},t)|\mathbf{Y}_{-I_{\bm{s}_{i}}}\big{)}$ and
$\pi\big{(}x(\bm{s}_{i},t)|\mathbf{Y}\big{)}$, where the latter is simply the
posterior predictive distribution given the complete data and the former is
the posterior predictive distribution with the data excluding the leave-out
set. Since there are several leave-out sets, we compute the mean standard
deviation (MSD) and the mean Kullback-Leibler divergence (MKLD) for all
testing points. We look into different sizes of the leave-out set by varying
the radius of the circle around the testing point. We considered four values
for the radius of the leave-out set: 60, 80, 125, and 150 km.

We compare three models: the stations-only model, the regression calibration
approach which is described in Equation (13), and the proposed data fusion
model. Figure 18 shows a comparison of the posterior prediction scores for the
temperature model. The proposed data fusion model has the smallest MDS scores
for all sizes of the leave-out set. Moreover, the MDS scores for the stations-
only model and the regression calibration approach increase as the size of the
leave-out set increases. The proposed data fusion model also has the smallest
RMSE, MAE, and MAPE. The prediction errors for the stations-only model and the
regression calibration model increase with the size of the leave-out set. The
MLKD and MSD for the proposed data fusion model is also the smallest among the
three models. Figure 18 shows that the proposed data fusion model outperforms
the other two approaches, and that the regression calibration approach fares
better than the stations-only model.

Figure 18: Comparison of LGOCV results for temperature from three models:
stations-only model, regression calibration model, and the proposed data
fusion model

Figure 35 in Appendix 10.5 shows a comparison of LGOCV results for the
relative humidity model. The results show that the proposed data fusion model
outperforms the stations-only model and the regression calibration approach
for all the posterior prediction scores considered. The MDS scores for the
proposed data fusion model decrease even as the size of the leave-out set
increases. The difference in the RMSE, MAE, and MAPE among the three
approaches is relatively small for the cases when the size of the leave-out
set is small, but the difference is substantially big with a large leave-out
set. The MKLD and MSD of the proposed data fusion model is also substantially
smaller compared to the other two models.

Finally, Figure 36 in Appendix 10.5 shows a comparison of LGOCV metrics for
the rainfall model. The MDS scores for the proposed data fusion model is
consistently smallest for the different sizes of the leave-out set. The data
fusion approach also has the smallest MAE. In terms of the RMSE and the MAPE,
there is no substantial difference in the scores among the three approaches
when the leave-out set is small, but for the largest leave-out set considered,
the proposed data fusion model outperforms the other two methods. Lastly, the
proposed data fusion model also has the smallest MSD and MKLD for all sizes of
the leave-out sets considered.

The LGOCV results show that the benefits from doing data fusion for
temperature and relative humidity is higher compared to that of rainfall. One
potential reason for this is that the quality of the GSM outcomes for rainfall
is lower compared to the other two meterological variables which is apparent
from Figure 2 and from the model results, particularly with the estimated
value of $\hat{\alpha}_{1}<1$ for the rainfall data fusion model. Nonetheless,
the LGOCV results show that the proposed data fusion model generally
outperforms the other two approaches for the three meteorological variables
considered.

##  6 Conclusions

Data fusion, which aims to combine information from different data sources,
has the potential benefit of improving the accuracy of model estimation and
the quality of model predictions while facing issues with varying qualities of
the data sources (bauer2015quiet, ; gettelman2022future, ; lawson2016handbook,
). This work addresses a data fusion problem motivated by meteorological data
in the Philippines. The proposed model is based on the Bayesian melding model,
which assumes a common latent process for the data sources, but we extend the
model specification from related work in the literature
moraga2017geostatistical ; villejo2023data ; zhong2023bayesian ;
forlani2020joint . In particular, we specify a time-varying random field for
the additive calibration bias for the outcomes from the numerical prediction
model which we call an error field, and a constant multiplicative bias
parameter. The proposed data fusion model defines a single and interpretable
latent process for the data outcomes, specifies a measurement error for all
data sources, gives flexibility in accounting for the different potential
sources of biases, allows for several spatially-misaligned data to inform
about the common latent process, and allows us to gauge the quality of the
different data sources. Although we assumed that the multiplicative bias
parameter $\alpha_{1}$ is constant in this paper, the model can be extended by
assuming that $\alpha_{1}$ is spatially-varying or time-varying. This adds
another layer of model complexity and also additional challenge in the
computation and model estimation, but this could be more relevant for other
applications.

The proposed data fusion model gave better predictions compared to the
stations-only model and the regression calibration approach based on the
leave-group-out cross-validation in the data application. We used the
approximation of the LGOCV in INLA as proposed in liu2022leave , but we
manually defined the leave-out sets in approximating the log predictive
densities and in computing the posterior prediction scores. Based on the mean
Dawid-Sebastiani score, the proposed data fusion model outperforms the other
two approaches. The proposed data fusion model also has the smallest mean KLD
between the predictive distribution from using the full data and that from
excluding the leave-out sets. We also considered other posterior predictions
scores such as the RMSE, MAE, MAPE, and the uncertainty (standard deviations)
in the posterior predictive distributions. The proposed data fusion model
generally outperforms the other two approaches, and that the discrepancy in
the scores among the three approaches increases as the size of the leave-out
set increases.

We conducted a simulation study in order to compare the performance of a data
fusion model and a stations-only model under different levels of sparsity of
the stations data and the specification of the priors of the model
hyperparameters. The results of the simulation study showed that the data
fusion model gave lower average squared-error of the estimated field and lower
average posterior uncertainty, especially when the stations data are sparse.
The data fusion model also gave lower mean Dawid-Sebastiani scores in the
observational data. Lastly, the data fusion model also has smaller relative
errors and smaller posterior uncertainty in the parameter estimates compared
to the stations-only model.

This work considers only two data sources for meteorological data in the
Philippines since these are the two primary data sources available in the
country. However, the current framework can be extended to incorporate
additional data sources such as satellite data. In this case, an additional
likelihood model for the new observation data can be introduced, which will
have its own set of bias parameters and an error field which can also be time-
varying. The proposed modelling approach allows for more than two data
sources, even spatially-misaligned data, to inform about the latent process of
interest.

This work was primarily motivated by the meteorological data problem in the
Philippines, but the proposed data fusion model is also applicable in
different contexts such as in air quality modelling. As an example, in the
United Kingdom, daily concentrations of air pollutants are collected by a
network of monitoring stations called the Automatic Urban and Rural Network
(AURN) lee2017rigorous . Additional data sources for air quality are outcomes
from a weather and chemistry transport model called the Air Quality Unified
Model (AQUM) which is maintained by the Met Office savage2013air , and
outcomes from another dispersion model called the Pollution Climate Mapping
(PCM) model which are run by Ricardo Energy & Environment forlani2020joint ;
DEFRA . The proposed data fusion model in this paper allows one to model the
AURN, AQUM, and PCM jointly in order to come up with good predictions of air
quality, which are important for generating good forecasts of the
concentration of air pollutants and for assessing the impact of these on
public health and safety. As noted by forlani2020joint , most of the data
fusion studies for air quality modelling in the literature only consider one
extra data source at a time. The proposed data fusion model in this work
allows us to use the three aforementioned data sources jointly, while also
allowing us to carefully account for their biases. Also, the proposed model
allows us to gauge the quality of the outcomes from the numerical models by
looking at the magnitude of the error field estimates and the multiplicative
bias parameter. Another important area where the proposed data fusion
framework is applicable is in species distributions modelling. With
technological advancements, there has been a rise in data collection
capabilities to study the natural world, and also an increase in the volume of
data from increased participation of the general public in data collection
august2015emerging ; belmont2024spatio . However, there is also a concern with
regards to the quality of these data, referred to as citizen science data, due
to systemic biases among others van2013opportunistic ; august2015emerging .
The appropriate use of citizen science data together with data from well-
planned surveys falls within the realm of data fusion. The proposed data
fusion framework in this work views these different data sources as realized
outcomes of the same latent process. The proposed data fusion model can then
be extended in this context by assuming some probabilistic structure of the
biases in the citizen science data while borrowing strength from the accuracy
of the outcomes from well-planned surveys.

The INLA and the SPDE approach was used for model inference because they
provide fast and reliable fitting of complex spatio-temporal models
rue2009approximate ; lindgren2011explicit . A potential challenge in the
computational aspect is that the proposed data fusion model involves a non-
linear predictor, particularly the product of the multiplicative bias
parameter $\alpha_{1}$ and a Gaussian field given by
$\bm{Z}_{t}\bm{\beta}_{t}+\bm{\xi}_{t}$. Although the proposed data fusion
model can be estimated using a data augmentation approach martins2013bayesian
; ruiz2012direct , this approach turned out to be less numerically stable.
Hence, we used the model averaging approach so that we can fit the data fusion
model conditional on $\alpha_{1}$. This is a viable approach since we have an
intuitive understanding of the plausible values of this bias parameter.
Another approach for fitting the model is to perform a linearization on the
non-linear predictors using a first-order Taylor approximation, and to
iteratively do this by looking for the optimal linearization point
serafini2023approximation ; bachl2019inlabru . Since the convergence of this
approach and the properties of the approximation depend on the non-linear
nature of the problem, it can also be computationally challenging for some
cases. The use of a model averaging approach successfully removed the
computational challenges due to non-linear predictors in the model, but the
linearized INLA approach as previously described is also a viable approach for
the problem and will be further explored in a future work.

An immediate future work is to use the predicted fields from the data fusion
model as an input in an epidemiological model in order to understand the link
between climate data and health outcomes. This introduces another layer of
spatial misalignment since typically data for the health outcomes are areal
while the predicted fields from the climate data fusion models are point-
referenced. Also, in this two-stage modeling framework, accounting for the
uncertainty in the data fusion model when fitting the health model should be
carefully considered. The problem of uncertainty propagation has been studied
in the context of health modeling blangiardo2016two ; lee2017rigorous ;
gryparis2009measurement  and which generally falls under the area of
measurement error models berry2002bayesian .

##  7 Competing interests

We declare that there is no conflict of interest regarding the publication of
this paper.

##  8 Author contributions statement

S.V.: Conceptualization, Methodology, Formal Analysis, Writing - Original
Draft, Writing - Review & Editing. J.I.: Conceptualization, Methodology,
Investigation, Writing - Review & Editing, Supervision. S.M.:
Conceptualization, Methodology, Investigation, Writing - Review & Editing,
Supervision. F.L.: Methodology, Writing - Review & Editing.

##  9 Acknowledgments

The authors thank the Philippine Atmospheric, Geophysical and Astronomical
Services Administration (PAGASA) office for generously sharing the
meteorological data and for providing additional model insights.

## References

  * [1] Nur Athen Mohd Hardy Abdullah, Nazri Che Dom, Siti Aekball Salleh, Hasber Salim, and Nopadol Precha.  The association between dengue case and climate: A systematic review and meta-analysis.  One Health, page 100452, 2022.
  * [2] A Adin, E Krainski, A Lenzi, Z Liu, J Martínez-Minaya, and H Rue.  Automatic cross-validation in structured models: Is it time to leave out leave-one-out?  arXiv preprint arXiv:2311.17100, 2023.
  * [3] Ali Arab, Monica C Jackson, and Cezar Kongoli.  Modelling the effects of weather and climate on malaria distributions in west africa.  Malaria journal, 13:1–9, 2014.
  * [4] Tom August, Martin Harvey, Paula Lightfoot, David Kilbey, Timos Papadopoulos, and Paul Jepson.  Emerging technologies for biological recording.  Biological Journal of the Linnean Society, 115(3):731–749, 2015.
  * [5] Fabian E Bachl, Finn Lindgren, David L Borchers, and Janine B Illian.  inlabru: an r package for bayesian spatial modelling from ecological survey data.  Methods in Ecology and Evolution, 10(6):760–766, 2019.
  * [6] Peter Bauer, Alan Thorpe, and Gilbert Brunet.  The quiet revolution of numerical weather prediction.  Nature, 525(7567):47–55, 2015.
  * [7] Hylke E Beck, Niklaus E Zimmermann, Tim R McVicar, Noemi Vergopolan, Alexis Berg, and Eric F Wood.  Present and future köppen-geiger climate classification maps at 1-km resolution.  Scientific data, 5(1):1–12, 2018.
  * [8] Jafet Belmont, Sara Martino, Janine Illian, and Håvard Rue.  Spatio-temporal occupancy models with inla.  arXiv preprint arXiv:2403.10680, 2024.
  * [9] Veronica J Berrocal, Alan E Gelfand, and David M Holland.  A spatio-temporal downscaler for output from numerical models.  Journal of agricultural, biological, and environmental statistics, 15:176–197, 2010.
  * [10] Scott M Berry, Raymond J Carroll, and David Ruppert.  Bayesian smoothing and regression splines for measurement error problems.  Journal of the American Statistical Association, 97(457):160–169, 2002.
  * [11] Roger Bivand and Gianfranco Piras.  Comparing implementations of estimation methods for spatial econometrics.  Journal of Statistical Software, 63:1–36, 2015.
  * [12] Marta Blangiardo, Francesco Finazzi, and Michela Cameletti.  Two-stage bayesian model to evaluate the effect of air pollution on chronic respiratory diseases using drug prescriptions.  Spatial and spatio-temporal epidemiology, 18:1–12, 2016.
  * [13] Cheng Chen, Qiuwen Chen, Gang Li, Mengnan He, Jianwei Dong, Hanlu Yan, Zhiyuan Wang, and Zheng Duan.  A novel multi-source data fusion method based on bayesian inference for accurate estimation of chlorophyll-a concentration over eutrophic lakes.  Environmental Modelling & Software, 141:105057, 2021.
  * [14] Lung-Chang Chien and Hwa-Lung Yu.  Impact of meteorological factors on the spatiotemporal patterns of dengue fever incidence.  Environment international, 73:46–56, 2014.
  * [15] A Philip Dawid and Paola Sebastiani.  Coherent dispersion criteria for optimal experimental design.  Annals of Statistics, pages 65–81, 1999.
  * [16] DEFRA.  Uk air information resource.  https://uk-air.defra.gov.uk/research/air-quality-modelling?view=modelling, 2024.  Accessed: 2024-04-04.
  * [17] MG Falk, RJ Denham, and KL Mengersen.  Estimating uncertainty in the revised universal soil loss equation via bayesian melding.  Journal of agricultural, biological, and environmental statistics, 15:20–37, 2010.
  * [18] Chiara Forlani, Samir Bhatt, Michela Cameletti, Elias Krainski, and Marta Blangiardo.  A joint bayesian space–time model to integrate spatially misaligned air pollution data in r-inla.  Environmetrics, 31(8):e2644, 2020.
  * [19] Montserrat Fuentes and Adrian E Raftery.  Model evaluation and spatial interpolation by bayesian combination of observations with outputs from numerical models.  Biometrics, 61(1):36–45, 2005.
  * [20] Geir-Arne Fuglstad, Daniel Simpson, Finn Lindgren, and Håvard Rue.  Constructing priors that penalize the complexity of gaussian random fields.  Journal of the American Statistical Association, 114(525):445–452, 2019.
  * [21] Andrew Gettelman, Alan J Geer, Richard M Forbes, Greg R Carmichael, Graham Feingold, Derek J Posselt, Graeme L Stephens, Susan C van den Heever, Adam C Varble, and Paquita Zuidema.  The future of earth system prediction: Advances in model-data fusion.  Science Advances, 8(14):eabn3488, 2022.
  * [22] Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  Journal of the American statistical Association, 102(477):359–378, 2007.
  * [23] Virgilio Gómez-Rubio, Roger S Bivand, and Håvard Rue.  Bayesian model averaging with the integrated nested laplace approximation.  Econometrics, 8(2):23, 2020.
  * [24] Richard Goody.  Principles of atmospheric physics and chemistry.  1995\.
  * [25] Sonja Greven, Francesca Dominici, and Scott Zeger.  An approach to the estimation of chronic air pollution effects using spatio-temporal information.  Journal of the American Statistical Association, 106(494):396–406, 2011.
  * [26] Alexandros Gryparis, Christopher J Paciorek, Ariana Zeka, Joel Schwartz, and Brent A Coull.  Measurement error caused by spatial misalignment in environmental epidemiology.  Biostatistics, 10(2):258–274, 2009.
  * [27] Darrel Hess and Dennis Tasa.  McKnight’s physical geography: a landscape appreciation.  Prentice Hall Upper Saddle River, NJ, 2011.
  * [28] I Gede Nyoman Mindra Jaya and Henk Folmer.  Spatiotemporal high-resolution prediction and mapping: methodology and application to dengue disease.  Journal of geographical systems, 24(4):527–581, 2022.
  * [29] Andrew B Lawson, Sudipto Banerjee, Robert P Haining, and María Dolores Ugarte.  Handbook of spatial epidemiology.  CRC press, 2016.
  * [30] A Lee, A Szpiro, SY Kim, and L Sheppard.  Impact of preferential sampling on exposure prediction and health effect inference in the context of air pollution epidemiology.  Environmetrics, 26(4):255–267, 2015.
  * [31] Duncan Lee, Sabyasachi Mukhopadhyay, Alastair Rushworth, and Sujit K Sahu.  A rigorous statistical framework for spatio-temporal pollution prediction and estimation of its long-term impact on health.  Biostatistics, 18(2):370–385, 2017.
  * [32] Finn Lindgren, Håvard Rue, and Johan Lindström.  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, 2011.
  * [33] Zhedong Liu and Haavard Rue.  Leave-group-out cross-validation for latent gaussian models.  arXiv preprint arXiv:2210.04482, 2022.
  * [34] Zhong Liu, Nhu D Le, and James V Zidek.  An empirical assessment of bayesian melding for mapping ozone pollution.  Environmetrics, 22(3):340–353, 2011.
  * [35] Thiago G Martins, Daniel Simpson, Finn Lindgren, and Håvard Rue.  Bayesian computing with inla: new features.  Computational Statistics & Data Analysis, 67:68–83, 2013.
  * [36] Paula Moraga, Susanna M Cramb, Kerrie L Mengersen, and Marcello Pagano.  A geostatistical model for combined analysis of point-level and area-level data using inla and spde.  Spatial Statistics, 21:27–41, 2017.
  * [37] Stefanie Muff, Andrea Riebler, Leonhard Held, Håvard Rue, and Philippe Saner.  Bayesian analysis of measurement error models using integrated nested laplace approximations.  Journal of the Royal Statistical Society Series C: Applied Statistics, 64(2):231–252, 2015.
  * [38] Suchithra Naish, Pat Dale, John S Mackenzie, John McBride, Kerrie Mengersen, and Shilu Tong.  Climate change and dengue: a critical and systematic review of quantitative modelling approaches.  BMC infectious diseases, 14(1):1–14, 2014.
  * [39] PAGASA.  Climate of the philippines.  https://www.pagasa.dost.gov.ph/information/climate-philippines, 2023.  Accessed: 2023-06-12.
  * [40] David Poole and Adrian E Raftery.  Inference for deterministic simulation models: the bayesian melding approach.  Journal of the American Statistical Association, 95(452):1244–1255, 2000.
  * [41] Håvard Rue, Sara Martino, and Nicolas Chopin.  Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations.  Journal of the royal statistical society: Series b (statistical methodology), 71(2):319–392, 2009.
  * [42] Ramiro Ruiz-Cárdenas, Elias T Krainski, and Håvard Rue.  Direct fitting of dynamic models using integrated nested laplace approximations.  Computational Statistics & Data Analysis, 56(6):1808–1828, 2012.
  * [43] Sujit K Sahu, Alan E Gelfand, and David M Holland.  Fusing point and areal level space–time data with application to wet deposition.  Journal of the Royal Statistical Society Series C: Applied Statistics, 59(1):77–103, 2010.
  * [44] Nathan J Sanders, Jean-Philippe Lessard, Matthew C Fitzpatrick, and Robert R Dunn.  Temperature, but not productivity or geometry, predicts elevational diversity gradients in ants across spatial grains.  Global Ecology and Biogeography, 16(5):640–649, 2007.
  * [45] NH Savage, P Agnew, LS Davis, C Ordóñez, R Thorpe, CE Johnson, FM O’Connor, and M Dalvi.  Air quality modelling using the met office unified model (aqum os24-26): model description and initial evaluation.  Geoscientific Model Development, 6(2):353–372, 2013.
  * [46] Alexandra M Schmidt and Alan E Gelfand.  A bayesian coregionalization approach for multivariate pollutant data.  Journal of Geophysical Research: Atmospheres, 108(D24), 2003.
  * [47] Francesco Serafini, Finn Lindgren, and Mark Naylor.  Approximation of bayesian hawkes process with inlabru.  Environmetrics, 34(5):e2798, 2023.
  * [48] Hana Ševčíková, Adrian E Raftery, and Paul A Waddell.  Uncertain benefits: Application of bayesian melding to the alaskan way viaduct in seattle.  Transportation Research Part A: Policy and Practice, 45(6):540–553, 2011.
  * [49] Daniel Simpson, Håvard Rue, Andrea Riebler, Thiago G Martins, and Sigrunn H Sørbye.  Penalising model component complexity: A principled, practical approach to constructing priors.  2017\.
  * [50] Arco J Van Strien, Chris AM Van Swaay, and Tim Termaat.  Opportunistic citizen science data of animal species produce reliable estimates of distribution trends if analysed with occupancy models.  Journal of Applied Ecology, 50(6):1450–1458, 2013.
  * [51] Stephen Jun Villejo, Janine B Illian, and Ben Swallow.  Data fusion in a two-stage spatio-temporal model using the inla-spde approach.  Spatial Statistics, 54:100744, 2023.
  * [52] Ruiman Zhong and Paula Moraga.  Bayesian hierarchical models for the combination of spatially misaligned data: A comparison of melding and downscaler approaches using inla and spde.  Journal of Agricultural, Biological and Environmental Statistics, pages 1–20, 2023.

##  10 Appendix

###  10.1 Simulation Results

(a) Plot of relative errors and posterior uncertainty in $\beta_{0}$

(b) Plot of relative errors and posterior uncertainty in $\beta_{1}$

Figure 19: Plot of average relative errors and average posterior uncertainty
for the fixed effects: (a) $\beta_{0}$ and (b) $\beta_{1}$ Figure 20: Plot of
average relative errors and average posterior uncertainty for $\rho_{\xi}$
Figure 21: Average model averaging weights for different $\alpha_{1}$ values
in fitting the proposed data fusion model with respect to the sparsity of the
stations data and the priors used. The correct value of $\alpha_{1}$ has the
highest weight.

###  10.2 Temperature model

Table 4: Marginal log-likelihoods conditional an $\alpha_{1}$ and the
corresponding BMA weights for the temperature data fusion model

$\alpha_{1}$ | $\log\pi(\mathbf{Y}|\alpha_{1})$ | $w_{k}$
---|---|---
0.5 | -978.2951 | 0.0000
0.6 | -914.7904 | 0.0000
0.7 | -849.6732 | 0.0000
0.8 | -778.4929 | 0.0000
0.9 | -697.5013 | 0.0001
1 | -688.1419 | 1.0000
1.1 | -811.9266 | 0.0000
1.2 | -899.7620 | 0.0000
1.3 | -949.7187 | 0.0000
1.4 | -2265.0736 | 0.0000
1.5 | -2329.8480 | 0.0000

Table 5: Posterior estimates of hyperparameters for the temperature model -
stations-only model versus proposed data fusion model

| Stations only | Data fusion
---|---|---
Parameter | Mean | SD | P2.5% | P97.5% | Mean | SD | P2.5% | P97.5%
$\sigma_{e_{1}}$ | 0.1781 | 0.0106 | 0.1581 | 0.1967 | 0.2429 | 0.0103 | 0.2241 | 0.2636
$\sigma_{e_{2}}$ | - | - | - | - | 0.0221 | 0.0068 | 0.0145 | 0.0430
$\rho_{1}$ | 621.8876 | 50.6769 | 528.7427 | 728.1444 | 764.7483 | 59.3223 | 655.5974 | 888.9635
$\sigma_{1}$ | 5.1208 | 0.6116 | 4.0336 | 6.4384 | 7.6904 | 0.8720 | 6.1331 | 9.5545
$\phi_{1}$ | 0.9920 | 0.0023 | 0.9879 | 0.9951 | 0.9979 | 0.0005 | 0.9968 | 0.9987
$\rho_{2}$ | - | - | - | - | 112.7677 | 8.0757 | 97.7734 | 129.5544
$\sigma_{2}$ | - | - | - | - | 0.6681 | 0.0662 | 0.5465 | 0.8065
$\phi_{2}$ | - | - | - | - | 0.9371 | 0.0137 | 0.9060 | 0.9596

Figure 22: Comparison of estimated spatial fields $\xi(\bm{s},t)$ for August
2019 between the stations-only model versus proposed data fusion model. The
estimated spatial fields are roughly similar.

(a)

(b)

(c)

Figure 23: Model results from the regression calibration approach for
temperature: (a) estimated temperature field (b) posterior standard deviation
(c) estimated spatial field Table 6: Posterior estimates for the temperature
model - regression calibration approach

Parameter | Mean | SD | P2.5 | P97.5
---|---|---|---|---
$\beta_{0}$ | 6.7283 | 5.7536 | -4.5485 | 18.0051
$\beta_{1}$, log(Elevation) | -0.4324 | 0.0958 | -0.6202 | -0.2447
$\beta_{2}$, Cool | -0.6835 | 0.1994 | -1.0743 | -0.2927
$\beta_{3}$, Climate Type | 2.1489 | 0.6177 | 0.9382 | 3.3595
$\beta_{4}$, GSM | 0.8293 | 0.1982 | 0.4409 | 1.2178
$\sigma_{e_{1}}$ | 0.1824 | 0.01078 | 0.1616 | 0.2017
$\rho_{1}$ | 625.5146 | 51.2126 | 531.2444 | 732.7520
$\sigma_{1}$ | 4.5273 | 0.5445 | 3.5349 | 5.6755
$\rho_{1}$ | 0.9898 | 0.0024 | 0.9842 | 0.9936

###  10.3 Relative humidity model

Table 7: Marginal log-likelihoods conditional an $\alpha_{1}$ and the
corresponding BMA weights for the relative humidity data fusion model

$\alpha_{1}$ | $\log\pi(\mathbf{Y}|\alpha_{1})$ | $w_{k}$
---|---|---
0.5 | 5955.3408 | 0.0000
0.6 | 5830.2520 | 0.0000
0.7 | 6076.5968 | 0.0000
0.8 | 6136.2539 | 0.0000
0.9 | 5320.3512 | 0.0000
1 | 6248.2818 | 1
1.1 | 6116.9061 | 0.0000
1.2 | 6101.2002 | 0.0000
1.3 | 6071.8142 | 0.0000
1.4 | 6044.1221 | 0.0000
1.5 | 5847.2576 | 0.0000

Table 8: Posterior estimates of hyperparameters for the relative humidity
model - stations-only model versus proposed data fusion model

| Stations only | Data fusion
---|---|---
Parameter | Mean | SD | P2.5% | P97.5% | Mean | SD | P2.5% | P97.5%
$\sigma_{e_{1}}$ | 0.0123 | 0.0008 | 0.0109 | 0.0137 | 0.0199 | 0.0010 | 0.0181 | 0.0218
$\sigma_{e_{2}}$ | - | - | - | - | 0.0029 | 0.0008 | 0.0016 | 0.0045
$\rho_{1}$ | 287.5769 | 26.9517 | 238.1271 | 344.1479 | 589.1127 | 67.9761 | 468.3741 | 735.5712
$\sigma_{1}$ | 0.0872 | 0.0078 | 0.0729 | 0.1034 | 0.1112 | 0.0141 | 0.0866 | 0.1419
$\phi_{1}$ | 0.9282 | 0.0121 | 0.9018 | 0.9491 | 0.9700 | 0.0079 | 0.9522 | 0.9829
$\rho_{2}$ | - | - | - | - | 117.2557 | 8.9556 | 100.4029 | 135.6345
$\sigma_{2}$ | - | - | - | - | 0.0404 | 0.0020 | 0.0366 | 0.0446
$\phi_{2}$ | - | - | - | - | 0.8553 | 0.0152 | 0.8235 | 0.8834

Figure 24: Estimated spatial fields for log relative humidity, August 2019 and
January 2020, using the stations-only model

(a)

(b)

Figure 25: Estimated spatial fields for log relative humidity (left) and
estimated error fields for the GSM log relative humidity data (right), August
2019 and January 2020, using the proposed data fusion model

(a)

(b)

(c)

Figure 26: Plot of observed relative humidity alues versus predicted values
using the proposed data fusion model for (a) weather stations and (b) GSM
data, and (c) calibrated GSM data Figure 27: Estimated relative humidity
fields for August 2019 and January 2020 using the regression calibration
approach

(a)

(b)

Figure 28: Estimated posterior uncertainty for log relative humidity field
(left) and estimated spatial field for relative humidity (right), August 2019
and January 2020, using the regression calibration approach Table 9: Posterior
estimates for the relative humidity model - regression calibration approach

Parameter | Mean | SD | P2.5 | P97.5
---|---|---|---|---
$\beta_{0}$ | 3.7515 | 1.2240 | 1.3526 | 6.1505
$\beta_{1}$, log(Temp) | 0.5608 | 0.0490 | 0.4647 | 0.6568
$\beta_{2}$, log(Temp)2 | -0.1714 | 0.0147 | -0.2002 | -0.1425
$\beta_{3}$, log(Elevation) | -0.0083 | 0.0030 | -0.0142 | -0.0024
$\beta_{4}$, Climate Type | 0.0252 | 0.0187 | -0.0116 | 0.0619
$\beta_{5}$, GSM | 0.1585 | 0.2775 | -0.3853 | 0.7023
$\sigma_{e_{1}}$ | 0.0124 | 0.001 | 0.0108 | 0.0139
$\rho_{1}$ | 288.0160 | 26.9687 | 238.8441 | 344.9359
$\sigma_{1}$ | 0.0874 | 0.0078 | 0.0731 | 0.1039
$\rho_{1}$ | 0.9287 | 0.0120 | 0.9024 | 0.9497

###  10.4 Rainfall

Table 10: Posterior estimates of hyperparameters for the rainfall model -
stations-only model versus proposed data fusion model

| Stations only | Data fusion
---|---|---
Parameter | Mean | SD | P2.5% | P97.5% | Mean | SD | P2.5% | P97.5%
$\sigma_{e_{1}}$ | 0.4817 | 0.0193 | 0.4459 | 0.5208 | 0.5006 | 0.0218 | 0.4662 | 0.5493
$\sigma_{e_{2}}$ | - | - | - | - | 0.2560 | 0.0139 | 0.2290 | 0.2840
$\rho_{1}$ | 614.6165 | 62.4983 | 501.6406 | 747.4644 | 584.2587 | 59.8174 | 480.0420 | 711.3533
$\sigma_{1}$ | 1.1130 | 0.0737 | 0.9759 | 1.2658 | 1.1066 | 0.0650 | 0.9780 | 1.2265
$\phi_{1}$ | 0.6014 | 0.0484 | 0.5014 | 0.6917 | 0.6913 | 0.0369 | 0.6107 | 0.7542
$\rho_{2}$ |  |  |  |  | 434.9348 | 66.3488 | 323.0664 | 588.1861
$\sigma_{2}$ | - | - | - | - | 0.9418 | 0.0953 | 0.7638 | 1.1157
$\phi_{2}$ | - | - | - | - | 0.8696 | 0.0225 | 0.8201 | 0.9081

(a)

(b)

Figure 29: Posterior uncertainty of the estimated log rainfall fields in
Figure 16. The posterior uncertainty in the estimated fields from using the
stations-only model is higher. Figure 30: Estimated spatial fields for log
rainfall, August 2019 and January 2020, using the stations-only model

(a)

(b)

Figure 31: Estimated spatial fields for log rainfall (left) and estimated
error fields for the GSM log rainfall data (right), August 2019 and January
2020, using the proposed data fusion model

(a)

(b)

(c)

Figure 32: Plot of observed log rainfall values versus predicted values using
the proposed data fusion model: (a) weather stations, (b) GSM data, (c)
calibrated GSM data Figure 33: Estimated log rainfall fields for August 2019
and January 2020 using the regression calibration approach

(a)

(b)

Figure 34: Estimated posterior uncertainty for log rainfall field (left) and
estimated spatial fields for log rainfall (right), August 2019 and January
2020, using the regression calibration approach Table 11: Posterior estimates
for the rainfall model - regression calibration approach

Parameter | Mean | SD | P2.5 | P97.5
---|---|---|---|---
$\beta_{0}$ | 6.2004 | 1.0181 | 4.2051 | 8.1958
$\beta_{1}$,  log(Temp) | 1.5734 | 0.5533 | 0.4891 | 2.6578
$\beta_{2}$,  log(Temp)2 | -0.5133 | 0.1647 | -0.8361 | -0.1905
$\beta_{3}$,  Season | 0.8063 | 0.3074 | 0.2038 | 1.4088
$\beta_{4}$,  Climate Type | 1.1589 | 0.1448 | 0.8752 | 1.4426
$\beta_{5}$,  Climate Type $\times$ Season | -0.8526 | 0.1612 | -1.1685 | -0.5367
$\beta_{6}$ | -0.0687 | 0.0367 | -0.1407 | 0.0033
$\sigma_{e_{1}}$ | 4.3337 | 0.3722 | 3.6447 | 5.1091
$\rho_{1}$ | 619.7900 | 63.1403 | 506.0104 | 754.3626
$\sigma_{1}$ | 1.1030 | 0.0727 | 0.9678 | 1.2539
$\rho_{1}$ | 0.5909 | 0.0494 | 0.4887 | 0.6829

###  10.5 LGOCV

Figure 35: Comparison of LGOCV results for relative humidity from three
models: stations-only model, regression calibration model, and the proposed
data fusion model Figure 36: Comparison of LGOCV results for rainfall from
three models: stations-only model, regression calibration model, and the
proposed data fusion model