Mixed models for biodiversity analysis - Modelli misti per l'analisi della biodiversità: The LME modelIl modello lineare ad effetti misti (Il modello LME)

The LME model

The Linear Mixed Effects model has the following general form: \[ \mathbf D_i = \mathbf X_i \boldsymbol{\beta} + \mathbf Z_i \mathbf b_i + \boldsymbol{\epsilon}_i \] where $\mathbf D_i$ contains $n_i$ values of a biodiversity index for the $i$-th community or $n_i$ repeated measures for the $i$-th sampling unit ($i=1,\ldots,I$). In this model both the fixed term $\mathbf X_i \boldsymbol{\beta}$ and the random term $\mathbf Z_i \mathbf b_i$ contain explanatory variables. The term $\mathbf Z_i \mathbf b_i$ is the community specific random effect and $\mathbf X_i \boldsymbol{\beta}$ is the overall (meta-community) fixed component. Because we have a fixed and a random component, we call the model a mixed effects model. $\mathbf X_i$ and $\mathbf Z_i$ are design matrices of dimension $n_i\times p$ and $n_i\times q$, respectively, where $n_i$ is the number of elements in $\mathbf D_i$ (the number of observations per community), $p$ the number of explanatory variables in $\mathbf X_i$, and $q$ the number of explanatory variables in $\mathbf Z_i$. Explanatory variables are all numeric and qualitative biotic and abiotic factors useful to characterize the biodiversity of the local communities. The matrices $\mathbf X_i$ and $\mathbf Z_i$ may, or may not, contain the same explanatory variables, depending on the type of model. Further model assumptions include: $$ \begin{eqnarray*} \mathbf b_i & \sim & N(\boldsymbol{0},\mathbf{D})\\ \boldsymbol{\epsilon}_i & \sim & N(\boldsymbol{0},\boldsymbol{\Sigma}_i) \end{eqnarray*} $$ $$ \mathbf b_1,\ldots, \mathbf b_I,\boldsymbol{\epsilon}_1,\ldots,\boldsymbol{\epsilon}_I\mbox{ independent} $$
A mixed effects model that does not contain any explanatory variable in the fixed term, but only an intercept $\beta_0$, is called a random effects model. For the $i$-th community (with $i=1,\ldots,I$), such a model is defined as follows: $$ \left( \begin{array}{c} D_{i,1} \\ D_{i,2} \\ \vdots \\ D_{i,n_i} \\ \end{array} \right) = \left( \begin{array}{c} 1\\ 1\\ \vdots \\ 1\\ \end{array} \right) \times \beta_0 + \left( \begin{array}{c} 1 \\ 1 \\ \vdots \\ 1 \\ \end{array} \right) \times b_{i0} + \left( \begin{array}{c} \epsilon_{i,1} \\ \epsilon_{i,2} \\ \vdots \\ \epsilon_{i,n_i} \\ \end{array} \right) $$ The biodiversity is modeled as a common intercept $\beta_0$ plus a community specific random term $b_{i0}$. For $i=1,\ldots,I$ we assume that
- the terms $b_{i0}$ are normally distributed with mean 0 and variance $d^2$,
- the $n_i$-dimensional error vectors $\boldsymbol{\epsilon}_i$ are normally distributed with covariance matrix $\boldsymbol{\Sigma}_i$.
The easiest definition of $\boldsymbol{\Sigma}_i$ is the following $$ \boldsymbol{\Sigma}_i=\sigma^2\left( \begin{array}{ccccc} 1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 1 \\ \end{array} \right) $$ that corresponds to independent and identically distributed error terms. In general, the elements of $\boldsymbol{\Sigma}_i$ do not depend on $i$, but this may not always be the case as there are various methods to incorporate heterogeneity into the model that influence the structure of $\boldsymbol{\Sigma}_i$.
Suppose that we want to model biodiversity as a linear function of a numeric explanatory variable $X$ and that the intercept of the linear function is allowed to change with the community. One option is to use the community as a random effect, including a community effect in the model and assuming that the variation of the intercept is normally distributed with a certain variance. A small variance means that differences per community (in terms of the intercept) are small, whereas a large variance allows for more variation. Such a mixed effects model is called random intercept model and is defined as follows: $$ \left( \begin{array}{c} D_{i,1} \\ D_{i,2} \\ \vdots \\ D_{i,n_i} \\ \end{array} \right) = \left( \begin{array}{cc} 1 & X_{i,1} \\ 1 & X_{i,2} \\ \vdots \\ 1 & X_{i,n_i} \\ \end{array} \right) \times \left( \begin{array}{c} \beta_0 \\ \beta_1 \end{array} \right) + \left( \begin{array}{c} 1 \\ 1 \\ \vdots \\ 1 \\ \end{array} \right) \times b_{i0} + \left( \begin{array}{c} \epsilon_{i,1} \\ \epsilon_{i,2} \\ \vdots \\ \epsilon_{i,n_i} \\ \end{array} \right) $$ In this case for $i=1,\ldots,I$ we assume that
- item the random effects $b_{i0}$ are normally distributed: $N(0, d^2)$,
- item the elements of the $n_i$-dimensional vectors $\boldsymbol\epsilon_i$ are independent and normally distributed with mean 0 and variance $\sigma^2$.
The random intercept model allows for a random shift around the intercept resulting in (community specific) fitted lines parallel to an average (meta-community) fitted line. The model implies one meta-community regression line with intercept $\beta_0$ and slope $\beta_1$ that is allowed to be shifted up, or down, for each community by the quantity $b_{i0}$ that is normally distributed with a certain variance $d^2$. If $d^2$ is large, the vertical shifts will be large. If $d^2 = 0$, all shifts are zero and coincide with the meta-community regression line.
Suppose that the relationship between biodiversity and the explanatory variable $X$ is different at each community, i.e. that the factor $X$ has a different influence on biodiversity at each community. The mixed effects model with a random intercept and a random slope can then be applied. Such a mixed effects model is called random coefficients model and is defined as follows: $$ \left( \begin{array}{c} D_{i,1} \\ D_{i,2} \\ \vdots \\ D_{i,n_i} \\ \end{array} \right) = \left( \begin{array}{cc} 1 & X_{i,1} \\ 1 & X_{i,2} \\ \vdots \\ 1 & X_{i,n_i} \\ \end{array} \right) \times \left( \begin{array}{c} \beta_0 \\ \beta_1 \end{array} \right) + \left( \begin{array}{cc} 1 & X_{i,1} \\ 1 & X_{i,2} \\ \vdots \\ 1 & X_{i,n_i} \\ \end{array} \right) \times \left( \begin{array}{c} b_{i0} \\ b_{i1} \end{array} \right) + \left( \begin{array}{c} \epsilon_{i,1} \\ \epsilon_{i,2} \\ \vdots \\ \epsilon_{i,n_i} \\ \end{array} \right) $$ In this case the random effects $\mathbf b_i$ are assumed to be distributed as: $$ \mathbf b_i = \left( \begin{array}{c} b_{i0} \\ b_{i1} \end{array} \right) \sim N(0, \mathbf D) \;\mbox{where}\; \mathbf D = \left( \begin{array}{cc} d_{11}^2 & d_{12} \\ d_{12} & d_{22}^2 \\ \end{array} \right) $$ The variance $d_{11}^2$ determines the amount of variation around the meta-community intercept $\beta_0$. The model also allows for random variation around the meta-community slope $\beta_1$ in a similar way as it does for the intercept. The variance $d_{22}^2$ determines the variation in slopes at the $I$ communities. Finally, there is the correlation between the random intercepts and slopes, $d_{12}$.
When applying mixed effects modeling, parameter estimation can be done by ML and REML estimation, standing for maximum likelihood and restricted expected maximum likelihood respectively. ML and REML estimators for the parameters in $\boldsymbol{\beta}$ are not identical. If the number of fixed explanatory variables is small relative to the number of observations, there are not big differences, but for models with many fixed terms, this may not be the case.
Model selection in mixed effects modeling is generally based on the Akaike Information Criterion (AIC), or the Bayesian Information Criterion (BIC). Both the AIC and BIC are the sum of two terms that measure the fit of the model to the data and the complexity of the model. In the formulae below $\theta$ is a set of $p$ parameters including elements of $\boldsymbol{\beta}$, $\mathbf{D}$ and $\boldsymbol{\Sigma}_i$ for $i=1,\ldots,I$, $L$ is either the ML or REML likelihood and $n^* = n$ for ML, while $n^* =n - p$ for REML. $$ \begin{eqnarray*} AIC & = & -2\times L(\theta) + 2\times p\\ BIC & = & -2\times L(\theta) + 2\times p\times \ln(n^*) \end{eqnarray*} $$ In general, the lowest the values of the two criteria, the better the estimated model. Notice that an AIC based on REML is not comparable with an AIC obtained by ML and the same holds for the BIC. A second approach to find the optimal model is via hypothesis testing, with three options: the t-statistic, the F-statistic, or the likelihood ratio test.
The validation of an estimated model can be done by normalized residuals based on the REML fit. Residuals should be plotted against fitted values to identify violations of homogeneity, indicated by differences in spread. An increase in spread for larger fitted values can be addressed in several ways: (i) apply a transformation, (ii) check whether the increase in spread is due to one of the explanatory variables, and (iii) apply generalized linear mixed model with a Poisson assumption (if the biodiversity is measured as the number of species). If the increase in spread is due to an explanatory variable, then this source of heterogeneity should be incorporated into the model. Residuals should also be plotted against each explanatory variable. If these plots show a pattern in the spread or a pattern in the residuals, one should consider adding more explanatory variables, interactions or quadratic terms. Histograms and q-q plots of the residuals can be used to verify normality. Some software packages also provide normality tests like the Shapiro-Wilks test.

Il modello lineare ad effetti misti (il modello LME)

Il modello lineare di effetti misti (LME) ha questa forma generale: \[ \mathbf D_i = \mathbf X_i \boldsymbol{\beta} + \mathbf Z_i \mathbf b_i + \boldsymbol{\epsilon}_i \] dove $\mathbf D_i$ contiene $n_i$ valori dell'indice di biodiversitá per l'$i$-esima comunitá o $n_i$ misure ripetute per $i$-esima unitá di campionamento ($i=1,\ldots,I$). In questo modello entrambi i termini fissi $\mathbf X_i \boldsymbol{\beta}$ e il termine casuale $\mathbf Z_i \mathbf b_i$ contiene le variabili esplicative. Il termine $\mathbf Z_i \mathbf b_i$ é l'effetto casuale specifico per la comunitá e $\mathbf X_i \boldsymbol{\beta}$ il componente fisso della meta-comunitá. Siccome sono presenti una componente fisso e una componente casuale, il modello si chiamerá modello ad effetti misti. $\mathbf X_i$ e $\mathbf Z_i$ sono matrici di design di dimensioni $n_i\times p$ e $n_i\times q$, rispettivamente, dove $n_i$ é il numero di elementi in $\mathbf D_i$ (numero di osservazioni per comunitá), $p$ il numero delle variabili esplicative in $\mathbf X_i$, e $q$ il numero delle variabili esplicative in $\mathbf Z_i$. Le variabili esplicative sono fattori numerici e qualitativi biotici e abiotici utili per caratterizzare la biodiversitá delle comunitá locali. Le matrici $\mathbf X_i$ e $\mathbf Z_i$ possono, o non possono, contenere le stesse variabili esplicative, dipendenti dal tipo del modello. Ulteriori assunzioni per il modello includono: $$ \begin{eqnarray*} \mathbf b_i & \sim & N(\boldsymbol{0},\mathbf{D})\\ \boldsymbol{\epsilon}_i & \sim & N(\boldsymbol{0},\boldsymbol{\Sigma}_i) \end{eqnarray*} $$ $$ \mathbf b_1,\ldots, \mathbf b_I,\boldsymbol{\epsilon}_1,\ldots,\boldsymbol{\epsilon}_I\mbox{ independent} $$
Un modello a effetti misti che non contiene alcuna variabile esplicativa nel valore fisso, ma soltanto un tercetto $\beta_0$, é detto modello a a effeti casuali (random effects model). Per la comunitá $i$-esima (con $i=1,\ldots,I$), un tale modello é definito come segue: $$ \left( \begin{array}{c} D_{i,1} \\ D_{i,2} \\ \vdots \\ D_{i,n_i} \\ \end{array} \right) = \left( \begin{array}{c} 1\\ 1\\ \vdots \\ 1\\ \end{array} \right) \times \beta_0 + \left( \begin{array}{c} 1 \\ 1 \\ \vdots \\ 1 \\ \end{array} \right) \times b_{i0} + \left( \begin{array}{c} \epsilon_{i,1} \\ \epsilon_{i,2} \\ \vdots \\ \epsilon_{i,n_i} \\ \end{array} \right) $$ La biodiversitá é modellata come una intercettazione comune $\beta_0$ piú un termine casuale specifico per la comunitá $b_{i0}$. Per $i=1,\ldots,I$ si assume che
- i termini $b_{i0}$ sono distributi normalmente coon una media 0 e varianza $d^2$,
- I vettori d'errore $n_i$-dimensionali$\boldsymbol{\epsilon}_i$ sono distributiti normalmente con una matrice di covarianza $\boldsymbol{\Sigma}_i$.
La definizione piú semplice di $\boldsymbol{\Sigma}_i$ é la seguente $$ \boldsymbol{\Sigma}_i=\sigma^2\left( \begin{array}{ccccc} 1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 1 \\ \end{array} \right) $$ che corrisponde al termine d'errore indipendente e distribuito indeticamente (i.i.d.). In generale, gli elementi di $\boldsymbol{\Sigma}_i$ non dipendono da $i$, ma non sempre é questo il caso, siccome vi sono altri metodi per incorporare l'eterogeneitá nel modello che influenza la struttura di $\boldsymbol{\Sigma}_i$.
Si supponga voler modellare biodiversitá come una funzione lineare di una variabile esplicativa numerica $X$ e che l'intercettazione della funzione lineare puó cambiare con la comunitá. Una opzione é quella di usare la comunitá come un effetto casuale, incluso un effetto di comunitá nel modello e assumendo che la variazione delle intercettazioni é distribuita normalmente con una certa varianza. Una varianza piccola significa che le differenze per comunitá (in termini di intercettazioni) sono piccole, mentre una ampia varianza permette variazioni maggiori. Un tale modello a effetti misti é detto modello a intercettazioni casuali ed é definito come: $$ \left( \begin{array}{c} D_{i,1} \\ D_{i,2} \\ \vdots \\ D_{i,n_i} \\ \end{array} \right) = \left( \begin{array}{cc} 1 & X_{i,1} \\ 1 & X_{i,2} \\ \vdots \\ 1 & X_{i,n_i} \\ \end{array} \right) \times \left( \begin{array}{c} \beta_0 \\ \beta_1 \end{array} \right) + \left( \begin{array}{c} 1 \\ 1 \\ \vdots \\ 1 \\ \end{array} \right) \times b_{i0} + \left( \begin{array}{c} \epsilon_{i,1} \\ \epsilon_{i,2} \\ \vdots \\ \epsilon_{i,n_i} \\ \end{array} \right) $$ In questo caso per $i=1,\ldots,I$ si assume che
- gli effetti casuali $b_{i0}$ sono distribuiti normalmente: $N(0, d^2)$,
- gli elementi del vettore $n_i$-dimensionale $\boldsymbol\epsilon_i$ sono indipendenti e distribuite normalmente con media 0 e varianza $\sigma^2$.
il modello delle intercettazioni casuali consente per spostamenti casuali intorno all'intercettazione risultante in (specifici per la comunitá) linee parallele a una linea media (meta-comunitá). Il modello comporta una linea di regressione con una intercetazzione $\beta_0$ ed una inclinazione $\beta_1$ che puó essere alzata, o abbasata, per tutte le comunitá di una quantiá $b_{i0}$ che é normalmente distribuita con una certa varianza $d^2$. Se $d^2$ é grande, gli spostamenti verticali sono ampi. Se $d^2 = 0$, tutti gli spostamenti sono zero e coincidono con la linea di regressione della meta-comunitá.
Supposto che le relazioni tra biodiversitá e la variabile esplicativa $X$ é differente per ogni comunitá , cioè il fattore $X$ ha una differente influenza sulla biodiversitá per ogni comunitá. Si puó,quindi, applicare il modello a effetti misti con intercettazioni casuali e inclinazioni casuali puó. Tali modelli a effetti misti sono detti modelli a coefficienti casuali e sono così definiti: $$ \left( \begin{array}{c} D_{i,1} \\ D_{i,2} \\ \vdots \\ D_{i,n_i} \\ \end{array} \right) = \left( \begin{array}{cc} 1 & X_{i,1} \\ 1 & X_{i,2} \\ \vdots \\ 1 & X_{i,n_i} \\ \end{array} \right) \times \left( \begin{array}{c} \beta_0 \\ \beta_1 \end{array} \right) + \left( \begin{array}{cc} 1 & X_{i,1} \\ 1 & X_{i,2} \\ \vdots \\ 1 & X_{i,n_i} \\ \end{array} \right) \times \left( \begin{array}{c} b_{i0} \\ b_{i1} \end{array} \right) + \left( \begin{array}{c} \epsilon_{i,1} \\ \epsilon_{i,2} \\ \vdots \\ \epsilon_{i,n_i} \\ \end{array} \right) $$ In questo caso gli effetti casuali $\mathbf b_i$ si assume siano distribuiti come: $$ \mathbf b_i = \left( \begin{array}{c} b_{i0} \\ b_{i1} \end{array} \right) \sim N(0, \mathbf D) \;\mbox{where}\; \mathbf D = \left( \begin{array}{cc} d_{11}^2 & d_{12} \\ d_{12} & d_{22}^2 \\ \end{array} \right) $$ La varianza $d_{11}^2$ determina la il valore di variazioni intorno alle intercettazioni della meta-comunitá $\beta_0$. Il modello consente anche variazioni casuali intorno alle inclinazioni della meta-comunitá in modo simile a quello per le intercettazioni. La varianza $d_{22}^2$ determina la variazione nelle inclinazioni per le comunitá $I$. Infine vi é una corelazione fra le intercettazioni casuali e le inclinazioni (slopes), $d_{12}$.
Quando si utilizza la modellazione ad effetti misti, la valutazione dei parametri puó essere effettuata con le stime ML e REML, ossia maximum likelihood (Massima verosimiglianza) e restricted expected maximum likelihood (Massima verosimiglianza ristretta) rispettivamente. Gli stimatori ML e REML per i parametri in $\boldsymbol{\beta}$ non sono identici. Se il numero delle variabili esplicative e piccolo in relazione al numero delle osservazioni, non vi sono grosse differenze, ma per modelli con molti termini fissi, questo potrebbe non essere il caso.
Le selezione del modello nella modellazione ad effetti misti generalmente é basato sul Criterio di Informazione di Akaike (Akaike Information Criterion - AIC), o il Criterio di Informazione di Bayesia (Bayesian Information Criterion - BIC). Entrambi i criteri, AIC e BIC, sono le somme di due termini che misurano l'adeguatezza del modello ai dati e la complessitá del modello. Nelle sottoscritte formule $\theta$ vi é un set di parametri $p$ che includono elementi di $\boldsymbol{\beta}$, $\mathbf{D}$ e $\boldsymbol{\Sigma}_i$ per $i=1,\ldots,I$, $L$ é la verosimiglianz ML o REML e $n^* = n$ per ML, mentre $n^* =n - p$ per REML. $$ \begin{eqnarray*} AIC & = & -2\times L(\theta) + 2\times p\\ BIC & = & -2\times L(\theta) + 2\times p\times \ln(n^*) \end{eqnarray*} $$ In generale, piú bassi sono i valori dei due criteri, migliore é il modello stimato. Si noti che un AIC basato su REML non é paragonabile con un AIC ottenuto da ML. Lo stesso dicasi per il BIC.
Un secondo approccio per determinare il modello ottimale é attraverso i test di ipotesi, con tre opzioni: la statistica-t, la statistica-F e il test del rapporto di verosomiglianza.
La validazione dei modelli stimati puó essere effettuata dai residui normalizzati basati sul REML. I residui dovrebbero essere disegnati insieme ai valori adattati in modo da identificare violazioni di omogeneitá, indicati dalle differenze nella diffusione. Un aumento nella variazione (spread) per grossi valori puó essere affrontato in diversi modi: (i) applicare una trasformazione, (ii) verificare se l'aumento della variazione dipende da una delle variabili esplicative, e (iii) applicare il modello misto lineare generalizzato con un'assunzione di Poisson (la biodiversitá é vista come numero di specie). Se l'aumento dello spread é causato da una variabile esplicativa, allora questa fonte di eterogenitá dovrebbe essere incorporata nel modello. I residui dovrebbero anche essere disegnate in relazione alle variabili esplicative. Se questi punti indicano uno schema nello spread o nei residui, si potrebbe pensare di aggiungere piú variabili esplicative, interazioni o termini quadratici. Istogrammi e q-q plots dei residui possono essere usati per verificare la normalitá. Alcuni pacchetti software forniscono test di normalitá come il test di Shapiro-Wilks.