# Adjusting Posteriors from Composite and Misspecified Likelihoods with
Application to Spatial Extremes in R-INLA

Silius M. Vandeskog
Department of Mathematics, The Norwegian University of Science and Technology
(NTNU)
and
Sara Martino
Department of Mathematics, The Norwegian University of Science and Technology
(NTNU)
and
Raphaël Huser
Statistics program, CEMSE Division, King Abdullah University of Science and
Technology (KAUST)

###### Abstract

Methods such as low-rank approximations, covariance tapering and composite
likelihoods are useful for speeding up inference in settings where the true
likelihood is unknown or intractable to compute. However, such methods can
lead to considerable model misspecification, which is particularly challenging
to deal with and can provide misleading results when performing Bayesian
inference. Adjustment methods for improving properties of posteriors based on
composite and/or misspecified likelihoods have been developed, but these
cannot be used in conjunction with certain computationally efficient inference
software, such as the R-INLA package. We develop a novel adjustment method
aimed at postprocessing posterior distributions to improve the properties of
their credible intervals. The adjustment method can be used in conjunction
with any software for Bayesian inference that allows for sampling from the
posterior. Several simulation studies demonstrate that the adjustment method
performs well in a variety of different settings. The method is also applied
for performing improved modelling of extreme hourly precipitation from high-
resolution radar data in Norway, using the spatial conditional extremes model
and R-INLA. The results show a clear improvement in performance for all model
fits after applying the posterior adjustment method.

Keywords: Model misspecification, Composite likelihoods, Bayesian inference,
Spatial conditional extremes

##  1 Introduction

Bayesian modelling of high-dimensional spatial and/or temporal data has
attracted much attention in recent decades, see, e.g., [11], [1] and [14] for
a non-extensive list of methods and applications. Likelihood-based inference
for such high-dimensional models is challenging and often computationally
intractable. When inference is intractable the modeller has two choices:
Either switch to a more efficient inference method or adjust the model into
one with similar properties but lower computational cost. The first choice is
often preferred, as typically the model was chosen for a good reason. Methods
that help speed up Bayesian inference include approximate Bayesian computation
(ABC) methods [35, 47], Max-and-Smooth methods [19, 23], improved Markov chain
Monte Carlo (MCMC) methods such as Hamiltonian Monte Carlo, with efficient
implementations in the Stan software [5], and the integrated nested Laplace
approximation [41, INLA;] with efficient implementations in the R-INLA package
for the R programming language [42]. INLA has become very popular thanks to
its combination of computational efficiency and high accuracy. Based on
deterministic numerical approximations of the posterior distribution, the
method can be applied for inference on latent Gaussian models, which are
popular and commonly used in a large variety of applications, including high-
dimensional spatial and/or temporal modelling [45, 33, 36].

When inference is intractable even with improved inference methods like the
aforementioned, it is necessary to change the model itself. As an example,
high-dimensional spatial and/or temporal phenomena are often modelled using
(possibly transformed) Gaussian processes, not necessarily because they are
believed to be exactly Gaussian, but because Gaussian processes are a rather
flexible class of models that are computationally efficient with appealing
theoretical properties [1, e.g.,]. Gaussian processes are even applied in
settings where they are known in advance to provide a misspecified
representation of the true data generating process, because alternative models
may be too complex or computationally intensive to fit [49, 9]. However,
inference for Gaussian processes typically requires computing the inverse of
the covariance matrix, whose cost scales cubicly with the model dimension.
Thus, in truly high-dimensional problems, even inference based on Gaussian
processes becomes intractable.

A common approach for further speeding up inference is to enforce (realistic
or artificial) assumptions of conditional or unconditional independence that
lead to simplified likelihood expressions. Covariance tapering, low-rank
estimation and composite likelihoods are all based on this approach.
Covariance tapering [24] relies on the working assumption that random
variables become independent beyond a specified distance. This leads to sparse
covariance matrices that are easier to invert. Low-rank estimation is based on
approximating random processes with linear combinations of lower-dimensional
processes, leading to low-rank covariance matrices that are easier to deal
with computationally. The list of low-rank approximation methods includes
predictive process models [2], fixed-rank Kriging [8] and models based on
stochastic partial differential equations [30, SPDE;]. Composite likelihoods
are functions composed of many lower-dimensional likelihoods [31], which may
completely remove the need to invert any high-dimensional covariance matrices.
The class of composite likelihoods includes the independence likelihood [6],
the pairwise likelihood [7], the block likelihood [12] and the Vecchia
approximation [51], among others. A composite likelihood must be treated
carefully as, in general, it does not correspond to the likelihood of a proper
statistical model. However, the theory on composite likelihoods is well
developed, and corresponding maximum composite likelihood estimators have been
shown to perform well both in theory and in practice [50].

A common challenge with all methods above is that they often introduce large
amounts of model misspecification through unrealistic distributional
assumptions. While composite likelihoods are not valid likelihood functions,
they have a lot in common with misspecified likelihoods, as shown in Section
2, and in this paper we consider composite likelihoods as a type of
misspecified likelihoods. Misspecified likelihoods complicate Bayesian
inference and can result in posterior distributions with poor frequency
properties, meaning that, as the sample size increases and the prior loses all
influence, $(1-\alpha)$-credible intervals do not, in general, coincide with
$(1-\alpha)$-confidence intervals [39, 27]. Adjustments for improving the
properties of posterior credible intervals have been proposed by, e.g., [6],
[38], [39], and [48]. Unfortunately, all such adjustments require specific
modifications of the likelihood function. Therefore, they usually cannot be
used in conjunction with optimised software for Bayesian inference, such as
R-INLA, where the user must choose from a predefined list of likelihood
functions that cannot be modified. In this paper, we develop a novel method
for improving the frequency properties of a posterior distribution, based on
postprocessing posterior samples. This method can be applied to a broad class
of computational software, including R-INLA.

We apply our method to the spatial conditional extremes model [52], and use it
to study spatial extreme precipitation from a high-resolution dataset.
Inference for this model traditionally relies on composite likelihoods, not
for computational efficiency, but for using more data for estimation. In
recent years, this modelling framework has gained much traction in the field
of extreme value theory, as it allows for a more flexible description of
extremal dependence than classical extreme value models [20]. The model
consists in choosing one spatial location and modelling surrounding data
conditional on extreme behaviour at that site. As extreme behaviour is, by
definition, rare, inference commonly relies on a composite likelihood that
combines separate models from multiple conditioning sites in order to provide
one global model fit that, on average, represents the data well for any given
conditioning site [18, 52, 40, 46]. A modification of the spatial conditional
extremes model implemented in R-INLA is proposed by [45] in order to perform
truly high-dimensional Bayesian inference using INLA and the SPDE approach.
The authors only employ single site conditioning, stating that their reason
for not providing a global model fit comes from the problems that occur when
performing Bayesian inference with a composite likelihood using R-INLA. We
demonstrate that our adjustment method solves these problems and allows for a
global Bayesian model fit with improved frequency properties. We also improve
upon parts of the model construction proposed by [52] and [45], and we
demonstrate a general methodology for implementing different variants of the
conditional extremes model in R-INLA using the recent cgeneric framework.

The remainder of the paper is organised as follows: In Section 2, we examine
problems that can occur when performing Bayesian inference based on composite
and/or misspecified likelihoods. Previous adjustment methods are presented,
and our novel method for post-processing posteriors is developed. The spatial
conditional extremes model and inferential methods based on composite
likelihoods are reviewed in Section 3. There, we also build upon [45] for
improving how such models can be implemented in R-INLA. Section 4 contains
several simulation studies for evaluating both the adjustment method and our
conditional extremes model implementation in R-INLA. Section 5 contains a case
study where the spatial conditional extremes model and the posterior
adjustment method are used for modelling extreme precipitation from high-
resolution radar data in Norway. Finally, we conclude in Section 6 with some
discussion and perspectives on future research.

##  2 Methodology

###  2.1 Inference with composite and/or misspecified likelihoods

Let $\bm{y}\in\mathbb{R}^{d}$ be a realisation of a random vector $\bm{Y}$
with distribution function $G$, and assume that one models the data $\bm{y}$
using the parametric family
$\mathcal{F}=\\{F(\cdot;\bm{\theta}):\bm{\theta}\in\Theta\\}$. If $G$ is not
contained in $\mathcal{F}$, we say that $\mathcal{F}$ is misspecified. Denote
the likelihood function corresponding to $F(\cdot;\bm{\theta})$ as
$L(\bm{\theta};\cdot)$. [31] defines composite likelihoods as functions
composed of likelihood functions of conditional or marginal events for
$\bm{y}$. Thus, given a set of events
$A_{1},A_{2},\ldots,A_{m},B_{1},B_{2},\ldots,B_{m}\subseteq\mathbb{R}^{d}$ and
nonnegative weights $\omega_{1},\omega_{2},\ldots,\omega_{m}\geq 0$, a
composite likelihood is defined as the weighted product

$L_{c}(\bm{\theta};\bm{y})=\prod_{i=1}^{m}L^{(i)}(\bm{\theta};\bm{y})^{\omega_{i}},\text{
where }L^{(i)}(\bm{\theta};\bm{y})=L(\bm{\theta};\left[\bm{y}\in
A_{i}\mid\bm{y}\in B_{i}\right]),$

Given $n$ independent realisations
$\mathcal{Y}=\\{\bm{y}_{1},\ldots\bm{y}_{n}\\}$ of $\bm{Y}$, the composite
log-likelihood equals

$\ell_{c}(\bm{\theta};\mathcal{Y})=\sum_{j=1}^{n}\ell_{c}(\bm{\theta};\bm{y}_{j})=\sum_{j=1}^{n}\sum_{i=1}^{m}\omega_{ij}\log
L^{(i)}(\bm{\theta};\bm{y}_{j})=\sum_{j=1}^{n}\sum_{i=1}^{m}\omega_{ij}\
\ell^{(i)}(\bm{\theta};\bm{y}_{j}),$ (1)

where we allow for different weights for different realisations of $\bm{Y}$.
Incorrectly interpreting the composite likelihood as a true likelihood is
tantamount to specifying a model for $\mathcal{Y}$ using a family of
distributions in which $[\bm{y}_{i}\in A_{j}\mid\bm{y}_{i}\in B_{j}]$ is
(wrongly) assumed to be independent of $[\bm{y}_{i}\in A_{k}\mid\bm{y}_{i}\in
B_{k}]$ for all $i$ and $j\neq k$. One can therefore consider a composite
likelihood as a type of misspecified model with unrealistic independence
assumptions. From (1) it is clear that any proper likelihood is a composite
likelihood with $m=1$ and $\omega_{1}=1$. Consequently, all the theory of
composite likelihoods trivially holds for valid likelihoods.

Using the theory of unbiased estimating equations one can show that the
maximum composite likelihood estimator,
$\hat{\bm{\theta}}_{c}=\arg\max_{\bm{\theta}}\ell_{c}(\bm{\theta};\mathcal{Y})$,
is asymptotically Gaussian distributed,

$\bm{\mathcal{I}}(\bm{\theta}^{*})^{1/2}\left(\hat{\bm{\theta}}_{c}-\bm{\theta}^{*}\right)\rightsquigarrow\mathcal{N}(\bm{0},\bm{I}),\text{
as }n\rightarrow\infty,$

under some weak regularity conditions [54], where $\bm{I}$ is the identity
matrix and $\bm{\theta}^{*}$ is the minimiser of the Kullback-Leibler distance
[28, KLD;] between the likelihood of the true distribution $G$ and the
composite likelihood $L_{c}$. Furthermore, $\bm{\mathcal{I}}(\bm{\theta})$ is
the so-called Godambe sandwich information matrix [16]:

$\begin{gathered}\bm{\mathcal{I}}(\bm{\theta})=\bm{H}(\bm{\theta})\bm{J}(\bm{\theta})^{-1}\bm{H}(\bm{\theta}),\\\
\bm{H}(\bm{\theta})=-\text{E}\left[\nabla_{\bm{\theta}}^{2}\ell_{c}(\bm{\theta};\mathcal{Y})\right],\quad\bm{J}(\bm{\theta})=\text{Cov}\left(\nabla_{\bm{\theta}}\ell_{c}(\bm{\theta};\mathcal{Y})\right),\end{gathered}$
(2)

where expectations are taken with respect to $G$. If the composite likelihood
is a valid and correctly specified likelihood then
$\bm{H}(\bm{\theta}^{*})=\bm{J}(\bm{\theta}^{*})$ and the Godambe sandwich
information is reduced to the Fisher information
$\bm{\mathcal{I}}(\bm{\theta}^{*})=\bm{H}(\bm{\theta}^{*})=\bm{J}(\bm{\theta}^{*})$.

Given a prior distribution $\pi(\bm{\theta})$ for $\bm{\theta}$, and a
composite likelihood $L_{c}(\bm{\theta};\cdot)$ such that
$\int_{\Theta}L_{c}(\bm{\theta};\mathcal{Y})\pi(\bm{\theta})\text{d}\bm{\theta}$
is finite, the posterior distribution of $[\bm{\theta}\mid\mathcal{Y}]$ is
defined as

$\pi_{c}(\bm{\theta}\mid\mathcal{Y})=\frac{L_{c}(\bm{\theta};\mathcal{Y})\pi(\bm{\theta})}{\int_{\Theta}L_{c}(\bm{\theta};\mathcal{Y})\pi(\bm{\theta})\text{d}\bm{\theta}}.$

Under appropriate regularity conditions, this distribution converges
asymptotically to a Gaussian distribution with mean $\bm{\theta}^{*}$ and
covariance matrix $\bm{H}(\bm{\theta}^{*})^{-1}$ [3, 27, 39]. Notice that the
asymptotic covariance of the posterior is different from the asymptotic
covariance of the maximum likelihood estimator, with equality only when the
composite likelihood is a valid and correctly specified likelihood. This
causes the posterior to lose its frequency properties when
$L_{c}(\bm{\theta};\cdot)$ is misspecified. Consequently, care must always be
taken when interpreting credible intervals that are based on composite
likelihoods or other types of clearly misspecified models.

Adjustments to the composite likelihood have been proposed for improving its
frequency properties. [38] and [39] propose to rescale the weights
$\omega_{ij}$ in (1) so that the composite and the true likelihood ratio
statistics achieve the same distribution asymptotically. [48] also rescale
$\omega_{ij}$, using bootstrap techniques to estimate the best rescaling
value. [6] and [39] propose an affine transformation to the argument of the
composite likelihood,

$L_{\text{adj}}(\bm{\theta};\mathcal{Y})=L_{c}\left(\bm{\theta}^{*}+\bm{C}(\bm{\theta}-\bm{\theta}^{*});\mathcal{Y}\right),$
(3)

where the matrix $\bm{C}$ is a solution of the equation
$\bm{C}^{T}\bm{H}(\bm{\theta}^{*})\bm{C}=\bm{\mathcal{I}}(\bm{\theta}^{*})$,
with $\bm{\mathcal{I}}(\bm{\theta})$ and $\bm{H}(\bm{\theta})$ as defined in
(2). All of the aforementioned methods are based on modifying the likelihood
function before performing inference. However, this option is not available
for users of inference software like R-INLA, where one must choose between a
preselected set of likelihood functions. Therefore, in Section 2.2, we develop
a novel adjustment method for posterior distributions based on composite
and/or misspecified likelihoods that can be used in conjunction with any
Bayesian inference program, including R-INLA.

###  2.2 Postprocessing posterior samples

Inspired by [6] and [39], we propose an affine transformation for samples from
the posterior distribution, aimed at improving the frequency properties of
credible intervals. Assume that the parameter $\bm{\theta}$ is sampled from a
posterior distribution based on a composite and/or misspecified likelihood. We
then define the adjusted version of $\bm{\theta}$ as

$\bm{\theta}_{\text{adj}}=\bm{\theta}^{*}+\bm{C}(\bm{\theta}-\bm{\theta}^{*}),$
(4)

where the matrix $\bm{C}$ is chosen such that the asymptotic distribution of
$\bm{\theta}_{\text{adj}}$ as $n\rightarrow\infty$ is Gaussian with mean
$\bm{\theta}^{*}$ and covariance matrix
$\bm{\mathcal{I}}(\bm{\theta}^{*})^{-1}$. To achieve this, $\bm{C}$ must solve
the equation
$\bm{C}\bm{H}(\bm{\theta}^{*})^{-1}\bm{C}^{T}=\bm{\mathcal{I}}(\bm{\theta}^{*})^{-1}.$
A solution of this equation is

$\bm{C}=\left(\bm{M}_{1}^{-1}\bm{M}_{2}\right)^{T},\text{ with
}\bm{M}_{1}^{T}\bm{M}_{1}=\bm{H}({\bm{\theta}}^{*})^{-1}\text{ and
}\bm{M}_{2}^{T}\bm{M}_{2}=\bm{\mathcal{I}}({\bm{\theta}}^{*})^{-1},$

where the matrix square roots $\bm{M}_{1}$ and $\bm{M}_{2}$ are computed
using, e.g., singular value decomposition. Notice that our $\bm{C}$ matrix is
different from the $\bm{C}$ matrix of [6] and [39] in (3), since we solve a
slightly different equation.

In order to estimate $\bm{C}$ one must first estimate the values of
$\bm{\theta}^{*}$, $\bm{H}(\bm{\theta}^{*})$ and $\bm{J}(\bm{\theta}^{*})$.
The mode of the likelihood, $\bm{\theta}^{*}$, can be estimated by the mode of
the posterior $\pi_{c}(\bm{\theta}\mid\mathcal{Y})$, denoted
$\hat{\bm{\theta}}^{*}$. The estimator $\hat{\bm{\theta}}^{*}$ provides a good
approximation to $\bm{\theta}^{*}$ if the prior is not overly informative and
the sample size is large enough, but other estimators may be more suitable if
this does not hold. The first step towards estimating
$\bm{J}(\bm{\theta}^{*})$ is to compute the gradients
$\nabla_{\bm{\theta}}\ell^{(i)}(\hat{\bm{\theta}}^{*};\bm{y}_{j})$ for
$i=1,2,\ldots m$ and $j=1,2,\ldots,n$. These can be computed analytically or
estimated using numerical derivation methods. Having computed all the
gradients, we estimate $\bm{J}(\bm{\theta}^{*})$ with

$\hat{\bm{J}}(\hat{\bm{\theta}}^{*})=\sum_{i,j}\omega_{ij}\nabla_{\bm{\theta}}\ell^{(i)}(\hat{\bm{\theta}}^{*};\bm{y}_{j})\sum_{(i^{\prime},j^{\prime})\in\Delta(i,j)}\omega_{i^{\prime}j^{\prime}}\left(\nabla_{\bm{\theta}}\ell^{(i^{\prime})}(\hat{\bm{\theta}}^{*};\bm{y}_{j^{\prime}})\right)^{T},$
(5)

where $\Delta(i,j)$ is the set of neighbours to $(i,j)$ in the sense that
$(i^{\prime},j^{\prime})\in\Delta(i,j)$ if and only if
$\nabla_{\bm{\theta}}\ell^{(i^{\prime})}(\cdot;\bm{y}_{j^{\prime}})$ is
correlated with $\nabla_{\bm{\theta}}\ell^{(i)}(\cdot;\bm{y}_{j})$. Summing
over all the non-correlated pairs of tuples would introduce unnecessary noise
that could cause the estimator to approximately equal zero [34]. In practice
one can therefore often compute $\hat{\bm{J}}(\hat{\bm{\theta}}^{*})$ using a
sliding window approach.

Estimation of $\bm{H}(\bm{\theta}^{*})$ is often easier than that of
$\bm{J}(\bm{\theta}^{*})$, since the former is a matrix of expected values,
whereas the latter is a covariance matrix. The law of large number implies
that, for $n$ and $m$ large enough, a good estimator for
$\bm{H}(\bm{\theta}^{*})$ is

$\hat{\bm{H}}(\hat{\bm{\theta}}^{*})=-\sum_{j=1}^{n}\sum_{i=1}^{m}\omega_{ij}\nabla^{2}_{\bm{\theta}}\ell^{(i)}(\hat{\bm{\theta}}^{*};\bm{y}_{j}).$
(6)

Thus, all we need for estimating $\bm{H}(\bm{\theta}^{*})$ is to compute the
Hessian of the composite log-likelihood terms at $\hat{\bm{\theta}}^{*}$. For
users of R-INLA, this is especially simple, as an estimator of
$\bm{H}(\hat{\bm{\theta}}^{*})^{-1}$ is computed and returned by the program
when fitting a model to data.

##  3 The spatial conditional extremes model

Having developed an adjustment method for posterior samples based on composite
and/or misspecified likelihoods, we aim to perform Bayesian inference for the
spatial conditional extremes model using a combination of R-INLA and a
composite likelihood. The model is described in Section 3.1, while Section 3.2
details how to perform inference based on single site conditioning, and how to
further improve it using a composite likelihood.

###  3.1 Model description

Statistical models for spatial extremes are crucial for properly assessing
risks in a variety of applications, see [20] for a non-extensive list of use-
cases. An important component of spatial extreme value theory is the
description and characterisation of a spatial process’ extremal dependence
properties. Two locations that have a positive limiting probability to
experience extremes simultaneously are denoted extremally dependent. If this
limiting probability is zero, the locations are denoted extremally
independent. Clearly, simultaneous extreme precipitation over a large area has
more dire consequences than extreme precipitation at one single point in
space. Therefore, correct estimation of extremal dependence properties is of
utmost importance when estimating environmental risks. Most of the classical
models within spatial extreme value theory are either based on max-stable
processes [10, e.g.,] or latent Gaussian models with conditional independence
assumptions [9, 37, e.g.,]. The latent Gaussian models often capture the
marginal structure of spatial processes well, and are computationally
efficient. However, due to common assumptions of conditional independence,
they fail to properly describe extremal dependence properties. Max-stable
processes are based on strong mathematical theory and allow for rich modelling
of extremal dependence, but in practice they are often too computationally
demanding and yield unsatisfactory models for extremal independence [20]. In
recent years the conditional extremes model [18, 17, 52] has gained much
traction, as it allows for a more flexible description of extremal dependence
by modelling the behaviour of a random vector conditional on one of its
components being extreme.

Let $Y(\bm{s})$ be a random process defined over space
$(\bm{s}\in\mathcal{S}\subset\mathbb{R}^{2})$ with Laplace margins. For this
random process, [52] assume the existence of standardising functions
$a(\bm{s};\bm{s}_{0},y_{0})$ and $b(\bm{s};\bm{s}_{0},y_{0})$ such that, for a
large enough threshold $t$,

$\left[Y(\bm{s})\mid
Y(\bm{s}_{0})=y_{0}>t\right]\overset{d}{=}a(\bm{s};\bm{s}_{0},y_{0})+b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0}),\quad\bm{s},\bm{s}_{0}\in\mathcal{S},$
(7)

where $Z(\bm{s};\bm{s}_{0})$ is a random process satisfying
$Z(\bm{s}_{0};\bm{s}_{0})=0$ almost surely, and
$a(\bm{s};\bm{s}_{0},y_{0})\leq y_{0}$, with equality when
$\bm{s}=\bm{s}_{0}$. It can be shown that $Y(\bm{s})$ and $Y(\bm{s}_{0})$ are
extremally dependent when $a(\bm{s};\bm{s}_{0},y_{0})=y_{0}$ and
$b(\bm{s};\bm{s}_{0},y_{0})=1$, while they are extremally independent when
$a(\bm{s};\bm{s}_{0},y_{0})<y_{0}$, with larger values of $a(\cdot)$ and
$b(\cdot)$ corresponding to larger degrees of subasymptotic dependence [18].
[52] provides some guidance on parametric functions for $a(\cdot)$ and
$b(\cdot)$ together with parametric distributions for $Z(\cdot)$ that cover a
large range of already existing models. To perform inference with high-
dimensional models, they further propose to model $Z(\cdot)$ as a random
process with a Gaussian copula and delta-Laplace marginal distributions. Their
proposed model for $Z(\cdot)$ has later seen usage by [43] and [40].

In order to perform Bayesian inference for the spatial conditional extremes
model with R-INLA, [45] cast (7) into the latent Gaussian modelling framework
by adding a Gaussian nugget effect and requiring $Z(\cdot)$ to be a Gaussian
random field. This gives the model

$\left[Y(\bm{s})\mid
Y(\bm{s}_{0})=y_{0}>t\right]\overset{d}{=}a(\bm{s};\bm{s}_{0},y_{0})+b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0})+\epsilon(\bm{s};\bm{s}_{0}),$
(8)

where $\epsilon(\bm{s};\bm{s}_{0})$ is Gaussian white noise with constant
precision, satisfying $\epsilon(\bm{s}_{0};\bm{s}_{0})=0$ almost surely. To
further speed up inference, they also assume that $Z(\cdot)$ has zero mean and
a Matérn covariance structure, so that it can be approximated by a sparse
Gaussian Markov random field using the SPDE approach of [30]. Based on the
guidelines of [52], they choose $b(\bm{s};\bm{s}_{0},y_{0})=y_{0}^{\beta}$
with $0\leq\beta<1$, whereas for $a(\cdot)$ they propose the semiparametric
spline model

$a(\bm{s};\bm{s}_{0},y_{0})=\alpha(\|\bm{s}-\bm{s}_{0}\|)y_{0}+\gamma(\|\bm{s}-\bm{s}_{0}\|),$

where $\alpha(\cdot)$ and $\gamma(\cdot)$ are spline functions with
$\alpha(0)=1$ and $\gamma(0)=0$.

Here, we adopt the model in (S1.1) and keep the somewhat strict requirement of
modelling $Z(\cdot)$ as a Gaussian random field using the SPDE approximation.
However, we propose some important modifications that allow for faster
inference and more flexible model building:

  * •

Using the recent cgeneric framework, which allows us to implement model
components that are not readily available in R-INLA using the C programming
language, we propose a general method for implementing the model (S1.1) in
R-INLA such that $a(\cdot)$ and $b(\cdot)$ can take any kind of form, both
parametric and semiparametric.

  * •

We explain why it is problematic to constrain $Z(\bm{s};\bm{s}_{0})$ as in
[52] and [45], and demonstrate an efficient constraining method that gives
better properties for $Z(\bm{s};\bm{s}_{0})$.

Both of these modifications are described in Section S1 in the supplementary
materials.

###  3.2 Inference

Inference for the conditional extremes model can be performed by choosing one
conditioning site $\bm{s}_{0}\in\mathcal{S}$ and a large threshold $t$. The
parameters of the standardising functions $a(\cdot)$ and $b(\cdot)$, the
residual field $Z(\cdot)$ and the nugget effect $\epsilon(\cdot)$ can then be
estimated using all observations where $Y(\bm{s}_{0})>t$. Since this approach
only allows us to use observations from times when $Y(\bm{s}_{0})$ is extreme,
it yields separate model fits for each conditioning site in $\mathcal{S}$.
Under stationarity, it is reasonable to believe that inference could be
strengthened by combining information from multiple conditioning sites into
one global model fit. If we assume a stationary model, in the sense that the
parameters of $a(\cdot)$, $b(\cdot)$, $Z(\cdot)$ and $\epsilon(\cdot)$ are
independent of the choice of $\bm{s}_{0}$, a global model fit can be achieved
by performing inference with the composite likelihood of [18] and [52]. Given
observations
$\mathcal{Y}=\\{y_{j}(\bm{s}_{i}):j=1,2,\ldots,n,i=1,2,\ldots,d\\}$ of the
random process at $n$ time points and $d$ locations, the composite log-
likelihood may be expressed as

$\ell_{c}(\bm{\theta};\mathcal{Y})=\sum_{i=1}^{d}\sum_{j=1}^{n}\ell(\bm{\theta};\bm{y}_{j}\mid
y_{j}(\bm{s}_{i}))I(y_{j}(\bm{s}_{i})>t).$ (9)

Here, $\ell(\cdot)$ is the log-likelihood of the conditional extremes model in
(S1.1), $\bm{y}_{j}$ is a $d$-dimensional vector containing all observations
from time point $j$, $I(\cdot)$ is an indicator function and $\bm{\theta}$
contains all the parameters of $a(\cdot)$, $b(\cdot)$, $Z(\cdot)$ and
$\epsilon(\cdot)$.

As described in Section 2, the downside of performing Bayesian inference with
a composite likelihood is that the resulting posterior distribution might have
poor frequency properties and become too focused or too wide. To the best of
our knowledge, our paper is the first attempt to perform Bayesian inference
for the spatial conditional extremes model based on the composite likelihood
in (9).

##  4 Simulation studies

We evaluate the performance of our adjustment method through several
simulation studies, considering different composite and/or misspecified
likelihoods described in Section 1. In Section 4.1, a Gaussian distribution is
fitted to univariate data from a Student’s $t$-distribution. In Section 4.2,
an SPDE approximation of low rank is applied for modelling data from a spatial
Gaussian random field. In Section 4.3, inference on a spatial Gaussian random
field is performed using a block composite likelihood. Finally, in Section
4.4, the composite likelihood approach for the spatial conditional extremes
model is considered. In each study we adjust the posteriors and compare
frequency properties between the unadjusted and the adjusted credible
intervals.

###  4.1 Univariate misspecified model

We here simulate data from a Student’s $t$-distribution with one degree of
freedom and estimate their median $\theta=0$ using a Gaussian distribution
with a fixed variance of 1. Since both the Student’s $t$-distribution and the
Gaussian distribution are symmetric, the KLD minimiser between the true
Student’s $t$-distribution and the misspecified Gaussian distribution is
$\theta^{*}=\theta=0$. However, due to the heavy tails of the Student’s
$t$-distribution and the small, fixed variance of the Gaussian model, the
posterior for $\theta$ will often be too focused. Thus, credible intervals
will be too thin, meaning that an estimated $(1-\alpha)$-credible interval
will contain $\theta^{*}$ considerably less often than $(1-\alpha)\times
100\%$ of the times, which should be the correct asymptotic coverage
probability when the effect of the prior is negligible. A correct adjustment
will therefore increase the variance of the posterior distribution and create
wider credible intervals.

To check this, we sample 100 observations from the Student’s $t$-distribution
and choose a Gaussian prior for $\theta$ with mean 0 and precision $0.001$.
Adjusted and unadjusted credible intervals are then created using R-INLA and
the adjustment method from Section 2.2. This is repeated 1000 times. Coverage
probabilities are then estimated by counting how many of the 1000 credible
intervals include $\theta^{*}$. Results are displayed in Table 1. As expected,
the unadjusted credible intervals are considerably too thin, with, e.g., the
$95\%$ intervals containing $\theta^{*}$ in only $13\%$ of the experiments,
whereas the adjusted coverage percentage is $98\%$. Similar improvements are
found for the $90\%$ and $99\%$ credible intervals. In practice, one would of
course never fix the variance to 1 when fitting a Gaussian distribution in
this context. However, this toy example serves as a simple demonstration for
how misspecified models can result in credible intervals with poor frequency
properties, and how these can be greatly improved by simply adjusting the
posterior, using our proposed method.

Table 1: Coverage percentage for the unadjusted and adjusted credible intervals for three different credible interval sizes. Aim | Unadjusted | Adjusted
---|---|---
$90\%$ | $11\%$ | $94\%$
$95\%$ | $13\%$ | $98\%$
$99\%$ | $15\%$ | $100\%$

###  4.2 Low-rank SPDE approximation

We now simulate data from a Gaussian random field $Z(\bm{s})$ with Matérn
covariance function

$\text{Cov}(Z(\bm{s}),Z(\bm{s}^{\prime}))=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}(\kappa\|\bm{s}-\bm{s}^{\prime}\|)^{\nu}K_{\nu}(\kappa\|\bm{s}-\bm{s}^{\prime}\|),$

with variance $\sigma^{2}$, range $\rho=\sqrt{8\nu}/\kappa>0$ for some
$\kappa>0$, and known smoothness parameter $\nu>0$. Here, $K_{\nu}$ is the
modified Bessel function of the second kind and order $\nu$. We also add a
Gaussian nugget effect with precision $\tau$. Parameter estimation is then
performed using an SPDE approximation of low rank, i.e., based on a coarse
mesh used to discretise the spatial domain. Such low-rank approximations are
typically unable to capture all the variability in the data, which means that
the nugget effect has to explain a large percentage of the variance, leading
to underestimation of the precision $\tau$. Thus, we expect the KLD minimiser
$\bm{\theta}^{*}$ to be different from the true parameters
$\bm{\theta}=(\tau,\rho,\sigma)^{T}$.

We sample 200 replicates of the spatial field on the domain
$\mathcal{S}=[0,25]\times[0,25]$ and observe it at 400 random locations. True
parameters are set to $\bm{\theta}=(\tau,\rho,\sigma)^{T}=(100,12,1)^{T}$ with
smoothness $\nu=1.5$. We assign $\tau$ a gamma prior with shape 1 and scale
$2\times 10^{4}$, while $\rho$ and $\sigma$ are given a joint penalised
complexity (PC) prior [44, 13], setting P$(\rho<12)=0.5$ and
P$(\sigma>1)=0.5$. We compute both unadjusted and adjusted credible intervals
for the three parameters. This is repeated 300 times, using the same
locations, but new observations each time.

To estimate the unknown KLD minimiser $\bm{\theta}^{*}$, we simulate $10^{4}$
replications of the Gaussian Matérn field and compute the maximum likelihood
estimator for the misspecified SPDE model. This gives
$\bm{\theta}^{*}=(\tau^{*},\rho^{*},\sigma^{*})\approx(13.0,14.5,1.2)^{T}$. As
expected, $\tau$ is severely underestimated, while $\rho$ and $\sigma$ are
slightly overestimated.

Table 2: Coverage percentages for unadjusted and adjusted credible intervals using the SPDE approach with a too coarse mesh. Aim | $\tau$ | $\tau_{\text{adj}}$ | $\rho$ | $\rho_{\text{adj}}$ | $\sigma$ | $\sigma_{\text{adj}}$
---|---|---|---|---|---|---
$90\%$ | $48\%$ | $93\%$ | $91\%$ | $90\%$ | $90\%$ | $90\%$
$95\%$ | $55\%$ | $97\%$ | $95\%$ | $95\%$ | $95\%$ | $96\%$
$99\%$ | $69\%$ | $99\%$ | $99\%$ | $98\%$ | $100\%$ | $99\%$

Table 2 displays estimated coverage probabilities detailing how often the
parameters of $\bm{\theta}^{*}$ are included in their respective credible
intervals. The adjustment of the posterior yields a considerable improvement
for $\tau$. The unadjusted frequency properties of $\rho$ and $\sigma$,
however, are already good, and our adjustment method does not deteriorate the
credible intervals for these parameters.

###  4.3 Block composite spatial Gaussian likelihood

In this other experiment, we test how the adjustment method works with a block
composite likelihood [12]. To do so we sample 100 realisations from a Gaussian
Matérn field with parameters and domain as in Section 4.2. This time we
generate samples directly using the SPDE approximation, so no model
misspecification is introduced when using it for inference. As in Section 4.2
we observe the fields at 400 random locations.

Inference is performed using a block composite likelihood, whereby the spatial
domain is divided into 25 disjoint blocks of size $5\times 5$, all wrongly
assumed to be independent of each other. This is in some sense a milder
misspecification than in Section 4.2, as it gives
$\bm{\theta}^{*}=\bm{\theta}$. On the other hand, such independence
assumptions typically lead to overconfidence in the posterior, meaning that
credible intervals will become too thin. A correct adjustment is therefore
expected to increase the variance of the posterior. Unadjusted and adjusted
credible intervals are created using R-INLA with the SPDE approach, where the
priors for $\bm{\theta}=(\rho,\sigma^{2},\tau)$ are the same as in Section
4.2. The experiment is repeated 300 times.

Table 3 displays estimated coverage probabilities for the model parameters. As
expected, the unadjusted credible intervals are often too thin, and the
adjustment method successfully increases the posterior variance to provide
improved coverage probabilities.

Table 3: Coverage percentages for unadjusted and adjusted credible intervals using the block composite likelihood. Aim | $\tau$ | $\tau_{\text{adj}}$ | $\rho$ | $\rho_{\text{adj}}$ | $\sigma$ | $\sigma_{\text{adj}}$
---|---|---|---|---|---|---
$90\%$ | $88\%$ | $89\%$ | $61\%$ | $90\%$ | $61\%$ | $91\%$
$95\%$ | $96\%$ | $94\%$ | $71\%$ | $93\%$ | $67\%$ | $96\%$
$99\%$ | $99\%$ | $100\%$ | $83\%$ | $99\%$ | $77\%$ | $99\%$

###  4.4 Spatial conditional extremes model

We evaluate our adjustment method when applied to the spatial conditional
extremes model (S1.1) with a global model fit based on the composite
likelihood in (9). In the simulation experiment we replicate the setup of the
case study in Section 5. The standardising functions $a(\cdot)$ and $b(\cdot)$
are defined as

$a(\bm{s};\bm{s}_{0},y_{0})=y_{0}\exp\left\\{-\left(d/\lambda\right)^{\kappa}\right\\},\quad
b(\bm{s};\bm{s}_{0},y_{0})=\sigma\sqrt{1-\exp\left(-2d/\rho_{b}\right)},$ (10)

with $d=\|\bm{s}-\bm{s}_{0}\|$ and parameters
$(\lambda,\kappa,\sigma,\rho_{b})$, all strictly positive. Furthermore, the
threshold $t$ is set equal to the $99.75\%$ quantile of the Laplace
distribution. The residual field, $Z(\cdot)$, is an SPDE approximation to a
Gaussian Matérn field with a fixed variance of $1$ and range parameter $\rho$,
created using the mesh seen in the left subplot of Figure 1. Note that, since
$b(\bm{s}_{0};\bm{s}_{0},y_{0})=0$, we do not need to enforce that
$Z(\bm{s}_{0};\bm{s}_{0})=0$. The nugget effect $\epsilon(\cdot)$ is defined
as Gaussian white noise with zero mean and precision $\tau$. In order to
produce realisations that mimic the real data application, parameter values
are set to
$\bm{\theta}=(\lambda,\kappa,\sigma,\rho_{b},\rho,\tau)^{T}\approx(19.1,0.6,1.9,4.6,13.0,23.1)^{T},$
where $\lambda$, $\rho_{b}$ and $\rho$ have unit km. Priors for $\bm{\theta}$
are displayed in Table 4. See Section 5 for more details and explanations.

Figure 1:  Left plot: Given a conditioning site $\bm{s}_{0}$ (displayed with ($\blacktriangle$)), locations used for inference are displayed as big black dots $(\bullet)$ and locations in $\mathcal{S}$ that are not used for inference are displayed as small black dots $(\cdot)$. The SPDE mesh is displayed in grey. Right plot: Elevation (m) map, over the Fosen area in central Norway. The study area $\mathcal{S}$ is located inside the black rectangle, and the Rissa radar is displayed using a triangle ($\blacktriangle$)).  Table 4: Prior distributions of the conditional extremes parameters. $\mathcal{N}(\mu,\sigma^{2})$ denotes the Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$. $\Gamma(a,b)$ denotes the gamma distribution with shape parameter $a$ and scale parameter $b$. The parameters $\rho$ and $\sigma$ are given the joint PC prior of [13] such that P$(\rho<60)=0.95$ and P$(\sigma>4)=0.05$. $\log(\lambda)\sim\mathcal{N}(4,3^{2})$, |  $\log(\rho_{b})\sim\mathcal{N}(\log(6),2^{2})$, |  $\rho\sim\text{PC}(60,0.95)$,
---|---|---
$\log(\kappa)\sim\mathcal{N}(-0.4,3^{2})$, |  $\tau\sim\Gamma(1,2\cdot 10^{4})$, |  $\sigma\sim\text{PC}(4,0.05)$,

To ensure that our model implementation works as it should we first choose one
location $\bm{s}_{0}\in\mathcal{S}$ and simulate from the conditional extremes
model using $\bm{s}_{0}$ as conditioning site. Parameters are then recovered
from the data using our implemented model in R-INLA. We find that the
implemented methodology accurately estimates the parameters. However, we also
find that R-INLA struggles considerably when initial values are too far away
from $\bm{\theta}$. Consequently, we choose to divide inference into two
steps: First we compute the maximum likelihood estimator using the R function
optim(). Then we use the maximum likelihood estimator as initial values for
R-INLA. This requires approximately the same amount of computational
resources, but yields considerably better results, than simply using R-INLA
without first computing the maximum likelihood estimator.

For evaluating the adjustment method, we sample from the distribution of

$\left[\left\\{Y(\bm{s}):\bm{s}\in\mathcal{S}\right\\}\mid\max_{\bm{s}\in\mathcal{S}}Y(\bm{s})>t\right],$
(11)

defined such that the spatial conditional extremes model holds for each
conditioning site $\bm{s}\in\mathcal{S}$. Different simulation algorithms for
sampling from this distribution have been developed by [26] and [52]. We use
the method of [52] as we find it to be more efficient for high-dimensional
spatial models where $a(\cdot)$ does not decay quickly to zero and $b(\cdot)$
is of the same magnitude as the threshold $t$. A considerable problem with the
distribution of (11), described by [18] and [32], is that it lacks self-
consistency, i.e., it is not consistent with the individual models based on
single conditioning sites (see Appendix A for a more detailed explanation). In
short, the consequence of this self-inconsistency is that
$\bm{\theta}^{*}\neq\bm{\theta}$ even when we fit the distribution of (11) to
data sampled from the distribution of (11) with parameters $\bm{\theta}$.
Thus, we must estimate $\bm{\theta}^{*}$ before we can evaluate frequency
properties. To estimate $\bm{\theta}^{*}$, we sample $5\times 10^{4}$
replications of (11) and compute the maximum composite likelihood estimator
using (9). This gives
$\hat{\bm{\theta}}^{*}\approx(15.5,0.8,1.9,1.5,13.1,20.1)^{T}\neq\bm{\theta}.$

Table 5: Coverage percentages for unadjusted and adjusted credible intervals based on the composite spatial conditional extremes likelihood. Aim | $\lambda$ | $\lambda_{\text{adj}}$ | $\kappa$ | $\kappa_{\text{adj}}$ | $\sigma$ | $\sigma_{\text{adj}}$ | $\rho_{b}$ | ${\rho_{b}}_{\text{adj}}$ | $\rho$ | $\rho_{\text{adj}}$ | $\tau$ | $\tau_{\text{adj}}$
---|---|---|---|---|---|---|---|---|---|---|---|---
$90\%$ | $27\%$ | $69\%$ | $54\%$ | $72\%$ | $32\%$ | $71\%$ | $78\%$ | $83\%$ | $36\%$ | $71\%$ | $66\%$ | $84\%$
$95\%$ | $33\%$ | $75\%$ | $63\%$ | $80\%$ | $40\%$ | $79\%$ | $85\%$ | $87\%$ | $40\%$ | $76\%$ | $76\%$ | $90\%$
$99\%$ | $43\%$ | $84\%$ | $78\%$ | $88\%$ | $49\%$ | $85\%$ | $92\%$ | $93\%$ | $54\%$ | $87\%$ | $87\%$ | $96\%$

Finally, we evaluate the adjustment method in a similar fashion as the
previous simulation studies. We sample 10 replications from the distribution
of (11) and fit the composite conditional extremes model to all observed
threshold exceedances. For the chosen design, 10 replications from the
distribution of (11) correspond to between 100 and 200 threshold exceedances
over all the chosen conditioning sites. The simulation experiment is repeated
300 times. Coverage percentages for the unadjusted and adjusted posteriors are
displayed in Table 5. Once again, the adjusted posterior has considerably
better frequency properties than the unadjusted posterior, though, the
adjusted posterior properties are not as good as in the previous simulation
studies. As demonstrated in Section S2 in the supplementary material, this is
most likely caused by the complications of simulating from a self-inconsistent
distribution. There, we perform a similar experiment where the global spatial
conditional extremes model is fitted to extremes from a spatial Gaussian
random process, and the adjusted posteriors once again perform as well as in
the three previous simulation studies.

##  5 Case study: Extreme precipitation in Norway

We apply our proposed methodology to the modelling of extreme hourly
precipitation in Norway. Data are presented in Section 5.1 and the inference
is described in Section 5.2. Results are presented and evaluated in Section
5.3.

###  5.1 Data

We consider $1\times 1$ km2 maps of mean hourly precipitation, produced by the
Norwegian Meteorological Institute by processing raw reflectivity data from
the weather radar located in Rissa ($63^{\circ}41^{\prime}26^{\prime\prime}$N,
$10^{\circ}12^{\prime}14^{\prime\prime}$E) in central Norway. Such maps are
available online (https://thredds.met.no), dating back to 1 January 2010. We
extract data from a rectangular domain, close to the Rissa radar, of size
$31\times 31$ km2. Denote the set of all grid points in the rectangular domain
as $\mathcal{S}$. We then have $|\mathcal{S}|=961$ unique locations containing
hourly precipitation estimates. A map containing $\mathcal{S}$ and the Rissa
radar is displayed in the right subplot of Figure 1. For each
$\bm{s}\in\mathcal{S}$, we extract all hourly observations from the summer
months (June, July and August) for the years 2010–2021. Removal of missing
data gives a total of 25512 observations at each location.

###  5.2 Modelling and inference

To fit the conditional extremes model to precipitation data, we must first
standardise the data to have Laplace marginals and choose parametric forms for
the functions $a(\cdot)$ and $b(\cdot)$. This is described in Section S3 in
the supplementary material, where we choose the parametric forms displayed in
(10). Furthermore, the threshold $t$ is set equal to the $99.75\%$ quantile of
the Laplace distribution, and the residual field $Z(\cdot)$ is assumed to have
a fixed variance of $1$ and range parameter $\rho$. Bayesian inference with
the spatial conditional extremes model is performed using the composite
likelihood in (9). Using all locations in $\mathcal{S}$ as conditioning sites
is computationally demanding, so we define a regular sub-grid of resolution
$4\times 4$ km2 and use all locations on the sub-grid as conditioning sites.
This yields a total of 49 conditioning sites, each containing between 15 to 42
available threshold exceedances. The posterior adjustment method is then
applied to the output of R-INLA for improving its frequency properties. In
order to compare the global model to models based on only one conditioning
site, we also create ten single site model fits, based on ten different
conditioning sites. The first conditioning site is chosen as the location on
the $4\times 4$ km2 resolution sub-grid that contains the most threshold
exceedances, while the remaining nine locations are randomly selected in
$\mathcal{S}$. The single site model fits are all based on valid likelihoods,
but they might still be misspecified, e.g., since the choice of modelling the
residual field as a Gaussian SPDE approximation was made solely for
computational reasons, meaning that it does not necessarily provide a perfect
fit to the actual precipitation data. Therefore, the posterior adjustment
method is also applied to all ten single site fits.

The preliminary study in Section S3 of the supplementary material shows that
the largest changes in $a(\cdot)$ and $b(\cdot)$ occur close to $\bm{s}_{0}$.
However, most of the locations in $\mathcal{S}$ are far away from
$\bm{s}_{0}$, where there is less information to be gained. To account for
this and give more weight to close-by locations, we discard some of the
observations far away from $\bm{s}_{0}$ during inference, which also speeds up
inference. The left subplot of Figure 1 shows an example of the locations used
to perform inference for 1 of the 49 chosen conditioning sites of the global
model fit. Similar patterns are used for all other conditioning sites, and
also for the single site model fits.

Using the SPDE approach to model $Z(\cdot)$ requires that we define a
triangular mesh. The mesh needs to be very fine in order to capture the rapid
changes in $b(\cdot)$ close to $\bm{s}_{0}$ (see Section S1 of the
supplementary material for a detailed discussion). Consequently, our mesh has
a node at each $\bm{s}\in\mathcal{S}$, allowing us to accurately capture the
changes in $b(\cdot)$ at every location on the grid. The mesh used for all
conditioning sites is displayed in Figure 1.

Priors for all model parameters ($\lambda,\kappa,\sigma$ and $\rho_{b}$ from
(10), the range $\rho$ in $Z(\bm{s};\bm{s}_{0})$ and $\tau$ in
$\epsilon(\bm{s};\bm{s}_{0})$) are described in Table 4. All priors are
relatively informative to avoid unrealistic behaviour, but have large enough
variance to ensure that they are not overly informative.

The posterior adjustment method is applied to the global and all the single
site model fits. As described in Section 2.2, estimates for $\bm{\theta}^{*}$
and $\bm{H}(\bm{\theta}^{*})$ are provided directly from R-INLA, while
$\bm{J}(\bm{\theta}^{*})$ is estimated by computing gradients of the log-
likelihood, using numerical derivation, and then computing the expression in
(5), using a sliding window approach with a window radius of $5$ hours. For
each posterior distribution, $\bm{\theta}^{*}$, $\bm{H}(\bm{\theta}^{*})$ and
$\bm{J}(\bm{\theta}^{*})$ are estimated using only the data that were
available during inference, meaning that different estimates are computed for
each model fit, and the estimates for the global model are based on more data
than the estimates for the single site models.

###  5.3 Results

Figure 2: Posterior distributions for a subset of the fitted models. The
displayed single site posteriors correspond to single site models nr. 2, 3, 7,
8 and 9.

We fit the global and the ten single site models to data. Posteriors are
displayed in Figure 2. For better visibility, not all posterior distributions
are represented, but those displayed are representative of the patterns
observed. The posterior from the unadjusted global model is considerably more
focused than all other posterior distributions for most parameters, whereas
the adjusted global posterior has a more reasonable width. Interestingly, the
opposite holds for $\rho_{b}$, where the unadjusted global posterior is wider
than all other posteriors, and the adjusted global posterior is the most
focused. The adjusted single site posteriors are somewhat similar to the
unadjusted single site posteriors for $\kappa$ and $\rho_{b}$, and
considerably wider than the unadjusted posteriors for all other parameters.
There is also much variability between the single site model fits, and the
global posteriors are mainly located somewhere in the middle of all the single
site posteriors.

The estimated values of $\rho_{b}$ are very small, yielding an estimated
function $b(\cdot)$ that is approximately constant for all
$\bm{s}\neq\bm{s}_{0}$. We also examine the posteriors of
$a(\bm{s};\bm{s}_{0},\cdot)$ and $b(\bm{s};\bm{s}_{0},\cdot)$ for different
distances $d=\|\bm{s}-\bm{s}_{0}\|$, and find similar patterns as those in
Figure 2: there is considerable variability between the single site
posteriors, and the global posteriors are located approximately in the middle
of the single site posteriors (results not shown).

We then perform model evaluation using the log-score [15, e.g.,]. For each of
the model fits, marginal likelihoods are estimated using Monte Carlo
estimation, i.e., given $n_{s}$ samples
$\bm{\theta}_{1},\ldots\bm{\theta}_{n_{s}}$ from the posterior of model
$\mathcal{M}$, the marginal likelihood is estimated as

$\hat{L}(\mathcal{M};\mathcal{Y})=\frac{1}{n_{s}}\sum_{i=1}^{n_{s}}L_{\mathcal{M}}(\bm{\theta}_{i};\mathcal{Y}),$

where $L_{\mathcal{M}}$ is the likelihood of model $\mathcal{M}$. The
estimated log-score $\hat{\ell}(\mathcal{M};\mathcal{Y})$ is then the
logarithm of $\hat{L}(\mathcal{M};\mathcal{Y})$. For two competing models,
$\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, we conclude that $\mathcal{M}_{1}$
provides a better fit to data $\mathcal{Y}$ if
$\hat{\ell}(\mathcal{M}_{1};\mathcal{Y})>\hat{\ell}(\mathcal{M}_{2};\mathcal{Y})$.

In order to evaluate out-of-sample performance, log-scores are estimated using
threshold exceedances from all the 912 conditioning sites in $\mathcal{S}$
that were not used for fitting the global model. To the best of our knowledge,
the theory behind the log-score is based on valid likelihoods, but the
conditional extremes composite likelihood is not a valid likelihood function.
Thus, we choose not to directly compare global log-scores based on the entire
test set. Instead, for each of the 912 conditioning sites, we estimate the
marginal likelihood of each model using only threshold exceedances from that
single site. Log-scores are then separately compared and ranked for each
conditioning site. We estimate the log-scores using $n_{s}=500$ samples.

Figure 3 displays histograms of the 912 rankings of all the model fits. The
posterior adjustment method improves the log-score of every single model, not
just for the global model based on a composite likelihood. This once again
demonstrates the effectiveness of the posterior adjustment method. Through its
usage we have both learned that there is considerable misspecification present
in all original model fits, and we have been able to improve them by
accounting for this misspecification.

Figure 3 also shows that the global model fit is outperformed by approximately
half of the single site model fits, even though the global model fit is based
on approximately $50$ times more data. This suggests that the distribution of
the precipitation data is not stationary in the sense that $a(\cdot)$,
$b(\cdot)$, $Z(\cdot)$ and $\epsilon(\cdot)$ only depend on the distance
$\|\bm{s}-\bm{s}_{0}\|$ between $\bm{s}$ and the conditioning site
$\bm{s}_{0}$. The value of the KLD minimiser $\bm{\theta}^{*}$ for the global
model is therefore the value that on average provides a best fit to data from
the 49 conditioning sites used for inference. However, this average is not
guaranteed to be close to the best fit for any one of the conditioning sites.
The adjustment method generally widens the global posterior to improve its
frequency properties, but this is different from changing the meaning of the
posterior from “what it is” to “what we want it to be” [53]. Note that the
global model might, however, still be useful for the given data set. If, e.g.,
the threshold $t$ was chosen so large that less than two exceedances were
available at each conditioning site, inference would be impossible for the
single site models, while the global model would still provide useful results.

Figure 3: Histograms displaying the log-score rankings of all 22 model fits,
at 912 different conditioning sites, where 1 is best and 22 is worst.
Unadjusted model fits are displayed in grey, while adjusted ones are displayed
in black.

##  6 Conclusion

A novel adjustment method is proposed for improving the frequency properties
of posterior distributions that are based on composite and/or misspecified
likelihoods. The method performs the required adjustment through
postprocessing, which makes it able to work in conjunction with a large
variety of efficient inference programs for Bayesian inference, such as
R-INLA. Several simulation studies are performed, and they all show
considerable improvements in model fit after applying the adjustment method.
The adjustment method is further applied to high-dimensional Bayesian
modelling of precipitation extremes with the spatial conditional extremes
model, using both single site likelihoods and a global composite likelihood.
An improved implementation of the spatial conditional extremes model in R-INLA
is also developed. The results once again show that the adjustment method
performs satisfactorily, with all adjusted model fits outperforming their
unadjusted counterparts. The global model fit is outperformed by several
single site model fits. This indicates that the distribution of extreme
precipitation might be nonstationary over the study domain.

Our simulation studies and the real data application show that the adjustment
method consistently improves model fits for Bayesian models based on composite
and/or misspecified likelihoods, without deteriorating the performance of
models that already have good frequency properties. Since all models for
natural phenomena are misspecified to a certain degree, the adjustment method
may allow for a simple way of improving model fits in a large variety of
settings, including spatial, temporal and spatio-temporal applications.

##  Appendix A Self-inconsistency of the conditional extremes model

As discussed by [4], if one defines a joint distribution by specifying a set
of conditional distributions, the resulting joint distribution may become
self-inconsistent. This means that joint probabilities may not be well
defined, because different ways of integrating over the conditional
distributions can yield different answers. Since the conditional extremes
distribution is defined through a set of conditional distributions, its joint
distribution may therefore become self-inconsistent. This self-inconsistency
is examined by [18] and [32], who show that it is impossible to achieve self-
consistency for finite thresholds $t$, given extremal independence, i.e.,
$a(\bm{s};\bm{s}_{0},y_{0})<y_{0}$, when the distribution of the residual
field $Z(\cdot)$ is smooth. The consequence of this self-inconsistency is that
probabilities for events where more than one random variable is extreme are
ill-defined. As an example, assume that the bivariate vector $(X,Y)^{T}$ can
be described by the conditional extremes model with threshold $t$. The
probability P$(X>t,Y>t)$ that both components exceed the threshold can, e.g.,
be computed as

$\int_{t}^{\infty}\int_{t}^{\infty}\pi(x\mid
y)\pi(y)\text{d}y\text{d}x,\quad\text{or}\quad\int_{t}^{\infty}\int_{t}^{\infty}\pi(y\mid
x)\pi(x)\text{d}x\text{d}y,$

where $\pi(\cdot)$ denotes different kinds of probability density functions.
Due to the lack of self-consistency between conditional models for $\pi(x\mid
y)$ and $\pi(y\mid x)$, these two integration paths might provide different
answers.

The lack of self-consistency causes problems when one attempts to simulate
from the global conditional extremes distribution (11), which in the bivariate
case is the joint distribution of $[(X,Y)\mid\max(X,Y)>t]$. In order to sample
from this distribution, one must first choose one specific integration path
and define it as the “correct one”. Then, one must sample in a way that
corresponds to this integration path. Consequently, two different simulation
algorithms for the global conditional extremes distribution may produce
samples with considerably different properties. As an example, for a
$d$-dimensional random vector $\bm{X}$, the simulation algorithm of [26] is
based on estimating probabilities as

$\text{P}(\bm{X}\in\cdot\mid\max(\bm{X})>t)=\sum_{i=1}^{d}\pi_{i}^{(1)}\frac{\text{P}(\bm{X}\in\cdot,X_{i}=\max(\bm{X})\mid
X_{i}>t)}{\text{P}(X_{i}=\max(\bm{X})\mid X_{i}>t)},$ (A.1)

where $\pi_{i}^{(1)}=\text{P}(X_{i}=\max(\bm{X})\mid\max(\bm{X})>t)$, while
the [52] algorithm is based on estimating probabilities as

$\text{P}(\bm{X}\in\cdot\mid\max(\bm{X})>t)=\frac{\sum_{K}\sum_{i=1}^{d}\pi_{i}^{(2)}\text{P}(\bm{X}\in\cdot\cap
R_{K}\mid X_{i}>t)/|K|}{\sum_{K}\sum_{i=1}^{d}\pi_{i}^{(2)}\text{P}(\bm{X}\in
R_{K}\mid X_{i}>t)/|K|},$ (A.2)

where $\pi_{i}^{(2)}=\text{P}(X_{i}>t)/\sum_{j=1}^{d}\text{P}(X_{j}>t)$, $K$
is an element in the power set of $\\{1,2,\ldots,d\\}$, and

$R_{K}=\\{\bm{x}=(x_{1},\ldots,x_{d})^{T}\in\mathbb{R}^{d}:x_{i}>t\ \forall\
i\in K,x_{j}\leq t\ \forall\ j\not\in K\\}.$

These two algorithms clearly correspond to two different integration paths,
and may therefore provide different results, when combined with different
conditional extremes models.

We compute P$(X>t,Y>t\mid\max(X,Y)>t)$ using the two integration paths (A.1)
and (A.2), for the symmetrical conditional extremes model where $X$ and $Y$
have exponential marginals and the conditional distribution of $[X\mid Y>t]$
is Gaussian with mean $\mu+\alpha Y$ and variance $\sigma^{2}Y^{2\beta}$, and
vice versa for $[Y\mid X>t]$. If $\mu=0$, $\sigma=1$, $\alpha=0.9$,
$\beta=0.8$ and $t=4$, then equation (A.1) gives a probability of $0.17$,
whereas equation (A.2) gives a probability of $0.37$. We also estimate the
probability using Monte Carlo simulation and find, as expected, that the [26]
algorithm agrees with (A.1), while the [52] algorithm agrees with (A.2). This
massive difference in computed probabilities could potentially have severe
consequences for, e.g., infrastructure built to withstand a certain return
level of hourly precipitation. As none of the possible integration paths are
neither more nor less correct than the others, it is our opinion that one
should be extremely careful when estimating return levels by sampling from the
conditional extremes model where the conditioning variable is unknown, i.e.,
$[\bm{X}\mid\max(X)>t]$. Instead, we believe that one should prefer to
estimate return levels by first choosing a specific conditioning
site/variable, e.g., $X_{i}$, and then sampling from the conditional extremes
distribution of $[\bm{X}\mid X_{i}>t]$, which is both proper and self-
consistent.

## Acknowledgements

The authors are grateful to Jordan Richards and Håvard Rue for many helpful
discussions.

Funding Raphaël Huser was partially supported by the King Abdullah University
of Science and Technology (KAUST) Office of Sponsored Research (OSR) under
Award No. OSR- CRG2020-4394.

Conflict of interest The authors report there are no competing interests to
declare.

Code and data availability The necessary code and data for achieving these
results are available online at
https://github.com/siliusmv/posteriorAdjustment.

## References

  * [1] Sudipto Banerjee, Bradley P. Carlin and Alan E. Gelfand  “Hierarchical Modeling and Analysis for Spatial Data”, Monographs on statistics and applied probability  Chapman & Hall/CRC, 2014  DOI: 10.1201/b17115
  * [2] Sudipto Banerjee, Alan E. Gelfand, Andrew O. Finley and Huiyan Sang  “Gaussian Predictive Process Models for Large Spatial Data Sets”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 70.4, 2008, pp. 825–848  DOI: 10.1111/j.1467-9868.2008.00663.x
  * [3] Robert H. Berk  “Limiting Behavior of Posterior Distributions When the Model Is Incorrect”  In _The Annals of Mathematical Statistics_ 37.1 Institute of Mathematical Statistics, 1966, pp. 51–58  DOI: 10.1214/aoms/1177699597
  * [4] Julian Besag  “Spatial Interaction and the Statistical Analysis of Lattice Systems”  In _Journal of the Royal Statistical Society: Series B (Methodological)_ 36.2, 1974, pp. 192–225  DOI: 10.1111/j.2517-6161.1974.tb00999.x
  * [5] Bob Carpenter et al.  “Stan: A Probabilistic Programming Language”  In _Journal of Statistical Software_ 76.1, 2017  DOI: 10.18637/jss.v076.i01
  * [6] Richard E. Chandler and Steven Bate  “Inference for Clustered Data using the Independence Loglikelihood”  In _Biometrika_ 94.1, 2007, pp. 167–183  DOI: 10.1093/biomet/asm015
  * [7] D.. Cox and N. Reid  “A Note on Pseudolikelihood Constructed from Marginal Densities”  In _Biometrika_ 91.3, 2004, pp. 729–737  DOI: 10.1093/biomet/91.3.729
  * [8] Noel Cressie and Gardar Johannesson  “Fixed Rank Kriging for Very Large Spatial Data Sets”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 70.1, 2008, pp. 209–226  DOI: 10.1111/j.1467-9868.2007.00633.x
  * [9] A.. Davison, S.. Padoan and M. Ribatet  “Statistical Modeling of Spatial Extremes”  In _Statistical Science_ 27.2 The Institute of Mathematical Statistics, 2012, pp. 161–186  DOI: 10.1214/11-STS376
  * [10] Anthony C. Davison, Raphaël Huser and Emeric Thibaud  “Spatial Extremes”  In _Handbook of environmental and ecological statistics_ ChapmanHall/CRC, 2019, pp. 711–744
  * [11] P.. Diggle, J.. Tawn and R.. Moyeed  “Model-based Geostatistics”  In _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ 47.3, 1998, pp. 299–350  DOI: 10.1111/1467-9876.00113
  * [12] Jo Eidsvik et al.  “Estimation and Prediction in Spatial Models With Block Composite Likelihoods”  In _Journal of Computational and Graphical Statistics_ 23.2 Taylor & Francis, 2014, pp. 295–315  DOI: 10.1080/10618600.2012.760460
  * [13] Geir-Arne Fuglstad, Daniel Simpson, Finn Lindgren and Håvard Rue  “Constructing Priors That Penalize the Complexity of Gaussian Random Fields”  In _Journal of the American Statistical Association_ 114.525 Taylor & Francis, 2019, pp. 445–452  DOI: 10.1080/01621459.2017.1415907
  * [14] Alan E. Gelfand, Peter J. Diggle, Montserrat Fuentes and Peter Guttorp  “Handbook of Spatial Statistics”  CRC press, 2010  DOI: 10.1201/9781420072884
  * [15] Tilmann Gneiting and Adrian E. Raftery  “Strictly Proper Scoring Rules, Prediction, and Estimation”  In _Journal of the American Statistical Association_ 102.477 Taylor & Francis, 2007, pp. 359–378  DOI: 10.1198/016214506000001437
  * [16] V.. Godambe  “An Optimum Property of Regular Maximum Likelihood Estimation”  In _The Annals of Mathematical Statistics_ 31.4 Institute of Mathematical Statistics, 1960, pp. 1208–1211  URL: http://www.jstor.org/stable/2237819
  * [17] Janet E. Heffernan and Sidney I. Resnick  “Limit Laws for Random Vectors With an Extreme Component”  In _The Annals of Applied Probability_ 17.2 The Institute of Mathematical Statistics, 2007, pp. 537–571  DOI: 10.1214/105051606000000835
  * [18] Janet E. Heffernan and Jonathan A. Tawn  “A Conditional Approach for Multivariate Extreme Values (with discussion)”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 66.3, 2004, pp. 497–546  DOI: 10.1111/j.1467-9868.2004.02050.x
  * [19] Birgir Hrafnkelsson et al.  “Max-and-Smooth: A Two-Step Approach for Approximate Bayesian Inference in Latent Gaussian Models”  In _Bayesian Analysis_ 16.2 International Society for Bayesian Analysis, 2021, pp. 611–638  DOI: 10.1214/20-BA1219
  * [20] Raphaël Huser and Jennifer L. Wadsworth  “Advances in Statistical Modeling of Spatial Extremes”  In _WIREs Computational Statistics_ 14.1, 2022, pp. e1537  DOI: 10.1002/wics.1537
  * [21] Rikke Ingebrigtsen, Finn Lindgren and Ingelin Steinsland  “Spatial Models With Explanatory Variables in the Dependence Structure”  In _Spatial Statistics_ 8, 2014, pp. 20–38  DOI: 10.1016/j.spasta.2013.06.002
  * [22] Rikke Ingebrigtsen, Finn Lindgren, Ingelin Steinsland and Sara Martino  “Estimation of a Non-Stationary Model for Annual Precipitation in Southern Norway Using Replicates of the Spatial Field”  In _Spatial Statistics_ 14, 2015, pp. 338–364  DOI: 10.1016/j.spasta.2015.07.003
  * [23] Árni V. Jóhannesson et al.  “Approximate Bayesian Inference for Analysis of Spatiotemporal Flood Frequency Data”  In _The Annals of Applied Statistics_ 16.2 Institute of Mathematical Statistics, 2022, pp. 905–935  DOI: 10.1214/21-AOAS1525
  * [24] Cari G. Kaufman, Mark J. Schervish and Douglas W. Nychka  “Covariance Tapering for Likelihood-Based Estimation in Large Spatial Data Sets”  In _Journal of the American Statistical Association_ 103.484 Taylor & Francis, 2008, pp. 1545–1555  DOI: 10.1198/016214508000000959
  * [25] Caroline Keef, Ioannis Papastathopoulos and Jonathan A. Tawn  “Estimation of the Conditional Distribution of a Multivariate Variable Given That One of Its Components Is Large: Additional Constraints for the Heffernan and Tawn Model”  In _Journal of Multivariate Analysis_ 115, 2013, pp. 396–404  DOI: 10.1016/j.jmva.2012.10.012
  * [26] Caroline Keef, Jonathan A. Tawn and Rob Lamb  “Estimating the Probability of Widespread Flood Events”  In _Environmetrics_ 24.1, 2013, pp. 13–21  DOI: 10.1002/env.2190
  * [27] B.J.K. Kleijn and A.W. Vaart  “The Bernstein-Von-Mises Theorem under Misspecification”  In _Electronic Journal of Statistics_ 6 Institute of Mathematical StatisticsBernoulli Society, 2012, pp. 354–381  DOI: 10.1214/12-EJS675
  * [28] S. Kullback and R.. Leibler  “On Information and Sufficiency”  In _The Annals of Mathematical Statistics_ 22.1 Institute of Mathematical Statistics, 1951, pp. 79–86  URL: http://www.jstor.org/stable/2236703
  * [29] Finn Lindgren and Håvard Rue  “Bayesian Spatial Modelling with R-INLA”  In _Journal of Statistical Software_ 63.19, 2015, pp. 1–25  DOI: 10.18637/jss.v063.i19
  * [30] 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”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 73.4, 2011, pp. 423–498  DOI: 10.1111/j.1467-9868.2011.00777.x
  * [31] Bruce G. Lindsay  “Composite Likelihood Methods”  In _Contemporary mathematics_ 80.1, 1988, pp. 221–239
  * [32] Y. Liu and J.. Tawn  “Self-Consistent Estimation of Conditional Multivariate Extreme Value Distributions”  In _Journal of Multivariate Analysis_ 127, 2014, pp. 19–35  DOI: 10.1016/j.jmva.2014.02.003
  * [33] Luigi Lombardo et al.  “Space-Time Landslide Predictive Modelling”  In _Earth-Science Reviews_ 209, 2020, pp. 103318  DOI: 10.1016/j.earscirev.2020.103318
  * [34] T. Lumley and P. Heagerty  “Weighted Empirical Adaptive Variance Estimators for Correlated Data Regression”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 61.2, 1999, pp. 459–477  DOI: 10.1111/1467-9868.00187
  * [35] Jean-Michel Marin, Pierre Pudlo, Christian P. Robert and Robin J. Ryder  “Approximate Bayesian Computational Methods”  In _Statistics and Computing_ 22.6, 2012, pp. 1167–1180  DOI: 10.1007/s11222-011-9288-2
  * [36] Thomas Opitz, Haakon Bakka, Raphaël Huser and Luigi Lombardo  “High-resolution Bayesian Mapping of Landslide Hazard with Unobserved Trigger Event”  In _The Annals of Applied Statistics_ 16.3 Institute of Mathematical Statistics, 2022, pp. 1653–1675  DOI: 10.1214/21-AOAS1561
  * [37] Thomas Opitz, Raphaël Huser, Haakon Bakka and Håvard Rue  “INLA Goes Extreme: Bayesian Tail Regression for the Estimation of High Spatio-Temporal Quantiles”  In _Extremes_ 21.3, 2018, pp. 441–462  DOI: 10.1007/s10687-018-0324-x
  * [38] Francesco Pauli, Walter Racugno and Laura Ventura  “Bayesian Composite Marginal Likelihoods”  In _Statistica Sinica_ 21.1 Institute of Statistical Science, Academia Sinica, 2011, pp. 149–164  URL: http://www.jstor.org/stable/24309266
  * [39] Mathieu Ribatet, Daniel Cooley and Anthony C. Davison  “Bayesian Inference From Composite Likelihoods, With an Application To Spatial Extremes”  In _Statistica Sinica_ 22.2 Institute of Statistical Science, Academia Sinica, 2012, pp. 813–845  URL: http://www.jstor.org/stable/24310036
  * [40] Jordan Richards, Jonathan A. Tawn and Simon Brown  “Modelling Extremes of Spatial Aggregates of Precipitation using Conditional Methods”  In _The Annals of Applied Statistics_ 16.4 Institute of Mathematical Statistics, 2022, pp. 2693–2713  DOI: 10.1214/22-AOAS1609
  * [41] Håvard Rue, Sara Martino and Nicolas Chopin  “Approximate Bayesian Inference for Latent Gaussian Models By Using Integrated Nested Laplace Approximations”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 71.2, 2009, pp. 319–392  DOI: 10.1111/j.1467-9868.2008.00700.x
  * [42] Håvard Rue et al.  “Bayesian Computing With INLA: a Review”  In _Annual Review of Statistics and Its Application_ 4.1, 2017, pp. 395–421  DOI: 10.1146/annurev-statistics-060116-054045
  * [43] Rob Shooter, Jonathan Tawn, Emma Ross and Philip Jonathan  “Basin-wide Spatial Conditional Extremes for Severe Ocean Storms”  In _Extremes_ 24.2, 2021, pp. 241–265  DOI: 10.1007/s10687-020-00389-w
  * [44] Daniel Simpson et al.  “Penalising Model Component Complexity: a Principled, Practical Approach To Constructing Priors”  In _Statistical Science_ 32.1 The Institute of Mathematical Statistics, 2017, pp. 1–28  DOI: 10.1214/16-STS576
  * [45] Emma S. Simpson, Thomas Opitz and Jennifer L. Wadsworth  “High-dimensional Modeling of Spatial and Spatio-Temporal Conditional Extremes using INLA and Gaussian Markov Random Fields” doi: 10.48550/arxiv.2011.04486  arXiv, 2020  DOI: 10.48550/arxiv.2011.04486
  * [46] Emma S. Simpson and Jennifer L. Wadsworth  “Conditional Modelling of Spatio-Temporal Extremes for Red Sea surface Temperatures”  In _Spatial Statistics_ 41, 2021, pp. 100482  DOI: 10.1016/j.spasta.2020.100482
  * [47] Scott A Sisson, Yanan Fan and Mark Beaumont  “Handbook of approximate Bayesian computation”  CRC Press, 2018
  * [48] Nicholas Syring and Ryan Martin  “Calibrating General Posterior Credible Regions”  In _Biometrika_ 106.2, 2018, pp. 479–486  DOI: 10.1093/biomet/asy054
  * [49] Silius M. Vandeskog, Sara Martino, Daniela Castro-Camilo and Håvard Rue  “Modelling Sub-daily Precipitation Extremes with the Blended Generalised Extreme Value Distribution”  In _Journal of Agricultural, Biological and Environmental Statistics_ , 2022  DOI: 10.1007/s13253-022-00500-7
  * [50] Cristiano Varin, Nancy Reid and David Firth  “An Overview of Composite Likelihood Methods”  In _Statistica Sinica_ 21.1 Institute of Statistical Science, Academia Sinica, 2011, pp. 5–42  URL: http://www.jstor.org/stable/24309261
  * [51] A.. Vecchia  “Estimation and Model Identification for Continuous Spatial Processes”  In _Journal of the Royal Statistical Society: Series B (Methodological)_ 50.2, 1988, pp. 297–312  DOI: 10.1111/j.2517-6161.1988.tb01729.x
  * [52] J.L. Wadsworth and J.A. Tawn  “Higher-dimensional Spatial Extremes via Single-site Conditioning”  In _Spatial Statistics_ 51, 2022, pp. 100677  DOI: 10.1016/j.spasta.2022.100677
  * [53] Stephen G. Walker  “Bayesian Inference With Misspecified Models”  In _Journal of Statistical Planning and Inference_ 143.10, 2013, pp. 1621–1633  DOI: 10.1016/j.jspi.2013.05.013
  * [54] Halbert White  “Maximum Likelihood Estimation of Misspecified Models”  In _Econometrica_ 50.1 [Wiley, Econometric Society], 1982, pp. 1–25  URL: http://www.jstor.org/stable/1912526

## Supplementary material

##  Appendix S1 Implementing the spatial conditional extremes model in R-INLA

We implement the conditional extremes model,

$\left[Y(\bm{s})\mid
Y(\bm{s}_{0})=y_{0}>t\right]\overset{d}{=}a(\bm{s};\bm{s}_{0},y_{0})+b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0})+\epsilon(\bm{s};\bm{s}_{0}),$
(S1.1)

in R-INLA. Here, $Z(\bm{s};\bm{s}_{0})$ is a constrained Gaussian Matérn field
and $\epsilon(\bm{s};\bm{s}_{0})$ is constrained Gaussian white noise. See
Section 3 for more information.

R-INLA performs inference on latent Gaussian models of the form

$\displaystyle\left[y_{i}\mid\bm{u},\bm{\theta}_{1}\right]$
$\displaystyle\overset{i.i.d.}{\sim}\mathcal{\pi}(y_{i}\mid\bm{\eta}(\bm{u}),\bm{\theta}_{1}),\
i=1,2,\ldots,n,$ $\displaystyle\left[\bm{u}\mid\bm{\theta}_{2}\right]$
$\displaystyle\sim\mathcal{N}(\bm{\mu}(\bm{\theta}_{2}),\bm{Q}^{-1}(\bm{\theta}_{2})),$
$\displaystyle\left(\bm{\theta}_{1},\bm{\theta}_{2}\right)^{T}$
$\displaystyle\sim\pi(\bm{\theta}_{1})\pi(\bm{\theta}_{2}),$

where $\bm{u}$ is a latent Gaussian field with mean
$\bm{\mu}(\bm{\theta}_{2})$ and precision matrix $\bm{Q}(\bm{\theta}_{2})$.
Hyperparameters $\bm{\theta}=(\bm{\theta}_{1},\bm{\theta}_{2})^{T}$, are
assigned priors $\pi(\bm{\theta}_{1})$ and $\pi(\bm{\theta}_{2})$.
Observations $\bm{y}=(y_{1},\ldots,y_{n})^{T}$ are linked to the latent field
through the linear predictor $\bm{\eta}(\bm{u})=\bm{A}\bm{u}$, where $\bm{A}$
is a known design matrix. This linear predictor defines the location parameter
of the likelihood, via a possibly non-linear link function. In the R-INLA
implementation, the observations are assumed to be conditionally independent
given $\bm{\eta}$ and $\bm{\theta}_{1}$, so that
$\pi(\bm{y}\mid\bm{\eta}(\bm{u}),\bm{\theta}_{1})=\prod_{i}\pi(y_{i}\mid\eta_{i}(\bm{u}),\bm{\theta}_{1})$.
As mentioned in Section 1, the likelihood must be chosen from a predefined set
of likelihood functions and the R-INLA user cannot modify it. The linear
predictor can be decomposed into $d\geq 1$ components,
$\bm{\eta}(\bm{u})=\bm{A}^{(1)}\bm{u}^{(1)}+\dots+\bm{A}^{(d)}\bm{u}^{(d)}$,
where each component represents, e.g., an intercept term, a linear combination
of regression coefficients, an SPDE component, etc. All of these components
must either be predefined in R-INLA or defined by the user, using the rgeneric
framework or the recently added cgeneric framework.

The model in (S1.1) corresponds to a latent Gaussian model where the
likelihood is Gaussian with variance $\theta_{1}$ and the linear predictor is
equal to
$a(\bm{s};\bm{s}_{0},y_{0})+b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0})$,
where $\bm{\theta}_{2}$ contains the parameters of $a(\cdot)$, $b(\cdot)$ and
$Z(\cdot)$. In our implementation we use the cgeneric framework to define the
linear predictor as a sum of the two components $a(\bm{s};\bm{s}_{0},y_{0})$
(see Section S1.2) and $b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0})$ (see
Section S1.1).

[45] employ the rgeneric model for implementing their version of the spatial
conditional extremes model in R-INLA. However, rgeneric is considerably slower
than cgeneric, which did not yet exist at the time of their work. The addition
of the cgeneric framework in R-INLA allows us to implement a wide selection of
forms for $a(\cdot)$ and $b(\cdot)$ at a lower computational cost.

###  S1.1 Implementing $b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0})$ in
R-INLA

In the following we describe how to implement the constrained, non-stationary
field $b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0})$ using the cgeneric
framework. We first describe how to implement a non constrained non-stationary
field and then how to include the constrain.

The stochastic partial differential equation (SPDE) approach of [30] creates a
Gaussian Markov random field $\hat{Z}(\bm{s})$ that is an approximation to a
Gaussian random field $Z(\bm{s})$ with Matérn covariance function:

$\text{Cov}(Z(\bm{s}),Z(\bm{s}^{\prime}))=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}(\kappa\|\bm{s}-\bm{s}^{\prime}\|)^{\nu}K_{\nu}(\kappa\|\bm{s}-\bm{s}^{\prime}\|),$

where $\sigma^{2}$ is the marginal variance, $\nu>0$ is the smoothness
parameter and $\rho=\sqrt{8\nu}/\kappa$ is the range parameter of $Z(\bm{s})$.
Furthermore, $K_{\nu}$ is the modified Bessel function of the second kind and
order $\nu$. The smoothness parameter $\nu$ is difficult to estimate from data
and is therefore often given a fixed value [29]. The SPDE approximation
$\hat{Z}(\bm{s})$ is constructed as a linear combination of Gaussian Markov
random variables on a triangular mesh:

$\hat{Z}(\bm{s})=\sum_{i=1}^{n}\phi_{i}(\bm{s})W_{i}\overset{d}{\approx}Z(\bm{s}),$

where $W_{1},\ldots,W_{n}$ are Gaussian Markov random variables, and
$\phi_{i},\ldots,\phi_{n}$ are piecewise linear basis functions. Helper
functions for creating the triangular mesh and the precision matrix of
$\bm{W}=(W_{1},W_{2},\ldots,W_{n})^{T}$ are available in R-INLA. In order to
approximate the non-stationary Gaussian random field
$b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s})$ with the SPDE approach, for any function
$b(\bm{s};\bm{s}_{0},y_{0})$, we modify the weights $W_{i}$ to get

$\hat{Z}_{b}(\bm{s};\bm{s}_{0},y_{0})=\sum_{i=1}^{n}\phi_{i}(\bm{s})b(\bm{s}_{i};\bm{s}_{0},y_{0})W_{i}\overset{d}{\approx}b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s}),$
(S1.2)

where $\bm{s}_{1},\ldots\bm{s}_{n}$ are the locations of the $n$ Gaussian
weights in the triangular mesh. This idea shares some similarities with the
approach of [21, 22] for implementing non-stationary SPDE fields.

Since $\hat{Z}_{b}(\cdot)$ is created as a linear combination of Gaussian
random variables, its variance equals

$\text{Var}\left(\hat{Z}_{b}(\bm{s};\bm{s}_{0},y_{0})\right)=\sum_{i,j=1}^{n}\phi_{i}(\bm{s})\phi_{j}(\bm{s})b(\bm{s}_{i};\bm{s}_{0},y_{0})b(\bm{s}_{j};\bm{s}_{0},y_{0})\text{Cov}\left(W_{i},W_{j}\right),$

which in general is unequal to

$\text{Var}\left(b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s})\right)=b(\bm{s};\bm{s}_{0},y_{0})^{2}\sigma^{2}.$

However, if $\bm{s}=\bm{s}_{i}$ for $i=1,2,\ldots,n$, then $\phi_{i}=1$ while
the other basis functions equal zero, resulting in a variance of

$\text{Var}\left(\hat{Z}_{b}(\bm{s};\bm{s}_{0},y_{0})\right)=b(\bm{s};\bm{s}_{0},y_{0})^{2}\text{Var}(W_{i}),$

which is much closer to the correct variance. On the other hand, if
$b(\bm{s};\bm{s}_{0},y_{0})$ changes quickly and $\bm{s}$ is far away from a
mesh node, then the variance of $\hat{Z}_{b}(\bm{s};\bm{s}_{0},y_{0})$ can be
considerably different from the variance of
$b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s})$. Thus, if the mesh used for defining
$\hat{Z}_{b}(\cdot)$ is dense and all data locations in $\mathcal{S}$ are
close to a mesh node, $\hat{Z}_{b}(\cdot)$ can provide a good approximation to
$b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s})$, but if the mesh is too coarse, the
approximation may deteriorate in performance. It is a commonly known problem
that using a too coarse mesh when approximating Gaussian random fields with
the SPDE approach is problematic. This is also demonstrated in the simulation
study in Section 4.2. However, the addition of a possible nonlinear function
$b(\cdot)$, which turns the Gaussian random field into a nonstationary random
field, makes the requirements on the mesh even stricter. Users of the SPDE
approach should therefore be aware that a mesh that is good enough for
approximating a stationary Gaussian field $Z(\bm{s})$ might not be good enough
for approximating the non-stationary Gaussian field
$b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s})$.

####  S1.1.1 Constraining the residual field

$\hat{Z}_{b}(\cdot)$ provides a good approximation to the unconstrained
Gaussian random field $b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s})$. However, in order
to implement the conditional extremes model in R-INLA, we need to approximate
the constrained field $b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0})$, where
$Z(\bm{s}_{0};\bm{s}_{0})=0$ almost surely. [52] describe two different
methods for turning an unconstrained Gaussian field $Z(\bm{s})$ into a
constrained field $Z(\bm{s};\bm{s}_{0})$. The first is to constrain the field
by conditioning: $Z(\bm{s};\bm{s}_{0})=[Z(\bm{s})\mid Z(\bm{s})=0]$, and the
second is to constrain by subtraction:
$Z(\bm{s};\bm{s}_{0})=Z(\bm{s})-Z(\bm{s}_{0})$. [52] use the first method
while [45] use the second method in their case study. We argue that
constraining by subtraction provides unrealistic dependence structures, and
should be avoided if other alternatives are available. A quick computation
shows that if $Z(\bm{s})$ is a stationary random process, then

$\lim_{c\rightarrow\infty}\text{Corr}(Z(c\bm{s};\bm{s}_{0}),Z(-c\bm{s};\bm{s}_{0}))=\frac{1}{2},$
(S1.3)

and if $Z(\bm{s})$ is isotropic with autocorrelation
$\text{Corr}(Z(\bm{s}),Z(\bm{s}^{\prime}))=r(\|\bm{s}-\bm{s}^{\prime}\|)$,
then

$\text{Corr}(Z(\bm{s}_{0}+\Delta\bm{s};\bm{s}_{0}),Z(\bm{s}_{0}-\Delta\bm{s};\bm{s}_{0}))=\frac{1-2r(\|\Delta\bm{s}\|)+r(2\|\Delta\bm{s}\|)}{2(1-r(\|\Delta\bm{s}\|)}.$
(S1.4)

The limit of (S1.4) as $\|\Delta\bm{s}\|\rightarrow 0$ is often negative, and
equals $0$ for the exponential autocorrelation function and $-1$ for the
Gaussian autocorrelation function. Thus, points that are infinitely far away
from each other are strongly correlated while points that are infinitesimally
close to each other might be negatively correlated or independent. Figure S1.1
displays the correlation structure of a Gaussian Matérn field that has been
constrained by subtraction.

Figure S1.1: Correlation between $Z(\bm{s}_{1};\bm{0})$ and
$Z(\bm{s}_{2};\bm{0})$, where $Z(\bm{s};\bm{s}_{0})=Z(\bm{s})-Z(\bm{s}_{0})$
and $Z(\bm{s})$ is a stationary Gaussian Matérn random field with range $\rho$
and smoothness $\nu=1.5$. Here, $\bm{s}_{1}=(x_{1},0)^{T}$ and
$\bm{s}_{2}=(x_{2},0)^{T}$.

This is clearly not a preferable correlation structure for most spatial
processes.

R-INLA contains an implementation for constraining a random field by
conditioning on linear combinations of itself [41]. One can therefore easily
constrain $\hat{Z}_{b}(\cdot)$ by conditioning using the extraconstr option in
R-INLA. Unfortunately, this conditioning adds an extra computational cost that
can become huge when dealing with high-dimensional spatial or spatio-temporal
models. In practice, when using R-INLA with large data sets, we experienced
that constraining by conditioning required more than ten times the amount of
computer memory and computation time than constraining by subtraction.
Constraining by conditioning therefore quickly turns infeasible for large data
sets.

A third method for constraining the residual field is used by [40]. They model
the variance of the residual field as a function of the distance
$\|\bm{s}-\bm{s}_{0}\|$ such that the variance is zero when
$\|\bm{s}-\bm{s}_{0}\|=0$. This can easily be implemented by modifying the
standardising function $b(\cdot)$ in $\hat{Z}_{b}(\cdot)$ such that
$b(\bm{s}_{0};\bm{s}_{0},y_{0})=0$. As an example, one can multiply
$b(\bm{s};\bm{s}_{0},y_{0})$ with the function

$\sqrt{1-\exp\left(-2\|\bm{s}-\bm{s}_{0}\|/\rho_{b}\right)},$

where $\rho_{b}>0$ can be either fixed or estimated. This expression is
inspired by the variance of a conditional bivariate Gaussian distribution with
exponential correlation function. In general, constraining
$\hat{Z}_{b}(\cdot)$ by requiring that $b(\bm{s}_{0};\bm{s}_{0},y_{0})=0$ is
almost as computationally efficient as constraining by subtraction, but
without any of the undesirable properties discussed above. In our opinion one
should therefore always prefer to constrain $\hat{Z}_{b}(\cdot)$ in this way
when modelling spatial conditional extremes with R-INLA.

###  S1.2 Implementing $a(\bm{s};\bm{s}_{0},y_{0})$ in R-INLA

Based on the R-INLA implementation of
$b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s};\bm{s}_{0})$ in Section S1.1, it might
seem reasonable to simply represent the entire conditional extremes model
using the SPDE approximation

$\hat{Z}_{a,b}(\bm{s};\bm{s}_{0},y_{0})=\sum_{i=1}^{n}\phi_{i}(\bm{s})\left(a(\bm{s}_{i};\bm{s}_{0},y_{0})+b(\bm{s}_{i};\bm{s}_{0},y_{0})W_{i}\right)\overset{d}{\approx}a(\bm{s};\bm{s}_{0},y_{0})+b(\bm{s};\bm{s}_{0},y_{0})Z(\bm{s}).$

However, this approximation has a mean equal to

$\text{E}\left[\hat{Z}_{a,b}(\bm{s};\bm{s}_{0},y_{0})\right]=\sum_{i,j=1}^{n}\phi_{i}(\bm{s})\phi_{j}(\bm{s})a(\bm{s}_{i};\bm{s}_{0},y_{0}),$

which in general is different from $a(\bm{s};\bm{s}_{0},y_{0})$. Unlike in the
previous section, there is a way of improving this approximation by
representing $a(\cdot)$ using a separate random field that is independent of
the SPDE approximation.

Given any deterministic vector $\bm{a}$, we claim that the Gaussian random
field with mean $\bm{a}$ and covariance matrix $\delta^{2}\bm{I}$, where
$\bm{I}$ is the identity matrix and $\delta^{2}$ is a fixed marginal variance
close to zero, provides a good approximation to $\bm{a}$. Thus, we can
approximate any standardising function $a(\bm{s};\bm{s}_{0},y_{0})$ using a
cgeneric model component with mean $a(\bm{s};\bm{s}_{0},y_{0})$ and diagonal
covariance matrix $\delta^{2}\bm{I}$, where we set $\delta^{2}=\exp(-15)$.

##  Appendix S2 A second simulation study with the spatial conditional
extremes distribution

In the simulation study in Section 4.4, we found that, even though the
adjustment method improved the frequency properties of the fitted conditional
extremes models, it did not improve them as much as in the other three
simulation studies in Section 4. In order to test our theory that this
discrepancy is caused by self-inconsistency that complicates the distribution
of the simulated data, we perform a similar simulation study where we fit the
spatial conditional extremes distribution to extremes of a Gaussian random
process, instead of simulated data from the global conditional extremes
distribution.

First, we define a rectangular grid, $\mathcal{S}$, with resolution $1\times
1$ and size $15\times 15$. Then, we simulate data in $\mathcal{S}$ using the
SPDE approximation with variance $\sigma^{2}=1$ and range $\rho=8$.
Standardising the data to have Laplace marginals is easy, since the marginal
distributions of the simulated data are known. We choose a threshold $t$ equal
to the $95\%$ quantile of the Laplace distribution, and simulate data until we
have observed $5\times 10^{4}$ realisations such that each transformed
realisation exceeds the threshold $t$ in at least one location in
$\mathcal{S}$. Then, we perform the pairwise parameter estimation procedure
described in Section S3.2 to examine the best forms for
$a(\bm{s};\bm{s}_{0},y_{0})$ and $b(\bm{s};\bm{s}_{0},y_{0})$ and to test
whether the chosen forms from Section 4.4 can provide a reasonable fit to the
transformed data.

The pairwise parameter estimation procedure shows that the parametric forms
from Section 4.4 provide a good enough fit to the data, and using maximum
likelihood estimation with the global composite likelihood, we find that the
Kullback-Leibler distance [28, KLD;] minimiser is

$\bm{\theta}^{*}=(\lambda^{*},\kappa^{*},\sigma^{*},\rho_{b}^{*},\rho^{*},\tau^{*})^{T}\approx(3.9,1.1,1.4,0.2,6.8,23.4)^{T}.$

Table 6: Coverage percentages for unadjusted and adjusted credible intervals based on the composite conditional extremes likelihood with fixed $\rho_{b}$. Aim | $\lambda$ | $\lambda_{\text{adj}}$ | $\kappa$ | $\kappa_{\text{adj}}$ | $\rho$ | $\rho_{\text{adj}}$ | $\sigma$ | $\sigma_{\text{adj}}$ | $\tau$ | $\tau_{\text{adj}}$
---|---|---|---|---|---|---|---|---|---|---
$90\%$ | $69\%$ | $91\%$ | $83\%$ | $88\%$ | $66\%$ | $87\%$ | $48\%$ | $88\%$ | $63\%$ | $87\%$
$95\%$ | $77\%$ | $95\%$ | $89\%$ | $94\%$ | $70\%$ | $92\%$ | $53\%$ | $94\%$ | $74\%$ | $94\%$
$99\%$ | $89\%$ | $98\%$ | $97\%$ | $99\%$ | $82\%$ | $98\%$ | $64\%$ | $97\%$ | $88\%$ | $97\%$

We evaluate the adjustment method by simulating 500 new realisation of the
transformed Gaussian random field and fitting the global conditional extremes
distribution to the available threshold exceedances with R-INLA. Then, the
output from R-INLA is adjusted using the posterior adjustment method. This is
repeated 300 times. The low value of $\rho_{b}^{*}$ makes it difficult to
estimate this parameter consistently, so we choose to fix the value of
$\rho_{b}$ to the estimated value of $\rho_{b}^{*}$ during inference, and only
compute posteriors for the five remaining parameters. Table 6 displays
estimated coverage probabilities for how often the parameters in
$\bm{\theta}^{*}$ are included in adjusted and unadjusted credible intervals.
There is a considerable improvement for the adjusted credible intervals, and
the overall frequency properties of the adjusted posteriors seem to be on par
with those of the adjusted posteriors in Section 4.1, 4.2 and 4.3. The only
differences between this simulation study and the one in Section 4.4 is that
we simulated extreme data in different ways, and that we did not compute the
posterior of $\rho_{b}$. However, we also repeat the simulation study from
Section 4.4 where $\rho_{b}$ is fixed and not estimated. The results are
similar to those in Section 4.4: The adjusted posterior displays improved
frequency properties, but these are not as good as in all other simulation
studies where we sample from self-consistent distributions. This implies that
the observed problems in Section 4.4 are caused by the difficulties that arise
when simulating data from a self-inconsistent distribution.

##  Appendix S3 Case study preliminaries

###  S3.1 Transforming the data to obtain Laplace marginals

The conditional extremes model in (S1.1) is defined for a random process with
Laplace margins. However, most natural processes are not Laplace distributed.
Thus, in order to perform inference with the conditional extremes model for a
random process $X(\bm{s})$, one must first standardise it to a random process
$Y(\bm{s})$ with Laplace margins. This is performed using the probability
integral transform [25]:

$Y(\bm{s})=\begin{cases}\log\left\\{2F_{X(\bm{s})}(X(\bm{s}))\right\\},&X(\bm{s})<F_{X(\bm{s})}(1/2)\\\
-\log\left\\{2\left[1-F_{X(\bm{s})}(X(\bm{s}))\right]\right\\},&X(\bm{s})\geq
F_{X(\bm{s})}(1/2),\end{cases}$

where $F_{X(\bm{s})}$ is the marginal distribution function of the random
variable $X(\bm{s})$. We estimate the marginal distribution functions as the
empirical distribution functions. However, independent standardisation of data
at each location can lead to an unrealistic lack of smoothness in the
transformed process $Y(\bm{s})$. Therefore, we apply a sliding window approach
for computing the empirical distribution functions, where the distribution at
location $\bm{s}$ is estimated as the empirical distribution function of
pooled data from all locations $\bm{s}^{\prime}$ such that
$\|\bm{s}-\bm{s}^{\prime}\|\leq r$ for some radius $r$. Based on exploratory
analysis we find $r=5$ km to yield a realistic degree of smoothness in the
estimated marginal distributions of $X(\bm{s})$ (results not shown).

A problem with modelling precipitation is that the empirical distribution has
a point mass in zero. This leads to $Y(\bm{s})$ having a truncated Laplace
distribution with a point mass, which can cause problems during inference. In
order for $Y(\bm{s})$ to follow a non-truncated Laplace distribution, we
choose to remove all zeros from the process $X(\bm{s})$ and only focus on
positive precipitation. This makes us unable to model the absence of
precipitation, which can lead to overestimation of return levels for spatially
aggregated precipitation. However, the main focus of this case study is to fit
a Bayesian version of the spatial conditional extremes model to data and to
improve the posterior with our adjustment method, and we do not focus on
applying the fitted model for estimating properties of the untransformed
process $X(\bm{s})$. We believe that our choices of removing all zeros and
estimating marginals using empirical distribution functions of the positive
precipitation values are acceptable given the aim of the current case study.
More sophisticated transformations of the marginal distributions might be
necessary if the aim is to predict or estimate properties of the untransformed
process $X(\bm{s})$.

###  S3.2 Choosing models for $a(\bm{s};\bm{s}_{0},y_{0})$ and
$b(\bm{s};\bm{s}_{0},y_{0})$

Due to the relatively small size of the domain $\mathcal{S}$, described in
Section 5.1, we assume that the distribution of hourly precipitation data is
approximately stationary over all locations in the sense that $a(\cdot)$,
$b(\cdot)$, $Z(\cdot)$ and $\epsilon(\cdot)$ from (S1.1) are the same for all
conditioning sites. We also assume that the data are isotropic and that
$a(\bm{s};\bm{s}_{0},\cdot)$ and $b(\bm{s};\bm{s}_{0},\cdot)$ only depend on
the distance $\|\bm{s}-\bm{s}_{0}\|$ between $\bm{s}$ and $\bm{s}_{0}$. Based
on the suggestions and results of [52, 45, 40], we aim to model the
precipitation data using the model in (S1.1) with standardising functions

$\begin{aligned} &a(\bm{s};\bm{s}_{0},y_{0})=\alpha(d)\cdot
y_{0}+\gamma(d),\\\ &b(\bm{s};\bm{s}_{0},y_{0})=\zeta(d)\cdot
y_{0}^{\beta(d)},\end{aligned}\quad d=\|\bm{s}-\bm{s}_{0}\|,$ (S3.1)

where $\alpha(d)$, $\gamma(d)$, $\zeta(d)$ and $\beta(d)$ are parametric
functions with the constraints that $\alpha(0)=1$, $\gamma(0)=0$ and
$\zeta(0)=0$. Furthermore, $\alpha(d)$, $\gamma(d)$ and $\beta(d)$ are all
assumed to converge to zero when the distance is so large that $Y(\bm{s})$ and
$Y(\bm{s}_{0})$ are independent. We also assume that the residual field
$Z(\bm{s};\bm{s}_{0})$ is standardised to have zero mean and a variance of
one.

Similarly to [40], we get a sense of the best shape of the functions in (S3.1)
for modelling our data, by first assuming that the residual field consists of
Gaussian white noise, meaning that all observations are conditionally
independent given $Y(\bm{s}_{0})$. Two locations, $\bm{s}_{0}$ and
$\bm{s}_{1}$, are then drawn randomly from $\mathcal{S}$, and the conditional
model

$\left[Y(\bm{s}_{1})\mid Y(\bm{s}_{0})=y_{0}>t\right]\overset{d}{=}\alpha\cdot
y_{0}+\gamma+\zeta\cdot y_{0}^{\beta}Z$ (S3.2)

is fitted to the observations from $\bm{s}_{0}$ and $\bm{s}_{1}$, where $Z$ is
a standardised Gaussian random variable. This pairwise model fitting is
repeated for 5000 random location pairs in $\mathcal{S}$. We set the threshold
$t$ equal to the $99.75\%$ quantile of the Laplace distribution. This
corresponds to an average of between one and three threshold exceedances per
year of hourly summer data. In column 4 of Figure S3.1, the 5000 parameter
estimates are displayed as functions of the distances between $\bm{s}_{1}$ and
$\bm{s}_{0}$. The parameter estimates are so full of noise that it is
difficult to find any clear patterns. Consequently, we repeat the entire
pairwise fitting procedure three more times while fixing either $\gamma=0$,
$\beta=0$ or both $\gamma=0$ and $\beta=0$. This allows clearer patterns to
emerge from the data (see the three remaining columns in Figure S3.1). Note
that the parameter estimates in Figure S3.1 have been found using constrained
optimisation, as unconstrained optimisation sometimes lead to estimates of
$\zeta$ larger than $10^{4}$ when both $\beta$ and $\zeta$ are estimated
simultaneously, and other similar problems. The constrained optimisation
removes outliers like these to make it possible to extract useful information
from the plots.

Figure S3.1: Parameter estimates from pairwise model fits using different
pairs of $\bm{s}_{0}$ and $\bm{s}_{1}$. Model 1 is model (S3.2) with fixed
$\beta=\gamma=0$. Model 2 is the same model with fixed $\gamma=0$. Model 3 is
the same model with fixed $\beta=0$. Model 4 is the same model with all four
parameters freely estimated.

In order to choose parametric forms for the functions in (S3.1), we first note
that it is possible for both $\gamma(d)$ and $\beta(d)$ to be constantly equal
to zero. However, $\alpha(d)$ must decrease from $\alpha(0)=1$ to
$\alpha(\infty)=0$, while $\zeta(d)$ must increase from $\zeta(0)=0$ to
$\zeta(d)>0$ for $d>0$. Therefore, we consider $\alpha(d)$ and $\zeta(d)$ to
be the most fundamental parts of (S3.1), and in an attempt to get rid of the
identifiability issues, we first choose parametric functions for $\alpha(d)$
and $\zeta(d)$ while we fix the functions $\gamma(d)$ and $\beta(d)$ equal to
zero. The typical parametric choice for $\alpha(d)$ is

$\alpha(d)=\exp\left\\{-\left(d/\lambda\right)^{\kappa}\right\\},$

with parameters $\lambda>0$ and $\kappa>0$ [52, 40]. This function corresponds
well to the trends seen in Figure S3.1, so we choose it for modelling
$\alpha(d)$. The pairwise estimates for $\zeta(d)$ shows that it is almost
constant for large distances, and that it approaches zero for short distances.
As discussed in Section S1.1.1, we therefore choose the parametric function

$\zeta(d)=\sigma\sqrt{1-\exp\left(-2d/\rho_{b}\right)},$

with parameters $\sigma>0$ and $\rho_{b}>0$. The parameters of $\alpha(d)$ and
$\zeta(d)$ are estimated using maximum likelihood estimation, where we still
assume that the residual field consists of Gaussian white noise. The
estimators $\hat{\alpha}(d)$ and $\hat{\zeta}(d)$ are displayed in the two
leftmost subplots of Figure S3.2. They seem to agree well with the pairwise
estimates from Figure S3.1.

Having estimated $\alpha(d)$ and $\zeta(d)$, we compute the residuals

$\left[\frac{Y(\bm{s}_{1})-\hat{\alpha}(\|\bm{s}_{1}-\bm{s}_{0}\|)Y(\bm{s}_{0})}{\hat{\zeta}\left(\|\bm{s}_{1}-\bm{s}_{0}\|\right)}\
\bigg{|}\ Y(\bm{s}_{0})>t\right],\\\ $ (S3.3)

for all 5000 pairs of $\bm{s}_{1}$ and $\bm{s}_{0}$. We compute the mean and
variance of these residuals, and find that the residual mean is close to zero,
while the residual variance is close to one. The pairwise estimation procedure
used for creating Figure S3.1 is then repeated on the residuals, in order to
examine trends in $\gamma(d)$ and $\beta(d)$. Pairwise estimates of $\gamma$
and $\beta$ are displayed in the two rightmost subplots of Figure S3.2. The
patterns in the pairwise estimates are weak, full of noise and centred around
zero. Based on these findings we choose to set $\gamma(d)=\beta(d)=0$ and we
conclude that the spatial conditional extremes model with

$\begin{aligned} &a(\bm{s};\bm{s}_{0},y_{0})=\alpha(d)\cdot
y_{0}=y_{0}\exp\left\\{-\left(d/\lambda\right)^{\kappa}\right\\},\\\
&b(\bm{s};\bm{s}_{0},y_{0})=\zeta(d)=\sigma\sqrt{1-\exp\left(-2d/\rho_{b}\right)},\end{aligned}\quad
d=\|\bm{s}-\bm{s}_{0}\|,$ (S3.4)

provides a good enough fit to the data. We note that, clearly, $\zeta(d)$ does
not provide a perfect fit to the trends of the pairwise estimates in Figure
S3.1, and a better approximation would likely decrease slightly with the
distance $d$. There are also still some weak trends remaining in the pairwise
estimates for $\beta$. However, as mentioned in Section S3.1, our aim is not
to estimate return levels or similar properties of the precipitation data, but
simply to fit a Bayesian version of the spatial conditional extremes model to
data and to improve the posterior using our adjustment method. The chosen
models for $a(\bm{s};\bm{s}_{0},y_{0})$ and $b(\bm{s};\bm{s}_{0},y_{0})$ are
therefore good enough for our intended goal, even though they might not be
good enough for, e.g., reliable estimation of return levels.

Figure S3.2: Left: Pairwise parameter estimates from column 1 of Figure S3.1,
in grey, together with the maximum likelihood estimators of $\alpha(d)$ and
$\zeta(d)$, displayed as black lines. Right: Pairwise parameter estimates of
$\beta$ and $\gamma$ based on fitting the pairwise conditional extremes model
(S3.2) to the residuals in (S3.3).