In previous vignettes, we show how to fit spatial Poisson mixed models for high-dimensional areal count data, how to use parallel or distributed computation strategies, and how to use the **bigDM** package to analyse high-dimensional spatio-temporal count data. Here, we describe how to use this package to fit order-free multivariate scalable Bayesian models to smooth mortality (or incidence) risks of several diseases simultaneously (Vicente et al., 2022).

Let us assume that the area of interest is divided into \(n\) contiguous small areas and data are available for \(J\) diseases. Let \(O_{ij}\) and \(E_{ij}\) denote the number of observed and expected cases respectively in the \(i\)-th small area (\(i=1, \ldots, n\)) and for the \(j\)-th disease (\(j=1, \ldots, J\)). Conditional on the relative risks \(r_{ij}\), the number of observed cases in the \(i\)-th area and the \(j\)-th disease is assumed to follow a Poisson distribution with mean \(\mu_{ij}=E_{ij} \cdot r_{ij}\), that is, \[\begin{eqnarray*} O_{ij}| r_{ij} & \sim & Poisson(\mu_{ij}=E_{ij} \cdot r_{ij}), \\ \log \mu_{ij} & = & \log E_{ij}+\log r_{ij}. \end{eqnarray*}\] Here \(E_{ij}\) is computed using indirect standardization as \(E_{ij}=\sum_{k}n_{ijk}\cdot m_{jk}\), where \(k\) is the age-group, \(n_{ijk}\) is the population at risk in area \(i\) and age-group \(k\) for the \(j\)th disease, and \(m_{jk}\) is the overall mortality (or incidence) rate of the disease \(j\) in the total area of study for the \(k\)-th age group. The log-risk is modelled as \[\begin{equation} \label{log.risk} \log r_{ij}=\alpha_j + \theta_{ij}, \end{equation}\] where \(\alpha_j\) is a disease-specific intercept and \(\theta_{ij}\) is the spatial effect of area \(i\) for the \(j\)-th disease.

Following the work by Botella-Rocamora et al. (2015), we rearrange the spatial effects into the matrix \(\Theta=\lbrace \theta_{ij} \rbrace\), for \(i=1, \ldots, n\); \(j=1, \ldots, J\) to better comprehend the dependence structure. The main advantage of the multivariate modelling is that dependence between the spatial patterns of the different diseases can be included in the model, so that latent associations between diseases can help to discover potential risk factors related to the phenomena under study. These unknown connections can be crucial to a better understanding of complex diseases such as cancer.

The potential association between the spatial patterns of the different diseases are included in the model considering the decomposition of \(\Theta\) as \[\begin{equation} \Theta = \Phi M, \tag{1} \end{equation}\] where \(\Phi\) and \(M\) deal with dependency within and between diseases respectively. We refer to Equation (1) as the M-model. In the following, we briefly describe the two components of the M-model.

The matrix \(\Phi\) is a matrix of order \(n \times K\) and it is composed of stochastically independent columns that are distributed following a spatially correlated distribution. Usually, as many spatial distributions as diseases are considered, that is, \(K=J\), although \(J\) and \(K\) may be different (see Corpas-Burgos et al. (2019), for a discussion). To deal with spatial dependence, different spatial priors have been considered in the literature, most of them based on the well known intrinsic conditional autoregressive (iCAR) prior (Besag et al., 1991). See the vignette bigDM: fitting spatial models for more details.

On the other hand, \(M\) is a \(K \times J\) nonsingular but arbitrary matrix and it is responsible for inducing dependence between the different columns of \(\Theta\), i.e., for inducing corrlation between the spatial patterns of the diseases. In Equation (1), the cells of \(M\) act as coefficients, so they can be considered as coefficients of the log-relative risks on the underlying patterns captured in \(\Phi\) and treated as fixed effects with a normal prior distribution with mean 0 and a large (and fixed) variance. Note that, assigning \(N(0,\sigma)\) priors to the cells of \(M\) is equivalent to assigning a Wishart prior to \(M'M\), i.e., \(M'M \sim Wishart(J, \sigma^{2} \mathbf{I}_J)\). See Botella-Rocamora et al. (2015) for further details.

Once the between-diseases dependencies are incorporated into the model, the resulting prior distributions for \(\mbox{vec} \left( \Theta \right)\) with Gaussian kernel has a precision matrix given by \[\begin{equation} \Omega_{\mbox{vec}(\Theta)} = \left(M^{-1} \otimes I_n \right) \: \mbox{Blockdiag}(\Omega_{1},\ldots,\Omega_{J}) \: \left(M^{-1} \otimes I_n \right)'. \end{equation}\] Recall that this precision matrix accounts for both within and between-disease dependencies. On the one hand, the \(\Omega_{1},\ldots,\Omega_{J}\) matrices are introduced to control the within-diseases spatial variability and the between-diseases variability is captured through the matrix \(M\). Note that if \(\Omega_{1} = \ldots = \Omega_{J}= \Omega_{w}\), the covariance structure is separable and can be expressed as \(\Omega_{\mbox{vec}(\Theta)}^{-1}=\Omega_{b}^{-1} \otimes \Omega_{w}^{-1}\), where \(\Omega_{b}^{-1}=M'M\) and \(\Omega_{w}^{-1}\) are the between- and within-disease covariance matrices, respectively. This M-model based framework includes both separable and non-separable covariance structures, and can accommodate different spatial dependency structures with different within-disease covariance matrices.

Several priors distributions are implemented in the `MCAR_INLA()`

function to deal with spatial dependence within-diseases:

`prior="intrinsic"`

for the M-model implementation of the intrinsic multivariate CAR latent effect.`prior="Leroux"`

for the M-model implementation of the Leroux et al. (1999) multivariate CAR latent effect.`prior="proper"`

for the M-model implementation of the proper multivariate CAR latent effect.`prior="iid"`

for the M-model implementation of spatially non-structured multivariate latent effect.

As for the spatial prior distributions for univariate (single disease) models, appropriate sum-to-zero constraints must be imposed to solve identifiability problems with the disease-specific intercepts. See Vicente et al. (2022) for details about prior distributions for model hyperparameters.

**Note**: The M-model implementation of these models using `R-INLA`

requires the use of at least \(J \times (J+1)/2\) hyperparameters. So, the results must be carefully checked, specially when using the LCAR or pCAR models.

In addition to enlarge the effective sample size and improving smoothing by borrowing information from the different responses, one of the main advantages of multivariate disease mapping models is that they take into account correlations between the spatial patterns of the different diseases \({\rho}=(\rho_{12},\ldots,\rho_{J-1,J})^{'}\), that is, they reveal connections between diseases. In addition, it also provides the diagonal elements of the between-disease covariance matrix (\(\sigma^2_j\)), hereafter referred to as variance parameters. In the case of separable covariance structures (the kronecker product of between and within disease covariance matrices) these parameters control the amount of smoothing within diseases.

We compute the marginal posterior estimates of these parameters by sampling from the approximated joint posterior for the model hyperparameters using the `inla.hyperpar.sample()`

function and computing kernel density estimates of the derived samples for the elements of the correlation matrix of the random effects. The results (summary statistics and posterior marginal densities) are contained in the `summary.cor`

/`summary.var`

and `marginals.cor`

/`marginals.var`

elements of the `inla`

model.

`MCAR_INLA`

functionAs in the `CAR_INLA()`

and `STCAR_INLA()`

functions, three main modelling approaches can be considered:

- the usual model with a global spatial random effect whose dependence structure is based on the whole neighbourhood graph of the areal units (
`model="global"`

argument), - a disjoint model based on a partition of the whole spatial domain where independent spatial
CAR models are simultaneously fitted in each partition (
`model="partition"`

and`k=0`

arguments), - a modelling approach where \(k\)-order neighbours are added to each partition to avoid border effects in the disjoint model (
`model="partition"`

and`k>0`

arguments).

For both the disjoint and \(k\)-order neighbour models, parallel or distributed computation strategies can be performed to speed up computations by using the **future** package (Bengtsson, 2021). See the vignette bigDM: parallel and distributed modelling for some examples and further details.

The data and its associated cartography file need to be specified into the `MCAR_INLA()`

function. These are some of the most relevant arguments of this function:

`carto`

: an object of class`sf`

or`SpatialPolygonsDataFrame`

that must contain at least the target variable of interest specified in the argument`ID.area`

.`data`

: an object of class`data.frame`

that must contain the target variables of interest specified in the arguments`ID.area`

,`ID.disease`

,`O`

and`E`

.`ID.area`

: name of the variable that contains the IDs of spatial areal units. The values of this variable must match those given in the carto and data variable.`ID.disease`

: name of the variable that contains the IDs of the diseases.`ID.group`

: name of the variable that contains the IDs of the spatial partition (grouping variable). Only required if`model="partition"`

.`O`

: name of the variable that contains the observed number of disease cases for each areal and time point.`E`

: name of the variable that contains either the expected number of disease cases or the population at risk for each areal unit and time point.`W`

: optional argument with the binary adjacency matrix of the spatial areal units. If`NULL`

(default), this object is computed from the`carto`

argument (two areas are considered as neighbours if they share a common border).`merge.strategy`

: one of either`"mixture"`

or`"original"`

(default), that specifies the merging strategy to compute posterior marginal estimates of relative risks. See`mergeINLA()`

function for further details.

The `Carto_SpainMUN`

object included in the **bigDM** package, contains the spatial polygons of the municipalities of continental Spain and simulated colorectal cancer mortality data (see the examples of the CAR_INLA function).

```
library(INLA)
#> Warning: package 'INLA' was built under R version 4.2.2
library(bigDM)
data(Carto_SpainMUN)
head(Carto_SpainMUN)
#> Simple feature collection with 6 features and 10 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 485318 ymin: 4727428 xmax: 543317 ymax: 4779153
#> Projected CRS: ETRS89 / UTM zone 30N
#> ID name lat long area
#> 1 01001 Alegria-Dulantzi 4742974 4742974 19913794
#> 2 01002 Amurrio 4765886 4765886 96145595
#> 3 01003 Aramaio 4766225 4766225 73338806
#> 4 01004 Artziniega 4775453 4775453 27506468
#> 5 01006 Arminon 4729884 4729884 10559721
#> 6 01008 Arrazua-Ubarrundia (San Martin de Ania) 4752714 4752714 57502811
#> perimeter obs exp SMR region geometry
#> 1 34372.11 2 3.0237149 0.6614380 Pais Vasco MULTIPOLYGON (((538259 4737...
#> 2 63352.32 28 20.8456682 1.3432047 Pais Vasco MULTIPOLYGON (((503520 4760...
#> 3 41430.46 6 3.7527301 1.5988360 Pais Vasco MULTIPOLYGON (((533286 4759...
#> 4 22605.22 3 3.2093191 0.9347777 Pais Vasco MULTIPOLYGON (((491260 4776...
#> 5 17847.35 0 0.4817391 0.0000000 Pais Vasco MULTIPOLYGON (((509851 4727...
#> 6 64968.81 2 1.9643891 1.0181282 Pais Vasco MULTIPOLYGON (((534678 4746...
```

In this vignette, simulated cancer mortality data for three diseases in the 7907 municipalities of mainland Spain (excluding Baleareas and Canary Islands, and the autonomous cities of Ceuta and Melilla) included in the object `Data_MultiCancer`

will be used for illustration (modified in order to preserve the confidentiality of the original data).

```
data(Data_MultiCancer)
str(Data_MultiCancer)
#> 'data.frame': 23721 obs. of 5 variables:
#> $ ID : chr "01001" "01002" "01003" "01004" ...
#> $ disease: int 1 1 1 1 1 1 1 1 1 1 ...
#> $ obs : int 6 52 8 8 0 5 7 8 4 3 ...
#> $ exp : num 6.615 42.634 7.431 6.355 0.934 ...
#> $ SMR : num 0.907 1.22 1.077 1.259 0 ...
```

Note that both objects contain a common identification variable of the areal units named as `ID`

.

We refer as *Global model* to the spatial multivariate model described in Equation (1), where the whole neighbourhood graph of the areal units is considered to define the adjacency matrix \(\textbf{W}\).

The Global model with an iCAR prior for the spatial random effects is fitted using the `MCAR_INLA()`

function as

```
Global <- MCAR_INLA(carto=Carto_SpainMUN, data=Data_MultiCancer, ID.area="ID", ID.disease="disease",
O="obs", E="exp", prior="intrinsic", model="global", strategy="gaussian")
#> STEP 1: Pre-processing data
#> STEP 2: Fitting global model with INLA (this may take a while...)
summary(Global)
#>
#> Call:
#> c("inla.core(formula = formula, family = family, contrasts = contrasts,
#> ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
#> scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
#> ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
#> verbose, ", " lincomb = lincomb, selection = selection, control.compute
#> = control.compute, ", " control.predictor = control.predictor,
#> control.family = control.family, ", " control.inla = control.inla,
#> control.fixed = control.fixed, ", " control.mode = control.mode,
#> control.expert = control.expert, ", " control.hazard = control.hazard,
#> control.lincomb = control.lincomb, ", " control.update =
#> control.update, control.lp.scale = control.lp.scale, ", "
#> control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
#> ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
#> num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
#> working.directory = working.directory, ", " silent = silent, inla.mode
#> = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
#> .parent.frame)")
#> Time used:
#> Pre = 0.914, Running = 155, Post = 4.77, Total = 161
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> I1 -0.180 0.006 -0.192 -0.180 -0.168 -0.180 0
#> I2 -0.092 0.007 -0.106 -0.092 -0.078 -0.092 0
#> I3 0.127 0.009 0.108 0.127 0.145 0.127 0
#>
#> Random effects:
#> Name Model
#> idx RGeneric2
#>
#> Model hyperparameters:
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> Theta1 for idx -0.719 0.022 -0.762 -0.720 -0.676 -0.720
#> Theta2 for idx -1.332 0.050 -1.427 -1.334 -1.232 -1.338
#> Theta3 for idx -1.493 0.070 -1.635 -1.492 -1.358 -1.487
#> Theta4 for idx 0.289 0.013 0.263 0.289 0.316 0.288
#> Theta5 for idx 0.159 0.015 0.129 0.160 0.189 0.160
#> Theta6 for idx -0.114 0.019 -0.152 -0.114 -0.076 -0.114
#>
#> Deviance Information Criterion (DIC) ...............: 81458.25
#> Deviance Information Criterion (DIC, saturated) ....: 27013.98
#> Effective number of parameters .....................: 2925.78
#>
#> Watanabe-Akaike information criterion (WAIC) ...: 81405.80
#> Effective number of parameters .................: 2403.79
#>
#> Marginal log-Likelihood: -41616.74
#> CPO, PIT is computed
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## Posterior estimates of between-disease correlations ##
Global$summary.cor
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> rho12 0.7380874 0.02461054 0.68659678 0.7390330 0.7839916 0.7419397
#> rho13 0.5321505 0.04197466 0.44762071 0.5328930 0.6108123 0.5359363
#> rho23 0.1369928 0.06039153 0.01855721 0.1377508 0.2560740 0.1422967
## Posterior estimates of variance parameters ##
Global$summary.var
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> var1 0.23744037 0.01078455 0.21705139 0.23712610 0.2593840 0.23652096
#> var2 0.15377438 0.01030375 0.13476698 0.15331540 0.1753713 0.15333300
#> var3 0.08990545 0.00933185 0.07261316 0.08964764 0.1086900 0.08881862
```

When the number of areas is very large, the M-model approach can be computationally very intensive. In this situation, the computational burden of these models is so high that they could be unfeasible for users with limited computing capacity. In addition, to fit a single model in the whole region could be not the best strategy as the degree of smoothing does not need to be same in the whole region. In contrast, the scalable Bayesian multivariate modelling approach described in Vicente et al. (2022) can be used to jointly smooth incidence or mortality risks of several diseases for high-dimensional areal count data. This proposal divides the spatial domain into \(D\) subregions, so that local multivariate spatial models can be fitted simultaneously (using parallel and/or distributed computation strategies) substantially reducing the computational time. In what follows, we show how to fit the **Disjoint** and **k-order neighbourhood models** extending the methodology described in Orozco-Acosta et al. (2021) for estimating spatial multivariate disease risks using the `MCAR_INLA()`

function.

A natural way to think of partitions is to consider subregions based on administrative subdivisions of the area of interest. For our example data in `Data_MultiCancer`

we propose to divide the data into the \(D=15\) Autonomous Regions of Spain (`region`

variable of the `Carto_SpainMUN`

object.

```
library(tmap)
tm_shape(Carto_SpainMUN) +
tm_polygons(col="region") +
tm_layout(legend.outside=TRUE)
```

In the code below, we show how to fit the *Disjoint model* with an iCAR prior for the spatial random effects and Gaussian approximation strategy using 4 local clusters (in parallel)

```
Disjoint <- MCAR_INLA(carto=Carto_SpainMUN, data=Data_MultiCancer,
ID.area="ID", ID.disease="disease", O="obs", E="exp", ID.group="region",
prior="intrinsic", model="partition", k=0, strategy="gaussian",
plan="cluster", workers=rep("localhost",4))
#> STEP 1: Pre-processing data
#> STEP 2: Fitting partition (k=0) model with INLA
#> + Model 1 of 15
#> + Model 2 of 15
#> + Model 3 of 15
#> + Model 4 of 15
#> + Model 5 of 15
#> + Model 6 of 15
#> + Model 7 of 15
#> + Model 8 of 15
#> + Model 9 of 15
#> + Model 10 of 15
#> + Model 11 of 15
#> + Model 12 of 15
#> + Model 13 of 15
#> + Model 14 of 15
#> + Model 15 of 15
#> STEP 3: Merging the results
summary(Disjoint)
#>
#> Call:
#> c("inla.core(formula = formula, family = family, contrasts = contrasts,
#> ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
#> scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
#> ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
#> verbose, ", " lincomb = lincomb, selection = selection, control.compute
#> = control.compute, ", " control.predictor = control.predictor,
#> control.family = control.family, ", " control.inla = control.inla,
#> control.fixed = control.fixed, ", " control.mode = control.mode,
#> control.expert = control.expert, ", " control.hazard = control.hazard,
#> control.lincomb = control.lincomb, ", " control.update =
#> control.update, control.lp.scale = control.lp.scale, ", "
#> control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
#> ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
#> num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
#> working.directory = working.directory, ", " silent = silent, inla.mode
#> = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
#> .parent.frame)")
#> Time used:
#> Running = 110, Merging = 55.2, Total = 165, NA = NA
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> I1 -0.139 0.005 -0.149 -0.139 -0.129 -0.140
#> I2 -0.100 0.006 -0.112 -0.100 -0.088 -0.100
#> I3 0.164 0.008 0.148 0.164 0.179 0.164
#>
#> Random effects:
#> Name Model
#> idx RGeneric2
#>
#> Deviance Information Criterion (DIC) ...............: 81545.45
#> Deviance Information Criterion (DIC, saturated) ....: 27101.18
#> Effective number of parameters .....................: 3198.50
#>
#> Watanabe-Akaike information criterion (WAIC) ...: 81477.35
#> Effective number of parameters .................: 2614.94
#>
#> is computed
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
## Posterior estimates of between-disease correlations ##
Disjoint$summary.cor
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> rho12 0.7474043 0.02391784 0.69910154 0.7479993 0.7923351 0.7509652
#> rho13 0.5523660 0.03909825 0.47345758 0.5534401 0.6266611 0.5543283
#> rho23 0.1678871 0.05452831 0.06102586 0.1679434 0.2732885 0.1613189
## Posterior estimates of variance parameters ##
Disjoint$summary.var
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> var1 0.2250720 0.01094214 0.20398209 0.2248963 0.2472827 0.2241454
#> var2 0.1538919 0.01054461 0.13398111 0.1535961 0.1752698 0.1528272
#> var3 0.1061530 0.00992664 0.08763015 0.1057148 0.1268397 0.1046273
```

** Computations are made in personal computer with a 3.41 GHz Intel Core i5-7500 processor and 32GB RAM using R-INLA stable version INLA_22.05.07.*

The result is an object of class `inla`

where the full domain log-risk is just the union of the posterior marginal estimates of each subregion, i.e., \(\log R =\left( \log R^{(1)'},\cdots,\log R^{(D)'} \right)'\). If the `save.models=TRUE`

argument is included, a list with all the inla submodels is saved in a temporary folder, that can be used as input argument for the `mergeINLA()`

function. In addition, the `compute.intercept=TRUE`

argument can be included to compute disease-specific overall log-risks \(\alpha_j\) by sampling from the joint posterior distribution of the linear predictors. Further details are given below.

As in the case of the scalable spatial and spatio-temporal models fitted with `CAR_INLA()`

and `STCAR_INLA()`

functions, respectively, **k-order neighbourhood models** can be defined to avoid the border effect of considering disjoint partitions. In this case, the entire spatial region \(\mathscr{D}\) is divided into a set of overlapping subregions and some small areas will belong to more than one of such subdivisions, i.e., \(\mathscr{D} = \bigcup_{d=1}^D \mathscr{D}_d\) but \(\mathscr{D}_i \cap \mathscr{D}_j \neq \emptyset\) for neighbouring subregions. As \(\sum_{d=1}^{D} I_{d}>n\), the final risk \(R=(R'_1,\ldots,R'_J)'\) with \(R'_j = (R_{1j},\ldots,R_{Ij})'\), \(j=1,\ldots,J\), is no longer the union of the posterior estimates obtained for each submodel as areas located in the borders of the spatial partition would have more than one estimated posterior distribution.

Two different merging strategies can be considered to obtain a unique posterior estimate of the relative risk for those areas in more than one submodel:

If the

`merge.strategy="mixture"`

argument is specified, mixture distributions of the estimated posterior probability density functions with weights proportional to the conditional predictive ordinates (CPOs) are computed. See Orozco-Acosta et al. (2021) for further details.If the

`merge.strategy="original"`

argument is specified (default option), the posterior marginal distributions of the relative risk estimated from its original partition is used. See Orozco-Acosta et al. (2023).

In the code below, we show how to fit the 1st-order neighbourhood model with an iCAR prior for the spatial random effects and Gaussian approximation strategy using 4 local clusters (in parallel):

```
order1 <- MCAR_INLA(carto=Carto_SpainMUN, data=Data_MultiCancer,
ID.area="ID", ID.disease="disease", O="obs", E="exp", ID.group="region",
prior="intrinsic", model="partition", k=1, strategy="gaussian",
plan="cluster", workers=rep("localhost",4))
#> STEP 1: Pre-processing data
#> STEP 2: Fitting partition (k=1) model with INLA
#> + Model 1 of 15
#> + Model 2 of 15
#> + Model 3 of 15
#> + Model 4 of 15
#> + Model 5 of 15
#> + Model 6 of 15
#> + Model 7 of 15
#> + Model 8 of 15
#> + Model 9 of 15
#> + Model 10 of 15
#> + Model 11 of 15
#> + Model 12 of 15
#> + Model 13 of 15
#> + Model 14 of 15
#> + Model 15 of 15
#> STEP 3: Merging the results
summary(order1)
#>
#> Call:
#> c("inla.core(formula = formula, family = family, contrasts = contrasts,
#> ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
#> scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
#> ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
#> verbose, ", " lincomb = lincomb, selection = selection, control.compute
#> = control.compute, ", " control.predictor = control.predictor,
#> control.family = control.family, ", " control.inla = control.inla,
#> control.fixed = control.fixed, ", " control.mode = control.mode,
#> control.expert = control.expert, ", " control.hazard = control.hazard,
#> control.lincomb = control.lincomb, ", " control.update =
#> control.update, control.lp.scale = control.lp.scale, ", "
#> control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
#> ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
#> num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
#> working.directory = working.directory, ", " silent = silent, inla.mode
#> = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
#> .parent.frame)")
#> Time used:
#> Running = 117, Merging = 74.4, Total = 192, NA = NA
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> I1 -0.132 0.005 -0.141 -0.132 -0.122 -0.132
#> I2 -0.085 0.006 -0.096 -0.085 -0.074 -0.085
#> I3 0.157 0.008 0.142 0.157 0.172 0.157
#>
#> Deviance Information Criterion (DIC) ...............: 81467.13
#> Deviance Information Criterion (DIC, saturated) ....: 27022.87
#> Effective number of parameters .....................: 3029.00
#>
#> Watanabe-Akaike information criterion (WAIC) ...: 81428.97
#> Effective number of parameters .................: 2504.57
#>
#> is computed
## Posterior estimates of between-disease correlations ##
order1$summary.cor
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> rho12 0.7454028 0.02304774 0.69925546 0.7458378 0.7891141 0.7477634
#> rho13 0.5366397 0.03786048 0.46063384 0.5375537 0.6089979 0.5413227
#> rho23 0.1481862 0.05493908 0.03991677 0.1484508 0.2541525 0.1498929
## Posterior estimates of variance parameters ##
order1$summary.var
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> var1 0.2179981 0.010271941 0.1982159 0.2178266 0.2388675 0.2177671
#> var2 0.1533085 0.010105308 0.1342173 0.1529948 0.1737581 0.1513422
#> var3 0.1049042 0.009510871 0.0870917 0.1045005 0.1246553 0.1033483
```

`mergeINLA`

functionThis function takes local models fitted for each subregion of the whole spatial domain and unifies them into a single `inla`

object. It is called by the main function `MCAR_INLA()`

, and is valid for both *Disjoint* and *k-order neighbourhood* models. In addition, approximations to model selection criteria such as the deviance information criterion (DIC) (Spiegelhalter et al., 2002) and Watanabe-Akaike information criterion (WAIC) (Watanabe, 2010) are also computed. See Vicente et al. (2022) for details on how to compute these measures for the multivariate scalable models described in this vignette.

**Computation of the disease-specific overall log-risk**If the

*Disjoint model*is fitted (`k=0`

argument), posterior marginal estimates of the disease-specific overall log-risk \(\alpha_j\) can be computed (`compute.intercept=TRUE`

argument) by sampling from the joint posterior distribution of the linear predictors using the`inla.posterior.sample()`

function of the**R-INLA**package. After joining the \(S\) samples (default to 1000) from each submodel, the following intercepts can be defined\[\begin{equation*} \tilde{\alpha}^s_j = \frac{1}{n}\sum_{i=1}^n \log{r_{ij}}, \quad \mbox{for} \quad s=1,\ldots,S \end{equation*}\] and then it is possible to compute the kernel density estimate of \(\tilde{\alpha_j}\).

**Computation of between-disease correlation and variance parameters**Partition models provide extra information about local relationships between the diseases in the subdivisions, as they provide local estimates of model’s parameters of interest: between-disease correlations \({\rho}=(\rho_{12},\ldots,\rho_{J-1,J})^{'}\), and variance parameters \({\sigma}^2=(\sigma^2_1,\ldots,\sigma^2_J)^{'}\) (diagonal elements of the between-disease covariance matrix).

In the

`mergeINLA()`

function, we implement an adaptation of the consensus Monte Carlo algorithm (Scott et al., 2016) to obtain global estimates of these parameters in the overall study domain from the marginal estimates of the partition models. Further details are given in Vicente et al. (2022).