# The Self-organized Criticality Behaviors of Two New Parameters in SGR
J1935+2154

Shuo Xiao,1,2 Shuang-Nan Zhang,3,4 Shao-Lin Xiong,3 Ping Wang,3 Xiu-Juan Li,5
Ai-Jun Dong,1,2 Qi-Jun Zhi,1,2 Di Li6,4
1School of Physics and Electronic Science, Guizhou Normal University, Guiyang
550001, People’s Republic of China
2Guizhou Provincial Key Laboratory of Radio Astronomy and Data Processing,
Guizhou Normal University, Guiyang 550001, People’s Republic of China
3Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China
4University of Chinese Academy of Sciences, Chinese Academy of Sciences,
Beijing 100049, China
5School of Cyber Science and Engineering, Qufu Normal University, Qufu 273165,
People’s Republic of China
6CAS Key Laboratory of FAST, NAOC, Chinese Academy of Sciences, Beijing
100101, China
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The minimum variation timescale (MVT) and spectral lag of hundreds of X-ray
bursts (XRBs) from soft gamma-ray repeater (SGR) J1935+2154 were analyzed in
detail for the first time in our recent work, which are important probes for
studying the physical mechanism and radiation region. In this work, we
investigate their differential and cumulative distributions carefully and find
that they follow power-law models. Besides, the distributions of fluctuations
in both parameters follow the Tsallis $q$-Gaussian distributions and the $q$
values are consistent for different scale intervals. Therefore, these results
indicate that both parameters are scale-invariant, which provides new
parameters for the study of self-organized criticality systems. Interestingly,
we find that the $q$ values for MVT and spectral lag are similar with duration
and fluence, respectively.

###### keywords:

stars: magnetars

††pubyear: 2022††pagerange: The Self-organized Criticality Behaviors of Two
New Parameters in SGR J1935+2154–LABEL:lastpage

##  1 Introduction

The phenomenon of self-similarities widely exists in temporal and spatial
scales in astronomical systems, including soft gamma repeaters (SGRs) (Cheng
et al., 1996; Göğüş et al., 1999; Chang et al., 2017; Wei et al., 2021; Sang &
Lin, 2022), X-ray flares of afterglows of gamma-ray burst (GRB) (Wang & Dai,
2013), fast radio bursts (FRBs) (Wang & Yu, 2017; Wei et al., 2021; Sang &
Lin, 2022), solar flares (Goodman et al., 2020), pulsar glitches (Melatos et
al., 2008), etc (see Aschwanden (2014) and Aschwanden et al. (2016) for a
review). Interestingly, by discovering similar distributions of released
energies for earthquakes and SGRs (Göğüş et al., 1999), it was suggested that
SGRs originate from starquakes in magnetars (Duncan & Thompson, 1992; Thompson
& Duncan, 1995).

The Self-organized criticality (SOC) in slowly driven nonlinear dissipative
systems was proposed by Bak et al. (1987), Aschwanden (2012) and Aschwanden et
al. (2016) to explain the phenomenon of self-similarities (i.e. power-law size
distributions), the spatial dimension and the classical diffusion of the
fractal-diffusive SOC system can be estimated through the power-law index for
the parameters’ length scales, durations, peak fluxes and fluences. On the
other hand, the distribution for the parameter fluctuations at different times
in a SOC system was found to follow a Tsallis $q$-Gaussian function, and the
$q$ values are constant for different scale intervals (i.e. scale invariance
of the avalanche size differences) (Caruso et al., 2007). The power-law index
of the parameter distribution and the $q$ value of the parameter difference
distribution can be inter-convertible (Caruso et al., 2007; Celikoglu et al.,
2010).

SGR J1935+2154 is one of the most active, and the only magnetar from which an
X-ray burst (XRB) and a FRB (a bright ms-long radio burst) have been observed
in association (Li et al., 2021; Bochenek et al., 2020; Mereghetti et al.,
2020; Younes et al., 2021; Ridnaia et al., 2021; Tavani et al., 2021)).
Recently Wei et al. (2021) and Sang & Lin (2022) reported that the XRBs
originated from a SOC system through analyzing the parameters the duration,
waiting time, fluence and flux, that is, these parameters are scale-invariant.
Similarly, the scale invariance of these parameters is also found in other
magnetars (e.g. SGR J1550–5418, Wang & Yu (2017); SGR 1806–20, Cheng et al.
(1996); 1900+14, Göğüş et al. (1999)) and astronomical sources.

Figure 1: The differential (left panel) and cumulative (right panel)
distributions of MVTs observed by GECAM, HXMT and Fermi/GBM. The fitted yellow
lines are obtained by MCMC with reduced-$\chi^{2}$ 1.9 and 0.6, respectively.

Figure 2: Left panel: the distribution of fluctuations (MVT) for $n$=1, the
yellow lines are fitted by a Tsallis $q$-Gaussian distribution ($q=2.36\pm
0.12$), and reduced-$\chi^{2}$ is 0.8. Right panel: the fitted q values for
different $n$, the yellow lines are fitted by a constant, which is $q=2.19\pm
0.01$ with reduced-$\chi^{2}$ is 0.9.

Figure 3: The differential (left panel) and cumulative (right panel)
distributions of spectral lags observed by GECAM, HXMT and Fermi/GBM. The
fitted yellow lines are obtained by MCMC, and reduced-$\chi^{2}$ are 2.3 and
0.9, respectively.

Figure 4: Left panel: the distribution of fluctuations (spectral lag) for
$n$=1, the yellow lines are fitted by a Tsallis $q$-Gaussian distribution
($q=2.94\pm 0.20$), and reduced-$\chi^{2}$ is 1.5. Right panel: the fitted q
values for different $n$, the yellow lines are fitted by a constant, which is
$q=2.86\pm 0.02$, and reduced-$\chi^{2}$ is 1.1.

In our recent work, we reported the small spectral lags (the distribution
peaks $\sim$ 1.3 ms) (Xiao et al., 2023b) and the minimum variation timescale
(MVT, the distribution peaks at $\sim$ 2 ms) (Xiao et al., 2023a) by analyzing
hundreds of XRBs from SGR J1935+2154 observed by GECAM, HXMT and GBM. The two
parameters are important probes to studying the physical mechanism and
radiation region of XRBs (Schmidt, 1978; Fenimore et al., 1993; Titarchuk et
al., 2007; Zhang et al., 2009; Huppenkothen et al., 2013; Ackermann et al.,
2014; Golkhou & Butler, 2014; Golkhou et al., 2015). However, to our
knowledge, whether these parameters are also scale-invariant has not been
studied, which is the focus of this work.

In this work, we present our sample selection and scale invariance analysis in
Section 2 and a discussion and summary are given in Section 3.

##  2 Scale invariance analysis

###  2.1 Sample selection and method

The bursts from SGR J1935+2154 and the values of MVT and spectral lag observed
by GECAM, HXMT and GBM are collected from our previous work (Xiao et al.,
2023a, b). Since all three satellites have large fields of view and high
sensitivity (Xiong, 2020; Zhang et al., 2020; Meegan et al., 2009), we have a
relatively comprehensive sample. For bursts observed jointly by multiple
satellites, we adopt the GBM measurements by default. In total, 669 bursts are
selected from July 2014 to January 2022 after removing those with “timing
glitches” (a known GBM hardware anomaly with dips and peaks in a light curve
(Briggs et al., 2013).) or saturation.

To investigate whether the parameters are scale-invariant or not, we adopt two
approaches like Refs. (Chang et al., 2017; Wei et al., 2021; Sang & Lin, 2022;
Li et al., 2023). The first is to investigate whether the differential and
cumulative distributions of the parameters follow power-law models
(Aschwanden, 2015) (see Equations 1 and 2) and the power-law index are
consistent within the uncertainties,

$\begin{split}n_{0}=N(1-\alpha_{\rm d})[(x_{2}+x_{0})^{(1-\alpha_{\rm
d})}-(x_{1}+x_{0})^{(1-\alpha_{\rm d})}]^{-1},\\\ N_{\rm diff}=n_{0}(x_{\rm
0d}+x)^{-\alpha_{\rm d}},\end{split}$ (1) $N_{\rm
cum}(>x)=1+(N-1)\times\frac{(x_{2}+x_{\rm 0c})^{1-\alpha_{\rm c}}-(x+x_{\rm
0c})^{1-\alpha_{c}}}{(x_{2}+x_{\rm 0c})^{1-\alpha_{c}}-(x_{1}+x_{\rm
0c})^{1-\alpha_{\rm c}}},$ (2)

where $n_{0}$, $\alpha_{d}$ and $\alpha_{c}$ are normalization constant, the
power-law index of differential and cumulative distributions, respectively.
$N$ is the number of bursts, $x_{0d}$ and $x_{0c}$ are constants due to the
threshold effects.

The second is to test whether the distributions of the differences of the
parameters $X_{n}$ follow the Tsallis $q$-Gaussian distribution (Tsallis,
1988; Tsallis et al., 1998) (Equation 4) and the $q$ values for different time
scales $n$ are consistent for different scale intervals,

$X_{n}=S_{i+n}-S_{i},$ (3)

where $S_{i}$ is the value of the parameter (i.e. MVT and spectral lag) of the
$i$th burst ordered according to time and $n$ is the temporal interval scale.
We obtain $x_{n}$ by dividing the difference $X_{n}$ by a scale factor
($\sigma_{X_{n}}$), which is the standard deviation of $X_{n}$.

$f(x_{n})=\alpha\left[1-\left(1-q\right)x_{n}^{2}/\beta\right]^{\frac{1}{1-q}},$
(4)

where $\alpha$ is the normalization factor, the $\beta$ and $q$ (the parameter
we are interested in) affect the width and sharpness of the peak,
respectively. We adopted the Markov Chain Monte Carlo (MCMC) method to fit the
data, and tests the goodness-of-fit by reduced-$\chi^{2}$ (i.e. close to 1
implies a good fit).

###  2.2 Minimum variation timescales

The MVT in 8-100 keV of each burst is calculated in Xiao et al. (2023a). The
left panel of Figure 1 shows the differential distribution, which can be
fitted fairly well (reduced-$\chi^{2}$=$\chi^{2}/dof$=37.4/19=1.9) by the
power-law model (Equation 1) with $\alpha_{\rm d}=2.56\pm 0.40$ and $x_{\rm
0d}=0.010\pm 0.003$. Note that less than $\sim$0.001 ms does not fit well due
to threshold effects (Aschwanden, 2015). The cumulative distribution of MVT is
shown in the right panel of Figure 1, which also can be fitted well
(reduced-$\chi^{2}$=$\chi^{2}/dof$=12.6/23=0.6) by the power-law model
(Equation 2) with $\alpha_{\rm c}=2.89\pm 0.35$ and $x_{\rm 0c}=0.014\pm
0.003$. Note that $\alpha_{\rm d}$ and $\alpha_{\rm c}$ are consistent within
the uncertainties, which indicates that the MVT parameter of XRBs from SGR
J1935+2154 exhibits SOC behavior.

The distributions of fluctuations for different time scales $n$ are
investigated, for example, the left panel of Figure 2 shows the distribution
for $n$=1, and it can be fitted well
(reduced-$\chi^{2}$=$\chi^{2}/dof$=58.1/76=0.8) by the Tsallis $q$-Gaussian
distribution (Equation 4) with $q=2.36\pm 0.12$. The q-values obtained by
fitting the distributions for different scales $n$ are shown in the right
panel of Figure 2, which are steady and can be fitted well with a constant
($q=2.19\pm 0.01$, with reduced-$\chi^{2}$=$\chi^{2}/dof$=83.3/98=0.9), which
also verifies that the MVT parameter of XRB from SGR J1935+2154 exhibits SOC
behavior.

###  2.3 Spectral lags

The spectral lag between light curves in the 10-20 and 60-100 keV energy bands
is calculated in Xiao et al. (2023b), and a positive lag is defined as that
low-energy photons follow high-energy photons. The differential and cumulative
distributions are shown in Figure 3, both follow power-law models well. The
result of the fit to the difference distribution is $\alpha_{\rm d}=1.54\pm
0.16$ and $x_{\rm 0d}=0.003\pm 0.001$ with
reduced-$\chi^{2}$=$\chi^{2}/dof$=43.9/19=2.3. For the cumulative
distribution, the result of the fit is $\alpha_{\rm c}=1.80\pm 0.11$ and
$x_{\rm 0c}=0.003\pm 0.001$, with
reduced-$\chi^{2}$=$\chi^{2}/dof$=11.5/13=0.9. Note that the $\alpha_{\rm d}$
and $\alpha_{\rm c}$ are also consistent within the uncertainties, which
indicates that the spectral lag parameter of XRB from SGR J1935+2154 exhibits
SOC behavior.

The left panel of Figure 4 shows the distribution of fluctuations for time
scale $n=1$, which can be fitted well
(reduced-$\chi^{2}$=$\chi^{2}/dof$=52.3/36=1.5) by the Tsallis $q$-Gaussian
distribution (Equation 4) with $q=2.94\pm 0.20$. The q-values obtained by
fitting the distributions for different scales $n$ are also consistent within
the uncertainties and fitted with a constant ($q=2.86\pm 0.02$,
reduced-$\chi^{2}$=$\chi^{2}/dof$=103.7/98=1.1), which also verifies that the
spectral lag parameter of XRB from SGR J1935+2154 exhibit SOC behavior.

##  3 Discussion and conclusion

In this work, we perform a detailed study of the statistical properties of the
two parameters, i.e., MVT and spectral lag, of hundreds of XRBs from SGR
J1935+2154 observed by GECAM, HXMT and GBM. With two methods, we report that
both parameters exhibit SOC behaviors, that is, the differential and
cumulative distributions of the parameters well follow power-law models and
the power-law index ($\alpha_{\rm d}$ and $\alpha_{\rm c}$) are consistent
within the uncertainties, as well as the distributions of fluctuations follow
the Tsallis $q$-Gaussian distribution and the $q$ values for different time
scales are consistent for different scale intervals.

Interestingly, the $q$ values (i.e. $q=2.19\pm 0.01$) for MVT are similar to
the result (i.e. $2.28\pm 0.15$) for the duration reported by Wei et al.
(2021). The duration of a burst only describes the total emission properties,
but MVT captures the information concerning individual pulses, that is, the
MVT is approximately equal to the rise time of the shortest pulse in a burst.
In a previous work, we did not find a significant correlation between duration
and MVT for XRBs from SGR J1935+2154 (Xiao et al., 2023a). The MVT of XRBs can
be used to estimate the emission region in pulsar-like models, as well as the
radius and Lorentz factor of the relativistic jet in GRB-like models (see
Zhang, 2020, 2022 for review). Therefore, the SOC behavior of MVT implies that
the radius of the radiation region may have a similar origin in pulsar-like
models.

On the other hand, although the spectral lag is a time-domain parameter, the
$q$ values (i.e. $q=2.86\pm 0.02$) for spectral lag are similar to the result
(i.e. $2.78\pm 0.12$) for fluence reported by Wei et al. (2021). The
relationship between fluence and spectral lag has been found in some transient
sources such as GRBs (Ukwatta et al., 2012), but there is a lack of studies
for XRBs. Further work is required to investigate this.

According to the study of Caruso et al. (2007) and Celikoglu et al. (2010),
there is an exact relation (i.e. $q=(\alpha+2)/2$) between the power-law index
($\alpha_{\rm d}$ or $\alpha_{\rm c}$) and the $q$ values based on the
assumption that the sizes of two events have no correlation. We find that this
relationship is also satisfied for MVT and spectral lag within the
uncertainties. However, MVT and spectral lag are different from other
parameters, such as duration, peak flux, and fluence or energy, which have
detailed theoretical explanations, e.g. the spatial dimension can be estimated
based on the Fractal-diffusive Avalanche Model (Aschwanden, 2012). Therefore,
we believe that the SOC behaviours of MVT and spectral lag will provide new
insights in future SOC studies of these phenomen.

## Acknowledgments

We acknowledge the public data from Fermi/GBM. This work is supported by the
National Natural Science Foundation of China (NOs. 12303043, 12273042 and
12273008) and the National Key R&D Program of China 2022YFF0711404. The
authors also thank supports from the Strategic Priority Research Program on
Space Science, the Chinese Academy of Sciences (Grant No. XDA15360300,
XDA15052700), the National SKA Program of China (Nos. 2022SKA0130100,
2022SKA0130104), the Natural Science and Technology Foundation of Guizhou
Province (No. [2023]024) and the Foundation of Guizhou Provincial Education
Department (No. KY (2020) 003).

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Ackermann et al. (2014) Ackermann M., et al., 2014, Science, 343, 42
  * Aschwanden (2012) Aschwanden M. J., 2012, Astronomy & Astrophysics, 539, A2
  * Aschwanden (2014) Aschwanden M. J., 2014, ApJ, 782, 54
  * Aschwanden (2015) Aschwanden M. J., 2015, The Astrophysical Journal, 814, 19
  * Aschwanden et al. (2016) Aschwanden M. J., et al., 2016, Space Science Reviews, 198, 47
  * Bak et al. (1987) Bak P., Tang C., Wiesenfeld K., 1987, Physical review letters, 59, 381
  * Bochenek et al. (2020) Bochenek C. D., Ravi V., Belov K. V., Hallinan G., Kocz J., Kulkarni S. R., McKenna D. L., 2020, Nature, 587, 59
  * Briggs et al. (2013) Briggs M. S., et al., 2013, Journal of Geophysical Research: Space Physics, 118, 3805
  * Caruso et al. (2007) Caruso F., Pluchino A., Latora V., Vinciguerra S., Rapisarda A., 2007, Physical Review E, 75, 055101
  * Celikoglu et al. (2010) Celikoglu A., Tirnakli U., Queirós S. M. D., 2010, Physical Review E, 82, 021124
  * Chang et al. (2017) Chang Z., Lin H.-N., Sang Y., Wang P., 2017, Chinese Physics C, 41, 065104
  * Cheng et al. (1996) Cheng B., Epstein R. I., Guyer R. A., Young A. C., 1996, Nature, 382, 518
  * Duncan & Thompson (1992) Duncan R. C., Thompson C., 1992, Astrophysical Journal, Part 2-Letters (ISSN 0004-637X), vol. 392, no. 1, June 10, 1992, p. L9-L13. Research supported by NSERC., 392, L9
  * Fenimore et al. (1993) Fenimore E., Epstein R., Ho C., 1993, Astronomy and Astrophysics Supplement Series, 97, 59
  * Göğüş et al. (1999) Göğüş E., Woods P. M., Kouveliotou C., Van Paradijs J., Briggs M. S., Duncan R. C., Thompson C., 1999, The Astrophysical Journal, 526, L93
  * Golkhou & Butler (2014) Golkhou V. Z., Butler N. R., 2014, The Astrophysical Journal, 787, 90
  * Golkhou et al. (2015) Golkhou V. Z., Butler N. R., Littlejohns O. M., 2015, The Astrophysical Journal, 811, 93
  * Goodman et al. (2020) Goodman M. L., Kwan C., Ayhan B., Shang E. L., 2020, Frontiers of Physics, 15, 1
  * Huppenkothen et al. (2013) Huppenkothen D., et al., 2013, The Astrophysical Journal, 768, 87
  * Li et al. (2021) Li C., et al., 2021, Nature Astronomy, 5, 378
  * Li et al. (2023) Li X.-J., Zhang W.-L., Yi S.-X., Yang Y.-P., Li J.-L., 2023, The Astrophysical Journal Supplement Series, 265, 56
  * Meegan et al. (2009) Meegan C., et al., 2009, ApJ, 702, 791
  * Melatos et al. (2008) Melatos A., Peralta C., Wyithe J., 2008, The Astrophysical Journal, 672, 1103
  * Mereghetti et al. (2020) Mereghetti S., et al., 2020, The Astrophysical Journal Letters, 898, L29
  * Ridnaia et al. (2021) Ridnaia A., et al., 2021, Nature Astronomy, 5, 372
  * Sang & Lin (2022) Sang Y., Lin H.-N., 2022, Monthly Notices of the Royal Astronomical Society, 510, 1801
  * Schmidt (1978) Schmidt W. K. H., 1978, Nature, 271, 525
  * Tavani et al. (2021) Tavani M., et al., 2021, Nature Astronomy, 5, 401
  * Thompson & Duncan (1995) Thompson C., Duncan R. C., 1995, Monthly Notices of the Royal Astronomical Society, 275, 255
  * Titarchuk et al. (2007) Titarchuk L., Shaposhnikov N., Arefiev V., 2007, ApJ, 660, 556
  * Tsallis (1988) Tsallis C., 1988, Journal of Statistical Physics, 52, 479
  * Tsallis et al. (1998) Tsallis C., Mendes R., Plastino A. R., 1998, Physica A Statistical Mechanics and its Applications, 261, 534
  * Ukwatta et al. (2012) Ukwatta T., et al., 2012, Monthly Notices of the Royal Astronomical Society, 419, 614
  * Wang & Dai (2013) Wang F., Dai Z., 2013, Nature Physics, 9, 465
  * Wang & Yu (2017) Wang F., Yu H., 2017, Journal of Cosmology and Astroparticle Physics, 2017, 023
  * Wei et al. (2021) Wei J.-J., Wu X.-F., Dai Z.-G., Wang F.-Y., Wang P., Li D., Zhang B., 2021, The Astrophysical Journal, 920, 153
  * Xiao et al. (2023a) Xiao S., et al., 2023a, The Astrophysical Journal Supplement Series, 268, 5
  * Xiao et al. (2023b) Xiao S., et al., 2023b, MNRAS, 521, 5308
  * Xiong (2020) Xiong S., 2020, SCIENTIA SINICA Physica, Mechanica & Astronomica, 50, 129501
  * Younes et al. (2021) Younes G., et al., 2021, Nature Astronomy, 5, 408
  * Zhang (2020) Zhang B., 2020, Nature, 587, 45
  * Zhang (2022) Zhang B., 2022, arXiv preprint arXiv:2212.03972
  * Zhang et al. (2009) Zhang B., et al., 2009, The Astrophysical Journal, 703, 1696
  * Zhang et al. (2020) Zhang S.-N., et al., 2020, SCIENCE CHINA Physics, Mechanics & Astronomy, 63, 1


# Automatic cross-validation in structured models:
Is it time to leave out leave-one-out?

Adin, A.1,2, Krainski, E.3, Lenzi, A.4, Liu, Z.5, Martínez-Minaya, J.6, and
Rue, H.3
1 Department of Statistics, Computer Sciences and Mathematics, Public
University of Navarre, Spain.
2 Institute for Advanced Materials and Mathematics, InaMat2, Public University
of Navarre, Spain.
3 Statistics Program, Computer, Electrical and Mathematical Sciences and
Engineering Division,
King Abdullah University of Science and Technology (KAUST).
4 School of Mathematics, University of Edinburgh, Scotland.
5 RIKEN Center for AI Project, Tokyo, Japan.
6 Department of Applied Statistics, Operations Research and Quality,
Universitat Politècnica de València, Spain.

###### Abstract

Standard techniques such as leave-one-out cross-validation (LOOCV) might not
be suitable for evaluating the predictive performance of models incorporating
structured random effects. In such cases, the correlation between the training
and test sets could have a notable impact on the model’s prediction error. To
overcome this issue, an automatic group construction procedure for leave-
group-out cross validation (LGOCV) has recently emerged as a valuable tool for
enhancing predictive performance measurement in structured models. The purpose
of this paper is (i) to compare LOOCV and LGOCV within structured models,
emphasizing model selection and predictive performance, and (ii) to provide
real data applications in spatial statistics using complex structured models
fitted with INLA, showcasing the utility of the automatic LGOCV method. First,
we briefly review the key aspects of the recently proposed LGOCV method for
automatic group construction in latent Gaussian models. We also demonstrate
the effectiveness of this method for selecting the model with the highest
predictive performance by simulating extrapolation tasks in both temporal and
spatial data analyses. Finally, we provide insights into the effectiveness of
the LGOCV method in modelling complex structured data, encompassing spatio-
temporal multivariate count data, spatial compositional data, and spatio-
temporal geospatial data.

Keywords: Cross-validation, Hierarchical models, INLA, Spatial statistics

##  1 Introduction

Cross-validation is a general approach used to evaluate the predictive
performance of a statistical model in the absence of new data, and it is
commonly used to compute score rules for model comparison or selection. The
fundamental concept underlying cross-validation is to split the observed data
into a training set (the sample data used for estimating the model’s
parameters) and a testing set (the set of data points used to compute
prediction errors based on the trained model) multiple times and use the
combined prediction error to estimate the predictive accuracy of a model
(Gelman et al.,, 1995; Hastie et al.,, 2009).

The $k$-fold cross-validation (KCV) is a widely used variant of cross-
validation, where data is first randomly partitioned into $k$ equally (or
nearly equally) sized mutually exclusive groups or folds, so that one fold is
used to test the performance of the model fitted on the remaining $k-1$ folds.
The procedure is repeated $k$ times, with each repetition employing a distinct
subset of the observed data as the validation set. However, a poorly chosen
value for $k$ may lead to a misrepresentative idea of the predictive
performance of the model. For instance, it could result in a score with a high
variance if the data used to fit the model give rise to highly variable test
error measures. Alternatively, it could lead to a high bias, such as an
overestimate of the performance of the model. Although no formal rule exists,
a value of $k$ equal to 5 or 10 is widely prevalent in the fields of
statistics and applied machine learning as a trade-off between computational
burden and bias reduction of the predictive error (Hastie et al.,, 2009; Kuhn
et al.,, 2013).

Leave-one-out cross-validation (LOOCV) represents limiting case of KCV,
wherein individual observations are systematically excluded, one at a time, to
test the model’s performance trained on the remaining data. That is, the data
is divided into $k=n$ folds (where $n$ is the size of the observed data) so
that each data point is sequentially used as a testing set. The main
advantages of LOOCV are that (i) the size of the training set closely
approximates the full data set and, (ii) it provides and approximately
unbiased estimate for the predictive error. While the exact computation of
LOOCV is a “simple” and straightforward strategy, as we only need to fit
models on all possible training sets and compute the utility function on the
corresponding testing set, this method can become practically unfeasible due
to the prohibitive computational cost of fitting the model as many times as
the number of observations in the data. For this reason, several methods have
been proposed to efficiently compute fast approximations of LOOCV. See, e.g.
Bürkner et al., (2021); Held et al., (2010); Liu and Rue, (2022) and
references therein.

Standard cross-validation techniques like KCV and LOOCV operate under the
assumption that data is sampled independently from a joint distribution, so
that it can be freely splitted into training sets that represent the observed
data and testing sets that represent the unobserved data (Arlot and Celisse,,
2010). However, it is well-known that these techniques might not be suitable
to measure the predictive performance of a model for correlated data, where
commonly hierarchical models with structured random effects are considered to
capture the temporal, spatial and/or hierarchical clustering dependence
structures within the data (Roberts et al.,, 2017; Rabinowicz and Rosset,,
2022). In such situations, the correlation existing between the training and
testing sets can have a notable impact on the model’s prediction error, and
consequently affect the decision-making process for model selection, since
standard cross-validation techniques are generally over-optimistic about the
predictive performance of structured models.

A common approach for enhancing independence in cross-validation techniques is
to adapt LOOCV to structured models by adjusting the training sets for each
testing set, a process frequently determined by the specific demands of
prediction task. For example, in the context of time series data (or other
intrinsically ordered data) the prediction task could encompass various
objectives such as interpolating a missing value in the observed data,
predicting the response in an independent replicate of the observed time
series, or forecasting the response in future time points. For the first case,
using LOOCV could be a suitable strategy because assessing the model’s
predictive performance by removing a single data point replicates the
interpolation task effectively. However, a different cross-validation scheme
should be used to evaluate the model’s ability when a time series model is
used for forecasting purposes (extrapolation task). A detailed comparison of
different cross-validation approaches for time series forecasting is described
in Bergmeir and Benítez, (2012).

Therefore, when handling structured models, it is essential to adapt the
cross-validation design to ensure the calculation of dependable scoring rules
for model comparison or selection. In practice, many real data analyses demand
the use of complex models with multiple correlation structures such as random
effects at different aggregation levels, making it very difficult to determine
the optimal design and adapting it to each of the testing points. These
problems are exacerbated when dealing with complex structures on large data
sets. To address both issues concurrently, Liu and Rue, (2022) introduced an
automatic group construction procedure for leave-group-out cross-validation
(LGOCV) to estimate the predictive performance of structured models by
deriving an efficient and accurate approximation of LGOCV for latent Gaussian
models within the INLA framework (Rue et al.,, 2009; Van Niekerk et al.,,
2023).

The aim of the current paper is two-fold: (i) to compare the behavior of
predictive measures based on leave-one-out and leave-group-out cross-
validation techniques within temporally or spatially structured models, and
(ii) to provide real data applications of complex structured models fitted
with INLA in the context of spatial statistics, while demonstrating the use of
the automatic LGOCV method. In addition, we conduct a simulation study to
evaluate the performance of the automatic group construction method for model
comparison (see Supplementary Material).

The rest of the paper is structured as follows. Section 2 briefly describes
the most relevant aspects of the LGOCV method for latent Gaussian models. In
Section 3, we compare the extrapolation performance of several structured
models using both temporal and spatial data. Our objective is to emphasize
disparities between model selection driven by conventional criteria for
Bayesian model comparison and predictive performance measures based on LOOCV
and the recently proposed automatic LGOCV method. We also validate the models’
predictive performance by imitating real extrapolation tasks. Section 4
presents three different real data applications in the field of spatial
statistics. Specifically, a joint modelling of pancreatic cancer mortality and
incidence data in England using spatio-temporal multivariate disease mapping
models (Section 4.1), the estimation of the membership proportions to the
different genetic clusters of Arabidopsis thaliana on the Iberian Peninsula
using a Logistic Normal Dirichlet Model for spatial compositional data
(Section 4.2), and modelling of wind speed data in the United Kingdom using
spatio-temporal models for continuous domains with non-separable covariance
structures (Section 4.3). Finally, Section 5 ends with some conclusions. The
data and R code to reproduce the examples, simulation studies, and real data
analysis presented in this paper are available at
https://github.com/spatialstatisticsupna/INLA_groupCV.

##  2 Leave-group-out cross-validation for latent Gaussian models

LGOCV takes each data point, represented by $y_{i}$, as a testing set, and
creates an index set $I_{i}$ specific to each point based on the prediction
task. Using all observed data except the data indexed by $I_{i}$, denoted by
$\boldsymbol{y}_{-I_{i}}$, it computes the posterior predictive density
$\pi(Y_{i}|\boldsymbol{y}_{-I_{i}})$ and validates the prediction using the
observed data point $y_{i}$ and scoring rules. However, when the model
structure is complex or the prediction task is unclear, constructing $I_{i}$
can be challenging. To overcome this, we can automatically construct $I_{i}$
by selecting the most informative data points to predict $y_{i}$. This
automatic construction is motivated by the fact that most extrapolation tasks
have less information available from the data used for prediction. While
constructing $I_{i}$ automatically can be challenging in general model
settings, it becomes straightforward when the model is represented as a latent
Gaussian model. Furthermore, using a latent Gaussian model allows for fast and
accurate approximation of posterior densities when excluding partial data.

A latent Gaussian model is composed of three layers that can be formulated by

$\begin{split}&y_{i}|\eta_{i},\boldsymbol{\theta}\sim\pi(y_{i}|\eta_{i},\boldsymbol{\theta}),\\\
&\boldsymbol{\eta}=\boldsymbol{A}\boldsymbol{f},\\\
&\boldsymbol{f}|\boldsymbol{\theta}\sim
N(0,\boldsymbol{P}_{\boldsymbol{f}}(\boldsymbol{\theta})),\\\
&\boldsymbol{\theta}\sim\pi(\boldsymbol{\theta}).\end{split}$

The observed data $\boldsymbol{y}$ is modelled conditionally independent given
the hyperparameter $\boldsymbol{\theta}$ and its corresponding linear
predictor $\eta_{i}$ using the likelihood function
$\pi(y_{i}|\eta_{i},\boldsymbol{\theta})$. The linear predictor
$\boldsymbol{\eta}$ is a linear combination of model effects $\boldsymbol{f}$,
which has a multivariate Gaussian prior with a zero mean and a precision
matrix $\boldsymbol{P}_{\boldsymbol{f}}(\boldsymbol{\theta})$ that depends on
$\boldsymbol{\theta}$. We assign a prior distribution to
$\boldsymbol{\theta}$. Using INLA method, we can have a Gaussian approximation
of the posterior density
$\pi(\boldsymbol{f}|\boldsymbol{\theta},\boldsymbol{y}_{-I_{i}})$. The
automatic groups are constructed using the correlation matrix of
$\boldsymbol{f}$, which can be derived from either the prior density or the
approximated posterior density at the mode of hyperparameter configuration.
See Liu and Rue, (2022) for further details.

We illustrate how to construct $I_{i}$ using a simple example. Suppose we have
correlation coefficients $\boldsymbol{C}_{i}$, which contain the correlation
of the linear predictor $\eta_{i}$ with all $\eta_{j}$, for $j=1,\ldots,n$. If
$\boldsymbol{C}_{i}=\\{1,1,0.9,0.9,0.8,0.8,-0.1,-0.1,0,0\\}$, we define
$\boldsymbol{L}_{i}$ as the sorted unique absolute values of
$\boldsymbol{C}_{i}$, giving us
$\boldsymbol{L}_{i}=\\{1.0,0.9,0.8,0.1,0.0\\}$. We then include the most
correlated data points based on the $m$ largest correlation levels in $I_{i}$.
For example, if we choose $m=3$, then $I_{i}=\\{1,2,3,4,5,6\\}$. We refer to
the parameter $m$ as the number of level sets, which is determined by the
degree of extrapolation implied by LGOCV. When the information about the
degree of extrapolation is missing, an empirical suggestion is to use $m=3$,
or we could compute a sequence of $m$ values to determine the predictive
performance.

##  3 Evaluating cross-validation measures in structured models

In the forthcoming examples, our focus is on analyzing the extrapolation
performance of structured models, in which the unobserved data is often
assumed less dependent on the observed data. For this end, we compare the
behaviour of different models in terms of (i) Bayesian model selection
criteria such as the deviance information criterion (DIC; Spiegelhalter et
al.,, 2002) and the Watanabe-Akaike information criterion (WAIC; Watanabe,,
2010), and (ii) predictive performance measures based on LOOCV and the
recently proposed automatic group construction procedure for LGOCV for latent
Gaussian models using the well-known INLA estimation technique. Specifically,
we compute commonly used utility functions such as the logarithmic score
(Gneiting and Raftery,, 2007)

$\mbox{LS}=-\frac{1}{n}\sum\limits_{i=1}^{n}\log\pi(Y_{i}=y_{i}|\boldsymbol{y}_{-I_{i}}),$

and the mean square prediction error

$\mbox{MSPE}=\frac{1}{n}\sum\limits_{i=1}^{n}\left(E[Y_{i}|\boldsymbol{y}_{-I_{i}}]-y_{i}\right)^{2},$

where $\pi(Y_{i}=y_{i}|\boldsymbol{y}_{-I_{i}})$ denotes the cross-validated
predictive density at the observed data $y_{i}$. Note that in the case of
LOOCV the training set $\boldsymbol{y}_{-I_{i}}$ reduces to ${\bf
y}_{-i}=(y_{1},\ldots,y_{i-1},y_{i+1},y_{n})^{{}^{\prime}}$. See Liu and Rue,
(2022) for details about the computation of
$E[Y_{i}|\boldsymbol{y}_{-I_{i}}]$.

In our first example, we compare the predictive performance exhibited by
different temporal models for long-term forecasting in time series data. Our
second example delves into the extrapolation performance of spatial disease
mapping models for areal count data, that is, predictions in non-observed
spatial units. We show that DIC, WAIC and leave-one-out cross-validation
(LOOCV) measures points to different models from those selected by using
utility (loss) functions based on the recently proposed automatic LGOCV. We
also validate the predictive performance of these models by imitating real
extrapolation tasks for both examples based on temporal and spatial data. All
computations are made using the recently implemented “compact” mode (Van
Niekerk et al.,, 2023) in R-INLA (stable) version INLA_23.09.09 of R-4.3.1.

###  3.1 Example 1: temporal models

In this example, we want to reproduce a real situation where the aim is to
perform long-term forecasting through the use of temporal models. Let us
consider a scenario where our interest lies in forecasting surface
temperatures using shortwave radiation as a predictor variable. We start by
simulating radiation data $X$ for a total of $n=2000$ time points using a
first-order autoregressive (AR1) model, where the autocorrelation coefficient
is set at $\rho=0.9$ and the marginal variance equals to one. Then, we
generate values of surface temperature $Y$ using the following linear model:

$Y=\beta_{0}+X+e,\quad\mbox{with}\quad e\sim N(0,0.1).$

Finally, let us assume that $\tilde{X}=X\cdot|X|$ corresponds to an imperfect
observation of the shortwave radiation $X$ (arising from potential
instrumental deficiencies, for instance).

Two different models are considered to estimate the surface temperature
$y_{t}$ based on observed radiation data $\tilde{x}_{t}$ for $t=1,\ldots,2000$
time points.

  * Linear model:

$\begin{split}&y_{t}|\eta_{t},\boldsymbol{\theta}\sim N(\eta_{t},\tau_{y}),\\\
&\eta_{t}=\beta_{0}+\beta\cdot\tilde{x}_{t},\end{split}$

where $\beta_{0}$ is an intercept and $\beta$ is the fixed effect coefficient
associated to the covariate $\tilde{x}_{t}$.

  * AR1 model:

$\begin{split}&y_{t}|\eta_{t},\boldsymbol{\theta}\sim N(\eta_{t},\tau_{y}),\\\
&\eta_{t}=\beta_{0}+u_{t},\\\ &u_{t}=\rho\cdot u_{t-1}+\epsilon_{t},\\\
&\epsilon_{t}\sim N(0,\tau_{u})\end{split}$

where the temporally structured random effect ${\bf
u}=(u_{1},\ldots,u_{2000})^{{}^{\prime}}$ with AR1 prior distribution is used
to model the surface temperature as a substitute of the observed predictor
$\tilde{x}_{t}$.

Figure 1 shows how the models capture the underlying dependence structure
between the observed radiation and temperature values. Clearly, the AR1 model
better captures the non-linear relationship between these variables. Posterior
median estimates for the temperature values ${\bf
y}_{t}=(y_{1},\ldots,y_{2000})^{{}^{\prime}}$ obtained with both models are
plotted in Figure 2.

Figure 1: Dependence structure captured between the observed shortwave
radiation $\tilde{X}$ and estimated surface temperature $\hat{Y}$. Figure 2:
True values (black dots) and posterior median estimates (colored lines) for
the surface temperature ${\bf y}_{t}=(y_{1},\ldots,y_{2000})^{{}^{\prime}}$.

Table 1 displays the values of model selection criteria and predictive
performance measures (logarithmic score and mean square prediction error) for
LOOCV and LGOCV with automatic groups construction (m=3), derived from either
the prior precision matrix or the posterior correlation matrix of each model.
See Liu and Rue, (2022) for further details. The groups derived from the AR1
model (as this is the true generating model) are used as a reference to make
predictive measures comparable. As expected, the AR1 model is pointed out as
the best model in terms of model fitting and predictive performance based on
the LOOCV approach. However, predictive measures based on LGOCV points to the
linear model. The reason for the difference in model choice is because LOOCV
technique in temporally structured models measures predictive performance by
interpolating missing values within the observed time series. In contrast,
LGOCV assesses the model’s capability when the prediction task is to forecast
the response variable for future time points.

Table 1: Model comparison for temporal data: model selection criteria
(posterior mean deviance $\bar{D}$, effective number of parameters $p_{D}$,
DIC and WAIC) and predictive performance measures for LOOCV and LGOCV with
automatic groups construction (m=3) based on prior and posterior correlations
for each temporal model.

| Model selection criteria | LOOCV | LGOCV
---|---|---|---
| "prior" | "posterior"
Model | $\bar{D}$ | $p_{D}$ | DIC | WAIC | LS | MSPE | LS | MSPE | LS | MSPE
Linear | 1500.0 | 3.0 | 1503.0 | 1507.7 | 0.376 | 0.125 | 0.383 | 0.126 | 0.383 | 0.126
AR1 | -439.5 | 1286.7 | 847.2 | 622.3 | 0.349 | 0.116 | 0.824 | 0.303 | 0.823 | 0.302

To validate these measures, we conduct the following simulation study. We
start by considering our training set as the time series encompassing the
period $t=\\{1,\ldots,1900\\}$ and compute long-term predictions for the next
100 time points using INLA. We repeat this process by shifting the time window
one step backward for a total of 500 data sets. That is, our last historic
data corresponds to the period $t=\\{1,\ldots,1400\\}$, so that predictions
are computed for future time points $\\{1401,\ldots,1500\\}$. We compute the
mean absolute prediction errors (MAPE) and root mean square prediction errors
(RMSPE) between the posterior median estimates and true values of temperature
for each data set, $s=1,\ldots,500$, as

$\mbox{MAPE}_{s}=\frac{1}{100}\sum\limits_{j=1}^{100}\left|\hat{y}_{t+j}-y_{t+j}\right|\qquad\mbox{and}\qquad\mbox{RMSPE}_{s}=\sqrt{\frac{1}{100}\sum\limits_{j=1}^{100}(\hat{y}_{t+j}-y_{t+j})^{2}}.$

Table 2 presents the average values of these metrics calculated across the 500
data sets. Our simulation study corroborates that the linear model outperforms
the AR1 model when the objective is to perform long-term forecasting. This is
something expected, as the linear model incorporates additional auxiliary
information derived from the observed covariate, resulting in more accurate
predictions. Conversely, the AR1 model converges toward the overall mean as
time advances, as illustrated in Figure 3.

Table 2: Average values over 500 data sets of MAPE and RMSPE for long-term predictions. Model | MAPE | RMSPE
---|---|---
Linear | 0.297 | 0.334
AR1 | 0.653 | 0.800
Figure 3: True values (black dots), posterior median estimates (blue line) and
forecasted values (red lines) obtained from linear and AR1 models for the
first data set of the simulation study.

###  3.2 Example 2: spatial models

The objective of this example is to validate the effectiveness of automatic
groups construction within LGOCV for extrapolation prediction tasks in disease
mapping models, that is, when the goal is to predict disease risks or rates in
unobserved regions using spatially structured models.

We use data concerning dowry deaths across the 70 districts of Uttar Pradesh
(the most populous state in India) during the year of 2011, a form of crime
against women which is very specific to India. For each district, we have
records of the number of observed and expected dowry deaths. The expected
deaths are calculated by considering women aged 15 to 49 years as the
population at risk. In addition, we also consider the following predictors as
potential areal risk factors: sex ratio (number of females per 1,000 males),
per capita income, and number of murders per 100,000 inhabitants. In Figure 4
we plot the maps of the standardized mortality ratio (SMR) and standardized
covariates at small area level. Similar data were analyzed by Vicente et al.,
2020a  and Adin et al., 2023b  to study the association between these
potential risk factors and dowry deaths in Uttar Pradesh during the period
2001-2011.

Figure 4: Maps of the standardized mortality ratio (SMR) of dowry deaths in
the districts of Uttar Prades and standardized potential risk factors.

Let $y_{i}$ and $E_{i}$ be the number of observed and expected deaths,
respectively. Conditional on the relative risks $r_{i}$, the number of
observed cases is assumed to follow a Poisson distribution

$\begin{split}&y_{i}|r_{i},\boldsymbol{\theta}\sim\mbox{Poisson}(\mu_{i}=E_{i}r_{i})\\\
&\log{\mu_{i}}=\log{E_{i}}+\log{r_{i}}\end{split}$

where $\log{E_{i}}$ is an offset and depending on the specification of
$\eta_{i}=\log{r_{i}}$ different models are defined. Here, different spatially
structured models are considered to estimate the relative risks of dowry
deaths for each small area $i=1,\ldots,70$.

  * CAR model:

$\eta_{i}=\beta_{0}+\xi_{i},$

where $\beta_{0}$ is an intercept representing the overall log-risk and
$\xi_{i}$ is a spatially structured random effect for which the so-called BYM2
conditional autoregressive (CAR) prior distribution (Riebler et al.,, 2016) is
assumed.

  * CAR model + covariates:

$\eta_{i}=\beta_{0}+x_{i}^{{}^{\prime}}{\mbox{\boldmath{$\beta$}}}+\xi_{i},$

where $x_{i}^{{}^{\prime}}$ is the vector of standardized covariates in the
$i$-th small area and
${\mbox{\boldmath{$\beta$}}}=(\beta_{1},\beta_{2},\beta_{3})^{{}^{\prime}}$ is
the vector of fixed effects coefficients.

  * P-spline model:

$\eta_{i}=\beta_{0}+f(s_{1i},s_{2i}),$

where $f(s_{1i},s_{2i})$ is a spatial smooth function that is approximated as
$f({\bf s}_{1},{\bf s}_{2})\approx{\bf
B}_{s}{\mbox{\boldmath{$\theta$}}}_{s}$, where ${\bf B}_{s}={\bf
B}_{2}\Box{\bf B}_{1}$ is the two-dimensional B-spline basis of dimension
$70\times k_{s}$ (with $k_{s}=k_{1}k_{2}$ depending on the number of knots and
the degree of the polynomials in the bases ${\bf B}_{1}$ and ${\bf B}_{2}$),
obtained from the row-wise Kronecker product of the marginal bases for
longitude ${\bf s}_{1}=(s_{1,1},\ldots,s_{1,70})^{{}^{\prime}}$ and latitude
${\bf s}_{2}=(s_{2,1},\ldots,s_{2,70})^{{}^{\prime}}$. To achieve smoothness,
the following prior distribution is assumed for the unknown coefficients
${\mbox{\boldmath{$\theta$}}}_{s}=(\theta_{1},\ldots,\theta_{k_{s}})^{{}^{\prime}}$

${\mbox{\boldmath{$\theta$}}}_{s}\sim N({\bf 0},{\bf P}_{s}^{-})$

where ${\bf P}_{s}$ is a spatial anisotropic precision matrix that induces
different amount of smoothing in longitude an latitude dimensions (see Ugarte
et al.,, 2017).

  * P-spline model + covariates:

$\eta_{i}=\beta_{0}+x_{i}^{{}^{\prime}}{\mbox{\boldmath{$\beta$}}}+f(s_{1i},s_{2i})$

The maps illustrating the posterior median estimates of $r_{i}=\exp(\eta_{i})$
are presented in Figure 5. In general, similar geographical distributions of
mortality risks related to dowry deaths are observed across the different
spatial models. Table 3 compares the values of model selection criteria and
predictive performance measures considering LOOCV and LGOCV with automatic
groups construction (m=3) based on posterior correlations for each spatial
model. According to DIC/WAIC measures, CAR models are preferred over P-spline
models in terms of model selection criteria. In addition, predictive measures
derived from LOOCV suggest that the incorporation of spatial covariates into
the linear predictor of the CAR/P-spline models does not enhance their
predictive performance. Considering these measure, the P-spline model with
covariates exhibits the poorest performance in both the trade-off between
model fit and complexity, and in terms of predictive accuracy. In contrast,
the predictive metrics obtained through LGOCV indicate that the CAR models
perform less effectively in predicting relative risks for unobserved small
areas. This is something expected considering that CAR models induce local
correlations among neighbouring areas, while P-spline models are able to
borrow spatial information for areas that are further apart. Similar
conclusions are obtained for higher values of $m$.

Figure 5: Maps of posterior median estimates of dowry death relative risks obtained from the different spatial models. Table 3: Model comparison for spatial data: model selection criteria (posterior mean deviance $\bar{D}$, effective number of parameters $p_{D}$, DIC and the WAIC) and predictive performance measures for LOOCV and LGOCV with automatic groups construction (m=3) based on posterior correlations for each spatial model. | Model selection criteria | LOOCV | LGOCV
---|---|---|---
Model | $\bar{D}$ | $p_{D}$ | DIC | WAIC | LS | MSPE | LS | MSPE
CAR | 433.2 | 47.5 | 480.7 | 473.4 | 3.604 | 90.529 | 3.611 | 89.480
CAR+covariates | 435.4 | 45.6 | 481.1 | 476.8 | 3.616 | 105.042 | 3.614 | 99.982
Pspline | 454.0 | 30.4 | 484.4 | 492.0 | 3.600 | 93.196 | 3.485 | 76.146
Pspline+covariates | 459.5 | 30.8 | 490.3 | 501.0 | 3.667 | 101.102 | 3.508 | 82.271

As in the previous example, we should notice that LOOCV is not an appropriate
approach to estimate the models predictive performance when the task is to
predict relative risks in unobserved group of areas (extrapolation task). This
concern is particularly relevant when employing spatially structured models,
where the linear predictor of each area is highly correlated with those of its
neighbouring regions. Once again, our focus lies in validating these measures
through a simulation study. For each small area, we remove its k-order
neighbouring areas ($k=1,2$ and 3) and compute posterior predictive
distributions for the relative risks on those areas under the different models
fitted with INLA. In Figure 6 we show a map with the 1st, 2nd and 3rd-order
neighbourhood areas for a randomly selected district of Uttar Pradesh. In
Table 4 we show average values of mean relative prediction errors (MRPE) and
relative root mean square prediction errors (RRMSPE) computed over each data
set. As expected, higher MRPE and RRMSPE values are obtained as the number of
removed areas increases, which is controlled by the parameter $k$. Our
simulation study corroborates that, as pointed out by the LGOCV approach, CAR
models exhibit inferior predictive performance compared to P-spline models in
the context of spatial extrapolation. In general, we observe that better
predictive count estimates are obtained when including auxiliary information
of area-level covariates as predictors in the models.

Figure 6: Maps with 1st-order (left), 2nd-order (center) and 3rd-order (right) neighbourhood areas for a selected district of Uttar Pradesh. Table 4: Average values over 70 data sets of MRPE and RRMSPE for predicted counts. | MRPE | RRMSPE
---|---|---
Model | k=1 | k=2 | k=3 | k=1 | k=2 | k=3
CAR | 0.273 | 0.325 | 0.370 | 0.354 | 0.442 | 0.521
CAR + covariates | 0.281 | 0.299 | 0.310 | 0.357 | 0.395 | 0.410
Pspline | 0.270 | 0.316 | 0.344 | 0.347 | 0.422 | 0.470
Pspline + covariates | 0.161 | 0.171 | 0.172 | 0.215 | 0.246 | 0.253

##  4 Applications

###  4.1 Joint modelling of pancreatic cancer mortality and incidence data

Bayesian hierarchical models are a powerful tool for analyzing area-level
incidence or mortality data, providing valuable insight into latent spatial,
temporal, and spatio-temporal patterns. These models typically rely on
generalized linear mixed models that incorporate spatially and temporally
structured random effects, allowing for the smoothing of disease risks or
rates by incorporating information from neighboring areas and time periods.
Further extensions of spatio-temporal disease mapping models aim to jointly
analyze several responses by the construction of multivariate proposals based
on Gaussian Markov random fields. An exhaustive review on the topic can be
found in the work by MacNab, (2018). Alternatively, shared-component models
have been widely used to investigate closely related phenomena or diseases
that share a recognized common risk factor (see Knorr-Held and Best,, 2001;
Held et al.,, 2005).

####  4.1.1 The models

In this section, we demonstrate the use of the automatic group construction
procedure for LGOCV when the objective is to jointly model incidence and
mortality data for a highly lethal disease such as pancreatic cancer using
different spatio-temporal models. Specifically, we conduct an analysis using
yearly count data for the male population across the 105 Clinical
Commissioning Groups (CCG) of England during the period 2001-2020 (official
data provided by the National Health Service in England available at
https://www.cancerdata.nhs.uk/incidence_and_mortality). In these types of
studies, we expect that using multivariate models will enhance the performance
of individual models by leveraging the strong correlation between incidence
and mortality data.

The setting for our study is the following. Let us define as $y_{it1}$ and
$y_{it2}$ the number of mortality and incidence cases, respectively, by
health-area (CCG) $i=1,\ldots,S=105$ and time period $t=1,\ldots,T=20$
(corresponding to the period 2001-2020). Conditional on the relative risks
$r_{itj}$ for mortality ($j=1$) and incidence ($j=2$) data, we assume that the
number of observed cases follow a Poisson distribution with mean
$\mu_{itj}=E_{itj}\cdot r_{itj}$, that is,

$\begin{split}&y_{itj}|r_{itj},\boldsymbol{\theta}\sim\mbox{Poisson}(\mu_{itj}=E_{itj}\cdot
r_{itj}),\\\ &\log{\mu_{itj}}=\log{E_{itj}}+\log{r_{itj}}.\end{split}$

Here $E_{itj}$ is computed using indirect standardization as
$E_{itj}=\sum_{k}n_{itjk}\cdot m_{jk}$, where $k$ denotes the age-groups
($<25$, $[25,50)$, $[50,60)$, $[60,70)$, $[70,80)$ and $\geq 80$), $n_{itjk}$
is the population at risk in the area $i$, time $t$ and age-group $k$, and
$m_{jk}$ is the overall mortality/incidence rate in England during the whole
study period for the $k$-th age-group. Then, the log-risk is modelled as

$\log r_{itj}=\alpha_{j}+\phi_{ij}+\gamma_{tj}+\delta_{itj},$ (1)

where $\alpha_{j}$ is an intercept representing the overall log-risk,
$\phi_{ij}$ and $\gamma_{tj}$ are the spatial and temporal main effects, and
$\delta_{itj}$ is the spatio-temporal interaction for mortality and incidence
data respectively.

Three main different disease mapping models are considered to jointly estimate
pancreatic cancer mortality and incidence data under the model formulation of
Equation (1). Firstly, independent spatio-temporal models (M1) are considered.
For comparison purposes with multivariate models, if we denote as ${\bf
r}=({\bf r}_{1}^{{}^{\prime}},{\bf r}_{2}^{{}^{\prime}})^{{}^{\prime}}$ to the
joined vector of mortality and incidence risks with ${\bf
r}_{j}=(r_{11j},\ldots,r_{S1j},\ldots,r_{1Tj},\ldots,r_{STj})^{{}^{\prime}}$,
we formulate the model in matrix form as

$\begin{pmatrix}\log{\bf r}_{1}\\\ \log{\bf r}_{2}\end{pmatrix}=({\bf
I}_{2}\otimes{\bf 1}_{T}\otimes{\bf 1}_{S})\begin{pmatrix}\alpha_{1}\\\
\alpha_{2}\end{pmatrix}+({\bf I}_{2}\otimes{\bf 1}_{T}\otimes{\bf
I}_{S})\begin{pmatrix}{\mbox{\boldmath{$\phi$}}}_{1}\\\
{\mbox{\boldmath{$\phi$}}}_{2}\end{pmatrix}+({\bf I}_{2}\otimes{\bf
I}_{T}\otimes{\bf 1}_{S})\begin{pmatrix}{\mbox{\boldmath{$\gamma$}}}_{1}\\\
{\mbox{\boldmath{$\gamma$}}}_{2}\end{pmatrix}+({\bf I}_{2}\otimes{\bf
I}_{T}\otimes{\bf I}_{S})\begin{pmatrix}{\mbox{\boldmath{$\delta$}}}_{1}\\\
{\mbox{\boldmath{$\delta$}}}_{2}\end{pmatrix},$ (2)

where ${\bf 1}_{S}$ and ${\bf 1}_{T}$ are column vectors of ones of length $S$
and $T$ respectively, ${\bf I}_{S}$ and ${\bf I}_{T}$ are identity matrices of
dimension $S\times S$ and $T\times T$ respectively,
${\mbox{\boldmath{$\phi$}}}_{j}=(\phi_{1j},\ldots,\phi_{Sj})^{{}^{\prime}}$ is
a spatial random effect with intrinsic CAR prior distribution,
${\mbox{\boldmath{$\gamma$}}}_{j}=(\gamma_{1j},\ldots,\gamma_{Tj})^{{}^{\prime}}$
is a temporally structured random effect that follows a first order RW prior
distribution, and
${\mbox{\boldmath{$\delta$}}}_{j}=(\delta_{11j},\ldots,\delta_{STj})^{{}^{\prime}}$
is a spatio-temporal random effect where the four different types of
interactions originally proposed by Knorr-Held, (2000). That is, for $j=1,2$

${\mbox{\boldmath{$\phi$}}}_{j}\sim N({\bf 0},[\tau_{\phi_{j}}{\bf
R}_{\phi}]^{-}),\qquad{\mbox{\boldmath{$\gamma$}}}_{j}\sim N({\bf
0},[\tau_{\gamma_{j}}{\bf
R}_{\gamma}]^{-}),\qquad\mbox{and}\qquad{\mbox{\boldmath{$\delta$}}}_{j}\sim
N({\bf 0},[\tau_{\delta_{j}}{\bf R}_{\delta}]^{-}),$

where ${\bf R}_{\phi}$, ${\bf R}_{\gamma}$ and ${\bf R}_{\delta}$ are
structure matrices for the corresponding spatial, temporal and spatio-temporal
random effects respectively, and $\tau_{\phi_{j}}$, $\tau_{\gamma_{j}}$ and
$\tau_{\delta_{j}}$ are precision parameters. See, for example, Goicoa et al.,
(2018) for details about identifiability issues in this type of spatio-
temporal models and its implementation using INLA.

Alternatively, two distinct multivariate models that enable the incorporation
of correlation structures between mortality and incidence data are considered:
a spatio-temporal extension of the so-called M-based model described in
Vicente et al., 2020b  (M2) and models including shared-component terms for
space and time similar to the ones described in Etxeberria et al., (2023)
(M3). Denote as
${\mbox{\boldmath{$\phi$}}}=({\mbox{\boldmath{$\phi$}}}_{1}^{{}^{\prime}},{\mbox{\boldmath{$\phi$}}}_{2}^{{}^{\prime}})^{{}^{\prime}}$
and
${\mbox{\boldmath{$\gamma$}}}=({\mbox{\boldmath{$\gamma$}}}_{1}^{{}^{\prime}},{\mbox{\boldmath{$\gamma$}}}_{2}^{{}^{\prime}})^{{}^{\prime}}$
the joined vectors of spatial and temporal random effects for mortality and
incidence risks, respectively. Then, the M-model formulation (M2) is identical
to the expression in Equation (2), assuming the following prior distribution

${\mbox{\boldmath{$\phi$}}}\sim N\left({\bf
0},{\mbox{\boldmath{$\Sigma$}}}_{\phi}^{-1}\otimes{\bf
R}_{\phi}^{-}\right),\qquad\mbox{and}\qquad{\mbox{\boldmath{$\gamma$}}}\sim
N\left({\bf 0},{\mbox{\boldmath{$\Sigma$}}}_{\gamma}^{-1}\otimes{\bf
R}_{\gamma}^{-}\right),$

where ${\mbox{\boldmath{$\Sigma$}}}_{\phi}$ and
${\mbox{\boldmath{$\Sigma$}}}_{\gamma}$ are the covariance matrices between
the spatial and temporal effects, respectively, of mortality and incidence
data. See Adin et al., 2023a  for details about model implementation in INLA.
Finally, we formulate our shared component model (M3) as

$\begin{split}&\log
r_{it1}=\alpha_{1}+c_{\phi}\phi_{i}+c_{\gamma}\gamma_{t}+\delta_{it1},\\\\[4.30554pt]
&\log
r_{it2}=\alpha_{2}+\dfrac{1}{c_{\phi}}\phi_{i}+\dfrac{1}{c_{\gamma}}\gamma_{t}+\delta_{it2},\end{split}$

where $\phi_{i}$ is a shared spatial effect with intrinsic CAR prior
distribution, $\gamma_{t}$ is a shared temporal effect with first order RW
prior distribution, and $\delta_{itj}$ is a specific spatio-temporal random
effect to model the interaction in mortality and incidence log-risks,
respectively. Finally, $c_{\phi}$ and $c_{\gamma}$ are scaling parameters to
allow a different “risk gradient” (on the log-scale) to be associated with the
spatial and temporal components, respectively, for both responses. A summary
with the main prior distributions for the three models considered in this
section has been included in Table 5.

Table 5:  Summary of the prior distributions for the spatial, temporal and
spatio-temporal random effects under the three models considered in Section
4.1.

| M1 | M2 | M3
---|---|---|---
| (Independent models) | (M-model) | (Shared-component model)
Spatial effect | ${\mbox{\boldmath{$\phi$}}}_{j}\sim N({\bf 0},[\tau_{\phi_{j}}{\bf R}_{\phi}]^{-}),\quad j=1,2$ | $\begin{pmatrix}{\mbox{\boldmath{$\phi$}}}_{1}\\\ {\mbox{\boldmath{$\phi$}}}_{2}\end{pmatrix}\sim N\left({\bf 0},{\mbox{\boldmath{$\Sigma$}}}_{\phi}^{-1}\otimes{\bf R}_{\phi}^{-}\right)$ | ${\mbox{\boldmath{$\phi$}}}\sim N({\bf 0},[\tau_{\phi}{\bf R}_{\phi}]^{-})$
Temporal effect | ${\mbox{\boldmath{$\gamma$}}}_{j}\sim N({\bf 0},[\tau_{\gamma_{j}}{\bf R}_{\gamma}]^{-}),\quad j=1,2$ | $\begin{pmatrix}{\mbox{\boldmath{$\gamma$}}}_{1}\\\ {\mbox{\boldmath{$\gamma$}}}_{2}\end{pmatrix}\sim N\left({\bf 0},{\mbox{\boldmath{$\Sigma$}}}_{\gamma}^{-1}\otimes{\bf R}_{\gamma}^{-}\right)$ | ${\mbox{\boldmath{$\gamma$}}}\sim N({\bf 0},[\tau_{\gamma}{\bf R}_{\gamma}]^{-})$
Spatio-temporal | ${\mbox{\boldmath{$\delta$}}}_{j}\sim N({\bf 0},[\tau_{\delta_{j}}{\bf R}_{\delta}]^{-})\quad j=1,2$
effect | $\mbox{Type I: }{\bf R}_{\delta}={\bf I}_{T}\otimes{\bf I}_{S},\qquad\mbox{Type II: }{\bf R}_{\delta}={\bf R}_{\gamma}\otimes{\bf I}_{S}$
$\mbox{Type III: }{\bf R}_{\delta}={\bf 1}_{T}\otimes{\bf
R}_{\phi},\qquad\mbox{Type IV: }{\bf R}_{\delta}={\bf R}_{\gamma}\otimes{\bf
R}_{\phi}$

To fit the models, improper uniform prior distributions are given to all the
standard deviations (square root inverse of precision parameters), and a
Gamma(10,10) distribution is considered for the scaling parameters $c_{\phi}$
and $c_{\gamma}$ of the shared-component model. Finally, a vague zero mean
normal distribution with a precision close to zero (0.001) is given to the
model intercepts.

####  4.1.2 Results

Table 6 presents a comparison of the different models used to estimate
pancreatic cancer mortality and incidence relative risks in terms of model
selection criteria and predictive performance measures. Both DIC/WAIC and
predictive measures using a LOOCV approach consistently indicate that Type II
interaction models should be selected. Among the considered models, the
shared-component model M3 exhibits a slightly better level of performance. To
compute comparable LS measures from the posterior predictive density functions
$\pi(Y_{itj}|\boldsymbol{y}_{-I_{itj}})$, we use the groups $I_{itj}$ derived
from Type II interaction models with number of level sets $m=5$. As shown in
Table 6, very close values are obtained for both multivariate modelling
approaches M2 and M3.

An interesting aspect of the automatic group construction approach, which
relies on the posterior correlation matrix of the linear predictor, is its
inherent interpretability. For each testing point (area $i$, time $t$ and
response $j$), we can examine its derived groups to identify which structured
effect of the model contributes to a stronger correlation in the risk
estimates. In the independent model (M1), the groups assigned to all the
testing points consistently correspond to adjacent years. In the M-model (M2),
while the temporally structured random effect remains as the dominant pattern
in the model, there are some testing points where the group of highly
correlated data points also includes spatial neighbours (areas that share a
common border). Finally, the groups derived from the shared-component model
(M3) reveal a strong correlation between both responses (mortality and
incidence data) for the majority of areas in the year 2011. This finding
suggest that the shared-component model would be the most suitable choice for
predicting pancreatic cancer mortality based on incidence data.

Table 6: Pancreatic cancer data: model selection criteria (DIC) and predictive performance measures (Logarithmic Score) with automatic groups construction for $m=3,5$, and $10$. | DIC | Logarithmic Score
---|---|---
| $m=3$ | $m=5$ | $m=10$
M1 - Type I | 26095.3 | 3.113 | 3.115 | 3.122
M1 - Type II | 26065.4 | 3.109 | 3.112 | 3.119
M1 - Type III | 26120.9 | 3.115 | 3.118 | 3.125
M1 - Type IV | 26074.1 | 3.110 | 3.112 | 3.119
M2 - Type I | 26000.0 | 3.099 | 3.101 | 3.104
M2 - Type II | 25961.7 | 3.094 | 3.093 | 3.094
M2 - Type III | 26018.8 | 3.101 | 3.102 | 3.106
M2 - Type IV | 25973.5 | 3.095 | 3.095 | 3.096
M3 - Type I | 25995.0 | 3.098 | 3.099 | 3.106
M3 - Type II | 25953.9 | 3.092 | 3.092 | 3.099
M3 - Type III | 26012.1 | 3.099 | 3.100 | 3.107
M3 - Type IV | 25964.7 | 3.093 | 3.093 | 3.100

###  4.2 Spatial Compositional Data: the case of Arabidopsis thaliana

Compositional Data Analysis has gained popularity for comprehending processes
wherein values correspond to $D$ distinct categories, with their sum
constituting a constant. These values typically represent proportions (or
percentages) summing up to one (or 100). The data stemming from such processes
are commonly referred as Compositional Data (CoDa). Without loss of
generality, we assume the constant to be one. A substantial body of literature
exists, delving into the application of CoDa analysis across various domains,
including Ecology (Kobal et al.,, 2017; Douma and Weedon,, 2019), Geology
(Buccianti and Grunsky,, 2014; Engle and Rowan,, 2014), Genomics (Tsilimigras
and Fodor,, 2016; Shi et al.,, 2016; Washburne et al.,, 2017; Creus Martí et
al.,, 2022), Environmental Sciences (Aguilera et al.,, 2021; Mota-Bertran et
al.,, 2022) or Medicine (Dumuid et al.,, 2018; Fairclough et al.,, 2018).

Recently, two approaches have been introduced to incorporate the CoDa analysis
within the framework of latent Gaussian models. One involves utilizing
Dirichlet regression (Martínez-Minaya et al.,, 2023), while the other proposes
the Logistic Normal Dirichlet Model (LNDM Martínez-Minaya and Rue,, 2023),
which mainly uses logistic-normal distribution with Dirichlet covariance
through the additive log-ratio transformation as likelihood promoting
interpretation as log-ratios. In this example, we focus on the latter,
illustrating the utility of an automatic group construction procedure for
LGOCV in a spatial CoDa problem: the case of Arabidopsis thaliana.

####  4.2.1 The data

We conduct our study using a dataset consisting of $N=301$ accessions of the
annual plant Arabidopsis thaliana located on the Iberian Peninsula. For each
accession, the genetic cluster (GC) membership proportions for the four
genetic clusters (GC1, GC2, GC3, GC4) are available as identified in Martínez-
Minaya et al., (2019). These probabilities, shown in Figure 7, sum up to one.
Our main interest lies in estimating the membership probability, which in this
specific context can be interpreted as the habitat suitability for each
genetic cluster. To achieve this, we employ LNDMs incorporating climate-
related covariates and spatially structured effects within the linear
predictor. Specifically, we employ two bioclimatic variables to define the
climatic aspect: annual mean temperature (BIO1) and annual precipitation
(BIO12). The dataset is available at http://doi.org/10.5281/zenodo.2552025
(Martínez-Minaya et al.,, 2019). Prior to analysis, climate covariates are
standardized.

####  4.2.2 The models

As mentioned earlier, four categories are employed in this problem: GC1, GC2,
GC3 and GC4. So, we deal with proportions in the simplex $\mathbb{S}^{4}$. To
construct the LNDM, the additive log-ratio transformation ($alr$) is required,
and we select GC4 as the reference category because it has the lowest variance
in its logarithm. We are thus dealing with a three dimensional normal
distribution with Dirichlet covariance,
$\mathcal{ND}(\boldsymbol{\mu},\boldsymbol{\Sigma})$. The data, as well as the
transformed data, are displayed in Figure 7.

Figure 7: CoDa of 301 accessions of the annual plant Arabidopsis thaliana
located on the Iberian Peninsula. A: Membership proportion to each of the 4
genetic clusters extracted from Martínez-Minaya et al., (2019). B:
$alr$-coordinates generated from the $alr$ transformation using GC4 as
reference.

The models employed to address this issue present the subsequent structure:

$\displaystyle\begin{array}[]{rcl}alr(\boldsymbol{Y})&\sim&\mathcal{ND}((\boldsymbol{\mu}^{(1)},\ldots,\boldsymbol{\mu}^{(3)}),\boldsymbol{\Sigma}),\\\\[4.30554pt]
\boldsymbol{\mu}^{(d)}&=&\boldsymbol{X}\boldsymbol{\beta}^{(d)}+\boldsymbol{\omega}^{(d)}\,,\
d=1,\ldots,3,\end{array}$ (5)

where $alr(\boldsymbol{Y})$ represents the $alr$ transformation applied to the
compositional variable. From now on, we will refer to the transformed variable
as $alr$-coordinates, and the vector
$\boldsymbol{\mu}^{(d)}=(\mu^{(d)}_{1},\ldots,\mu^{(d)}_{301})^{{}^{\prime}}$
is the mean of the $d$-th $alr$-coordinate. The design matrix $\boldsymbol{X}$
of dimension $301\times 3$ consists of ones in its first column and climate
covariate values in the remaining columns. The spatial random effect for each
$d$-th $alr$-coordinate is denoted by $\boldsymbol{\omega}^{(d)}$,
characterized by a Matérn covariance structure. Specifically,
$\boldsymbol{\omega}^{(d)}$ follows a normal distribution
$\mathcal{N}(0,\boldsymbol{Q}^{-1}(\sigma_{\boldsymbol{\omega}},\phi))$, where
$\sigma_{\boldsymbol{\omega}}$ represents the standard deviation of the
spatial effect, and $\phi$ its range. The parameter vector
$\boldsymbol{\beta}=(\beta^{(d)}_{0},\beta^{(d)}_{1},\beta^{(d)}_{2})^{{}^{\prime}}$
corresponds to fixed effects. The latent field comprises both the fixed
effects parameters and realizations of the random field. This configuration
can be succinctly presented as:

$\boldsymbol{\mathcal{X}}=\\{\boldsymbol{\beta}^{(d)},\boldsymbol{\omega}^{(d)}:d=1,\ldots,3\\}.$

In contrast, the vector
$\boldsymbol{\theta}_{1}=\\{\sigma_{d}^{2},\gamma:d=1,\ldots,3\\}$ encompasses
hyperparameters linked to the likelihood. On the other hand,
$\boldsymbol{\theta}_{2}=\\{\sigma_{\boldsymbol{\omega}},\phi\\}$ comprises
hyperparameters pertaining to the spatial random effect. These combined
elements collectively form the set of hyperparameters. Gaussian prior
distributions are assumed for the fixed effects, while PC-priors are employed
for the hyperparameters, following Simpson et al., (2017). Based on the model
structure defined in Equation (5) where a spatial effect is included, we focus
on the eight structures presented in Martínez-Minaya and Rue, (2023):

  * •

Type I: share the same parameters for fixed effects, and do not include
spatial random effects.

  * •

Type II: have different parameters for fixed effects, and do not include
spatial random effects.

  * •

Type III: share the same parameters for fixed effects, and share the same
spatial effect.

  * •

Type IV: different parameters for fixed effects, and share the same spatial
effect.

  * •

Type V: share the same parameters for fixed effects, and the spatial effects
between linear predictors are proportional.

  * •

Type VI: different parameters for fixed effects, and the spatial effects
between linear predictors are proportional.

  * •

Type VII: share the same parameters for fixed effects, and different
realisations of the spatial effect for each linear predictor.

  * •

Type VIII: different parameters for fixed effects, and different realisations
of the spatial effect for each linear predictor.

Note that in Type V and Type VI structures, realisations of the spatial field
are the same, but a proportionality hyperparameter is added in two of the
three linear predictors, denoted as $\alpha^{(1)}$ and $\alpha^{(2)}$. On the
other hand, in Type VII and Type VIII structures, although realisations of
random effects are different, they share the same hyperparameters. Table 7
provides an overview of the distinct structures, their associated latent
fields, and the accompanying set of hyperparameters.

Table 7: Different structures included in the model in an additive way with
their corresponding latent field and the set hyperparameters to be estimated.

Models | Predictor | Latent Field ($\boldsymbol{\mathcal{X}}$) | Hyperparameters ($\boldsymbol{\theta}$)
---|---|---|---
Type I | $\boldsymbol{X}\boldsymbol{\beta}$ | $\\{\beta_{0},\ldots,\beta_{M}\\}$ | $\\{\sigma_{d}^{2},\gamma\\}$
Type II | $\boldsymbol{X}\boldsymbol{\beta}^{(d)}$ | $\\{\beta_{0}^{(d)},\ldots,\beta_{M}^{(d)}\\}$ | $\\{\sigma_{d}^{2},\gamma\\}$
Type III | $\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\omega}$ | $\\{\beta_{0},\beta_{1},\beta_{2},\omega_{1},\ldots,\omega_{N}\\}$ | $\\{\sigma_{d}^{2},\gamma,\sigma_{\boldsymbol{\omega}},\phi\\}$
Type IV | $\boldsymbol{X}\boldsymbol{\beta}^{(d)}+\boldsymbol{\omega}$ | $\\{\beta_{0}^{(d)},\beta_{1}^{(d)},\beta_{2}^{(d)},\omega_{1},\ldots,\omega_{N}\\}$ | $\\{\sigma_{d}^{2},\gamma,\sigma_{\boldsymbol{\omega}},\phi\\}$
Type V | $\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\omega}^{*(d)}$ | $\\{\beta_{0},\beta_{1},\beta_{2},\omega_{1},\ldots,\omega_{N}\\}$ | $\\{\sigma_{d}^{2},\gamma,\sigma_{\boldsymbol{\omega}},\phi,\alpha^{(1)},\alpha^{(2)}\\}$
Type VI | $\boldsymbol{X}\boldsymbol{\beta}^{(d)}+\boldsymbol{\omega}^{*(d)}$ | $\\{\beta_{0}^{(d)},\beta_{1}^{(d)},\beta_{2}^{(d)},\omega_{1},\ldots,\omega_{N}\\}$ | $\\{\sigma_{d}^{2},\gamma,\sigma_{\boldsymbol{\omega}},\phi,\alpha^{(1)},\alpha^{(2)}\\}$
Type VII | $\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\omega}^{(d)}$ | $\\{\beta_{0},\beta_{1},\beta_{2},\omega_{1}^{(d)},\ldots,\omega_{N}^{(d)}\\}$ | $\\{\sigma_{d}^{2},\gamma,\sigma_{\boldsymbol{\omega}},\phi\\}$
Type VIII | $\boldsymbol{X}\boldsymbol{\beta}^{(d)}+\boldsymbol{\omega}^{(d)}$ | $\\{\beta_{0}^{(d)},\beta_{1}^{(d)},\beta_{2}^{(d)},\omega_{1}^{(d)},\ldots,\omega_{N}^{(d)}\\}$ | $\\{\sigma_{d}^{2},\gamma,\sigma_{\boldsymbol{\omega}},\phi\\}$

####  4.2.3 LGOCV in Compositional Data

As highlighted in Martínez-Minaya and Rue, (2023), excluding a category from a
CoDa point during cross-validation procedures may not be practical due to the
inherent constraint imposed by CoDa, where the sum of its components must
equal one. This implies that the remaining categories hold valuable
information about the category we intend to exclude. Consequently, the
remaining log-ratio coordinates offer insights into the category we have
removed during cross-validation. In this manner, the concept of ‘friendship’
emerges. Accordingly, we can assert that the first $alr$-coordinate of $i$-th
data point is associated with the second $alr$-coordinate of $i$-th data
point, thereby contributing pertinent information.

Hence, to perform cross-validation for the $i$-th data point and the $d$-th
$alr$-coordinate, it becomes necessary to exclude the values associated with
all $alr$-coordinates for that data point and for the corresponding group. To
do this, the group $I_{i}^{(d)}$ must be constructed. If it is determined that
the $i$-th data point in its $d$-th $alr$-coordinate is associated with the
data point indexed as $i+1$ in its $d$-th $alr$-coordinate, the data point
$i+1$ will also be included in the set $I_{i}^{(d)}$ in all of its
$alr$-coordinates. It should be noted that $I_{i}^{(d)}$ could be different
from $I_{i}^{(d^{*})}$, being $d\neq d^{*}$, as even though we are referring
to the same compositional point $i$, we are modelling each $alr$-coordinate.
So we need to compute LS measures for CoDa in the following way:

LS $\displaystyle=$
$\displaystyle-\frac{1}{N\cdot(D-1)}\sum\limits_{d=1}^{D-1}\sum\limits_{i=1}^{N}\log\left(\pi(\text{alr}(\boldsymbol{y})_{i}^{(d)}\mid\text{alr}(\boldsymbol{y})_{-I_{i}^{(d)}}^{\bullet})\right),$

where $alr(\boldsymbol{y})_{i}^{(d)}$ is the observed vector for the $i$-th
data point and the $d$-th $alr$-coordinate, and
$alr(\boldsymbol{y})_{-i}^{\bullet}$ represents the observed data in
$alr$-coordinates excluding the $I_{i}^{(d)}$ data points with its
corresponding $D-1$ $alr$-coordinates.

####  4.2.4 Results

Table 8 presents a comparison of various models used for estimating the
membership proportions to the different genetic clusters of Arabidopsis
thaliana. As in the previous section, the evaluation is based on model
selection criteria and predictive performance measures. To compute comparable
LS measures from the posterior predictive density functions, we use the groups
derived from Type VIII model. We use values of $m=3,5$ and $10$ for
constructing the groups and compute performance measures based on LGOCV.
Although, Type VII model seems to have a good predictive performance, DIC/WAIC
and predictive measures based on LGOCV approaches consistently suggest that
Type VIII model is the most suitable choice.

Table 8: Arabidopsis thaliana data: model selection criteria (DIC) and predictive performance measures (Logarithmic Score) with automatic groups construction for $m=3,5$, and $10$. Type | DIC | Logarithmic Score
---|---|---
$m=3$ | $m=5$ | $m=10$
I | 3353.403 | 1.894 | 1.897 | 1.904
II | 3294.855 | 1.869 | 1.873 | 1.883
III | 3202.082 | 1.784 | 1.780 | 1.795
IV | 3145.198 | 1.757 | 1.756 | 1.778
V | 3059.359 | 1.470 | 1.579 | 1.634
VI | 3002.910 | 1.396 | 1.526 | 1.588
VII | 2754.844 | 1.427 | 1.513 | 1.567
VIII | 2741.440 | 1.415 | 1.513 | 1.568

In Figure 8, we depict automatic groups formed from the posterior correlations
in the Type VIII model. Specifically, this representation pertains to the
three $alr$-coordinates associated with the $88$-th data point when $m=10$. By
employing the ‘friendship’ feature, we made the assumption that these three
observations share a bond, meaning that if we intend to create groups for the
first $alr$-coordinate of observation 88, we must not only include that
$alr$-coordinate of observation $i=88$ in $I_{88}^{(1)}$, but also the second
and third $alr$-coordinates for that data point. As shown in Figure 8, it
becomes evident that the groups established for data point $i=88$ differ for
each $alr$-coordinate, i.e., $I_{88}^{(1)}\neq I_{88}^{(2)}\neq I_{88}^{(3)}$.
This observation underscores not only our assumption of distinct realizations
of the spatial effect for each $alr$-coordinate but also varying effects of
the climatic variables. This is due to the fact that posterior correlations
consider all components of the linear predictor.

Figure 8: Automatically constructed groups (with $m=10$) under Type VIII model
for data point $i=88$ located in the region of Madrid. The blue circle
represents the focal point, while the crosses indicate the points that belong
to the constructed group $I_{88}^{(d)}$, $d=1,\ldots,3$.

###  4.3 Space-time models for the United Kingdom wind speed data

In what follows, we describe how to use the procedure outlined in Section 2 to
perform cross-validation on space-time models for wind speed data in the
United Kingdom (UK). We consider a dataset from a regularly updated 20th and
21st-century database from the European Climate Assessment & Dataset (ECA&D)
project (Klein Tank et al.,, 2002). We select daily wind speed measurements,
in meters per second (m/s), from July 2021 at 189 stations in Ireland and the
UK. The locations of the stations are presented as coloured dots in Figure 9.
The observations at each location are displayed with different colors in this
figure as grouped time series over July 2021. The groups are constructed
according to the spatial grid overlayed in dashed gray line. Similar patterns
can be observed at nearby locations.

Figure 9: Map displaying the station locations with the data presented as
grouped time series (respect to the spatial grid overlayed in dashed gray
lines) over the period of July 2021.

####  4.3.1 The models

We compare four models with different characterizations of the space-time
dynamics in the UK wind speed data. Consider the process $Y(\mathbf{s},t)$,
which represents the wind speed at location $\mathbf{s}$ and at time $t$. We
model observations $y(\mathbf{s}_{i},t_{i})$ from this process collected at
site $\mathbf{s}_{i}$ and time $t_{i}$ with a Gaussian distribution
N($\mu(\mathbf{s}_{i},t_{i})$, $\sigma^{2}_{e}$), where $\sigma^{2}_{e}$ is a
common variance parameter and the mean function is given by

$\mu(\mathbf{s},t)=\beta_{0}+\beta_{1}x(\mathbf{s})+r(\mathbf{s})+u(\mathbf{s},t),$
(6)

where $\beta_{0}$ is the intercept, which is constant in space and time, $x$
is the elevation (in meters), and $\beta_{1}$ represents the associated
regression coefficient. The $r(.)$ term is a spatial field, which accommodate
regional persistent spatial variations and is modeled with a Matérn covariance
using the stochastic partial differential equation (SPDE) approach described
in Lindgren et al., (2011) with smoothness parameter fixed to $\alpha=2$,
range $\rho_{r}$, and marginal standard error $\sigma_{r}$, which need to be
estimated. Space-time interactions are captured through the space-time random
field $u(.)$, defined using the models proposed in Lindgren et al., (2023).
These models are defined as solutions to a diffusion-based family of equations
that induces physically realistic processes, such as wind. Specifically, it is
the solution to the following SPDE

$\gamma_{e}(\gamma_{s}^{2}-\Delta)^{\alpha_{e}/2}\left(\gamma_{t}\frac{\partial}{\partial
t}+(\gamma_{s}^{2}-\Delta)^{\alpha_{s}}\right)^{\alpha_{t}}\
u(s,t)=\mathcal{W}(s,t).$ (7)

The power parameters $\alpha_{t}$, $\alpha_{s}$ and $\alpha_{e}$ in Equation
(7) specify the smoothness of $u(.)$, whereas $\gamma_{s}$, $\gamma_{t}$, and
$\gamma_{e}$ govern the marginal properties and are related to dynamics of the
process. We consider here four cases of the models proposed in Lindgren et
al., (2023): Model $M_{A}$ (with $\alpha_{t}=1$, $\alpha_{s}=0$, and
$\alpha_{e}=2$); Model $M_{B}$ (with $\alpha_{t}=1$, $\alpha_{s}=2$, and
$\alpha_{e}=1$); Model $M_{C}$ (with $\alpha_{t}=2$, $\alpha_{s}=0$, and
$\alpha_{e}=2$); and Model $M_{D}$ (with $\alpha_{t}=2$, $\alpha_{s}=2$, and
$\alpha_{e}=0$). When $\alpha_{s}=0$ (Models $M_{A}$ and $M_{C}$), the space-
time SPDE in Equation (7) generates a precision matrix that can be factorized
into two precision matrices, one for the purely spatial model and another for
the purely temporal model. While separability is broadly used due to a
substantial reduction of the computation time when working with large data
sets with complex dependencies, it is usually too simplistic to model natural
phenomena. Due to recent developments, with the SPDE approach both, separable
or non-separable, have similar computation time.

For fitting the models, we considered a flat prior for $\beta_{0}$ and a
Gaussian with zero mean and variance $0.01$ for $\beta_{1}$. The penalized
complexity prior proposed in Fuglstad et al., (2018), are used for the spatial
and temporal range parameters of $r$ and $u$, as well as for all the variance
parameters $\sigma^{2}_{e}$, $\sigma^{2}_{r}$ (marginal variance of $r$), and
$\sigma^{2}_{u}$ (marginal variance of $u$).

####  4.3.2 Results

We use the methods described in Section 2 to perform LGOCV on the four models
presented above. For all the models, we set the groups to those generated from
model $M_{A}$. Figure 10 illustrates the groups for $m=10$, at time $1$ for
seven locations from different parts of the spatial domain. Even though $m$ is
fixed to 10, the total number of observations in each group is always larger
than $10$ (see the last row of the summary table on the top left of this
figure). This is due to several observation pairs having very similar
correlation values. For instance, the group for station ‘id’ = 5 is of size 17
and for ‘id’ = 131 it is of size 50, which is the number we set as the maximum
group size. In the top left table, for each station, the removed neighbours
are classified into self (same location and time), time (same location but
different times), spatial (same time but different locations), or space-time
(different location and different time). Stations located in isolated places
such as ‘id’ = 131 tend to select more neighbours in time. Indeed, only three
spatial neighbors were selected and the entire time series of this station was
left out, that is fifteen observations before and fifteen after the selected
time. On the other hand, stations that have several several neighbors close by
such as ‘id’ = 18 tend to mostly include spatial and space-time neighbours.

Figure 10: Location of the stations and examples of seven observations at time
16 with their respective automatic selected groups from the LGOCV method.

Finally, Table 9 shows the DIC for the different models and the LS values for
$m=3$, $m=5$, and $m=10$. Model $M_{B}$ provides the best fit with the lowest
DIC. This model also performed the best for prediction purposes, as shown by
the LS values for all level sets. Model $M_{D}$ has the second lowest DIC but
the worst LS, indicating that it is too smooth to predict the highly variable
wind trajectories (see Figure 9).

Table 9: wind speed data: model selection criteria (DIC) and predictive performance measures (Logarithmic Score) with automatic groups construction for $m=3$, $m=5$, and $m=10$. Model for $u$ | DIC | Logarithmic Score
---|---|---
m=3 | m=5 | m=10
$M_{A}$ | 3499.78 | 6202.79 | 6692.26 | 7461.16
$M_{B}$ | 1622.16 | 6054.93 | 6563.31 | 7413.71
$M_{C}$ | 3433.40 | 6214.94 | 6678.89 | 7435.29
$M_{D}$ | 2078.38 | 7123.03 | 7587.39 | 8342.42

##  5 Discussion

Cross-validation techniques are at the core of any prediction task as they
reveal the adequacy of a statistical model to unobserved data. Among the
several cross-validation methods developed over the years, LOOCV stands out as
the prevailing choice in numerous practical applications due to its simplicity
and straightforward implementation. Although widely used, this method is only
optimal under the assumption of independent and identically distributed (iid)
data, an assumption rarely verified in many practical applications. In
contrast, the newly developed LGOCV method handles non iid data with multiple
correlation structures by automatically excluding a determined set of data
points from the training set, expanding beyond the practice of solely
excluding individual observations. The R-INLA package provides a convenient
implementation of LGOCV for laten Gaussian models that is highly efficient and
readily accessible.

In this paper, we showed that the automatic group construction procedure for
LGOCV provides a robust measure of predictive performance in many settings
where dependencies in space and/or time are present. Our simulation studies
and real data analysis underscore that a common practice in modelling
dependent data often leads to erroneous model selection and interpretations.
In contrast, the LGOCV method, integrated into the R-INLA package, stands as
an improved alternative to LOOCV.

## Acknowledgments

This research has been supported by projects
PID2020-113125RB-I00/MCIN/AEI/10.13039/501100011033 (AA) and
PID2020-115882RB-I00 (JM-M).

## References

  * (1) Adin, A., Goicoa, and Ugarte, M. D. (2023a).  Multivariate Disease Mapping Models to Uncover Hidden Relationships Between Different Cancer Sites.  In Larriba, Y., editor, Statistical Methods at the Forefront of Biomedical Advances, pages 1–20. Springer International Publishing.
  * (2) Adin, A., Goicoa, T., Hodges, J. S., Schnell, P., and Ugarte, M. D. (2023b).  Alleviating confounding in spatio-temporal areal models with an application on crimes against women in India.  Statistical Modelling, 23(1):9–30.
  * Aguilera et al., (2021) Aguilera, A., Bautista, F., Gutiérrez-Ruiz, M., Ceniceros-Gómez, A. E., Cejudo, R., and Goguitchaichvili, A. (2021).  Heavy metal pollution of street dust in the largest city of Mexico, sources and health risk assessment.  Environmental Monitoring and Assessment, 193(4):1–16.
  * Arlot and Celisse, (2010) Arlot, S. and Celisse, A. (2010).  A survey of cross-validation procedures for model selection.  Statistics Surveys, 4:40–79.
  * Bergmeir and Benítez, (2012) Bergmeir, C. and Benítez, J. M. (2012).  On the use of cross-validation for time series predictor evaluation.  Information Sciences, 191:192–213.
  * Buccianti and Grunsky, (2014) Buccianti, A. and Grunsky, E. (2014).  Compositional data analysis in geochemistry: Are we sure to see what really occurs during natural processes?  Journal of Geochemical Exploration, 141:1–5.
  * Bürkner et al., (2021) Bürkner, P.-C., Gabry, J., and Vehtari, A. (2021).  Efficient leave-one-out cross-validation for Bayesian non-factorized normal and Student-t models.  Computational Statistics, 36(2):1243–1261.
  * Creus Martí et al., (2022) Creus Martí, I., Moya, A., and Santonja, F. (2022).  Bayesian Hierarchical Compositional Models for Analysing Longitudinal Abundance Data from Microbiome Studies.  Complexity, 2022.
  * Douma and Weedon, (2019) Douma, J. C. and Weedon, J. T. (2019).  Analysing continuous proportions in Ecology and Evolution: A practical introduction to beta and Dirichlet regression.  Methods in Ecology and Evolution, 10(9):1412–1430.
  * Dumuid et al., (2018) Dumuid, D., Stanford, T. E., Martin-Fernández, J.-A., Pedišić, Ž., Maher, C. A., Lewis, L. K., Hron, K., Katzmarzyk, P. T., Chaput, J.-P., Fogelholm, M., et al. (2018).  Compositional data analysis for physical activity, sedentary time and sleep research.  Statistical Methods in Medical Research, 27(12):3726–3738.
  * Engle and Rowan, (2014) Engle, M. A. and Rowan, E. L. (2014).  Geochemical evolution of produced waters from hydraulic fracturing of the Marcellus Shale, Northern Appalachian Basin: A multivariate compositional data analysis approach.  International Journal of Coal Geology, 126:45–56.
  * Etxeberria et al., (2023) Etxeberria, J., Goicoa, T., and Ugarte, M. D. (2023).  Using mortality to predict incidence for rare and lethal cancers in very small areas.  Biometrical Journal, 65(3):2200017.
  * Fairclough et al., (2018) Fairclough, S. J., Dumuid, D., Mackintosh, K. A., Stone, G., Dagger, R., Stratton, G., Davies, I., and Boddy, L. M. (2018).  Adiposity, fitness, health-related quality of life and the reallocation of time between children’s school day activity behaviours: A compositional data analysis.  Preventive Medicine Reports, 11:254–261.
  * Fuglstad et al., (2018) Fuglstad, G. A., Simpson, D., Lindgren, F., and Rue, H. (2018).  Constructing priors that penalize the complexity of Gaussian random fields.  Journal of the American Statistical Association, 114(525):445–452.
  * Gelman et al., (1995) Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (1995).  Bayesian data analysis.  Chapman and Hall/CRC.
  * Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association, 102(477):359–378.
  * Goicoa et al., (2018) Goicoa, T., Adin, A., Ugarte, M. D., and Hodges, J. S. (2018).  In spatio-temporal disease mapping models, identifiability constraints affect PQL and INLA results.  Stochastic Environmental Research and Risk Assessment, 32(3):749–770.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., Friedman, J. H., and Friedman, J. H. (2009).  The elements of statistical learning: data mining, inference, and prediction, volume 2.  Springer.
  * Held et al., (2005) Held, L., Natário, I., Fenton, S. E., Rue, H., and Becker, N. (2005).  Towards joint disease mapping.  Statistical Methods in Medical Research, 14(1):61–82.
  * Held et al., (2010) Held, L., Schrodle, B., and Rue, H. (2010).  Posterior and Cross-validatory Predictive Checks: A Comparison of MCMC and INLA.  In Statistical Modelling and Regression Structures, pages 111–131. Springer.
  * Klein Tank et al., (2002) Klein Tank, A., Wijngaard, J., Können, G., Böhm, R., Demarée, G., Gocheva, A., Mileta, M., Pashiardis, S., Hejkrlik, L., Kern-Hansen, C., et al. (2002).  Daily dataset of 20th-century surface air temperature and precipitation series for the European Climate Assessment.  International Journal of Climatology: A Journal of the Royal Meteorological Society, 22(12):1441–1453.
  * Knorr-Held, (2000) Knorr-Held, L. (2000).  Bayesian modelling of inseparable space-time variation in disease risk.  Statistics in Medicine, 19(17-18):2555–2567.
  * Knorr-Held and Best, (2001) Knorr-Held, L. and Best, N. G. (2001).  A shared component model for detecting joint and selective clustering of two diseases.  Journal of the Royal Statistical Society Series A: Statistics in Society, 164(1):73–85.
  * Kobal et al., (2017) Kobal, M., Kastelec, D., and Eler, K. (2017).  Temporal changes of forest species composition studied by compositional data approach.  iForest-Biogeosciences and Forestry, 10(4):729–738.
  * Kuhn et al., (2013) Kuhn, M., Johnson, K., et al. (2013).  Applied predictive modeling, volume 26.  Springer.
  * Lindgren et al., (2023) Lindgren, F., Bakka, H., Bolin, D., Krainski, E., and Rue, H. (2023).  A diffusion-based spatio-temporal extension of Gaussian Matérn fields.  arXiv:2006.04917.
  * Lindgren et al., (2011) Lindgren, F., Rue, H., and Lindström, J. (2011).  An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 73(4):423–498.
  * Liu and Rue, (2022) Liu, Z. and Rue, H. (2022).  Leave-group-out cross-validation for latent Gaussian models.  arXiv:2210.04482.
  * MacNab, (2018) MacNab, Y. C. (2018).  Some recent work on multivariate Gaussian Markov random fields.  Test, 27(3):497–541.
  * Martínez-Minaya et al., (2019) Martínez-Minaya, J., Conesa, D., Fortin, M.-J., Alonso-Blanco, C., Picó, F. X., and Marcer, A. (2019).  A hierarchical bayesian beta regression approach to study the effects of geographical genetic structure and spatial autocorrelation on species distribution range shifts.  Molecular Ecology Resources, 19(4):929–943.
  * Martínez-Minaya et al., (2023) Martínez-Minaya, J., Lindgren, F., López-Quílez, A., Simpson, D., and Conesa, D. (2023).  The Integrated Nested Laplace Approximation for Fitting Dirichlet Regression Models.  Journal of Computational and Graphical Statistics, 32(3):805–823.
  * Martínez-Minaya and Rue, (2023) Martínez-Minaya, J. and Rue, H. (2023).  A flexible Bayesian tool for CoDa mixed models: logistic-normal distribution with Dirichlet covariance.  arXiv:2308.13928.
  * Mota-Bertran et al., (2022) Mota-Bertran, A., Saez, M., and Coenders, G. (2022).  Compositional and Bayesian inference analysis of the concentrations of air pollutants in Catalonia, Spain.  Environmental Research, 204:112388.
  * Rabinowicz and Rosset, (2022) Rabinowicz, A. and Rosset, S. (2022).  Cross-validation for correlated data.  Journal of the American Statistical Association, 117(538):718–731.
  * Riebler et al., (2016) Riebler, A., Sørbye, S. H., Simpson, D., and Rue, H. (2016).  An intuitive Bayesian spatial model for disease mapping that accounts for scaling.  Statistical Methods in Medical Research, 25(4):1145–1165.
  * Roberts et al., (2017) Roberts, D. R., Bahn, V., Ciuti, S., Boyce, M. S., Elith, J., Guillera-Arroita, G., Hauenstein, S., Lahoz-Monfort, J. J., Schröder, B., Thuiller, W., et al. (2017).  Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure.  Ecography, 40(8):913–929.
  * Rue et al., (2009) Rue, H., Martino, S., and Chopin, N. (2009).  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2):319–392.
  * Shi et al., (2016) Shi, P., Zhang, A., Li, H., et al. (2016).  Regression analysis for microbiome compositional data.  The Annals of Applied Statistics, 10(2):1019–1040.
  * Simpson et al., (2017) Simpson, D., Rue, H., Riebler, A., Martins, T. G., and Sørbye, S. H. (2017).  Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors.  Statistical Science, 32(1):1–28.
  * Spiegelhalter et al., (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P., and Van Der Linde, A. (2002).  Bayesian measures of model complexity and fit.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4):583–639.
  * Tsilimigras and Fodor, (2016) Tsilimigras, M. C. and Fodor, A. A. (2016).  Compositional data analysis of the microbiome: fundamentals, tools, and challenges.  Annals of Epidemiology, 26(5):330–335.
  * Ugarte et al., (2017) Ugarte, M. D., Adin, A., and Goicoa, T. (2017).  One-dimensional, two-dimensional, and three dimensional B-splines to specify space–time interactions in Bayesian disease mapping: Model fitting and model identifiability.  Spatial Statistics, 22:451–468.
  * Van Niekerk et al., (2023) Van Niekerk, J., Krainski, E., Rustand, D., and Rue, H. (2023).  A new avenue for bayesian inference with inla.  Computational Statistics & Data Analysis, 181:107692.
  * (44) Vicente, G., Goicoa, T., Fernandez-Rasines, P., and Ugarte, M. D. (2020a).  Crime against women in India: unveiling spatial patterns and temporal trends of dowry deaths in the districts of Uttar Pradesh.  Journal of the Royal Statistical Society Series A: Statistics in Society, 183(2):655–679.
  * (45) Vicente, G., Goicoa, T., and Ugarte, M. D. (2020b).  Bayesian inference in multivariate spatio-temporal areal models using INLA: analysis of gender-based violence in small areas.  Stochastic Environmental Research and Risk Assessment, 34:1421–1440.
  * Washburne et al., (2017) Washburne, A. D., Silverman, J. D., Leff, J. W., Bennett, D. J., Darcy, J. L., Mukherjee, S., Fierer, N., and David, L. A. (2017).  Phylogenetic factorization of compositional data yields lineage-level associations in microbiome datasets.  PeerJ, 5:e2969.
  * Watanabe, (2010) Watanabe, S. (2010).  Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory.  Journal of Machine Learning Research, 11(116):3571–3594.

##  A Appendix: additional figures

### Joint modelling of pancreatic cancer mortality and incidence data

Figure A1 shows the maps with the posterior median estimates of the shared
spatial random effects $\exp(\phi_{i})$ (top left) and the posterior
exceedence probabilities $P(\exp(\phi_{i})>1|\boldsymbol{y})$ (top right)
obtained under the M3-TypeII model. Dark blue areas in the maps indicate
significantly high area-specific risks during the period 2001-2020 in
comparison to England as a whole. Posterior median estimates and 95% credible
intervals of the temporal shared random effect $\exp(\gamma_{t})$ is also
plotted (bottom). Figure A2 display the maps with posterior median estimates
of pancreatic cancer mortality (top) and incidence (bottom) risks.

Figure A1: Top: posterior median estimates $\exp(\phi_{i})$ and posterior
exceedence probabilities $P(\exp(\phi_{i})>1|\boldsymbol{y})$ for the shared
spatial random effect (M3). Bottom: Posterior median estimates and 95%
credible intervals of the temporal shared random effect (M3 - TypeII).

Figure A2: Posterior median estimates of pancreatic cancer mortality (top) and
incidence risks (bottom) for the shared-component model (M3) with Type II
interaction effect.

### Spatial Compositional Data: the case of Arabidopsis thaliana

Figure A3 shows the posterior estimates (mean and standard deviation) of the
spatial effects $\boldsymbol{\omega}^{(d)},d=1,2,3$ for each $alr$-coordinate.
We observe how the posterior spatial pattern shifts for each $alr$-coordinate,
a crucial factor for obtaining more flexible predictions.

Figure A3: Posterior mean (A) and standard deviation (B) of the spatial
effects by $alr$-coordinate under Type VIII structure model.

### Space-time models for the United Kingdom wind speed data

Figure A4: Posterior mean of the spatial effect $r(\mathbf{s})$ from the
fitted model B.  Figure A5: Posterior mean of the space-time effect
$u(\mathbf{s},t)$ from the fitted model B.