\section{Introduction}
Massively Multiplayer Online Role-Playing Game (MMORPG) is one of the major genres of online games. Users create their characters and catch monsters, hunt on scenarios designed by the game designer; they perform quests to strengthen items, thereby raising experience levels and thus raising the level. The game provides a variety of contents for frequent user interaction activities such as collecting and producing items, a battle between users (PvP), party play, exchange of items, etc.

Research on the acquisition and maintenance strategy of game users directly related to the profit of the game company has become the most crucial concern in the game industry. Besides, since the cost of acquiring new customers is higher than the cost of maintaining existing customers \cite{b1}, research on securing and managing existing customers has been actively studied in various fields in conjunction with the needs of companies \cite{b2, b3}. Users in the game can lose interest in the game over time and then eventually leave the game. These churn not only lead to a decline in revenues of game companies, but they also harm other users and can lead to their departure \cite{b4}. As a result, game companies try to extend the life cycle of game users by providing ongoing events and rewards.

In this study, we confirm whether churn rate decreases when game companies provide events and promotions through analyzing AION's actual dataset. AION is one of the largest MMORPGs in east Asia, developed and serviced by NCSOFT, one of the most famous game companies in Korea. Besides, we extract a few characteristics such as the user's character level, item purchasing frequency, and game-playing time band among the logs in which are recorded the user's various activities; and we confirm how the events have different effects per extracted characteristics.

\section{Contribution}
\begin{itemize}
\item This paper is the first case to focus on engineering analysis with a business perspective in online game research.
\item We discover the event is effective, especially for old users.
\item We discover that the event affects different influences according to user characteristics. Thus, game companies know which events are the most effective for maximizing their profits.
\end{itemize}

\section{Related works}
User churn issues directly related to the profitability of an online service company are the eternal assignment to be solved in all online services industry. Online games have also made many attempts to solve these problems \cite{b5, b6}. Borbora and Srivastava proposed a method of comparing and analyzing the user’s lifecycle using a distance-based classification system called wClusterDist as a method for distinguishing between churn users and users who continue to play game \cite{b7}, Kawale \textit{et al.} analyzed the relationship between social activity and user churn and discovered that social activities are influencing user churn from the game \cite{b4}. Seo \textit{et al.} proposed a churn prediction model by measuring social activities such as co-play between guild members \cite{b8}. Kim \textit{et al.} presented the study results that the difference in the churn prediction performance varies according to the observation period and the churn prediction period through actual data of three games. According to the results of these researches, it is crucial to set an observation period and the churn prediction period according to game characteristics\cite{b9}. Kim \textit{et al.} analyzed bad effect from bot users in game. They revealed that low-level users (e.g., newbie users) can be easily affected by malicious bot users' game plays \cite{b10}.

These studies can be categorized into survey-based, Open API-based, and server-log based on the data collection method, as shown in Table \ref{tab1}.

\begin{table}[htbp]
\caption{Data collection methods for user churn analysis and related games and research}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\textbf{Category} & \textbf{Related online game} & \textbf{Reference} \\
\hline \hline
Survey & World of Warcraft & \cite{b5} \\
\hline
Open API & BATTLEGROUNDS & \cite{b12} \\
\hline
Server-Log based & AION & \cite{b8, b10, b11, b13, b16}\\
\hline
\end{tabular}
\label{tab1}
\end{center}
\end{table}

The first method, the survey-based method, is used to survey game users to extract the characteristics that want to verify and whether the user leaves the games. Typically, there is research on World of Warcraft game users, and this research investigated how user’s churn rate varies according to game motivation and demographic information. The research results show that users with strong motivation tend to have low churn rates, while those with strong social motivation tend to have higher churn rates\cite{b5}.

The second method, Open-API (Application Programming Interface) is used to collect data provided by game companies. A typical example is the BATTLEGROUND game. There is an advantage in that a lot of log data about a game can be obtained by only registering without having a direct relationship with a game company.

Finally, this is a method of estimating churn through the measure of the user’s behavior and churn based on analysis of stored logs in the game server. This method can have the highest prediction performance as much as it can utilize the most information. Recently, several game companies have provided a sample log to host a contest for churn prediction analysis \cite{b12}.

Moon \textit{et al.} developed a model which can quickly estimate a user's prediction with 90\% of accuracy (4 weeks earlier than real churn occurs) by analyzing social activities such as quest, party-play and trading \cite{b11}.

Lee \textit{et al.} proposed a model that optimizes profit by predicting churn of loyal customers who provided a lot of profits for a long time in AION, it confirmed that this method is more cost-effective than performing the churn prediction from the entire users after applying to the actual game environment \cite{b14}.

Oh \textit{et al.} proposed a churn prediction method based on various user types. They classified churn users into 5 types (light user, light-contents enjoying, heavy-contents enjoying, item making, and hunting with party-play type) by game users' play style. They found 98.3\% of churn users are in the `light user' type \cite{b16}.

\section{Methodology}
In this section, we define the churn criterion used to analyze the event effect and describe the classification method of the user type used in this study, the event reference date which is the basis of the event effect measurement, and the methodology for analyzing the churn change by user type according to an event.
\subsection{Defining churn criterion}
A reasonable churn criterion must be established for accurate churn user classification. Defining appropriate measures that do not overestimate or underestimated churn users is essential to create a practical abstraction classifier model \cite{b15}. In this study, we refer to the method used by Yang \textit{et al.} \cite{b17} and define ‘a user who has no login for N days’ as a churn user. In this experiment, users who are not login continuously from 2 days to 6 days after login based on the start date of analyzed data are defined as churn user. A churn rate is calculated as follows.

Churn rate = (the number of users who have not logged in consecutively for N days) / (the number of cumulative users logged in per day)

\subsection{Classification of user type}
AION provides various analysis features on game logs. Among them, we selected the character’s level, item purchase frequency, and game-playing time band as features which are the most frequently used in churn prediction research as shown in Table \ref{tab2}.

\begin{table}[htbp]
\caption{Feature selection for analyzing user type}
\begin{center}
\begin{tabular}{|p{2.4cm}||p{5.2cm}| }
\hline
\textbf{Feature} & \textbf{Description} \\
\hline \hline
Character’s Level & Highest level of character before and after the event \\
\hline
Item purchase frequency & Number of purchases in store \\
\hline
Game-playing time band & Time band with the largest number of activity logs for two weeks \\
\hline

\end{tabular}
\label{tab2}
\end{center}
\end{table}

The character’s level correlates with game-playing time and loyalty to the game. Well-designed events can satisfy engaged users’ satisfaction, which increases game-playing time. Item purchase frequency is the most critical indicator of the game company’s profits. In general, user’s satisfaction and commitment to the game can increase item purchase. To check the increased number of items purchased after the event is a routine process to measure the event’s success by game companies. Game-playing time band is selected to identify the targeting users of the promotion. For example, we can classify the target user’s group (i.e., hard-core user and light-user) by analyzing game-playing time.

\begin{table}[htbp]
\caption{Criteria for classification of user type}
\begin{center}
\begin{tabular}{|p{2cm}||p{1.7cm}|p{1.7cm}|p{1.7cm}| }
\hline
\textbf{Feature}&\multicolumn{3}{|c|}{\textbf{Classification Criteria}} \\
\cline{2-4}
\hline \hline
Character’s Level & High (61-80) & Mid (31-60) & Low (1-30) \\
\hline
Item purchase frequency & High (Over 10 times) & Mid (6 to 10 times) & Low (0 to 5 times) \\
\hline
Game-playing time band & Dawn (0-6) & Day (6-19) & Night (19-24)\\
\hline
\end{tabular}
\label{tab3}
\end{center}
\end{table}

\subsection{Selection of event reference date}
In the AION game, several types of events are played in succession from two to four weeks. Besides, since more than 60 events are held every year, it is difficult to measure the effects of individual events because there are no independent events to be performed at a specific time.

In this study, the event effects were analyzed by comparing the data before and after on 3 February 2016, when four events occurred simultaneously with the mid-day of the experiment data as shown in Table \ref{tab4}. All event information is retrieved from AION’s official web site \cite{b18}.

\begin{figure*}[]
\includegraphics[width=\textwidth,height=6cm]{daily_churn_rate.png}
  \caption{Churn rate per day during the observation period. It shows the churn rate gradually decreases from the event start time.}
  \label{fig2}
\end{figure*}

\begin{table}[htbp]
\caption{Event reference date using in the experiment}
\begin{center}
\begin{tabular}{|p{2.4cm}|p{1.7cm}|p{1.5cm}|p{1.5cm}| }
\hline
\textbf{Event name}& \textbf{type} & \textbf{Start date of event} & \textbf{Event duration} \\
\hline \hline
Stigma          & Reinforcement  & 3 Feb, 2016 & 14 days \\
\hline
Sullenday      & Ballentine Day & 3 Feb, 2016  & 21 days \\
\hline
Strongest legions & Battle         & 3 Feb, 2016 & 24 days \\
\hline
Dr. Edishung      & Potion making  & 3 Feb, 2016 & 28 days \\
\hline
\end{tabular}
\label{tab4}
\end{center}
\end{table}

Each event's description is summarized as follows.
\begin{itemize}
\item \textbf{Stigma}: Stigma does not disappear even if a user fails to strengthen a stigma in the event period. Stigma usually vanishes into the air if a user fails to reinforce a stigma.

\item \textbf{Sullenday}: A user can purchase a Lucky Box with various items at a low price.

\item \textbf{Strongest legions}: Only level 46 or higher users can attend battles with legion team members, and they can get a reward as much as they win.

\item \textbf{Dr. Edishung}: If users purchase the ticket, they get essential items for level up and get a `40\% discount coupon' which can buy other tickets for 90 days or 900 hours.
\end{itemize}

\subsection{Analysis of churn by user type according to an event}
We propose two-step analysis with Method I and II as shown in Figure\ref{fig1}

\begin{figure}[h]
\centerline{\includegraphics[width=0.47\textwidth]{method.png}}
\caption{This image shows validation method of event effects. First, we verify the event effects by comparing numbers of login users before and after events. Second, we validate the event affects different influence per user type.}
\label{fig1}
\end{figure}

\section{Experiment}
\subsection{Definition of churn criterion}
In this study, we used about 54 million actual log data from AION provided by NCSOFT, the most famous online game companies in Korea. This data contains about 20,000 real user game activity information for two weeks from 27 January 2016 to 9 February 2016. To calculate the percentage of churn users during the two weeks, we extracted all the information about the users who did not login for two to six days in two weeks.

Figure \ref{fig2} shows the percentage of churn users within the observation period. For each event, it heats the peak value one day before the event date, and then churn rates gradually decrease from the event start time. This result indicates that the event is useful to prevent churn in online games.

\begin{figure}[]
\centerline{\includegraphics[width=0.47\textwidth]{comp_login_user.png}}
\caption{Comparison of the number of login users between before and after the event.}
\label{fig3}
\end{figure}

\subsection{Classification of user type}
Table \ref{tab5} summarizes the data classified according to the user type classification by analyzing the log data. First, AION’s maximum level is up to 80. The AION, which passed more than seven years since its first release, shows that more than 50 percent of users have a level of 61 or higher. Second, we confirmed that the number of items purchased per user is less than five times per week. Third, we know that about 64\% of gameplay users are playing games during the daytime.

\begin{table}[htbp]
\caption{This table shows user distribution by classes}
\begin{center}
\begin{tabular}{|p{2cm}||p{1.7cm}|p{1.7cm}|p{1.7cm}| }
\hline
\textbf{Feature}&\multicolumn{3}{|c|}{\textbf{User Distribution by Classes}} \\
\cline{2-4}
\hline
Character’s Level & High 50\% & Mid 45\% & Low 5\% \\
\hline
Item purchase frequency & High 13.5\%  & Mid 14.5\% & Low 72\% \\
\hline
Game-playing time band & Dawn 18\% & Day 64\% & Night 18\%\\
\hline
\end{tabular}
\label{tab5}
\end{center}
\end{table}

\subsection{Comparison of the number of active users between before and after the event}
Figure \ref{fig3} shows the number of active users before and after the event. Before the event, the average number of active users per week is 7,500, and the number of cumulative users is 12,737. After the event, the average number of active users 8,860 and the number of cumulative users is 18,159, which is increased by 18\% and 19\% respectively.

This result shows that the event is an effective way to increase the number of active users. However, it is difficult to confirm whether the increased user is an existing user, a new user, or a returning user.

\subsection{Distribution of user-level before and after the event}
Figure \ref{fig4} shows the level distribution of users before and after the event in quartiles. Before the event, the average level of login users was 60, but the average level after the event was lowered to 58, and the lowest level user was also lowered from 44 to 22 after the event. This graph shows that the event has a more significant impact on users at lower levels than average.

\begin{figure}[h]
\centerline{\includegraphics[width=0.47\textwidth]{dist.png}}
\caption{Distribution of user-level before and after the event; The blue bar shows the distribution of user-level before the event while the orange bar shows the distribution of user-level after the event. We know that the login of low-level users increases after starting the event.}
\label{fig4}
\end{figure}

\subsection{Event effect analysis by character’s level class}
Figure \ref{fig5} shows the analysis of the event effect by a character’s level class. At the event, the number of logins for the entire class level user has increased, especially at level 31 to 60, with middle-class user connections are growing to about 44\%.

\begin{figure}[h]
\centerline{\includegraphics[width=0.47\textwidth]{class.png}}
\caption{Event effect analysis by character’s level class; It shows mid-level user’s login is the most increased}
\label{fig5}
\end{figure}

\subsection{Event effect analysis by game-playing time bands}
Figure \ref{fig6} shows the effect of the time-based event on the game. It shows that dawn and daytime users increased while nighttime users decreased.

\subsection{Event effect analysis by item purchase frequency}
Figure \ref{fig7} shows the change in the item purchase frequency before and after the event. Contrary to expectations, purchase times after the event decrease by about 10\%, while the number of buyers increased by nearly 150\% due to event gifts.
The number of users who buy game items more than ten times per week reduces significantly while the number of purchases under five times per week increases drastically.

If AION game logs are including the detailed purchase log, it is possible to check the effect of the event in monetary units. However, since the analyzed logs do not include the purchase item price, only the number of simple purchases and the number of purchase users were analyzed.

Figure \ref{fig5} shows that mid-level user’s login is the most increased, Figure \ref{fig6} shows that the event influences on daytime users relatively, and Figure \ref{fig7} shows that the number of purchase decreases. If we analyze Figure \ref{fig5}, \ref{fig6} and \ref{fig7} collectively, we can find that users who play games during the daytime have a lot of mid-level and do not use the money for buying game items.
If game companies use the result of analyzing various event effects by the user's characteristics, the most profitable, effective marketing strategy can be formulated.

\begin{figure}[h]
\centerline{\includegraphics[width=0.47\textwidth]{timezone.png}}
\caption{Event effect analysis by game-playing time band; It shows that the event influence on daytime users relatively.}
\label{fig6}
\end{figure}

\begin{figure}[h]
\centerline{\includegraphics[width=0.47\textwidth]{buycount.png}}
\caption{Purchase pattern change before and after the event; It shows the number of purchases over 10 times decrease.}
\label{fig7}
\end{figure}

\subsection{Event effect analysis by game-playing time bands}
Figure \ref{fig6} shows the effect of the time-based event on the game. It shows that dawn and daytime users increased while nighttime users decreased.

\section{Conclusion and Future work}
In this study, we used the log data of AION, which is currently in service to check that game events are useful to churn prevention and how long time takes to generate effects to users in MMORPG. Besides, we confirmed that event effects are different depending on the character’s level, item purchase frequency and game-playing time band.

The most unusual aspect of the analysis is that the number of purchase activities in the game decreased even though the event has the effect of preventing the user churn and increasing the number of return users. This phenomenon means that additional analysis is needed internally for this fact in the game company to check whether the event results in substantial monetary gain.

We find the increase of new customers is minimal while the event increases the number of mid-level users who usually play in the daytime through the experiment. This finding shows that the event is more effective in preventing churn users rather than introducing new users.

In this study, we find that the users who purchased more than ten times during the event period and who enjoyed playing in the nighttime decreases while the users who purchased less than five times and who enjoyed playing in the daytime increases. This analysis shows that the game companies need to focus on events for the daytime user to retain users while focusing on events for the night user to increase revenues.

These results of these analyses are very instrumental and practical not only to game engineer but to game marketers. Thus, this study can be effectively used for marketing strategies and plans that can lead to an increase in return users as well as decreasing churn users.

Furthermore, the event effect analysis proposed in this study will be able to have a significant synergy effect through the target marketing if it is analyzed together with the classification of the user characteristics type in the game conducted in other studies.

In the future, we will verify the universality of the proposed method by testing with the extended period of AION dataset. Then, we will research the similarity and differences of event effects by analyzing other genre of games such as FPS (First Person Shooter games), arcade and so forth. Finally, we will study the how events in PC and mobile games have different effects.


\section{Introduction}
\label{sec:intro}

The integrated nested Laplace approximation \citep[INLA,][]{INLA} provides an
alternative to traditional Markov chain Monte Carlo
\citep[MCMC,][]{Gilksetal:1996} for Bayesian inference. The INLA methodology
focuses on estimating the posterior marginals of the model parameters instead
of their joint posterior distribution. INLA is implemented in the \pkg{INLA}
package for the \proglang{R} programming language, that provides a simple way to
fit models via the \code{inla()} function, which works in a similar way as
other functions to fit regression models such as \code{glm()} or \code{gam()}.

The \pkg{INLA} package implements several likelihoods, priors and latent
effects that can be used to build models. It is also capable of fitting models
with several likelihoods, which can be useful for multivariate modeling.
However, multivariate spatial models are not included. The \pkg{INLAMSM}
package adds a number of multivariate latent effects that implement well-known
multivariate spatial models for areal data. By fitting these models with INLA,
instead of MCMC, computing times should be reduced.

Lattice data are made of  $I$ areas, usually related to some
administrative boundaries, where data are collected. We will assume
that values of $K$ variables are obtained from each area. In addition,
an adjacency structure is defined by considering that areas with a shared
boundary are neighbors. The analysis of this type of datasets is often
made by resorting to multivariate regression models.

This paper is organized as follows. The remainder of this introduction includes
a review of current software for multivariate spatial modeling.
Section~\ref{sec:models} gives a description of the different multivariate
spatial models implemented in the \pkg{INLAMSM} package.   Two detailed
applications of the \pkg{INLAMSM} package and the corresponding results can be
founded in Section~\ref{sec:examples}. Finally, a summary and discussion
are provided in Section~\ref{sec:discussion}.

\medskip
The \proglang{R} programming language provides a number of standalone packages for
multivariate spatial analysis.  In the specific case of analyzing spatial point
patterns, some \proglang{R} packages are available, such as \pkg{spatstat}
\citep{spatstat} and \pkg{spatialkernel} \citep{spatialkernel}. Package
\pkg{spatstat} is able to model multitype point patterns as well as handling a
good deal of the models and methods for the analysis of spatial and
spatio-temporal point patterns such as estimators of the space-time
inhomogeneous K-function and pair correlation function.  Package
\pkg{spatialkernel} performs edge-corrected kernel density estimation and
binary kernel regression estimation for multivariate spatial point patterns.

\medskip
Regarding geostatistical data, \pkg{gstat} \citep{GSTAT,gstatRpackage} provides
functions to fit both univariate and multivariate models using different
types of kriging. Package \pkg{spBayes} \citep{spBayes} is able to fit a wide
variety of Gaussian spatial process models for univariate as well as
multivariate point-referenced data using efficient MCMC algorithms.

\medskip
Finally, in the case of areal data, \proglang{R} package \pkg{CARBayes}
\citep{CARBayes2013} offers the possibility of fitting a wide class of CAR
models using MCMC methods. This package provides a number of univariate and
multivariate likelihoods, and it also includes a multivariate CAR model with an
inverse-Wishart and a CAR prior proposed by \citet{Lerouxetal:2000} to
estimate the variability between the different variables and the spatial
variation, respectively.

\medskip
In addition to these packages, the \proglang{BUGS} language
\citep{BUGS,WINBUGS} is a flexible framework for the implementation of
multivariate spatial models for lattice data. Package \pkg{R2WinBUGS}
\citep{R2WINBUGS} provides a set of functions to call the \pkg{WinBUGS}
software from \proglang{R}.  Recently, \proglang{Stan} \citep{Stan:2017}
develops a flexible language for Bayesian inference that could also be used to
fit (multivariate) spatial models \citep{Morrisetal:2019}. The \pkg{Rstan}
package \citep{Rstan} is a convenient interface between \proglang{R} and
\proglang{Stan} to fit models with ease.  Both \proglang{BUGS} and
\proglang{Stan} rely on MCMC algorithms for model fitting and inference.

\medskip
A computationally efficient alternative to software based on MCMC algorithms
is described in \cite{INLA}. This has been implemented in an \proglang{R} package
called \pkg{INLA}, which is often referred to as \pkg{R-INLA} to distinguish it
from the main INLA methodology. INLA advantages include an easy way of
implementing hierarchical models and short computing times to
fit spatial or spatio-temporal models
\citep{INLABayesianSpatial,INLAspatialreview}. This software is designed for
Bayesian inference on latent Gaussian Markov random field models which
include (generalized) linear mixed and spatial and spatio-temporal models.
\pkg{INLA} can deal with lattice data as well as geostatistical data by means of
an approximation to continuous spatial processes based on stochastic partial
differential equations \cite[SPDE,][]{SPDE}.  This approach can be used to fit
log-Gaussian Cox processes \citep{Simpsonetal:2016} for spatial and spatio-temporal point patterns as
well \citep{SPDE}.  Multivariate models can be fit by considering different
likelihoods with shared terms. See, for example, \citet{SPDE},
\citet{INLAVirgilioAdvanced} and \citet{INLAbook:2020}.

\medskip
Recently, some \proglang{R} packages have been developed on top of \pkg{R-INLA}.
For instance, \pkg{inlabru} \citep{inlabru} uses \pkg{INLA}
to fit spatial density surfaces and estimating
abundance in a spatial point process, gridded and georeferenced context.

\medskip
As our review has shown, there are only a few packages that can fit
multivariate spatial models in a Bayesian framework. Furthermore, in the case
of analyzing multivariate areal data using Bayesian inference, there are
options based on MCMC and INLA. INLA is a widely used alternative to MCMC and
often provides shorter computing times when fitting spatial models.  However,
there are not a set of functions to fit multivariate spatial models for lattice
data using \pkg{INLA}. Developing these latent effects is the main goal of this
paper, which will enable users to include multivariate spatial effects in their
models to consider different ways of measuring spatial cross-correlation
\citep{SainCressie:2002}.

Some of the models implemented in the \pkg{INLAMSM} package have been described
by several authors \citep{GelfandVounatsou:2003,CarlinBanerjee:2003}.
Furthermore, other authors
\citep{Jinetal:2007,MartinezBeneito:2013,MacNab:2018} provide different reviews
of multivariate spatial models that describe the models implemented and many
others that could be added to the package in the future.  Furthermore,
\citet{bookMSM:2019} provide a recent review of multivariate spatial models in
the context of disease mapping.

\section{Models}
\label{sec:models}

\pkg{INLAMSM} implements different functions which correspond to different
multivariate spatial latent effects. These differ in structure, complexity and
number of hyperparameters.  However, all of them are defined in a multivariate
spatial context in which $i = 1, \ldots, I$ represents the spatial areas and $k
= 1, \ldots, K$ is used to index the variables measured in region $i$.

Random effects can be represented using a matrix $\Theta$ with entries
$\theta_{ik},\ i=1,\ldots,I,\ k=1,\ldots,K$. Hence, the $k$-th column of
$\Theta$ represent the spatial random effects associated to variable $k$.

A particular application of these multivariate spatial models is disease
mapping, as described in the examples in Section~\ref{sec:examples}.  In this
context, the number of the observed cases of $k$-th disease in the $i$-th
spatial area, $Y_{ik}$, is modeled as a Poisson random variable:

$$
Y_{ik} \sim Po(E_{ik} \cdot R_{ik})
$$
\noindent
where $E_{ik}$ and $R_{ik}$ represent the expected cases and the relative risk
for the $i$-th spatial area and the $k$-th disease, respectively. Then, the
logarithm of the relative risk is the sum of two terms:

$$
\log(R_{ik})= a_k + \theta_{ik}
$$
\noindent
where $a_k$ is the intercept of the $k$-th disease and $\theta_{ik}$ is the term
that models the spatial variability. Note that other covariates can be
included in the linear predictor on the right hand side of the previous
equation.

\medskip
Following \citet{MartinezBeneito:2013}, the variability of a multivariate
spatial latent effect can be divided in two terms: variability between the
values of the variables measured in the same area and the spatial variability
corresponding to the values of different areas for a particular variable. Both
sources of variability are modeled with their respective variance matrices, that have
different number of hyperparameters.

\medskip
Finally, before explaining any model, some useful specifications are defined.
Let $\Theta_{i\cdot}$ ($i=1, ...,I$) denote the $i$-th row of matrix
$\Theta$ and $\Theta_{\cdot j}$ ($j=1, ...,J$) denote the $j$th column of
matrix $\Theta$. In addition, on a matrix $A=[A_{\cdot 1}: \cdots : A_{\cdot
J}] $, operator $vec(\cdot)$ is defined as

$$
vec(A) = (A_{\cdot 1}^{\top}, \ldots, A_{\cdot J}^{\top})^{\top} .
$$

Note that all latent effects will be defined using the \code{rgeneric} latent
effect, which defines the latent effect using the representation of the Gaussian Markov random
field. This is described in the \pkg{INLA} documentation, that can be accessed
with \code{inla.doc("rgeneric")}, and also in Chapter 11 of
\citet{INLAbook:2020}.  This representation is essentially a multivariate
Gaussian vector with a sparse precision matrix.  Hence, $vec(\Theta)$ will be
defined using a Normal distribution with zero mean and a highly structured
variance matrix $\Sigma$, i.e.,

$$
vec(\Theta) \sim N(0, \Sigma) .
$$

Note that \pkg{INLA} works with the precision and that this is why
$\Sigma^{-1}$ will be required instead of the variance matrix $\Sigma$.
Furthermore, $\Sigma$ will be created so that it accounts for the association
between the different variables within areas as well as association between
areas for the same variable or disease, as explained above.  The new latent
effects included in the \pkg{INLAMSM} can be used as a guidance to develop
other types of multivariate spatial latent effects.

\subsection{Independent intrinsic MCAR}

One of the easiest options to model the variability within- and
between-variables is using an intrinsic CAR distribution \citep{Besagetal:1991}
for the spatial structure and a diagonal matrix for the covariance matrix
between variables, respectively. Specifically, $\Theta$ is modeled as

$$
vec(\Theta) \sim N\left(0, \Lambda^{-1} \otimes (D-W)^{-1}\right)
$$
\noindent
where $D=diag(n_1,...,n_I)$ is a diagonal matrix with values $n_i$  (the
number of neighbors of region $i$) and $W$ is an adjacency matrix with entries
$W_{ij}$. Each entry $W_{ij}$ is equal to 1 if units $i$ and $j$ are neighbors
and 0 otherwise.

Let $\tau_k$ be the marginal precision of the $k$-th variable. When variables
are independent of each other, the between-variables precision matrix $\Lambda$
is a $k \times k$ diagonal matrix defined as

\[
  \Lambda=
  \left[ {\begin{array}{cccccc}
   \tau_1 &   0    & \cdots & \cdots &    0   \\
   0      & \tau_2 &   0    &        & \vdots \\
   \vdots &   0    & \ddots & \ddots & \vdots \\
   \vdots &        & \ddots & \ddots &    0   \\
   0      & \cdots & \cdots &    0   &  \tau_K\\
  \end{array} } \right]
.
\]

Following \citet{Jinetal:2007}, matrix $\Lambda$ can be regarded as a
(precision) matrix for the association between variables within the same area
and (precision) matrix $(D-W)^{-1}$ to account for the association
between areas for a given variable.  In this case it is assumed that the
different variables are independent and this is why $\Lambda$ is a diagonal
matrix, but in other latent random effects introduced below $\Lambda$ will be a
dense matrix to model association between the different variables.

Thus, this model has $K$ hyperparameters, $ \{\tau_k \}^{K}_{k=1}$, equal to
the total number of variables. This model is computationally fast
because the number of hyperparameters is the lowest among the multivariate
spatial models implemented in the \pkg{INLAMSM} package, as it will be seen
below.  Therefore, this model can be used  as a baseline to compare a
na\"ive model with independent spatial terms with more complex
alternatives.

Hyperparameters are internally represented in \pkg{INLA} using the log-scale
so that they are not bounded, and the vector
of hyperparameters is $(\log(\tau_1), \ldots, \log(\tau_K))$.  Prior
distributions are assigned to the standard deviations instead of the precisions
\citep{Gelman:2006}. In particular, a uniform improper distribution between $0$
and $+\infty$ is assigned to each standard deviation:

$$
\sigma_k = \frac{1}{\sqrt{\tau_k}} \sim Un(0, +\infty ); \,\,\,\,\,\, k=1, \ldots, K.
$$

Note that \pkg{INLA} will report the results in the internal scale but we have
included several functions in the package to transform from the log-precision
into some other more meaningful scale. See Section~\ref{sec:examples} for
details. Furthermore, \pkg{INLA} provides a set of functions that can be
used to make this transformation \citep[see, for example,][Chapter 2]{INLAbook:2020}.

Finally, given that this is an improper distribution, it will be
necessary to add further constraints when fitting the model, similarly to the
one dimensional intrinsic CAR distribution \citep{BesagKooperberg:1995}. In this
case, a sum-to-zero constraint will be added to the random effects associated
to each variable.  This implies that a variable-specific intercept should be
added to the model as well.

\subsection{Independent proper MCAR}

A proper CAR distribution can be used to model the within-variables variability instead
of the intrinsic version. Here, a common spatial autocorrelation parameter
$\alpha$ for all the variables is introduced. Thus $\Theta$ is modeled as

$$
vec(\Theta) \sim N\left(0, \Lambda^{-1} \otimes (D- \alpha \cdot W)^{-1}\right)
$$
\noindent
where matrices $D$ and $W$ are defined similarly as above.  The matrix used to
model the between-variables variability keeps the same structure, thus
$\Lambda$ is defined as a diagonal matrix with entries $(\tau_1, \ldots,
\tau_K)$.

The proper CAR specification described above is not implemented in the
\pkg{INLA} package.  In our opinion this is an additional contribution of our
\pkg{INLAMSM} package. It is worth noting that the way in which the proper CAR
distribution is implemented here differs from the proper latent random effects
implemented in the \pkg{INLA} package as defined in latent effects
\code{besagproper} and \code{besagproper2}. The \code{besagproper} effect is
based on adding a constant to the diagonal entries of the precision matrix of
the intrinsic CAR distribution, while the \code{besagproper2} effect is
actually the model proposed by \cite{Lerouxetal:2000}.

Now, this model has $K+1$ hyperparameters which are the spatial autocorrelation
parameter, $\alpha$, and the $K$ precisions, $ \{\tau_k \}^{K}_{k=1}$.
Internally, the hyperparameters are $(\alpha^{*},\log(\tau_1), \ldots,
\log(\tau_K)).$ Hyperparameter $\alpha^{*}$ is defined by transforming $\alpha$
as follows:

$$
\alpha^{*} = logit(\frac{\alpha - \alpha_{min}}{\alpha_{max} - \alpha_{min}})
$$
\noindent
where $\alpha_{min}$ and $\alpha_{max}$ define the bounds of the domain of
$\alpha$ \citep{Sunetal:1999}. However, the proper CAR distribution has been
aduced to show weird artifacts when the spatial correlation parameter takes
negative values \citep{Wall:2004}. For this reason, it is common to assume
$\alpha_{min}=0$ and $\alpha_{max}=1$.

As in the previous model, uniform prior distributions are set on the standard
deviations, and a uniform prior on $\alpha$ is considered:

$$
\sigma_k =
\frac{1}{\sqrt{\tau_k}} \sim Un(0, +\infty); \,\,\,\,\,\, k=1, \ldots, K ,
$$

$$
\alpha \sim Un(\alpha_{min}, \alpha_{max}) .
$$

Once again, the independent proper MCAR distribution is a na\"ive
implementation of independent spatial patterns as the independent intrinsic
MCAR. Nevertheless, this model could also be used as a benchmark for comparing
multivariate proper CAR models assuming dependence between variables (as
described below).

\subsection{Improper MCAR model}

A diagonal precision matrix $\Lambda^{-1}$ is a na\"ive way to model the
between-variables variability which leads to short computing times. However,
setting the off-diagonal elements of $\Lambda$ to zero assumes that variables
are independent of each other.  The off-diagonal elements model dependence
between every pair of variables, which could be useful when looking for similar
spatial patterns of different variables.

Therefore, a dense precision matrix and an intrinsic conditional autoregressive
distribution can be chosen to model the between- and within-variables
variability, respectively. This can be regarded as a generalization of the
univariate intrinsic conditional autoregressive model. In particular, $\Theta$
is modeled similarly as the independent IMCAR:

$$
vec(\Theta) \sim N\left(0, \Lambda^{-1} \otimes (D-W)^{-1}\right)
$$
\noindent
but now $\Lambda$ is a dense, symmetric and positive-definite matrix. The
parameterization of $\Lambda$ follows that of other latent effects implemented
in \pkg{INLA} and the hyperparameters are the precisions of the variables and
the correlations between any pair of variables. Hence, instead of dealing
with the structure of $\Lambda$, $\Lambda^{-1}$ is defined as follows:

\[
  \Lambda^{-1}=
  \left[ {\begin{array}{cccccc}
             1/\tau_1          & \rho_{12}/\sqrt{\tau_1\tau_2} &           \cdots              &               \cdots                  &      \rho_{1K}/\sqrt{\tau_1\tau_K}    \\
 \rho_{12}/\sqrt{\tau_1\tau_2} &            1/\tau_2           & \rho_{23}/\sqrt{\tau_2\tau_3} &                                       &                  \vdots               \\
              \vdots           & \rho_{23}/\sqrt{\tau_2\tau_3} &           \ddots              &               \ddots                  &                  \vdots               \\
              \vdots           &                               &           \ddots              &               \ddots                  & \rho_{(K-1)K}/\sqrt{\tau_{K-1}\tau_K} \\
 \rho_{1K}/\sqrt{\tau_1\tau_K} &            \cdots             &           \cdots              & \rho_{(K-1)K}/\sqrt{\tau_{K-1}\tau_K} &                 1/\tau_K              \\
  \end{array} } \right]
\]

\noindent

Here, $\rho_{jk}$ is the correlation coefficient between variables $j$ and $k$
and $\tau_k$ is the marginal precision of the $k$-th variable. Therefore, the
set of hyperparameters contains $K$ marginal precisions and $K*(K-1)/2$
correlation parameters. In this case, a joint prior distribution is
considered for $\Lambda^{-1}$ instead of setting a prior distribution
for each hyperparameter. Thus, $\Lambda^{-1}$ follows a Wishart distribution as

$$\Lambda^{-1} \sim Wishart_K(r, R^{-1})
$$
\noindent
where $r$ is the number of degrees of freedom and $R^{-1}$ is a fixed symmetric
positive definite matrix of size $K \times K$. In our implementation, $r$ is
equal to $K$ \citep[following ][]{CarlinBanerjee:2003} and $R$ is the $K\times
K$ identity matrix.

The vector of hyperparameters contains now the precisions plus the correlations
in the lower triangular matrix of $\Lambda^{-1}$ (columnwise). However, these
parameters are re-scaled so that the log-precisions are used and the correlation
parameters are transformed using

$$
\rho^{*} = logit((\rho + 1) / 2),
$$
\noindent
for any correlation $\rho$. Hence, the vector of hyperparameters in the
internal scale is now defined as $(\log(\tau_1), \ldots, \log(\tau_K), \rho^{*}_{21}, \rho^{*}_{31}, \ldots)$.

As stated above, given that this distribution is improper, a sum-to-zero
constraint on the effects for each variable will be added when fitting the
model. This also means that variable-specific intercepts are required in the
model.

\subsection{Proper MCAR model}

A proper CAR distribution can be used instead of the intrinsic CAR in order to
model the within-variables variability. This alternative corresponds to the
multivariate generalization of the univariate proper conditional autoregressive
model. As in the previous case, a dense $\Lambda$ matrix is used to model the
between-variables variability. Specifically, $\Theta$ is modeled as

$$
vec(\Theta) \sim N\left(0, \Lambda^{-1} \otimes (D-\alpha \cdot W)^{-1}\right)
$$
\noindent
where $\alpha$ is a common spatial autocorrelation parameter and $\Lambda$ is a
dense symmetric positive-definite matrix with $\Lambda^{-1}$ defined as in the previous model.

In this case, the set of hyperparameters comprehends one common spatial
autocorrelation parameter $\alpha$, $K$ marginal precisions, $\{\tau_k
\}^{K}_{k=1}$, and $K*(K-1)/2$ correlation parameters, $
\{\rho_{jk}\}^{K}_{j,k=1}$.  Again, a Wishart distribution is
considered as a joint prior distribution for the $\Lambda^{-1}$ matrix, while a
uniform prior is consider for $\alpha$, i.e.,

$$
\alpha \sim Un(\alpha_{min}, \alpha_{max}) .
$$

In this case, the vector of hyperparameters contains the spatial
autocorrelation parameter $\alpha$, the precisions and the correlations in the
lower triangular elements of $\Lambda^{-1}$ (columnwise). As in previous
models, all hyperparameters are transformed to take values in the $(-\infty,
+\infty)$ interval.  Hence, the vector of hyperparameters is $(\alpha^{*},
\log(\tau_1),\ldots,\log(\tau_K), \rho^{*}_{21}, \rho^{*}_{31}, \ldots)$.

\subsection{M-model}

\citet{BotellaRocamoraetal:2015} describe a unifying modeling framework for
multivariate disease mapping when the number of diseases is potentially large.
Here, spatial effects are linear combinations of proper CAR spatial effects. In
particular, we will consider $K$ underlying proper CAR spatial effects defined by

$$ \phi_k \sim N\left(0, (D - \alpha_k W )^{-1}\right),\ k=1,\ldots,K.
$$
\noindent
Here, $\phi_k$ is a vector of length $I$.

The value of the spatial random effect $\Theta_{\cdot j}$ for variable $j$ is
defined as

$$
\Theta_{\cdot j} = \phi_1 m_{1j} + \ldots \phi_J m_{Jj} .
$$
\noindent
Hence, matrix $M$ with entries $m_{ij}$ defines the loadings of
the different underlying CAR spatial effects for each disease or variable.

The distribution of these random effects is given by

$$
vec(\Theta) \sim N\left(0, (M^{\top} \otimes I)
  diag((\Sigma_w)_1, \ldots, (\Sigma_w)_K) (M \otimes I)\right) .
$$
\noindent
Here, matrices $(\Sigma_w)_k$ are the variance matrices of the $K$
underlying proper CAR spatial effects, i.e.,

$$
(\Sigma_w)_k = (D - \alpha_k W)^{-1},\ k=1,\ldots,K.
$$

In addition, \citet{BotellaRocamoraetal:2015} show that for the separable case, with $\alpha_1 = \ldots = \alpha_K$, the between-variables
variance matrix is $M^{\top} M$. The prior of this model is
on $M^{\top} M$, and it follows a Wishart with parameters $K$ and
$\tau I$. Parameter $\tau$ is a fixed precision which is set to
$0.001$, but it could also be considered as another hyperparameter to be
estimated \citep{CorpasBurgosetal:2019}.

The vector of hyperparameters in this model is made of the $K$ autocorrelation
parameters (conveniently transformed) followed by the columns of matrix
$M$, for which no transformation is required.

\subsection{Summary of the models}
\label{sec:summary}

To summarize the different models implemented and their associated functions,
Table~\ref{tab:models} displays some basic information about them.
Note how all functions names follow a similar pattern for consistency,
and that the latent models differ in the number of hyperparameters. These
depend on the number of variables $K$. Furthermore, there are two different
functions to define the new latent effects: the one that implements the
latent effects using the \code{rgeneric} framework and a wrapper function.
The wrapper function is intended to provide a simpler way to define
the latent effect and avoid calling function \code{inla.rgeneric.define}
\citep[see, for example, Chapter 11 in][]{INLAbook:2020}. This will be further
discussed in the examples in Section~\ref{sec:examples} below.

\begin{table} [h]
\centering
\begin{scriptsize}
\begin{tabular}{|c|c|c|c|}
\hline
Latent effect & Wrapper function & \code{rgeneric} function & \# hyperparameters \\
\hline
  Independent IMCAR  &  \code{inla.INDIMCAR.model} & \code{inla.rgeneric.indep.IMCAR.model}  &  $K$   \\
  Independent PMCAR  &  \code{inla.INDMCAR.model} & \code{inla.rgeneric.indep.MCAR.model}   &  $K+1$   \\
      IMCAR     &  \code{inla.IMCAR.model} &    \code{inla.rgeneric.IMCAR.model}      &  $K*(K+1)/2$   \\
      PMCAR     &  \code{inla.MCAR.model} &    \code{inla.rgeneric.MCAR.model}     &  $K*(K+1)/2 + 1$   \\
M-model         & \code{inla.Mmodel.model} & \code{inla.rgeneric.Mmodel.model}         &  $K * K$ \\
\hline
\end{tabular}
\end{scriptsize}
\caption{Summary of the latent effects implemented in the \pkg{INLAMSM}
 package.}
\label{tab:models}
\end{table}

\section{Examples}
\label{sec:examples}

In this section, two examples are developed with the \pkg{INLAMSM} package. The
first one is on the well-known North Carolina sudden infant death syndrome (SIDS) data
\citep{CressieRead:1985}, which is used to show how the different models
implemented in the \pkg{INLAMSM} package are fit. This example is also intended
to show that multivariate spatial models can not only be used to model several
diseases but also the same disease across different time periods so that
dependence across different periods can be estimated. The second example is
based on the study of three causes of death in Comunidad Valenciana (Spain), a
region that comprises 540 municipalities and provides a more challenging
dataset to fit spatial models than the North Carolina SIDS data (with only 100
areas).  In this second example, the focus in on investigating dependence among
the different diseases.

\subsection{North Carolina SIDS data}

In the first example, the North Carolina SIDS data \citep{CressieRead:1985}
have been considered to test the methods implemented in the \pkg{INLAMSM}
package.  This dataset includes counts of number of live births and number of
deaths from sudden infant death syndrome for the 100 counties of North
Carolina for two time periods: from July 1st, 1974 to Jun 30th, 1978 and from
July 1st, 1979 to June 30th, 1984.  This dataset is available in a shapefile in
\proglang{R} package \pkg{spData} \citep{spData}, which will be loaded using
package \pkg{rgdal} \citep{rgdal}.

\begin{Schunk}
\begin{Sinput}
R> library("rgdal")
R> #Load SIDS data
R> nc.sids <- readOGR(system.file("shapes/sids.shp", package = "spData")[1],
+    verbose = FALSE)
R> proj4string(nc.sids) <- CRS("+proj=longlat +ellps=clrk66")
\end{Sinput}
\end{Schunk}

Next, the adjacency matrix of the 100 counties in North Carolina
is obtained with function \code{poly2nb} in the \pkg{spdep} package
 \citep{AppliedVirgilio,spdep} and converted into a sparse matrix of type \code{Matrix}
\citep{Matrix}.

\begin{Schunk}
\begin{Sinput}
R> library("spdep")
R> #Compute adjacency matrix, as nb object 'adj' and sparse matrix 'W'
R> adj <- poly2nb(nc.sids)
R> W <- as(nb2mat(adj, style = "B"), "Matrix")
\end{Sinput}
\end{Schunk}

The model that will be fit in this case is

$$
O_{ik} \sim Po(\mu_{ik});\
\log(\mu_{ik}) = \log(E_{ik}) + a_k + \theta_{ik},\ i=1,\ldots,100,\ k=1,2
$$
\noindent
Here, $O_{ik}$ is the number of SIDS cases in county $i$ and time period $k$,
$E_{ik}$ the expected counts, $a_k$ a period-specific intercept and
$\theta_{ik}$ the multivariate spatial effect, which is defined using one of
the models implemented in the \pkg{INLAMSM} package.
Note that other covariates and effects could be included in the linear
predictor as well.

The expected number of cases $E_{ik}$ are computed by multiplying the
period-specific mortality rate $r_k$ by the number of births $N_{ik}$:

$$
E_{ik} = r_k N_{ik};\ r_k = \frac{\sum_{i=1}^{100} O_{ik}}{\sum_{i=1}^{100} N_{ik}},\ i=1,\ldots,I,\ k=1,2 .
$$

In the next lines of \proglang{R} code
the expected counts for both time periods are computed, as well as the
standardized mortality ratio (SMR), $O_{ik} / E_{ik}$, and the proportion of
non-white births, which could be used as a covariate for both time periods.

\begin{Schunk}
\begin{Sinput}
R> # First time period
R> # Compute expected cases
R> r74 <- sum(nc.sids$SID74) / sum(nc.sids$BIR74)
R> nc.sids$EXP74 <- r74 * nc.sids$BIR74
R> # SMR
R> nc.sids$SMR74 <- nc.sids$SID74 / nc.sids$EXP74
R> # Proportion of non-white births
R> nc.sids$NWPROP74 <- nc.sids$NWBIR74 / nc.sids$BIR74
R> # Second time period
R> # Compute expected cases
R> r79 <- sum(nc.sids$SID79) / sum(nc.sids$BIR79)
R> nc.sids$EXP79 <- r79 * nc.sids$BIR79
R> # SMR
R> nc.sids$SMR79 <- nc.sids$SID79 / nc.sids$EXP79
R> # Proportion of non-white births
R> nc.sids$NWPROP79 <- nc.sids$NWBIR79 / nc.sids$BIR79
\end{Sinput}
\end{Schunk}

In order to prepare the data to fit the models with \pkg{INLA}, a new object is
created by appending the data from the first time period to the data from the
second one. In addition, to the observed, expected and proportion of non-white
births, and index is created to identify the counties and time periods.

\begin{Schunk}
\begin{Sinput}
R> d <- data.frame(OBS = c(nc.sids$SID74, nc.sids$SID79),
+    PERIOD = c(rep("74", 100), rep("79", 100)),
+    NWPROP = c(nc.sids$NWPROP74, nc.sids$NWPROP79),
+    EXP = c(nc.sids$EXP74, nc.sids$EXP79))
R> # County-period index
R> d$idx <- 1:length(d$OBS)
\end{Sinput}
\end{Schunk}

Now, we will fit the models. This requires a two step process:

\begin{enumerate}

\item Define the multivariate spatial latent effect using one of the functions
in the \pkg{INLAMSM} package \citep{INLAMSM}. This is done via the
\code{inla.rgeneric.define} function, that will take the function of the
multivariate spatial model to be included plus any other required arguments,
such as the adjacency matrix, the number of time periods, etc.  Alternatively,
the associated wrapper function can be used (see below).

\item Fit the model with \pkg{INLA} \citep{INLApkg} using a formula that includes
the newly defined latent effect.

\end{enumerate}

For example, for the independent IMCAR model the latent effect will be
defined as follows:

\begin{Schunk}
\begin{Sinput}
R> library("INLAMSM")
R> library("INLA")
R> # Number of variables (i.e., periods)
R> k <- 2
R> # Define bivariate latent effect
R> model.indimcar <- inla.rgeneric.define(inla.rgeneric.indep.IMCAR.model,
+    list(k = k, W = W))
\end{Sinput}
\end{Schunk}

The previous definition of the latent effect uses a call to function
\code{inla.rgeneric.define}, which takes function
\code{inla.rgeneric.indep.IMCAR.model} that defines the new latent effect plus
any other required arguments (i.e., a named list with the number of diseases
\code{k} and adjacency matrix \code{W}).  However, as explained in
Section~\ref{sec:summary} there are wrapper functions associated to each new
latent effect, that will be used from now on. For this particular case, the
definition of the new latent effect using a wrapper function is:

\begin{Schunk}
\begin{Sinput}
R> model.indimcar <- inla.INDIMCAR.model(k = k, W = W)
\end{Sinput}
\end{Schunk}

As explained in Section~\ref{sec:models}, given that this model includes
intrinsic CAR specifications, it is necessary to include some constraints so
that the random effects for each disease sum up to zero. Constraints in
\pkg{INLA} are added using equation $\mathbf{A} \mathbf{x}^{\top} =
\mathbf{e}^{\top}$, where matrix $\mathbf{A}$ and vector $\mathbf{e}$ are used
to define linear constraints on the (multivariate) vector of random effects
$\mathbf{x}$. These two elements are defined below:

\begin{Schunk}
\begin{Sinput}
R> A <- kronecker(Diagonal(k, 1), Matrix(1, ncol = nrow(W), nrow = 1))
R> e  = rep(0, k)
\end{Sinput}
\end{Schunk}

\noindent
The constraints are specified using argument \code{extraconstr} in the call to
the \code{f} function that defines the latent effect (see below).

The next step fits the model with \pkg{INLA}:

\begin{Schunk}
\begin{Sinput}
R> #Fit model
R> IIMCAR <- inla(OBS ~ 0 + PERIOD + f(idx, model = model.indimcar,
+      extraconstr = list(A = as.matrix(A), e = e)),
+    data = d,
+    E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}

Note that linear constraints have been added so that random effects associated
to each disease sum up to zero, as explained above.  Furthermore, argument
\code{control.predictor} has been set to compute the marginals of the linear
predictor (i.e., by setting \code{compute} equal to \code{TRUE}) and argument
\code{control.compute} has been set to compute two model selection criteria. In
particular, the deviance information criterion
\citep[DIC]{spiegelhalteretal:2002} and Watanabe-Akaike information criterion
\citep[WAIC]{Watanabe:2013} will be computed. See below for a discussion on how
to use these criteria for model choice.

This model can be summarized as follows:

\begin{Schunk}
\begin{Sinput}
R> summary(IIMCAR)
\end{Sinput}
\begin{Soutput}
Call:
   c("inla(formula = OBS ~ 0 + PERIOD + f(idx, model =
   model.indimcar, ", " extraconstr = list(A = as.matrix(A), e =
   e)), family = \"poisson\", ", " data = d, E = EXP,
   control.compute = list(dic = TRUE, waic = TRUE), ", "
   control.predictor = list(compute = TRUE))")
Time used:
    Pre = 1.21, Running = 3.51, Post = 0.1, Total = 4.82
Fixed effects:
           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
PERIOD74 -0.071 0.056     -0.184   -0.070      0.036 -0.067   0
PERIOD79 -0.017 0.046     -0.110   -0.016      0.071 -0.015   0

Random effects:
  Name	  Model
    idx RGeneric2

Model hyperparameters:
                mean    sd 0.025quant 0.5quant 0.975quant  mode
Theta1 for idx 0.774 0.339      0.129    0.764       1.46 0.731
Theta2 for idx 1.429 0.408      0.664    1.414       2.27 1.361

Expected number of effective parameters(stdev): 65.92(8.94)
Number of equivalent replicates : 3.03

Deviance Information Criterion (DIC) ...............: 909.59
Deviance Information Criterion (DIC, saturated) ....: 271.91
Effective number of parameters .....................: 67.35

Watanabe-Akaike information criterion (WAIC) ...: 910.58
Effective number of parameters .................: 53.61

Marginal log-Likelihood:  -484.51
Posterior marginals for the linear predictor and
 the fitted values are computed
\end{Soutput}
\end{Schunk}

Note that hyperparameters \code{Theta1} and \code{Theta2} correspond to the
log-precisions of the latent effect, i.e., the two independent ICAR spatial
effects. Because the two effects are independent of each other, this
model is equivalent to fitting two models using an intrinsic CAR model
(implemented as the \code{besag} latent effect in \pkg{INLA}).

Next, the fitted values of the relative risks of the IIMCAR model can be
compared with the SMR for both time periods as seen in Figure~\ref{fig:ncsids}.
Package \pkg{RColorBrewer} \citep{RColorBrewer} has been used to set the colors
in the palette. Note how the fitted values show smoothed spatial patterns as
compared with the raw SMR values.

\begin{figure}
\centering
\includegraphics{Palmi-PeralesF-INLAMSM-011}
\caption{Standardized mortality ratio and posterior means of the fitted values with the IIMCAR model for both time periods for the North Carolina SIDS data.}
\label{fig:ncsids}
\end{figure}

The other models implemented in the \pkg{INLAMSM} package can be fit in a
similar way, as listed below. For example, the independent proper MCAR model
requires the values of $\alpha_{min}$ and $\alpha_{max}$ to be passed in
addition to \code{k} and \code{W} when defining the latent effect. This is
exemplified below. We have constrained $\alpha$ to take just positive values
since the PCAR distribution has been said to have a counterintuitive
performance when this parameter takes negative values \citep{Wall:2004}.
Anyway, our implementation of the PCAR distribution admits, if wanted, the
whole range of admissible values for this parameter \citep{Sunetal:1999}.

\begin{Schunk}
\begin{Sinput}
R> # Independent Proper MCAR model
R> # Define range for the autocorrelation parameter
R> alpha.min <- 0
R> alpha.max <- 1
R> model.indmcar <- inla.INDMCAR.model(k = k, W = W, alpha.min = alpha.min,
+    alpha.max = alpha.max)
R> #Fit model
R> IPMCAR <- inla(OBS ~ 0 + PERIOD + f(idx, model = model.indmcar), data = d,
+    E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}

The remainder of the models in the \pkg{INLAMSM} package are defined and
fit as:

\begin{Schunk}
\begin{Sinput}
R> # IMCAR model
R> model.imcar <- inla.IMCAR.model(k = k, W = W)
R> #Fit model
R> IMCAR <- inla(OBS ~ 0 + PERIOD + f(idx, model = model.imcar,
+      extraconstr = list(A = as.matrix(A), e = e)),
+    data = d,
+    E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}

As discussed above, linear constraints have been added when fitting this model
too.  This way, the intercepts in the model can be identified as well.

\begin{Schunk}
\begin{Sinput}
R> # Proper MCAR model
R> model.mcar <- inla.MCAR.model(k = k, W = W, alpha.min = alpha.min,
+    alpha.max = alpha.max)
R> #Fit model
R> PMCAR <- inla(OBS ~ 0 + PERIOD + f(idx, model = model.mcar), data = d,
+    E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}

\begin{Schunk}
\begin{Sinput}
R> # M-model
R> model.mmodel <- inla.Mmodel.model(k = k, W = W, alpha.min = alpha.min,
+    alpha.max = alpha.max)
R> # Fit model
R> Mmodel <- inla(OBS ~ 0 + PERIOD + f(idx, model = model.mmodel), data = d,
+    E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}

Once we have fit all models, we can compare point estimates of the relative
risks. Posterior means of the relative risks from different models are
displayed in Figure~\ref{fig:ncsidsmsm}.

\begin{Schunk}
\begin{Sinput}
R> # Add results to nc.sids
R> nc.sids$IIMCAR74 <- IIMCAR$summary.fitted[1:100, "mean"]
R> nc.sids$IIMCAR79 <- IIMCAR$summary.fitted[100 + 1:100, "mean"]
R> nc.sids$IPMCAR74 <- IPMCAR$summary.fitted[1:100, "mean"]
R> nc.sids$IPMCAR79 <- IPMCAR$summary.fitted[100 + 1:100, "mean"]
R> nc.sids$IMCAR74 <- IMCAR$summary.fitted[1:100, "mean"]
R> nc.sids$IMCAR79 <- IMCAR$summary.fitted[100 + 1:100, "mean"]
R> nc.sids$PMCAR74 <- PMCAR$summary.fitted[1:100, "mean"]
R> nc.sids$PMCAR79 <- PMCAR$summary.fitted[100 + 1:100, "mean"]
R> nc.sids$Mmodel74 <- Mmodel$summary.fitted[1:100, "mean"]
R> nc.sids$Mmodel79 <- Mmodel$summary.fitted[100 + 1:100, "mean"]
\end{Sinput}
\end{Schunk}

\begin{figure}
\centering
\includegraphics{Palmi-PeralesF-INLAMSM-017}
\caption{Posterior means of the relative risks for the different multivariate spatial effects fit to the North Carolina SIDS data.}
\label{fig:ncsidsmsm}
\end{figure}

\subsubsection{Including covariates in the model}

The proportion of non-white births shows a very similar pattern to that of the
relative risk. These proportions for the time periods are shown in
Figure~\ref{fig:nwprop}, which can be compared with the estimates of the
relative risk in Figure~\ref{fig:ncsids} to appreciate the similar patterns.
Note the similar patterns between the proportion of non-white births and the
standardized mortality ratio for both time periods.  For this reason, several
authors \citep[see, for example,][]{CressieRead:1985} have mentioned the
importance of including this covariate in the model. Hence, the same models are
fitted now including the covariate so that a different coefficient is estimated
for each time period.

\begin{figure}
\centering
\includegraphics{Palmi-PeralesF-INLAMSM-018}
\caption{Proportion of non-white births per county for both time periods.}
\label{fig:nwprop}
\end{figure}

In order to have the covariate in the model so that a different coefficient is
estimated  for each time period, the covariate must be structured in a
two-column as follows:

\begin{Schunk}
\begin{Sinput}
R> # Number of areas
R> n <- nrow(W)
R> NWPROP <- matrix(NA, ncol = 2, nrow = 2 * n)
R> NWPROP[1:n, 1] <- nc.sids$NWPROP74
R> NWPROP[n + 1:n, 2] <- nc.sids$NWPROP79
\end{Sinput}
\end{Schunk}

Next, data used to fit the previous models are converted into a list so that
the covariate can be added:

\begin{Schunk}
\begin{Sinput}
R> d <- as.list(d)
R> d$NWPROP <- NWPROP
\end{Sinput}
\end{Schunk}

Finally, models are fit again:

\begin{Schunk}
\begin{Sinput}
R> # Independent ICAR model
R> IIMCAR2 <- inla(OBS ~ 0 + PERIOD + NWPROP + f(idx, model = model.indimcar,
+      extraconstr = list(A = as.matrix(A), e = e)),
+    data = d,
+    E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
R> # Independent proper CAR model
R> IPMCAR2 <- inla(OBS ~ 0 + PERIOD + NWPROP + f(idx, model = model.indmcar),
+    data = d,
+    E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
R> # Instrinsic MCAR model
R> IMCAR2 <- inla(OBS ~ 0 + PERIOD + NWPROP + f(idx, model = model.imcar,
+      extraconstr = list(A = as.matrix(A), e = e)),
+    data = d,
+    E = EXP, family = "poisson", control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
R> # Proper MCAR model
R> PMCAR2 <- inla(OBS ~ 0 + PERIOD + NWPROP + f(idx, model = model.mcar),
+    data = d, E = EXP, family = "poisson",
+    control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}
\begin{Schunk}
\begin{Sinput}
R> # M-model
R> Mmodel2 <- inla(OBS ~ 0 + PERIOD + NWPROP + f(idx, model = model.mmodel),
+    data = d, E = EXP, family = "poisson",
+    control.predictor = list(compute = TRUE),
+    control.compute = list(dic = TRUE, waic = TRUE))
\end{Sinput}
\end{Schunk}

Note that, in a similar manner that sum-to-zero constraints have been imposed
for intrinsic CAR distributions, further constraints could be also imposed
above assuming the spatial random effects to be orthogonal to the
covariate. In this manner, confounding between the covariate and spatial random
effects could be alleviated \citep{HodgesReich:2010}.

In order to show the estimates of the coefficients of the covariates, the
summary of the independent IMCAR model with covariates is shown below:

\begin{Schunk}
\begin{Sinput}
R> summary(IIMCAR2)
\end{Sinput}
\begin{Soutput}
Call:
   c("inla(formula = OBS ~ 0 + PERIOD + NWPROP + f(idx, model =
   model.indimcar, ", " extraconstr = list(A = as.matrix(A), e =
   e)), family = \"poisson\", ", " data = d, E = EXP,
   control.compute = list(dic = TRUE, waic = TRUE), ", "
   control.predictor = list(compute = TRUE))")
Time used:
    Pre = 1.37, Running = 4.22, Post = 0.0786, Total = 5.67
Fixed effects:
           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
PERIOD74 -0.687 0.130     -0.949   -0.684     -0.438 -0.679   0
PERIOD79 -0.211 0.124     -0.458   -0.210      0.031 -0.208   0
NWPROP1   1.978 0.351      1.297    1.974      2.681  1.966   0
NWPROP2   0.622 0.365     -0.103    0.623      1.336  0.627   0

Random effects:
  Name	  Model
    idx RGeneric2

Model hyperparameters:
               mean    sd 0.025quant 0.5quant 0.975quant mode
Theta1 for idx 1.91 0.587      0.867     1.87       3.16 1.72
Theta2 for idx 1.54 0.437      0.733     1.53       2.44 1.46

Expected number of effective parameters(stdev): 49.28(9.58)
Number of equivalent replicates : 4.06

Deviance Information Criterion (DIC) ...............: 899.49
Deviance Information Criterion (DIC, saturated) ....: 261.81
Effective number of parameters .....................: 51.51

Watanabe-Akaike information criterion (WAIC) ...: 904.67
Effective number of parameters .................: 46.14

Marginal log-Likelihood:  -479.10
Posterior marginals for the linear predictor and
 the fitted values are computed
\end{Soutput}
\end{Schunk}

These results show the positive association between the relative risk and the
proportion of non-white births. Other models have similar estimates of the
coefficients of this  covariate. In addition, Figure~\ref{fig:ncsidsmsm2} shows
the posterior means of the relative risks estimated with these models.  They
can be compared to the ones in Figure~\ref{fig:ncsidsmsm} to notice that the
estimates are usually very similar.

\begin{figure}
\centering
\includegraphics{Palmi-PeralesF-INLAMSM-027}
\caption{Posterior means of the relative risks from the different multivariate
spatial effects fit to the North Carolina SIDS data. These models include the
proportion of non-white births as a covariate.}
\label{fig:ncsidsmsm2}
\end{figure}

\subsubsection{Model selection}

Given the large number of models fit to the data it is necessary to select the
best model among all of them. The deviance information criterion
\citep[DIC]{spiegelhalteretal:2002} and Watanabe-Akaike information criterion
\citep[WAIC]{Watanabe:2013} have been computed when the models were fit and
they are shown in Table~\ref{tab:NCSIDSIC}.

\begin{table}[h!]
\centering
\begin{tabular}{c|c|c|c|c}
Model & \multicolumn{2}{c}{DIC} & \multicolumn{2}{|c}{WAIC} \\
\hline
 & \multicolumn{2}{c}{Covariate} & \multicolumn{2}{|c}{Covariate} \\
\cline{2-5}
 & No & Yes & No & Yes\\
\cline{2-5}
\hline
Independent IMCAR &
909.59 & 899.49 &
910.58 & 904.67 \\
Independent PMCAR &
910.59 & 897.87 &
903.75 & 895.61 \\
IMCAR &
904.26 & 894.96 &
905.96 & 897.38 \\
PMCAR &
905.47 & 893.8 &
898.34 & 888.76 \\
M-model &
911.47 & 893.3 &
917.71 & 898.17 \\
\end{tabular}
\caption{Information criteria computed for model selection.}
\label{tab:NCSIDSIC}
\end{table}

According to the values obtained of the DIC and WAIC, it is clear that
the covariate should be included in the model because it produces a decrease
in the DIC and WAIC of about 10 for all models. Secondly, models with
independent effects
in general do not perform as well as the other models. This means that
there is some within-area association between the different variables that
needs to be accounted for.

\subsection{Mortality in Comunidad Valenciana (Spain)}

The next example is based on simulated data of the mortality by cirrhosis, lung
and oral cancer in Comunidad Valenciana (Spain).  This dataset has been
obtained from \citet{bookMSM:2019} and it has been generated to mimic the
spatial pattern of the real data, that cannot be provided due to
confidentiality constraints. The original files are available at
\texttt{http://github.com/MigueBeneito/DisMapBook}.

Here, the number of deaths by these three causes are available at the
municipality level in Comunidad Valencia (Spain), as well as the expected
number of cases that have been computed using internal standardization.
Hence, the aim now is to estimate the spatial pattern of the different
diseases as well as their possible correlations.

Given that in the previous example we have already described how to fit
all the new models in the \proglang{R} package, we will focus now on the models
that include a term to model the covariance for several diseases. That is,
only models IMCAR, PMCAR and M-model will be fit.

Comunidad Valenciana data are available in the package as two separate objects:
\code{CV} and \code{CV.nb}. \code{CV} is a \code{SpatialPolygonsDataFrame}
\citep{AppliedVirgilio} that contains the data and boundaries of the
municipalities in Comunidad Valenciana, while \code{CV.nb} is an \code{nb}
object \citep{spdep} with the neighborhood structure. Both objects can be
loaded as:

\begin{Schunk}
\begin{Sinput}
R> #Load data
R> data(CV)
\end{Sinput}
\end{Schunk}

The adjacency matrix is computed using a sparse representation \citep{Matrix}
as follows:

\begin{Schunk}
\begin{Sinput}
R> #Compute sparse adjacency matrix W
R> W <- as(nb2mat(CV.nb, style = "B"), "Matrix")
\end{Sinput}
\end{Schunk}

Next, a \code{data.frame} is created with the observed and expected data,
together with an index variable to be passed to the definition of the latent
effects and a disease-specific intercept:

\begin{Schunk}
\begin{Sinput}
R> #Data
R> d <- data.frame(OBS = c(CV$Obs.Cirrhosis, CV$Obs.Lung, CV$Obs.Oral),
+    EXP = c(CV$Exp.Cirrhosis, CV$Exp.Lung, CV$Exp.Oral)
+  )
R> # Add disease specific intercept
R> d$Intercept <- rep(c("Cirrhosis", "Lung", "Oral"), each = nrow(W))
R> # Index for latent effect
R> d$idx <- 1:length(d$OBS)
\end{Sinput}
\end{Schunk}

In order to fit the different models, the following parameters
are defined:

\begin{Schunk}
\begin{Sinput}
R> #Number of diseases
R> k <- 3
R> # Range of autocorrelation parameter
R> alpha.min <- 0
R> alpha.max <- 1
\end{Sinput}
\end{Schunk}

Similarly to the SIDS example, some of the multivariate spatial
effects require additional constraints. The matrices required
are defined below:

\begin{Schunk}
\begin{Sinput}
R> A <- kronecker(Diagonal(k, 1), Matrix(1, ncol = nrow(W), nrow = 1))
R> e = rep(0, k)
\end{Sinput}
\end{Schunk}

Then, the latent effects are defined and the models are fit:

\begin{Schunk}
\begin{Sinput}
R> # Define latent IMCAR model
R> model <- inla.IMCAR.model(k = k, W = W)
R> # FIT IMCAR model
R> IMCAR.cval <- inla(OBS ~ 0 + Intercept + f(idx, model = model,
+      extraconstr = list(A = as.matrix(A), e = e)),
+    data = d, E = EXP, family = "poisson",
+    control.compute = list(config = TRUE, dic = TRUE, waic = TRUE),
+    # Increase starting values of marginal precisions for model fitting
+    control.mode = list(theta = c(1, 1, 1, 0, 0, 0), restart = TRUE),
+    # Required to obtain more accurate gradients for model fitting
+    control.inla = list(h = 0.001),
+    control.predictor = list(compute = TRUE))
R> # Define latent PMCAR model
R> model <- inla.MCAR.model(k = k, W = W, alpha.min = alpha.min,
+    alpha.max = alpha.max)
R> # Fit PMCAR model
R> PMCAR.cval <-  inla(OBS ~ 0 + Intercept + f(idx, model = model), data = d,
+    E = EXP, family = "poisson",
+    control.compute = list(config = TRUE, dic = TRUE, waic = TRUE),
+    control.predictor = list(compute = TRUE))
R> # Define latent M-model
R> model <- inla.Mmodel.model(k = k, W = W, alpha.min = alpha.min,
+    alpha.max = alpha.max)
R> # Fit M-model
R> Mmodel.cval <- inla(OBS ~ 0 + Intercept + f(idx, model = model), data = d,
+    E = EXP, family = "poisson",
+    control.compute = list(config = TRUE, dic = TRUE, waic = TRUE),
+    control.predictor = list(compute = TRUE)
+  )
\end{Sinput}
\end{Schunk}

Table~\ref{tab:times} shows the computing times for the models fit to both
examples. All models have been run on a Mac OS X computer with an Intel Core i5
processor (2,7 GHz), 4 cores and 16GB of RAM. Now the models take longer to run
than in the previous example because there are three diseases and about 5 times
more areas (as Comunidad Valenciana has 540 municipalities in the cartography
that we have used). \citet{BotellaRocamoraetal:2015} report times of about 16
minutes to fit the M-model with \pkg{WinBUGS} on this same dataset.
Hence, \pkg{INLA} can fit the same models in a fraction of the
time.

In addition, Table~\ref{tab:times} also shows the computing times of
models fit to the Comunidad Valenciana dataset using two diseases (cirrhosis
and lung cancer). These models have been computed as follows:

\begin{Schunk}
\begin{Sinput}
R> # New dataset
R> d2 <- subset(d,  Intercept != "Oral")
R> d2$idx <- 1:nrow(d2)
R> k <- 2
R> # Linear constraint for the new models
R> A <- kronecker(Diagonal(k, 1), Matrix(1, ncol = nrow(W), nrow = 1))
R> e = rep(0, k)
R> model <- inla.IMCAR.model(k = k, W = W)
R> # FIT IMCAR model
R> IMCAR.cval2 <- inla(OBS ~ 0 + Intercept + f(idx, model = model,
+      extraconstr = list(A = as.matrix(A), e = e)),
+    data = d2, E = EXP, family = "poisson",
+    control.compute = list(config = TRUE, dic = TRUE, waic = TRUE),
+    control.predictor = list(compute = TRUE))
R> # Define latent PMCAR model
R> model <- inla.MCAR.model(k = k, W = W, alpha.min = alpha.min,
+    alpha.max = alpha.max)
R> # Fit PMCAR model
R> PMCAR.cval2 <-  inla(OBS ~ 0 + Intercept + f(idx, model = model),
+    data = d2, E = EXP, family = "poisson",
+    control.compute = list(config = TRUE, dic = TRUE, waic = TRUE),
+    control.predictor = list(compute = TRUE))
R> # Define latent M-model
R> model <- inla.Mmodel.model(k = k, W = W, alpha.min = alpha.min,
+    alpha.max = alpha.max)
R> # Fit M-model
R> Mmodel.cval2 <- inla(OBS ~ 0 + Intercept + f(idx, model = model),
+    data = d2, E = EXP, family = "poisson",
+    control.compute = list(config = TRUE, dic = TRUE, waic = TRUE),
+    control.predictor = list(compute = TRUE)
+  )
\end{Sinput}
\end{Schunk}

\begin{table}[h]
\begin{scriptsize}
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
 &  &  & \multicolumn{5}{c|}{Model} \\
\cline{4-8}
Dataset & \# diseases & \# Areas & Ind. IMCAR & Ind. PMCAR & IMCAR & PMCAR & M-model\\
\hline
NC SIDS &  2 & 100 &
  4.82 &
  3.72 &
  8.88 &
  9.19 &
  15.64\\
C. Valenciana &  2 & 540 & -- & -- &
  10.5 &
  16.92 &
  62.41\\
C. Valenciana &  3 & 540 & -- & -- &
  36.51 &
  45.67 &
  319.46\\
\hline
\end{tabular}
\end{scriptsize}
\caption{Computing times (in seconds) for the models fit to the two examples.}
\label{tab:times}
\end{table}

The fitted values (i.e., posterior means of the relative risks) can be added to
the map in the \code{SpatialPolygonsDataFrame} object as follows:

\begin{Schunk}
\begin{Sinput}
R> n <- nrow(W)
R> CV$IMCAR.CIR <- IMCAR.cval$summary.fitted[1:n, "mean"]
R> CV$IMCAR.LUN <- IMCAR.cval$summary.fitted[n + 1:n, "mean"]
R> CV$IMCAR.ORA <- IMCAR.cval$summary.fitted[2 * n + 1:n, "mean"]
R> CV$PMCAR.CIR <- PMCAR.cval$summary.fitted[1:n, "mean"]
R> CV$PMCAR.LUN <- PMCAR.cval$summary.fitted[n + 1:n, "mean"]
R> CV$PMCAR.ORA <- PMCAR.cval$summary.fitted[2 * n + 1:n, "mean"]
R> CV$Mmodel.CIR <- Mmodel.cval$summary.fitted[1:n, "mean"]
R> CV$Mmodel.LUN <- Mmodel.cval$summary.fitted[n + 1:n, "mean"]
R> CV$Mmodel.ORA <- Mmodel.cval$summary.fitted[2 * n + 1:n, "mean"]
\end{Sinput}
\end{Schunk}

\begin{figure}
\centering
\includegraphics{Palmi-PeralesF-INLAMSM-036}
\caption{Posterior means of the relative risks of cirrhosis, lung cancer and oral cancer in Comunidad Valenciana (Spain).}
\label{fig:maps}
\end{figure}

The maps in Figure~\ref{fig:maps} show the different posterior means of the
relative risks from the models fitted for the different causes of death. In
general, the three models produce similar point estimates of the relative
risks.

The \pkg{INLAMSM} package includes a few functions to transform the marginals
and summary statistics of the  model hyperparameters in the internal scale into
the original scale in the model. For the models fitted, this
back-transformation  can be obtained as follows:

\begin{Schunk}
\begin{Sinput}
R> hyper.imcar <- inla.MCAR.transform(IMCAR.cval, 3)
R> hyper.pmcar <- inla.MCAR.transform(PMCAR.cval, 3, model = "PMCAR",
+    alpha.min = alpha.min, alpha.max = alpha.max)
R> hyper.mmodel <- inla.Mmodel.transform(Mmodel.cval, 3,
+    alpha.min = alpha.min, alpha.max = alpha.max)
\end{Sinput}
\end{Schunk}

These will provide a transformation of the model parameters so that they are
not in the internal scale anymore and make inference easier.  Spatial
autocorrelation parameters are transformed to be in the range between
$\alpha_{min}$ and $\alpha_{max}$, and log-precisions are transformed to be
variances. Correlation hyperparameters are transformed to be between $-1$ and
$1$.

Hence, summary statistics for the models are:

\begin{Schunk}
\begin{Sinput}
R> # IMCAR
R> hyper.imcar$summary.hyperpar
\end{Sinput}
\begin{Soutput}
                mean     sd quant0.025 quant0.25 quant0.5 quant0.75
Theta1 for idx 0.243 0.0396     0.1752     0.214    0.239     0.267
Theta2 for idx 0.127 0.0162     0.0971     0.115    0.126     0.137
Theta3 for idx 0.226 0.0461     0.1490     0.193    0.222     0.254
Theta4 for idx 0.469 0.0829     0.2976     0.415    0.472     0.527
Theta5 for idx 0.563 0.0964     0.3542     0.502    0.571     0.632
Theta6 for idx 0.705 0.0633     0.5657     0.666    0.711     0.751
               quant0.975
Theta1 for idx      0.330
Theta2 for idx      0.160
Theta3 for idx      0.329
Theta4 for idx      0.622
Theta5 for idx      0.730
Theta6 for idx      0.813
\end{Soutput}
\begin{Sinput}
R> # PMCAR
R> hyper.pmcar$summary.hyperpar
\end{Sinput}
\begin{Soutput}
                mean      sd quant0.025 quant0.25 quant0.5 quant0.75
Theta1 for idx 0.989 0.00634      0.973     0.986    0.991     0.994
Theta2 for idx 0.262 0.04355      0.190     0.231    0.257     0.288
Theta3 for idx 0.139 0.01804      0.106     0.126    0.138     0.150
Theta4 for idx 0.254 0.05207      0.169     0.217    0.248     0.285
Theta5 for idx 0.458 0.08174      0.291     0.403    0.460     0.515
Theta6 for idx 0.545 0.09598      0.340     0.483    0.551     0.613
Theta7 for idx 0.688 0.06531      0.547     0.646    0.692     0.734
               quant0.975
Theta1 for idx      0.998
Theta2 for idx      0.361
Theta3 for idx      0.176
Theta4 for idx      0.372
Theta5 for idx      0.610
Theta6 for idx      0.715
Theta7 for idx      0.802
\end{Soutput}
\begin{Sinput}
R> #M-model
R> hyper.mmodel$summary.hyperpar
\end{Sinput}
\begin{Soutput}
                  mean      sd quant0.025 quant0.25 quant0.5
Theta1 for idx  0.2659 0.18179     0.0250    0.1169   0.2296
Theta2 for idx  0.9963 0.00335     0.9874    0.9950   0.9972
Theta3 for idx  0.9900 0.01186     0.9573    0.9875   0.9938
Theta4 for idx  0.5367 0.05211     0.4358    0.5011   0.5359
Theta5 for idx  0.1836 0.03697     0.1084    0.1590   0.1846
Theta6 for idx  0.1566 0.04985     0.0610    0.1224   0.1554
Theta7 for idx  0.0373 0.04748    -0.0546    0.0049   0.0366
Theta8 for idx  0.2946 0.02608     0.2423    0.2771   0.2951
Theta9 for idx  0.1820 0.04260     0.1007    0.1527   0.1809
Theta10 for idx 0.1292 0.08541    -0.0398    0.0715   0.1295
Theta11 for idx 0.2395 0.04345     0.1544    0.2100   0.2393
Theta12 for idx 0.3893 0.04152     0.3068    0.3613   0.3896
                quant0.75 quant0.975
Theta1 for idx      0.382      0.685
Theta2 for idx      0.999      1.000
Theta3 for idx      0.997      0.999
Theta4 for idx      0.571      0.640
Theta5 for idx      0.209      0.253
Theta6 for idx      0.190      0.256
Theta7 for idx      0.069      0.131
Theta8 for idx      0.312      0.344
Theta9 for idx      0.210      0.268
Theta10 for idx     0.187      0.295
Theta11 for idx     0.269      0.325
Theta12 for idx     0.417      0.470
\end{Soutput}
\end{Schunk}

Note that the first hyperparameter in the PMCAR model is the spatial
autocorrelation, which is very close to one. All the other parameters are the
variances and correlation parameters.

In order to recover the variance matrix, the off-diagonal entries need to be
computed. Note that these depend on three parameters (i.e., correlation and
marginal variances) and that computing these entries involves multivariate
inference. As the joint posterior
distribution of these three parameters needs to be estimated sampling will
be used.

For approximate multivariate posterior inference, \pkg{INLA} can draw samples
from the (approximate) joint posterior of the hyperparameters using function
\code{inla.posterior.sample}. This sampling method is based on the internal
representation of the model, which is based on values of the ensemble
of the hyperparametes and their posterior log-densities \citep[see, for example, ][for details]{INLAbook:2020}.

The internal representation of the model stores different ensembles of values
of the hyperparameters $\{\gamma_g\}_{g=1}^G$ and associated values of the
log-posterior density. Instead of sampling with \code{inla.posterior.sample} we
will re-scale the log-posterior densities to obtain weights associated to
$\{\gamma_g\}_{g=1}^G$. By avoiding sampling we obtain short computing times
here. These weights are the posterior probabilities of the ensembles of values
of the hyperparameters and they can be used to compute posterior quantities of
interest. Furthermore, posterior estimates of transformations of the
hyperparameters can be computed as well. Note that, because the ensemble of
values is available, multivariate inference on the hyperparameters is also
possible. We will rely on this fact to estimate some of the quantities of
interest such as the posterior mean of the between-diseases variance matrix.

For the IMCAR and PMCAR models, the posterior means of the entries of the
between-diseases variance matrix are:

\begin{Schunk}
\begin{Sinput}
R> # IMCAR
R> hyper.imcar$VAR.m
\end{Sinput}
\begin{Soutput}
       [,1]   [,2]  [,3]
[1,] 0.2440 0.0811 0.130
[2,] 0.0811 0.1273 0.118
[3,] 0.1298 0.1178 0.224
\end{Soutput}
\begin{Sinput}
R> # PMCAR
R> hyper.pmcar$VAR.m
\end{Sinput}
\begin{Soutput}
       [,1]   [,2]  [,3]
[1,] 0.2653 0.0874 0.139
[2,] 0.0874 0.1410 0.129
[3,] 0.1395 0.1286 0.255
\end{Soutput}
\end{Schunk}

For the M-model, the variance of the between diseases variance is
given by $M^{\top}M$:

\begin{Schunk}
\begin{Sinput}
R> # M-model
R> hyper.mmodel$VAR.m
\end{Sinput}
\begin{Soutput}
      [,1]  [,2]  [,3]
[1,] 0.353 0.104 0.174
[2,] 0.104 0.127 0.147
[3,] 0.174 0.147 0.226
\end{Soutput}
\end{Schunk}

Point estimates of the variance matrix are very similar for the PMCAR and the
M-model. In addition, all models seem to point to a higher correlation between
lung cancer (2nd disease in the model)  and oral cancer (3rd disease) than
cirrhosis (1st disease) with any of the two other diseases. This makes sense as
lung and oral cancer are known to be highly correlated
\citep{BotellaRocamoraetal:2015}.

Finally, the different model selection criteria computed for the models
(DIC and WAIC) can be compared to choose the best model. The values of these
criteria are shown in Table~\ref{tab:CVALIC}. In this case, both criteria
point to the M-model as the best model in this example.

\begin{table}[h!]
\centering
\begin{tabular}{c|c|c|c|c}
Model & DIC & WAIC\\
\hline
IMCAR &
7407.79 & 7396.8 \\
PMCAR &
7410.04 & 7388.57 \\
M-model &
7384.47 & 7366.15 \\
\end{tabular}
\caption{Information criteria computed for model selection.}
\label{tab:CVALIC}
\end{table}

\section{Discussion}
\label{sec:discussion}

The \pkg{INLAMSM} package builds on top of the \pkg{INLA} package and implements
a number of multivariate spatial latent effects. Hence, this package allows
an easy and simple definition of these multivariate effects to be used within
a formula term to fit multivariate spatial models to lattice data.

Implementation of new latent effects for multivariate data is straightforward
with the \code{rgeneric} latent effect included in the \pkg{INLA} package.
This only requires the specification of the latent effect as a GMRF,
which means that the mean, precision matrix and priors for the hyperparameters
need to be provided. Once the model is implemented, it is easy to include
it in the model formula to fit multivariate models with \pkg{INLA}.

We find important to mention that although all the models in \pkg{INLAMSM} have
been defined through multivariate ensembles of spatial processes, \pkg{INLAMSM}
could be also used for reproducing some multivariate spatial models defined
through univariate or multivariate conditional distributions, which induce
cross-correlations between the different spatial patterns
\citep{SainCressie:2002}.  As shown, a close connection can be drawn between
both approaches \citep{MartinezBeneito:2020} and some important classes of
conditional spatial multivariate models can be reformulated as (marginal)
multivariate models, as those reproduced in \pkg{INLAMSM}. Therefore, this
package could be used as a starting point for multivariate models formulated as
ensembles of conditional distributions

In addition, the latent models implemented in the package can be used as
templates to implement new multivariate spatial models for lattice data.  In
the future, we hope to increase the number of multivariate spatial models
in the package to include other different types of spatial cross-correlation.

In the examples provided to illustrate the use of the package we have
considered a small dataset to fit all possible models. Times required
to fit the models are short. The second example deals with a region with
a larger number of areas which shows that our package can be used
together with \pkg{INLA} to fit multivariate models.

Despite our focus on multivariate spatial models for disease mapping, it is
worth mentioning that the multivariate models implemented in the \pkg{INLAMSM}
package can be used in other contexts.  Furthermore, these multivariate spatial
models can also be used to build temporal and spatio-temporal models by using a
temporal adjacency matrix.

\section*{Acknowledgments}

This work has been supported by grants PPIC-2014-001-P and SBPLY/17/180501/000491, funded by Consejer\'ia de Educaci\'on, Cultura y Deportes (Junta de Comunidades de Castilla-La Mancha, Spain) and FEDER, grant MTM2016-77501-P, funded by Ministerio de Econom\'ia y Competitividad (Spain) and a grant to support
research groups by the University of Castilla-La Mancha (Spain).

F. Palm\'i-Perales has been supported by a Ph.D. scholarship awarded by the University of Castilla-La Mancha (Spain).