PC priors for residual correlation parameters in one-factor mixed models

TL;DR

This paper introduces penalized complexity priors for residual correlation parameters in one-factor Bayesian mixed models, enabling better shrinkage towards independence and consistent scaling across models, demonstrated through a community ecology case study.

Contribution

It develops a unified PC prior framework for correlation parameters in one-factor mixed models, facilitating practical prior specification and model comparison.

Findings

PC priors effectively shrink towards independence

Prior scaling depends only on variance explained by groups

Application shows improved model comparison in ecology

Abstract

Lack of independence in the residuals from linear regression motivates the use of random effect models in many applied fields. We start from the one-way anova model and extend it to a general class of one-factor Bayesian mixed models, discussing several correlation structures for the within group residuals. All the considered group models are parametrized in terms of a single correlation (hyper-)parameter, controlling the shrinkage towards the case of independent residuals (iid). We derive a penalized complexity (PC) prior for the correlation parameter of a generic group model. This prior has desirable properties from a practical point of view: i) it ensures appropriate shrinkage to the iid case; ii) it depends on a scaling parameter whose choice only requires a prior guess on the proportion of total variance explained by the grouping factor; iii) it is defined on a distance scale common to all group models, thus the scaling parameter can be chosen in the same manner regardless the adopted group model. We show the benefit of using these PC priors in a case study in community ecology where different group models are compared.

Tables

Table 1. Table 1 : Comparison between group models ℳ 1 subscript ℳ 1 \mathcal{M}_{1} , ℳ 2 subscript ℳ 2 \mathcal{M}_{2} and ℳ 3 subscript ℳ 3 \mathcal{M}_{3} . Results are shown for different prior settings, i.e. prior median ICC equal to 0.1 0.1 0.1 , 0.5 0.5 0.5 and 0.9 0.9 0.9 . For the OU process, the displayed ρ 𝜌 \rho values refers to the correlation at one meter distance between locations along the transects.

prior setting 1: median ICC $=0.1$
		$\rho$
grouping factor	group model	0.025q	mean	0.975q	log.mlik
campaign	exch	0.152	0.3	0.492	-392.522
transect	exch	0.171	0.282	0.418	-401.053
transect	ou	0.036	0.109	0.207	-412.727
prior setting 2: median ICC $=0.5$
		$\rho$
grouping factor	group model	0.025q	mean	0.975q	log.mlik
campaign	exch	0.176	0.331	0.52	-392.848
transect	exch	0.184	0.3	0.439	-401.147
transect	ou	0.040	0.116	0.216	-413.779
prior setting 3: median ICC $=0.9$
		$\rho$
grouping factor	group model	0.025q	mean	0.975q	log.mlik
campaign	exch	0.191	0.349	0.538	-393.208
transect	exch	0.186	0.303	0.443	-401.639
transect	ou	0.040	0.117	0.218	-414.439

Equations

$$\mathbf{Y}=\mathbf{X}\mathbf{\beta }+\mathbf{Z}\mathbf{b}+\mathbf{ϵ},;\mathbf{b}\sim \mathcal{N}(0,{\mathbf{\Sigma }}_b);\mathbf{ϵ}\sim \mathcal{N}(0,\mathbf{I}{\sigma }_e^2)$$

$$y_{ij}$$

$$b_j$$

$${ϵ}_{ij}$$

$$y_{ij}=\alpha +{\mathbf{x}}_{ij}^{\text{T}}\mathbf{\beta }+{\theta }_{ij}i=1,\dots ,m_{j}j=1,\dots ,n,$$

$$Corr({\theta }_{ij},{\theta }_{hj})=\frac{{\sigma }_b^2}{{\sigma }_b^2+{\sigma }_{ϵ}^2}i,h\in \{1,\dots ,m_j\}\forall j.$$

$${\mathbf{\theta }}_j$$

$${\mathbf{R}}_j(\rho )=\left[\begin{array}{ccccc}1& \rho & \cdots & \cdots & \rho \\ \rho & 1& \ddots & & ⋮\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ ⋮& & \ddots & 1& \rho \\ \rho & \cdots & \cdots & \rho & 1\end{array}\right].$$

$$y_{ij}$$

$$b_{1j}$$

$${ϵ}_{ij}$$

$$Corr({\theta }_{ij},{\theta }_{hj})=\frac{{\sigma }_b^2}{{\sigma }_b^2+{\sigma }_{ϵ}^2}{\tilde{\rho }}^{|i-h|}i,h\in \{1,\dots ,m_j\}\forall j.$$

$${\mathbf{R}}_j(\rho )=\left[\begin{array}{ccccccc}1& \rho & {\rho }^2& \cdots & \cdots & {\rho }^{m-1}& {\rho }^m\\ \rho & 1& \rho & {\rho }^2& \dots & {\rho }^{m-2}& {\rho }^{m-1}\\ {\rho }^2& \rho & 1& \rho & {\rho }^2& \cdots & {\rho }^{m-2}\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ {\rho }^{m-2}& \cdots & {\rho }^2& \rho & 1& \rho & {\rho }^2\\ {\rho }^{m-1}& {\rho }^{m-2}& \dots & {\rho }^2& \rho & 1& \rho \\ {\rho }^m& {\rho }^{m-1}& \cdots & \cdots & {\rho }^2& \rho & 1\end{array}\right],$$

$$Corr({\theta }_{ij},{\theta }_{hj})=\frac{{\sigma }_b^2}{{\sigma }_b^2+{\sigma }_{ϵ}^2}\exp (-{\delta }_{ih}\varphi ),\varphi >0.$$

$${\pi }_0(\mathbf{\theta })=\mathcal{N}(0,{\mathbf{C}}_0);{\mathbf{C}}_0={\mathbf{I}}_M,$$

$$\pi (\mathbf{\theta })=\mathcal{N}(0,\mathbf{C});\mathbf{C}=\text{diag}\left\{{\mathbf{R}}_1(\rho ),\dots ,{\mathbf{R}}_n(\rho )\right\},$$

$$\text{KLD}(\pi ||{\pi }_0)=\int \pi (\mathbf{\theta })\log \left(\frac{\pi (\mathbf{\theta })}{{\pi }_0(\mathbf{\theta })}\right)d\mathbf{\theta }.$$

$$\text{KLD}(\pi ||{\pi }_0)=\frac{1}{2}\left[\text{trace}\left({\mathbf{C}}_0^{-1}\mathbf{C}\right)-M-\log \left(\frac{|\mathbf{C}|}{|{\mathbf{C}}_0|}\right)\right],$$

$$\text{KLD}(\pi ||{\pi }_0)=-\frac{{\sum }_{j=1}^n\log \left(|{\mathbf{R}}_j(\rho )|\right)}{2}.$$

$$d(\rho )=\sqrt{-\sum _{j=1}^n\log \left(|{\mathbf{R}}_j(\rho )|\right)}.$$

$$\pi (d(\rho ))=\lambda \exp \left(-\lambda d(\rho )\right),\lambda >0.$$

$$\pi (\rho )$$

$$d(\rho )=\sqrt{-n\log \left((1+(m-1)\rho )(1-\rho {)}^{m-1}\right)}0\le \rho <1.$$

$$\pi (\rho )=\frac{m-1}{2}\left(\frac{1}{1-\rho }-\frac{1}{1+(m-1)\rho }\right)\frac{{\lambda }^{\prime }}{\sqrt{-\log \left(|\mathbf{R}(\rho )|\right)}}\exp (-{\lambda }^{\prime }\sqrt{-\log \left(|\mathbf{R}(\rho )|\right)}),$$

$$d(\rho )=\sqrt{n(1-m)\log (1-{\rho }^2)}0\le \rho <1.$$

$$\pi (\rho )=\frac{\rho (m-1)}{1-{\rho }^2}\frac{{\lambda }^{\prime }}{\sqrt{-\log \left(|\mathbf{R}(\rho )|\right)}}\exp (-{\lambda }^{\prime }\sqrt{-\log \left(|\mathbf{R}(\rho )|\right)}),$$

$$d(\varphi )=\sqrt{n(1-m)\log (1-\exp (-2\varphi ))}\varphi >0.$$

$$\pi (\varphi )=\frac{(m-1)\exp (-2\varphi )}{1-\exp (-2\varphi )}\frac{{\lambda }^{\prime }}{\sqrt{-\log \left(|\mathbf{R}(\varphi )|\right)}}\exp (-{\lambda }^{\prime }\sqrt{-\log \left(|\mathbf{R}(\varphi )|\right)}),$$

$$\mathbb{P}(\rho <U)={\int }_0^U\lambda \exp (-\lambda d(\rho ))\left|\frac{\partial d(\rho )}{\partial \rho }\right|=1-\exp (-\lambda d(U))=a.$$

$$\lambda =-\log (1-a)/d(U),$$

$$\text{BF}({\mathcal{M}}_1;{\mathcal{M}}_2)$$

$$\pi (d(\rho ))=\lambda \exp (-\lambda d(\rho ))\lambda >0.$$

Full text

Keywords: Bayesian mixed models, Group model; One-way anova; INLA; Intra-class correlation; Within group residuals.

Address for Correspondence: Massimo Ventrucci, massimo.ventrucci@unibo.it

1 Introduction

1.1 Mixed models in community ecology

The understanding of factors determining the distribution of organisms is a striking goal of community ecology (Heino,, 2013) and the key for forecasting the future trajectories of communities (Wisz et al.,, 2013). Modelling the organization and the evolution of natural communities is not an easy task, since the abundance and distribution of organisms depend on environmental drivers as well as on single taxon features and interactions among different taxa (Ovaskainen et al.,, 2017; Wisz et al.,, 2013). Several authors report high levels of unexplained variation, after considering the effect of environmental variables (Lamouroux et al.,, 2004). This residual variation is often ascribed to biotic (intra- and inter-specific) interactions, including both negative (like competition, predation and parasitism) and positive interactions (like mutualism and commensalism) which can play a crucial role in shaping communities.

A study in community ecology typically consists of observations of species abundance (biotic response variable) and environmental covariates (abiotic factors) collected at different locations and/or time points, through several sampling campaigns. At an early stage of the analysis ecologists’ goal is to quantify the effect of the observed covariates, often assuming a linear regression model. Two alternative assumptions on the residuals are under examination:

• (a)

residuals are iid;

• (b)

residuals are correlated following some dependence structure.

Assumption (a) means that the linear regression model including the observed covariates is the true model; in other words, what is not explained by covariates is iid noise. Assumption (b) means that residuals are not iid, thus there is unobserved heterogeneity in the data. It is crucial for ecologists to investigate the correlation pattern in the residuals, as this may reflect the presence of unobserved processes playing a role in shaping the underlying communities.

Mixed models are the most used tools for evaluating the two assumptions above, and have been extensively used in analyzing ecological data (Zuur et al.,, 2009). The popularity of mixed model in ecology is probably due to the fact that the effect of the observed covariates and unobserved processes can be neatly separated in the model. In its general formulation, a linear mixed model for a Gaussian response $\bm{Y}$ and covariates $\bm{X}$ is expressed as

[TABLE]

where $\bm{\beta}$ are denoted as fixed effects and $\bm{b}$ as random effects, i.e. random variables with Gaussian distribution conditional on one or more variance (hyper-)parameter. The usual interpretation in ecology is that the $\bm{\beta}$’s account for variability in the data explained by observed abiotic factors, while the $\bm{b}$’s account for sources of variability in $\bm{Y}$ driven by unobserved abiotic or biotic factors (Warton et al.,, 2015).

1.2 One-way anova

Let us consider one-way anova, which is the simplest mixed model case. Assume data are grouped according to the levels of a grouping factor, with $y_{ij}$ being the response at unit $i=1,\ldots,m_{j}$ within group $j=1,\ldots,n$. The one-way anova model is:

[TABLE]

where $b_{j}$ and $\epsilon_{ij}$ are assumed as independent. The $b_{j}$’s are random effects quantifying group-specific deviations from the global intercept $\alpha$. It is important to note that i) when $\sigma_{b}^{2}=0$ model (1) corresponds to $y_{ij}=\alpha+\boldsymbol{x}_{ij}^{\textsf{T}}\bm{\beta}+\epsilon_{ij}$, where only covariates matter and the rest is iid variation; ii) the random effects $b_{j}$’s and $\epsilon_{ij}$’s compete to capture the total residual variability, i.e. the variance of the terms $y_{ij}-\alpha-\boldsymbol{x}_{ij}^{\textsf{T}}\bm{\beta}$.

By reparametrizing model (1) - see Section 3 - it can be shown that the $b_{j}$’s induce correlation in residuals belonging to the same group; in this case, within group residuals are exchangeable. A possible ecological explanation for such type of unobserved heterogeneity is that members of the same group interact with each other or share some common features.

1.3 Aim of the work

Model (1) is an example of what in this paper is denoted as group model, i.e. a model for the within group residual correlation structure. In particular, model (1) assumes an exchangeable group model. The first aim of this paper is to illustrate how group models beyond the exchangeable case can be constructed by extending formulation (1). The focus will be on models with one grouping factor, i.e. one-factor mixed models. Comparing different group models can provide information on the main sources of heterogeneity in the data. The need for reliable model comparison tools requires to build sensible priors for the correlation parameters in group models, which is the second purpose of this paper.

The usual practice in Bayesian analysis is to select independent prior distributions for the variance components $\sigma_{b}^{2}$ and $\sigma_{\epsilon}^{2}$ of model (1). The literature on priors for variance parameters is vast, see for instance Gelman, (2006), and it is outside the scope of this work to give a comprehensive review of it. We note that common choices like conjugate priors lead to overfitting, as shown by several papers (Frühwirth-Schnatter and Wagner,, 2010, 2011; Simpson et al.,, 2017).

The priors we propose exploit two pieces of prior knowledge, one regarding the structure of the model and the other concerning its variance components. The first piece of information is about the base model, i.e. - applying the definition in Simpson et al., (2017) - the simplest model in the class of group models. Coherently with the assumptions (a) and (b) aforementioned, an ecologist would find that the natural base model for (1) is the one where $\sigma_{b}^{2}=0$, corresponding to assumption (a). We stress the fact that this information is certain and for free, as it simply reflects the fact that a mixed model is an extension of linear regression. The second piece of information - not for free and uncertain - is about the relative weight of the two variance components of the mixed model. The idea is that while an ecologist may have no opinion on the range of plausible values for $\sigma_{b}^{2}$ and $\sigma_{\epsilon}^{2}$, he/she may have an intuition on the relative importance of correlated residuals (controlled by $\sigma_{b}^{2}$) versus iid residuals (controlled by $\sigma_{\epsilon}^{2}$). This is essentially asking for a prior guess on the proportion of total (residual) variance explained by the grouping factor.

In order to make good use of the knowledge about the base model we advocate priors that avoid overfitting by construction, i.e. priors that always give a chance to the iid base model to arise in the posterior, unless the data do require a more flexible one. Parsimony is a reasonable leading principle in modelling ecological data, where unobserved heterogeneity plays a major role. It is therefore desirable to have flexible models that are able to shrink to simple ones. For building such priors we apply the Penalized complexity prior framework by Simpson et al., (2017).

In order to exploit the second piece of information, we abandon the idea of selecting independent priors for the variance components. Instead, we exploit a common reparameterization of model (1) and define a prior on the intra-class correlation (ICC) paramater, i.e. the proportion of total residual variance explained by the grouping factor. The resulting PC prior depends on a scaling parameter that can intuitively be elicited based on a prior opinion on the ICC. Importantly, this scaling can be used in general for any group model because PC priors are invariant over reparameterization.

The plan of the paper is as follows. In Section 2 the motivating example is described. In Section 3 group models are presented, distinguishing between the cases of exchangeable and structured residuals. The PC prior for the correlation parameter in a generic group model is derived in Section 4, with additional results regarding the balanced design case. An illustration of the proposed group models and PC priors is given in Section 5, where we also describe a strategy to perform group model comparison by means of the Bayes factor. The paper closes with a discussion in Section 6.

2 Motivating example

The motivating example for this work concerns a study on macroinvertebrate communities from data collected in six sampling campaigns carried out in three different streams, tributaries of the Po River (Northern Italy): Nure Stream, Parma Stream and Enza Stream. For each river a sampling area was selected and sampled twice, once in summer and once in winter. The spatial design included $50$ random points in each area, aligned along several transects. At each point, abundance of macroinvertebrates (response) and environmental covariates such as flow velocity, water depth, substrate composition and benthic organic matter were recorded.

The application goal is to investigate the role of the environmental covariates and the presence of small scale processes within macroinvertebrate communities. In Section 5 we propose different group models for the residuals, as an exploratory analysis to understand the main sources of unobserved heterogeneity. If we consider campaign as grouping factor, data are grouped in $n=6$ group/campaigns having $m=50$ observations. If we consider transect as grouping factor, data are grouped in $n=38$ groups/transects, each having a varying number (between $7$ and $10$) of irregularly-spaced observations.

3 Group models

We distinguish between two broad classes of models for the within group residuals: exchangeable residuals and structured residuals. For the latter, we discuss in detail two group models, namely the autoregressive of order 1 (AR1) process and its continuous version known as the Ornstein Uhlenbeck (OU) process; these are particularly useful for the case study considered in this paper, but are also relevant in general applications within ecology where observations are often taken at different time points. Other group models can be constructed following the same idea illustrated here.

3.1 Exchangeable residuals

Let us reparametrize (1) as

[TABLE]

where $\theta_{ij}=b_{j}+\epsilon_{ij}$ is the residual at unit $i$ within group $j$. It follows that

[TABLE]

This means that the distribution of the residuals within group $j$, $\bm{\theta}_{j}=(\theta_{1j},\ldots,\theta_{mj})^{T}$, is unchanged under permutation of the indexes $1,\ldots,m$, leading to the exchangeable model

[TABLE]

where $\sigma^{2}=\sigma^{2}_{b}+\sigma_{\epsilon}^{2}$ is the total variance and the correlation matrix is

[TABLE]

Following the constrained interpretation (Verbeke and Molenberghs,, 2003), where the variance components are restricted to be non negative, the correlation parameter $\rho=\sigma^{2}_{b}/(\sigma^{2}_{b}+\sigma_{\epsilon}^{2})$ is the proportion of total (residual) variance explained by the grouping factor, i.e. the ICC. We use notation $\boldsymbol{R}_{j}(\rho)$ to emphasize that the correlation matrix depends on $\rho$. Note that if design is balanced then $m_{j}=m,\forall j$, hence $\bm{R}_{j}(\rho)=\bm{R}(\rho),\forall j$.

Matrix (3) implies that residuals within each group are mutually correlated, with correlation parameter equal to $\rho$. We note that for any $\sigma^{2}>0$, $\rho=0$ identifies the linear regression model, thus $\rho$ is the parameter responsible for the shrinkage towards the base model. In Section 4 we introduce a PC prior for $\rho$ depending on a scaling parameter that can be specified in a very intuitive manner, using a prior statement about the ICC.

3.2 Structured residuals

3.2.1 Autoregressive of order one

The AR1 process is often used to model correlation over time, when observations are taken at regularly-spaced time points (e.g. days, weeks, etc). Assume the model

[TABLE]

which defines an AR1 process on $b_{ij},i=1,\ldots,m$, with $\sigma_{I}^{2}$, $\tilde{\rho}$ and $\sigma_{b}^{2}=\frac{\sigma_{I}^{2}}{1-\tilde{\rho}^{2}}$ being, respectively, the innovation variance, the lag-one correlation and the marginal variance of the process.

Analogously to the exchangeable case, by reparameterizing $\theta_{ij}=b_{ij}+\epsilon_{ij}$, it follows that:

[TABLE]

In compact notation the group model can be rewritten as in Eq. (2) with total variance $\sigma^{2}=\sigma^{2}_{b}+\sigma^{2}_{\epsilon}$ and correlation matrix,

[TABLE]

where the correlation parameter at generic lag $|i-h|$ is equal to $\rho^{|i-h|}=\frac{\sigma^{2}_{b}}{\sigma^{2}_{b}+\sigma^{2}_{\epsilon}}\tilde{\rho}^{|i-h|}$.

This group model implies that within group residuals are correlated, with correlation structure driven by the ordering of the observations. We note that, analogously to the exchangeable case, $\rho=0$ identifies the linear regression model, for any $\sigma^{2}>0$. In Section 4 we will describe a PC prior for $\rho$ whose scaling parameter can be chosen in the same intuitive way as in the exchangeable case.

3.2.2 Ornstein Uhlenbeck

The OU process is the continuous version of the AR1 and is appropriate when the observations are not equally-spaced. The model is conceptually the same as the AR1 but reparameterized to account for the distances between locations. Let us assume $\delta_{ih}$ is the distance between observations $i$ and $h$, the correlation is

[TABLE]

This model can also be rewritten as in Eq. (2) with total variance $\sigma^{2}=\sigma^{2}_{b}+\sigma^{2}_{\epsilon}$ and correlation matrix as in (4), with correlation between $i$ and $h$ equal to $\rho_{i,h}=\frac{\sigma^{2}_{b}}{\sigma^{2}_{b}+\sigma^{2}_{\epsilon}}\exp(-\delta_{ih}\phi)$, thus the only difference w.r.t. the AR1 case is that the correlation parameter is expressed in log scale, $\phi=-\log(\rho)$. The OU process has good computational properties due to its sparse tridiagonal precision matrix, whose structure is detailed in Finley et al., (2009). We use the OU process to model residual correlation when observations are randomly located along transects, as it is the case in our motivating example.

Note that, for any $\sigma^{2}>0$, the linear regression model is achieved in the limit for $\phi\to\infty$. In Section 4 we describe a PC prior for $\phi$ whose parameters can be specified in the same intuitive way as in the exchangeable case.

4 PC priors for group models

In each of the group models presented above we have two (hyper-)parameters that need to be assigned a prior, the marginal variance $\sigma^{2}$ and the correlation $\rho$. What is relevant for us is to derive the PC prior for $\rho$, as this is the only parameter responsible for the shrinkage to the iid base model.

A PC prior is defined in Simpson et al., (2017) as an exponential distribution on a distance scale, measuring model complexity w.r.t. the base model. Simpson et al., (2017) argue that, because the mode of the exponential is at distance [math] (i.e. at the base model), PC priors guard against overfitting by construction as they always give non-zero probability to a neighbourhood of the base model. This argument is based on an informal definition of an overfitting prior as a prior density that is zero at the base model (Simpson et al.,, 2017); the idea is that an overfitting prior may drag the posterior away from the base model even when the latter is the true model; on the contrary a prior that contracts to the base model prevents overfitting by default. Several papers have confirmed tendency of PC priors to prevent overfitting via simulation studies (Fuglstad et al.,, 2018; Klein and Kneib,, 2016; Ventrucci and Rue,, 2016).

Another property of PC priors which is exploited in the context of this paper is that they are invariant over reparameterization, as they are defined on a distance scale (instead of the original scale of the parameter), then translated in the scale of the original parameter by the change of variable rule.

4.1 The PC prior for $\rho$

With no loss of generality we present the steps to construct the PC prior for $\rho$ assuming $\sigma^{2}=1$; for a full discussion of the principles underpinnning the PC prior framework see Simpson et al., (2017).

Let $\bm{\theta}=\left(\bm{\theta}_{1},\ldots,\bm{\theta}_{n}\right)^{\textsf{T}}$ denote the vector of residuals from all groups, $\pi(\bm{\theta})$ the flexible model and $\pi_{0}(\bm{\theta})$ the base model. The base model corresponds to $\rho=0$, as

[TABLE]

where $\bm{C}_{0}$ is the identity matrix of dimension $M=\sum_{j=1}^{n}m_{j}$. The flexible model is when $\rho>0$, hence

[TABLE]

where $\bm{C}$ is a block diagonal matrix containing all the within group correlation matrices.

The first step is computation of the Kullback Leibler divergence (KLD, Kullback and Leibler, (1951)) between the flexible and the base model,

[TABLE]

In our case $\pi$ and $\pi_{0}$ are zero-mean multivariate normal densities of dimension $M$ with covariances $\bm{C}$ and $\bm{C}_{0}$, respectively, thus the KLD simplifies to:

[TABLE]

where notation $|\cdot|$ indicates the matrix determinant. Given that $\bm{C}$ has block diagonal structure, we obtain:

[TABLE]

For mathematical convenience, the distance from the base model is expressed as $\sqrt{2{\text{KLD}}(\pi_{1}||\pi_{0})}$, which gives

[TABLE]

From Eq. (5) we see that the distance from the base model is a function of $\rho$ which takes values in the interval $[0,\infty)$; it is [math] at the base model ($\rho=0$) and goes to $\infty$ as $\rho\to 1$.

The next step requires to specify an exponential distribution on $d(\rho)$ with rate $\lambda$,

[TABLE]

Here $\lambda$ plays the role of a scaling parameter, controlling the degree of penalty for deviating from the base model. The larger $\lambda$, the stronger the penalty for deviating from the base model, at prior. Finally, the PC prior for $\rho$ is derived by the change of variable rule:

[TABLE]

The PC prior in Eq. (7) depends on $\rho$ through the determinant of the within group $j$ correlation matrix $|\bm{R}_{j}(\rho)|$ and its derivative $\frac{\partial|\bm{R}_{j}(\rho)|}{\partial\rho}$. Therefore, one can derive analytically (or compute numerically) the PC prior for different group models by just plugging-in the determinant and its derivative in (7). In order to implement the PC prior in our model we need to choose $\lambda$ in Eq. (7); we postpone the discussion on how to choose it to Section 4.3.

4.2 The balanced design case

If design is balanced, PC priors for all the group models of Section 3 can be derived analytically. In the unbalanced case, closed form expressions are more involved and will not be presented here; a practical solution for unbalanced designs is to evaluate the PC prior numerically. Below we report the PC priors for the correlation parameter in each group model of Section 4, for the balanced case. See Appendix A for the mathematical details.

4.2.1 Exchangeable residuals

In the exchangeable case, $\rho$ is the within group correlation. The distance is

[TABLE]

The PC prior is

[TABLE]

where $|\bm{R}(\rho)|=(1+(m-1)\rho)(1-\rho)^{m-1}$ and $\lambda^{\prime}=\lambda\sqrt{n}$.

4.2.2 Structured residuals: AR1

In the AR1 case, $\rho$ is the lag-one correlation parameter. The distance is

[TABLE]

The PC prior is

[TABLE]

where $|\bm{R}(\rho)|=(1-\rho^{2})^{m-1}$ and $\lambda^{\prime}=\lambda\sqrt{n}$.

4.2.3 Structured residuals: OU

In the OU case, $\phi=-\log(\rho)$ is the lag-one correlation parameter expressed in the log scale. The distance is

[TABLE]

The PC prior is

[TABLE]

where $|\bm{R}(\phi)|=(1+\exp(-2\phi))^{m-1}$ and $\lambda^{\prime}=\lambda\sqrt{n}$.

4.3 Choice of $\lambda$

The degree of informativeness of the PC prior in Eq. (6) can be managed through $\lambda$, that defines the “prior distance” from the base model. Simpson et al., (2017) proposed to select $\lambda$ through the following rule: set values $U$ and $a$ such that $\mathbb{P}(\rho<U)=a$. By working out the cumulative distribution function of $\pi(\rho)$ we have,

[TABLE]

Solving Eq. (11) for $\lambda$ we obtain

[TABLE]

hence, for the PC priors in Section 4.2, $\lambda^{\prime}=-(\sqrt{n}\log(1-a))/d(U)$. Note that $\lambda$ in Eq. (12) can be computed in a generic group model (e.g. exchangeable, AR1, OU) by just plugging-in the associated distance function evaluated at $U$.

Scaling the PC prior (i.e., defining $\lambda$) following this rule becomes very intuitive: if $a$ is small, $U$ can be thought of as a lower bound on $\rho$; if $a=0.5$, then $U$ is the prior median for $\rho$. In the exchangeable case, because of the interpretation of $\rho$ as the ICC, the user could simply set $U$ and $a$ according to a prior statement on the proportion of total variance explained by the grouping factor; e.g., for $a=0.5$, choice of $\lambda$ translates into eliciting the “median proportion of variance explained by the grouping factor”. In presence of weak (or no) information, a sensible strategy is to set the median ICC to $0.5$, so to be exactly half way between the two opposite scenarios: residuals are iid ($\rho=0$) vs residuals are completely predicted by the grouping factor ($\rho=1$). Figure 1 shows the PC prior in Eq. (8) for different choices of the median ICC. It can be seen that the smaller the median ICC set at prior, the stronger the penalty for deviating from the iid base model.

5 Application

We illustrate an application of the proposed models with the data described in Section 2. In order to investigate presence of small scale processes within macroinvertebrate communities we propose comparison of different group models by means of the Bayes factor (Kass and Raftery,, 1995). Sørbye and Rue, (2018) note that, in case the models under comparison are extensions of the same base model, a convenient approach would be to set, for both models, the same prior on the distance from the base model; in this way, the resulting Bayes factor will only respond to the models under comparison, while being insensitive to the choice of priors for the model hyper-parameters. In the remaining of this section we provide an example showing the advantages of PC priors connected to Bayes factors and then discuss the results for our case study.

5.1 Comparing different group models using the Bayes factor

Consider the reparametrized mixed model $y_{ij}=\alpha+\bm{x}_{ij}^{\textsf{T}}\bm{\beta}+\theta_{ij}$ as in Section 3, where $\bm{\theta}_{j}$ is the vector of residuals within group $j$. Let us compare the following models for $\bm{\theta}_{j}$:

• •

$\mathcal{M}_{1}$: $\bm{\theta}^{(1)}_{j}$ are exchangeable residuals like in Eq. (2); this model has hyper-parameters $\{\sigma^{2},\rho\}$: $0\leq\rho<1$ is the intraclass correlation and $\sigma^{2}$ is the variance;

• •

$\mathcal{M}_{2}$: $\bm{\theta}^{(2)}_{j}$ are structured residuals following an OU process like in Section 3.2.2; this model has hyper-parameters $\{\sigma^{2},\phi\}$: $\phi>0$ is the lag-one correlation in log scale and $\sigma^{2}$ is the variance.

The Bayes factor (BF) quantifies the strength of evidence of $\mathcal{M}_{1}$ compared to $\mathcal{M}_{2}$ by the ratio of the marginal likelihoods (considering $\alpha$ and $\bm{\beta}$ as known, without loss of generality):

[TABLE]

From Eq. (13) we can see that BF depends on the data $\bm{y}$ (and their likelihood) and choices made by the user regarding the models to compare, i.e. $\pi(\bm{\theta}^{(1)}|\sigma^{2},\rho)$ and $\pi(\bm{\theta}^{(2)}|\sigma^{2},\phi)$, and the prior on the hyper-parameters of such models, $\pi(\sigma^{2})$, $\pi(\rho)$ and $\pi(\phi)$.

In general, the effect of the prior cannot be separated from the effect of the model when the BF is used. Quantifying evidence of the alternative models while neutralizing the impact of the priors on hyper-parameters like $\rho$ and $\phi$ is desirable; especially because practitioners are interested in learning about model’s goodness of fit and often have no prior knowledge to inform the prior on $\rho$ and $\phi$. To achieve this goal, the same degree of uncertainty would need to be encoded in the priors on $\{\sigma^{2},\rho\}$ (for $\mathcal{M}_{1}$) and the priors on $\{\sigma^{2},\phi\}$ (for $\mathcal{M}_{2}$). Regarding $\sigma^{2}$, this can be done directly by choosing the same density $\pi(\sigma^{2})$ for both $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$. The hyper-parameter $\rho$ and $\phi$ live in different spaces and have different interpretations, hence to encode the same degree of uncertainty is a much more difficult task. The PC priors proposed in this paper provide a solution to this issue as they are defined on a distance scale $d(\cdot)$, common to both $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$. Following Sørbye and Rue, (2018) we assume PC priors on $d(\rho)$ (for $\mathcal{M}_{1}$) and $d(\phi)$ (for $\mathcal{M}_{2}$) with equal rate $\lambda$; in this way, $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ have the same distance, at prior, from the iid base model. This strategy is reminiscent of the concept of compatible priors (Dawid and Lauritzen,, 2001), i.e. priors that share similarities across the compared models, which were proposed to lessen the influence of priors on Bayes factors in an objective Bayes perspective.

Figure 2 helps in clarifying the benefit of using the same $\lambda$ for different group models. The top panel displays the PC prior on $\rho$ (exchangeable), the central panel shows the PC prior on $\phi$ (OU) and the bottom panel displays both PC priors in the common distance scale. The parameter $\lambda$ is set to have a median ICC equal to $0.5$. While the prior densities materialize differently in the original scales, they are the same in the distance scale (dotted red line and black solid lines are superimposed). In our opinion, this shows that specifying group models having the same degree of complexity w.r.t a common base model (or same median ICC) is, in general, a non trivial task that can easily be addressed using PC priors.

5.2 Results

We consider the model $y_{ij}=\alpha+\bm{x}_{ij}^{\textsf{T}}\bm{\beta}+\theta_{ij}$, where $y$ is the log-abundance of macroinvertebrates, covariates $\bm{x}$ include benthic organic matter (BOM), water depth (P), flow velocity (V), substrate composition (SUB) and season (winter, summer). As a first model we assume residuals $\bm{\theta}_{j}$ are iid. Visual inspection of residuals for this model indicates substantial structure (results not shown here). In order to understand the nature of such unobserved heterogeneity, we focus on several group models for the residuals, analysing the relevance of grouping factors like campaign and transect. We consider three group models for $\bm{\theta}_{j}$:

• •

$\mathcal{M}_{1}$: exchangeable residuals within campaign (grouping factor campaign, group model exch)

• •

$\mathcal{M}_{2}$: exchangeable residuals within transect (grouping factor transect, group model exch)

• •

$\mathcal{M}_{3}$: serially correlated residuals within transect (grouping factor transect, group model ou)

All these models offer a great improvement w.r.t the iid residual model in terms of Bayes factor (figures not shown here) and, perhaps, all three types of structure may be worthy to be included in the final model. However, our aim is not to find the best model for these data, but to illustrate the use of the proposed group models in a real case study. Comparing the above models is beneficial to generate hypotheses on the main sources of unobserved heterogeneity characterizing the ecological community under study.

The model comparison summaries are reported in Table 1, focusing on three different prior settings regarding the median ICC. From top to bottom, the prior on the median ICC is $0.1$ (more importance assigned to iid residuals), $0.5$ (equal weight to iid and grouping factor) and $0.9$ (more weight assigned to the grouping factor). Within each setting, the three group models have the same prior on the distance scale, thus the marginal likelihoods on the last column represent a fair model comparison tool. Importantly, the Bayes factor here is adopted to fairly compare models that differ as regards to both the grouping factor (e.g. campaign vs transect) and the group model (e.g. exch vs ou). Bayes factors can be evaluated in the log scale, by taking the difference of the log marginal likelihoods of the compared models.

All the models were implemented in INLA (Rue et al.,, 2009); in Web Appendix 1 an example of the R-INLA code to fit the exchangeable group model is provided. Regarding the total variance, we adopt the Gumbel$(1/2,\psi)$ type-2 distribution on the precision parameter $1/\sigma^{2}$ (which is the PC prior for the precision of a Gaussian random effect; see Simpson et al., (2017), Sec. 3.3, for a proposal on how to choose the scaling parameter $\psi$).

5.2.1 Exchangeability within campaign vs exchangeability within transect

We first compare $\mathcal{M}_{1}$ against $\mathcal{M}_{2}$, to assess whether residuals are more correlated within campaigns than within transects. From Table 1 we see there is clear evidence that the most relevant grouping factor is campaign; the difference of the log marginal likelihoods is around $9$ units ($401-392$) in favour of exchangeability within campaign, translating into a Bayes factor of around $8100$, i.e. very strong evidence according to the categories proposed by Kass and Raftery, (1995). Finally, the posterior mean for the correlation within campaign/transect is roughly the same, i.e. around $0.3$. All these results are stable for varying priors on the median ICC.

5.2.2 Serial correlation within transect vs exchangeability within transect

We have seen that there is more correlation in the residuals belonging to the same campaign, rather than to the same transect. This does not mean that correlation along transect is not the case. The goal of the study is to assess presence of small scale interactions between organisms. Even if transect has proved to be less important than campaign, it is worth to investigate serial correlation along transects, the latter being a surrogate of small scale interactions between organisms. We then compare $\mathcal{M}_{2}$ against $\mathcal{M}_{3}$, where the latter is an OU process on the transect. From Table 1 we see there is clearly more evidence in favour of the model implying exchangeability, than the one implying serial correlation along transects; the difference of the log marginal likelihoods is around $11$ units ($413-401$), translating into a Bayes factor of around $60000$, i.e. very strong evidence. The posterior for $\rho$ in the OU case refers to the correlation at one meter distance between locations along the transects; this correlation is around $0.11$. All results are stable for varying priors on the median ICC.

In conclusion, the presence of small scale interactions remains an open question that should be investigated further. More generally, conclusive evidence for spatially structured unobserved heterogeneity needs to be investigated via models including all relevant random effects simultaneously. We stress that the model comparisons presented above are intended as an exploratory tool, at a further stage an ecologist might consider as a first model the one including exchangeable random effects within campaign, then adding complexity on top of it.

In Figure 3 we explore covariates effects for the considered group models and the iid case too. All models track the responses of macroinvertebrates to environmental covariates roughly in the same way, except for the covariate substrate composition (SUB). It can be seen that SUB is found to be significant under the iid model and the OU within transect group model, but it is not under an exchangeable group model accounting for within campaign/transect correlation. This points out the importance of model selection for ecologists, in order to avoid type I errors and misinterpretation of the evidence in the data.

6 Discussion

Starting from the the one-way anova case, which assumes within group residuals to be exchangeable, we presented group models encoding different assumptions on the residual structure. All these models are an extension of the linear regression base model (iid residuals). Importantly, the generic group model was parametrized so that only one correlation (hyper-)parameter ($\rho$) is responsible for the shrinkage towards the iid case. We then derived an intuitive PC prior for the correlation parameter which is built from general principles and can therefore be applied in general to different types of one-factor mixed models.

In community ecology studies, residual correlation can be linked to the effect of unaccounted abiotic factors or unobserved biotic processes, like interactions among organisms (e.g. competition, predation etc). Different structures in the residuals match different ecological interpretations about the nature of the unobserved heterogeneity, hence reliable model comparison tools are needed to compare alternative models on the residuals. Saville and Herring, (2009) proposed approximated Bayes factors for testing random effects in linear mixed models, exploiting closed-form solutions derived using conjugate priors for the variance components; in contrast to their work, our paper focuses on the restricted class of one-factor mixed models but embraces group models with generic covariance, without being restricted to conjugate priors.

In our opinion assuming a PC prior for the correlation parameter $\rho$ gives several advantages to the user/ecologist. First, it ensures that the group model shrinks to the linear regression case, which avoids overfitting and is coherent with the way ecologists think about the random effect component, that is as an additional assumption required only when covariates are not enough to explain variability in the data. Second, the PC prior for $\rho$ is easy-to-elicit given a prior statement on the median ICC. The ICC represents a highly intuitive scale to quantify the distance from the base model. Third, since PC priors are invariant over reparametrization, such user-defined-scaling in terms of the ICC can be applied in general to any group model, regardless of the interpretation and the scale of $\rho$ in the group model itself. Fourth, we show in the application that there is clear advantage in using these PC priors in a model comparison setting: if the same $\lambda$ is used for all the compared group models, the impact of the prior in the compared marginal likelihoods is the same. For this reason, the Bayes factor becomes a convenient tool to compare alternative mixed models, varying according to both the grouping factor and the group model.

The proposed approach can be extended to the generalized case of a response variable belonging to the exponential family. By including an iid Gaussian term in the linear predictor (reflecting a measurement error), derivation of the PC prior follows straightforwardly. As future work, other group models will be considered like the Matérn covariance, with the spatial range playing the role of $\rho$. The PC prior for the spatial range can be derived numerically, by calculating the matrix determinant and its derivative in (7); to gain computational efficiency, it is convenient to work with sparse precision matrices exploiting the SPDE approach by Lindgren et al., (2011). Finally, an important extension of this work would be to consider more than one factor. This requires working out joint PC priors for more than one correlation parameter, in a group model with nested random effects.

Acknowledgements

Daniela Cocchi and Massimo Ventrucci are supported by the PRIN 2015 grant project n.20154X8K23 (EPHASTAT) founded by the Italian Ministry for Education, University and Research. Gemma Burgazzi is supported by the project PRIN NOACQUA – responses of communities and ecosystem processes in intermittent rivers a National Relevant Project funded by the Italian Ministry of Education and University (PRIN 2015, Prot. 201572HW8F). The authors thank Maria Franco Villoria and Hȧvard Rue for the stimulating comments received about this work.

Appendix A Proofs of results in Section 4.2

Recall the definition of PC prior as an exponential distribution on the distance $d(\rho)$, with rate parameter $\lambda$,

[TABLE]

If design is balanced then $m_{j}=m,\forall j=1,\ldots,n$; recall that $n$ is the number of groups while $m$ is the number of within group observations. In this case, the distance function in Eq. (5) simplifies to

[TABLE]

Fixing $\lambda=\lambda^{\prime}/\sqrt{n}$, the PC prior for $\rho$ results (by the change of variable rule)

[TABLE]

Below, the PC priors in Eq. (8), (9) and (10) are derived. In each case, the proof is completed by deriving the analytical expression for the term $\frac{1}{2}|\bm{R}(\rho)|^{-1}\left|\frac{\partial|\bm{R}(\rho)|}{\partial\rho}\right|$ and plugging it in (14).

Exchangeable

Proof of Eq. (8).

Let us consider the compound symmetric matrix $\bm{R}(\rho)$ as in (3), where subscript $j$ is removed as we are working under a balanced design. Riebler et al., (2012) showed that

[TABLE]

hence the distance function is equal to $d(\rho)=\sqrt{-n\log\left\{(1+(m-1)\rho)(1-\rho)^{m-1}\right\}}$. The derivative term in (14) is

[TABLE]

After some algebraic steps, we obtain

[TABLE]

which completes the proof. ∎

Autoregressive of order one

Proof of Eq. (9).

The PC prior for the lag-one correlation of an AR1 is derived by Sørbye and Rue, (2017). Here we extend it to group models having within group correlation matrix $\bm{R}(\rho)$ as in (4). It can be shown that

[TABLE]

where $|\bm{P}|=1-\rho^{2}$. Thus the determinant of the AR1 correlation matrix is

[TABLE]

hence the distance function is equal to $d(\rho)=\sqrt{-n(m-1)\log(1-\rho^{2})}$. The derivative term in (14) is

[TABLE]

After some algebraic steps, we obtain

[TABLE]

which completes the proof. ∎

Ornstein Uhlenbeck

Proof of Eq. (10).

This proof follows straightforwardly from the AR1 case, by recognizing that $\phi=-\log(\rho)$, hence $\rho=\exp(-\phi)$. In this case, the determinant is

[TABLE]

and the distance function is equal to $d(\phi)=\sqrt{-n(m-1)\log(1-\exp(-2\phi))}$. The derivative term in (14) is

[TABLE]

After some algebraic steps, we obtain

[TABLE]

which completes the proof. ∎

Appendix B Example code for an exchangeable group model

Assume data vector y and covariate vector x of length $N=mn$. Observations are clustered in $n$ groups (indexed by $j=1,\ldots,n$), the $j^{th}$ group having $m_{j}$ observations. Assume the group model for the residuals $\bm{\theta}$

[TABLE]

where $\tau=1/\sigma^{2}$ (prec in the code) and $\bm{R}_{j}$ is the exchangeable correlation matrix, with correlation $\rho$ (rho in the code). The grouping factor is group = rep(1:n,each=m) (although the code below works for a general unbalanced design). First, we compute the PC prior for $\rho$, scaling it based on a median ICC equal to 0.5.

In R-INLA it is possible to implement the group model either using the control.group feature or building it “manually” with inla.rgeneric.define(). For illustrative purposes, we show the rgeneric option; the code below can be adapted to other group models. The group model is specified inside the function myrgeneric.exch, by coding the precision matrix (Q), the log-prior (log.prior) and the normalizing constant (log.norm.const). The normalizing constant of $\pi(\bm{\theta}|\tau,\rho)$, for the exchangeable case, is (assuming a generic unbalanced design)

[TABLE]

Finally, the INLA call: