library(tibble)
library(dplyr)
library(ggplot2)

Introducció

Els models mixts són són una expansió dels models lineals que permet tenir en compte la possible falta d’independència entre les observacions d’un estudi estadístic. Sovint pot passar que les característiques del sistema que es vol estudiar o del disseny experimental o mostral impliquen incompliment del supòsit d’independència, com en el cas de mesures repetides sobre el mateix individu, la mateixa comunitat o ecosistema, etc. El supòsti d’independència és bàsics en els models LM, GLM i GAM.

Els models mixts tenen en compte aquesta no independència dels errors quan modelen l’estructura de la covariància, introduïda per agrupacions de dades.

Es diuen així perquè inclouen dos tipus de variables explicatives que fan referència a efectes fixos i efectes aleatoris. Els efectes fixos comprenen totes la variables contínues i els factors que hagin estat definits a priori.

Un factor fix és un factor d’interès real. Les condicions d’aquest factor no es trien arbitràriament, sinó que són les que creieu que haurien de tenir un efecte sistemàtic sobre la vostra variable dependent.

Els factors aleatoris, en canvi, són factors que creieu que poden tenir un paper, però que no són d’interès per se. En incloure’ls al vostre model, voleu controlar la seva variància no sistemàtica/aleatòria.

Definició de factor fix i factor aleatori

Primer de tot, qualsevol factor de 4 o menys nivells ha de ser considerat factor fix.

Els efectes del factor són allò que determina si un factor és fix o aleatori:

Els efectes del factor són informatius si ens aporten informació de cara a l’objecte del nostre estudi. El factor és per tant fix.
Els efectes no són informatius, son circumstancials respecte a l’objecte del nostre estudi. El factor és aleatori.

Formulació teòrica del model d’efectes mixts

Com sabem, el model de Regressió lineal múltiple suposa considerar una variable de resposta \(\small Y\) sobre què influeixen \(\small k\) variables explicatives (o independents o covariables), de la forma \[ \small Y= \beta_0+\beta_1X_1+...+\beta_kX_k+e \] on \(\small e\) és una variable d’error amb distribució normal \(\small N(0,\sigma^2)\), en una mostra aleatòria d’\(\small n\) observacions.

si denominem \(\small \textbf{y}=(y_1,...,y_n)^t\) el vector de les observacions de la variable de resposta, integrem tots els coeficients en un vector \(\small \textbf β=(\beta_0,...,\beta_k)^t\) i denominem \(\small \textbf{X}\) la matriu \(\small n\times(k+1)\)¹ de les observacions de les variables explicatives (matriu del disseny), el model de regressió lineal es pot expressar així \[ \small \textbf y = \textbf{Xβ}+\bf e \] on \(\small \textbf e=(e_1,..., e_n)^t\) es el vector d’errors, amb distribució normal multivariant n-dimensional de vector de mitjanes \(\small \textbf 0=(0,...0)^t\) i matriu de covariàncies \(\small \textbf I\sigma^2\) (sent \(\bf I\) la matriu d’identitat).

En el Model d’efectes mixts generalitzarem l’anterior, explicant \(\small\bf y\) de la següent manera: \[ \small \textbf y = \bf Xβ+Zb+e \]

on el vector \(\small \bf β\) representa els efectes fixos, però on hem afegit el vector \(\small \bf b\), que representa els efectes aleatoris, amb distribució normal multivariant, p-dimesional de vector de mitjanes \(\small \textbf 0=(0,...,0)^t\) i matriu de covariàncies \(\small \bf \psi_\theta\), on \(\small\bf Z\) és la matriu del model per als efectes aleatoris i on \(\small \textbf e= (e_1,...,e_n)^t\) és el vector d’errors, ara amb distribució normal multivariant, n-dimensional amb vector de mitjanes \(\small \textbf 0=(0,...,0)^t\) i matriu de covariàncies \(\small \textbf Λσ^2\) i on els vectors \(\small \bf b\) i \(\small \bf e\) són independents.

Referències

L. Cayuela, M. De la Cruz, Análisis de datos ecológicos en R, 2022.
S.N. Wood, Generalized Additive Models, - An Introduction with R, cap. 2, 2017.
A. García Pérez, Estadística aplicada avanzada con R, cap. 17, 2025.

Models mixts amb R: `nlme` i `nlme4`

EnR hi ha diferent paquets que poden executar models mixtos. Tot i així, els més coneguts i utilitzats són nlme i lme4. Ambdós disposen funcions per ajustar models lineals mixts, però només lme4 pot ajustar models GLMM.

Aquí, només utilitzarem lme4 que és més complet.

La funció `lmer()`

En la funció lmer(), la fórmula del model s’especifica en un sol argument que inclogui efectes fixos i aleatoris. El efectes fisos s’especifiquen primers a la dreta del tilde. Seguidament es col·loquen els efectes aleatoris anteposant el símbol + i entre parèntesis. R identificarà com a efectes aleatoris tot allò que vingui a la dreta de la barra vertical |.

Model amb un factor aleatori

\[ \small \tt lmer(y\sim 1 +(1|a)) \]

Model amb un factor aleatori amb efecte sobre un factor fix o covariable

\[ \small \tt lmer(y\sim x +(x|a)) \] o, que és equivalent \[ \small \tt lmer(y\sim x +(1+x|a)) \]

Model amb més d’un factor aleatori

Si hi ha més d’un factor aleatori, s’han d’identificar separadament entre parèntesi \[ \small \tt lmer(y\sim x + (x|a)+(1|b)+(1|c)) \]

Model amb factors anidats

\[ \small \tt lmer(y\sim 1 +(1|a/b/c)) \] o bé \[ \small \tt lmer(y\sim 1 +(1|a:b:c)) \]

Entendre els models lineals mixts (LMM) amb un exemple

Utilitzarem dades sobre bentos marí² procedents de nou zones intermareals holandeses. Les dades van ser recollides pel RIKZ neerlandès l’estiu de 2002. A cada zona es van recollir 5 mostres de la macrofauna i de variables abiòtiques, seguint sempre el mateix esquema.

El dataset `RIKZ`

Les dades són al dataset RIKZ de la llibreria ADER

library(ADER)
data(RIKZ)

str(RIKZ)

## 'data.frame':    45 obs. of  89 variables:
##  $ Sample       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ C1           : num  4 0 0 0 1 0 0 0 0 0 ...
##  $ P1           : num  0 0 0 0 0 0 0 0 0 2 ...
##  $ P2           : num  0 1 3 0 0 0 0 0 0 0 ...
##  $ P3           : num  0 0 0 0 0 0 1 0 21 2 ...
##  $ P4           : num  0 0 0 0 0 0 1 0 11 0 ...
##  $ P5           : num  0 0 0 0 0 1 0 0 3 0 ...
##  $ P6           : num  1 0 0 0 0 0 0 0 0 0 ...
##  $ P7           : num  0 0 0 0 0 0 0 0 34 7 ...
##  $ P8           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ P9           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ P10          : num  0 2 1 0 1 0 0 0 0 2 ...
##  $ P11          : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ P12          : num  0 0 0 0 0 1 3 0 14 0 ...
##  $ P13          : num  0 0 0 0 0 0 0 5 16 0 ...
##  $ P14          : num  0 0 0 0 0 889 21 10 49 0 ...
##  $ P15          : num  4 0 0 0 0 0 0 0 1 0 ...
##  $ P16          : num  32 22 16 79 15 0 3 7 1 0 ...
##  $ P17          : num  0 0 0 0 0 0 0 0 0 16 ...
##  $ P18          : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ P19          : num  0 0 1 0 0 0 0 0 1 0 ...
##  $ P20          : num  0 0 0 0 0 0 0 0 0 2 ...
##  $ P21          : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ P22          : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ P24          : num  0 0 0 92 0 0 0 0 0 0 ...
##  $ P25          : num  50 12 7 1 0 0 0 0 0 0 ...
##  $ N1           : num  0 0 0 0 0 0 0 1 0 0 ...
##  $ CR1          : num  2 0 0 0 0 0 0 0 0 0 ...
##  $ CR2          : num  0 0 0 2 0 0 0 0 0 0 ...
##  $ CR3          : num  0 1 2 0 25 1 0 0 0 0 ...
##  $ CR4          : num  0 8 25 0 13 0 0 0 0 0 ...
##  $ CR5          : num  0 0 0 1 1 0 175 163 0 0 ...
##  $ CR6          : num  0 0 0 1 0 0 0 0 0 0 ...
##  $ CR7          : num  0 0 1 0 0 0 0 0 0 0 ...
##  $ CR8          : num  0 0 0 0 0 0 0 0 0 1 ...
##  $ CR9          : num  0 0 0 1 4 0 0 0 0 0 ...
##  $ CR10         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ CR11         : num  0 0 0 0 0 0 0 0 2 0 ...
##  $ CR12         : num  0 0 0 0 0 25 35 4 17 0 ...
##  $ CR13         : num  0 0 0 0 0 0 0 0 1 0 ...
##  $ CR14         : num  30 3 2 4 2 0 0 1 1 1 ...
##  $ CR15         : num  9 0 0 15 4 0 1 1 0 0 ...
##  $ CR16         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ CR17         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ CR18         : num  0 0 2 0 0 0 0 0 0 1 ...
##  $ CR19         : num  8 0 0 0 0 0 0 0 0 0 ...
##  $ CR20         : num  0 0 2 0 0 0 0 0 0 0 ...
##  $ CR21         : num  0 1 7 0 1 0 0 0 0 0 ...
##  $ CR22         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ CR23         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ CR24         : num  0 0 0 0 0 0 0 0 0 1 ...
##  $ CR25         : num  0 1 0 0 0 0 0 0 0 0 ...
##  $ CR26         : num  2 0 0 1 0 0 0 0 0 0 ...
##  $ CR27         : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ CR28         : num  0 0 1 0 0 0 0 0 0 0 ...
##  $ M1           : num  0 0 0 0 0 0 0 0 2 1 ...
##  $ M2           : num  0 0 0 0 0 0 0 0 2 0 ...
##  $ M3           : num  0 0 0 0 0 0 0 0 0 1 ...
##  $ M4           : num  0 0 0 0 0 0 0 0 0 1 ...
##  $ M5           : num  0 0 0 0 0 0 0 0 0 1 ...
##  $ M6           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ M7           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ M8           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ M9           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ M10          : num  0 0 0 0 0 0 0 0 16 6 ...
##  $ M11          : num  0 0 0 0 0 0 1 0 18 0 ...
##  $ M12          : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ M14          : num  0 1 0 0 0 0 0 0 0 1 ...
##  $ M15          : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ M16          : num  0 0 0 0 0 0 0 0 0 2 ...
##  $ M17          : num  0 0 0 0 0 0 0 0 1 0 ...
##  $ I1           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ I2           : num  0 0 0 0 0 1 0 0 0 0 ...
##  $ I3           : num  1 0 0 2 0 24 0 0 0 0 ...
##  $ I4           : num  0 0 0 0 0 2 0 0 0 0 ...
##  $ I5           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ week         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ angle1       : int  32 62 65 55 23 129 126 52 26 143 ...
##  $ angle2       : int  96 96 96 96 96 89 89 89 89 89 ...
##  $ exposure     : int  10 10 10 10 10 8 8 8 8 8 ...
##  $ salinity     : num  29.4 29.4 29.4 29.4 29.4 29.6 29.6 29.6 29.6 29.6 ...
##  $ temperature  : num  17.5 17.5 17.5 17.5 17.5 20.8 20.8 20.8 20.8 20.8 ...
##  $ NAP          : num  0.045 -1.036 -1.336 0.616 -0.684 ...
##  $ penetrability: num  254 227 237 249 252 ...
##  $ grainsize    : num  222 200 194 221 202 ...
##  $ humus        : num  0.05 0.3 0.1 0.15 0.05 0.1 0.1 0.1 0.15 0 ...
##  $ chalk        : num  2.05 2.5 3.45 1.6 2.45 2.5 1.85 1.7 2.3 2.6 ...
##  $ sorting1     : num  69.8 59 59.2 67.8 57.8 ...
##  $ Beach        : int  1 1 1 1 1 2 2 2 2 2 ...

Com es pot veure, es tracta d’un dataset amb 89 variables: la primera columna indica la mostra les columnes de 2 a 76 són recomptes d’individus d’espècies (la capçalera indica alguna codificació de l’espècie). Les columnes de 77 a 89 representen valors de les covariables.

A nosaltres interessarà NAP, el nivell intermareal, una mesura de l’elevació de cada mostra respecte al nivell mitjà de la marea. En aquest exemple investigarem si existeix alguna relació entre riquesa d’espècies i el NAP.

Preparació de les dades

La riquesa d’espècies en aquest cas és en forma de simple recompte. Per tant el model a ajustar més adequat seria per tant un GLM amb distribució de Poisson. Tanmateix, amb l’ànim de simplificar les coses en aquesta fase de l’exemple, emprarem un model de regressió lineal amb distribució d’errors normal.

Prèviament, serà necessari retocar algun element en l’estructura de les dades, perquè sigui viable la nostra anàlisi.

Les dades sobre espècies venen com a recomptes d’individus (abundància) que hem de convertir en presència (0,1) per tal de poder fer sumes que ens donin el valor de la riquesa.

Com que ens agrada tidyverse, en farem servir les funcions, en compte del paquet base de R:

rikz <- RIKZ %>% 
  mutate(across(2:76, ~ ifelse(.> 0,1,0)),
         Beach = factor(Beach)) %>% 
  rowwise() %>% 
  mutate(riquesa= sum(c_across(2:76))) %>% 
  select(1,90,89,2:88)

D’aquesta manera hem afegit la variable riquesa d’espècies. També hem transformat la variable Beach de integer a factor, perquè ho R ho pugui entendre.

Aproximació amb model lineal

Cas 0: independència entre mostres

Si assumíssim que les mostres fossin totalment independents entre si, podríem expressar el model així: \[ \small Riquesa_i=\alpha+\beta_{}NAP_i+e_i~~~~~~e_i\sim N(0,\sigma^2) \] que quedaria representat per un sola recta

ml1 <- lm(riquesa ~ NAP, data = rikz)

rikz %>% 
  ggplot(aes(x=NAP, y= riquesa)) +
  geom_point() +
  geom_abline(intercept = ml1$coefficients[1], slope = ml1$coefficients[2])+
  labs(x= "Nivell intemareal (m)",
       y= "Riquesa de bentos") +
  theme_bw()

El model té 3 paràmetres desconeguts (els dos coeficients i la variància residual) i que assumeix que la relació riquesa i NAP sigui la mateixa per a totes les plathes (Beach). Però ja sabem que podria no ser així (veure Nota núm. 2 d’estadística aplicada a les ciències ambientals). Podríem tenir 3 opcions:

Cas 1: La pendent és la mateixa, però l’ordenada a l’origen és diferent.
Cas 2: L’ordenada a lorigen és la mateica i canvia la pendent.
Cas 3: Tant la pendent com l’ordenada a l’origen són diferents.

Cas 1: pendent constant

El model seria \[ \small Riquesa_{ij}=\alpha_{j}+\beta_{}NAP_{ij}+e_{ij}~~~~~~e_{ij}\sim N(0,\sigma^2) \]

Le prediccions quedarien representades per a 9 rectes paral·leles:

ml2 <- lm(riquesa ~ Beach + NAP, data = rikz)

new_data <- expand.grid(
  NAP = seq(min(rikz$NAP), max(rikz$NAP), by= 0.1),
  Beach = unique(rikz$Beach)
)
new_data$pred <- predict(ml2, newdata= new_data)

rikz %>% 
  ggplot(aes(x=NAP, y= riquesa, colour = Beach)) +
  geom_point() +
  geom_line(data = new_data, aes(x=NAP, y=pred, colour = Beach))+
  geom_vline(xintercept = 0, lty=3)+
  labs(x= "Nivell intermareal (m)",
       y= "Riquesa de bentos") +
  theme_bw() +
  theme(legend.position = "none")

Aquí representem els models lineals, suposant que, si bé la relació entre riquesa i NAP no variï entre platges, la riquesa mitjana, sí, pot variar.

Cas 2: intercept constant

En aquest cas el model seria \[ \small Riquesa_{ij}=\alpha_{}+\beta_{j}NAP_{ij}+e_{ij}~~~~~~e_{ij}\sim N(0,\sigma^2) \]

ml3 <- lm(riquesa ~ NAP + Beach:NAP, data = rikz)

new_data$pred <- predict(ml3, newdata= new_data)

rikz %>% 
  ggplot(aes(x=NAP, y= riquesa, colour = Beach)) +
  geom_point() +
  geom_line(data = new_data, aes(x=NAP, y=pred, colour = Beach))+
  geom_vline(xintercept = 0, lty=3)+
  labs(x= "Nivell intermareal (m)",
       y= "Riquesa de bentos") +
  theme_bw() +
  theme(legend.position = "none")

Aquí en canvi, representem els models lineals, suposant que la relació entre riquesa i NAP variï entre platges, mentre que la riquesa mitjana es manté fixa.

Cas 3: intercept i pendent variables

En aquest cas el model seria \[ \small Riquesa_{ij}=\alpha_{j}+\beta_{j}NAP_{ij}+e_{ij}~~~~~~e_{ij}\sim N(0,\sigma^2) \]

ml4 <- lm(riquesa ~ Beach*NAP, data = rikz)

new_data$pred <- predict(ml4, newdata= new_data)


rikz %>% 
  ggplot(aes(x=NAP, y= riquesa, colour = Beach)) +
  geom_point() +
  geom_line(data = new_data, aes(x=NAP, y=pred, colour = Beach))+
  geom_vline(xintercept = 0, lty=3)+
  labs(x= "Nivell intermareal (m)",
       y= "Riquesa de bentos") +
  theme_bw() +
  theme(legend.position = "none")

Aquest tercer model ja assumeix la variabilitat entre platges de la relació entre variables i del valor de l’ordenada a l’origen.

Conclusions

El model d’una sola recta és el més senzill i els darrer el més complex. El nombre de paràmetres a estimar així com el valor de \(\small R^2\) en els 4 casos seria

Cas	Model	Paràmetres	\(\small R^2\)
Cas 0	riquesa ~ NAP	3	0.3245413
Cas 1	riquesa ~ Beach + NAP	11	0.7024972
Cas 2	riquesa ~ NAP + Beach:NAP	11	0.3652364
Cas 3	riquesa ~ Beach*NAP	19	0.8493863

Clara ment el model més complex és el que més explica la variabilitat, però és el model més costós.

D’altra banda, si comprovem la significació dels diversos termes implicats en el model més complex emprant les sumes de quadrats de tipus I (seqüencial),

anova(ml4)

## Analysis of Variance Table
## 
## Response: riquesa
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Beach      8 543.24  67.906 11.0500 9.017e-07 ***
## NAP        1 230.66 230.658 37.5341 1.517e-06 ***
## Beach:NAP  8 161.82  20.227  3.2915  0.009434 ** 
## Residuals 27 165.92   6.145                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

veiem que tots els termes, interacció Beach-NAP inclosa, són significatius.

Però nosaltres estem interessats en l’efecte entre riquesa i NAP, i no tant en l’efecte platja. Tanmateix, eliminant la platja de l’anàlisi, diferències lligades a aquest terme passarien a formar part de la variància residual. A més a més deixaríem de tenir en compte que les observacions no són independents entre si, que podria afectar a l’error estàndard del model en cascada els valors p. Això podria comportar mals resultats com la detecció d’una relació riquesa-NAP que en realitat sigui inexistent. Per tant, cal introduir el efecte platja al model: això implica perdre 16 graus de llibertat amb pèrdua de potència també pot comportar no detectar relacions significatives entre variables (error de tipus II).

L’alternativa a aquest cúmul de dubtes pot ser la dels models mixts que, a més a més, poden aportar altres avantatges:

Si considerem la platja com a factor fix, les nostres conclusions es circumscriuen a aquelles platges concretes.
Però, si considerem la platja com un factor aleatori, les nostres conclusions es referiran a totes les platges de les quals, aquestes nou, en són una mostra. Per tant, d’aquesta manera podrem predir la relació entre riquesa i NAP en un sentit més general, sense haver de limitar-nos a les nou platges de la nostra mostra.

Aproximació amb model lineal mixt

Model amb `Beach` com a factor aleatori

Descripció del model

Si reanalitzem el primer dels models anteriors, el que consideraba la existència de diferents ordenades a l’origen de la recta de regressió entre riquesa i NAP. Si considerem la platja com un factor aleatori, el model quedaria així \[ \small Riquesa_{ij}=\alpha+\beta_{}NAP_{ij}+a_j+e_{ij}~~~~~~e_{ij}\sim N(0,\sigma^2) \] L’índex \(_j\), representant les platges, pren valor d’1 a 9, mentre \(_i\), representant les mostres dins de cada platja, pren valors d’1 a a 5. En aquest model assumim que només hi ha una única recta de regressió amb una sola ordenada a l’origen i una sola pendent. L’ordenada a l’origen \(\small\alpha\) i la pendent \(\small \beta\) són els paràmetres o efectes fixos del model que a més a més inclou un efecte aleatori \(\small a_j\) que afegeix a l’ordenada a l’origen la variació particular de cada una de les platges. Els paràmetres són quatre en aquest cas: \(\small\alpha\), \(\small \beta\), la variància residual \(\small\sigma^2\) i la variància de l’ordenada a l’origen \(\small\sigma_a^2\). Per tant, un model equivalent al model lineal, però amb 4 paràmetres en comptes de 11, perquè, en comptes de tenir 9 ordenades a l’origen, estimades per a cada una de les rectes de regressió, tenim 9 valors estimats \(\small a_1,...,a_9\) que assumim que procedeixen d’una distribució normal i el que s’estima en el model és la variància en la distribució.

Ajustament del model amb `lmer()`

lmer1 <- lme4::lmer(riquesa ~ NAP + (1|Beach), data = rikz)
summary(lmer1)

## Linear mixed model fit by REML ['lmerMod']
## Formula: riquesa ~ NAP + (1 | Beach)
##    Data: rikz
## 
## REML criterion at convergence: 239.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4227 -0.4848 -0.1576  0.2519  3.9794 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Beach    (Intercept) 8.668    2.944   
##  Residual             9.362    3.060   
## Number of obs: 45, groups:  Beach, 9
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.5819     1.0958   6.007
## NAP          -2.5684     0.4947  -5.192
## 
## Correlation of Fixed Effects:
##     (Intr)
## NAP -0.157

Imprimint el summary,

obtenim en primer lloc la variància dels efectes aleatoris (Random effects): la de l’ordenada en origen \(\small\sigma_a^2=8.668\) (recordem que no s’estimen valors individuals dels efectes aleatoris, sinó la variància) i la variància residual \(\small\sigma^2=9.362\).
en segon lloc, pel que fa els efectes fixos, es faciliten els valors de \(\small\hat\alpha=6.582\) i de \(\small\hat\beta=-2.568\).
A baix hi ha una taula que mostra la correlació entre els coeficients. En aquest cas tenim dos coeficients, per tant, una correlació.

Què ens diu el model en resum? En resum el model estima una respost general de la forma \(\small 6.582-2.568NAP_{ij}\). Per a cada platja, la constant augmenta i disminueix per un valor aleatori que segueix una distribució \(\small N(\mu=0,\sigma_a^2=8.668)\). La variància residual, és a dir, l’error que podem afegir a cada una de les nostres prediccions és \(\small \sigma^2= 9.362\).

Realitzarem un test ANOVA amb Anova() de car, ja que anova() no genera valors de significació per a objectes lmer.

car::Anova(lmer1)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: riquesa
##      Chisq Df Pr(>Chisq)    
## NAP 26.952  1  2.085e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Com podem observar, l valor de significació del valor p per a una distribució \(\small \chi^2\) és infinitesimal (2.085e-07), que indica que existeix un cert efecte del nivell intermareal en la riquesa del bentos.

Model amb `Beach` i `NAP` com a factors aleatoris

Descripció del model

Continuant amb el procediment que hem seguit per la regressió lineal, podem contemplar la possibilitat que hi hagi, a més de diverses ordenades a l’origen, també diverses pendents per a cada una de les platges.

En aquest cas, tindríem el model següent \[ \small Riquesa_{ij}=\alpha+a_j+\beta_{}NAP_{ij}+b_jNAP_{ij}+e_{ij}~~~~~~e_{ij}\sim N(0,\sigma^2) \] \[ \small \begin {pmatrix} a_j\\ b_j \end{pmatrix}\sim N \bigg( \textbf{0}= \begin {bmatrix} 0\\ 0\\ \end {bmatrix}, ~\textbf Λσ^2=\begin {bmatrix} \sigma_a^2 & \rho\sigma_a\sigma_b\\ \rho\sigma_a\sigma_b & \sigma_b^2\\ \end {bmatrix} \bigg) \]

El model és igual que l’anterior, excepte pel que fa al terme \(b_jNAP_{ij}\) que permet la variació aleatòria de la pendent a cada una de les platges. Com hem vist al principi, els efectes aleatoris procedeixen d’una distribució normal multivariada amb mitjana \(\small \bf 0\), i matriu de variàncies-covariàncies \(\small \bf \Lambda\sigma^2\). El model és molt semblant al que hem analitzat abans al Cas 3. Tanmateix amb molt menys paràmetres 5 (aquell en tenia 19!), dos per als efectes fixos (\(\small \alpha,\beta\)) i 3 per a les variàncies dels termes aleatoris (\(\small \sigma^2, \sigma_a^2,\sigma_b^2\)).

Ajustament del model amb `lmer()`

lmer2 <- lme4::lmer(riquesa ~ NAP + (NAP|Beach), data = rikz)
summary(lmer2)

## Linear mixed model fit by REML ['lmerMod']
## Formula: riquesa ~ NAP + (NAP | Beach)
##    Data: rikz
## 
## REML criterion at convergence: 232.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8212 -0.3410 -0.1674  0.1925  3.0397 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Beach    (Intercept) 12.600   3.550         
##           NAP          2.942   1.715    -0.99
##  Residual              7.307   2.703         
## Number of obs: 45, groups:  Beach, 9
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.5884     1.2649   5.208
## NAP          -2.8301     0.7231  -3.914
## 
## Correlation of Fixed Effects:
##     (Intr)
## NAP -0.820
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00267151 (tol = 0.002, component 1)
##   See ?lme4::convergence and ?lme4::troubleshooting.

La part del resultat dels efectes aleatoris mostra que la variància atribuïble a l’ordenada a l’origen és \(\small \sigma_a^2=12.6\) i la variància atribuïda a la pendent és \(\small \sigma_b^2=2.942\) i finalment la variància residual és \(\small \sigma^2=7.307\). Hi ha per tant una disminució de la variància residual respecte al model anterior. Pel que fa als efectes fixos els valors estimats són \(\small \hat\alpha=6.5884\) i \(\small \hat\beta=-2.8301\).

La correlació entre efectes aleatoris és \(\small \rho=-0.99\), negativa, com la correlació en els efectes fixos. Aquestes correlacions tan fortes (properes als valos -1 i 1) poden generar algun problema d’ajustament: me això el missatge “Model failed to converge...”, que ja veurem com resoldre.

El valors dels efectes aleatoris estimats per a cada platja es poden obtenir amb coef()

coef(lmer2)

## $Beach
##   (Intercept)       NAP
## 1    8.419607 -3.658230
## 2   12.362712 -5.540119
## 3    3.806803 -1.505151
## 4    3.562397 -1.385380
## 5   11.200323 -5.134689
## 6    4.425965 -1.776055
## 7    4.082558 -1.643322
## 8    5.099701 -2.107085
## 9    6.335511 -2.721025
## 
## attr(,"class")
## [1] "coef.mer"

No es tracta de paràmetres del model, sinó modes condicionades (denominades també mitjanes posteriors).

Els podem representar com a diferències amb els seus paràmetres fixos corresponents. Per això podem emprar la funció ranef() de lme4 que ens donarà el valor d’aquesta diferència, condval, incloent-hi l’error estàndard, condsd. I així, calculant els intervals de confiança condicionals, representar-los amb dos panells, amb ggplot2:

as.data.frame(lme4::ranef(lmer2, condVar = TRUE)) %>% 
  mutate(lower= condval-1.96*condsd,
         upper= condval+1.96*condsd) %>% 
  ggplot(aes(x= condval, y=grp, xmin=lower,xmax=upper)) +
  facet_wrap(~ term) +
  geom_point() +
  geom_segment(aes(x = lower, y = grp, 
                   xend = upper, yend = grp)) +
  labs(y= "Beach",
       x= "Dif. entre mitjanes posteriors i termes fixos",
       title = stringr::str_wrap("Diferències en els efectes aleatòris (modes condicinades), respecte als valors mitjans dels efectes fixos per a cada nivell del factor Beach, amb IC del 95%",60)) +
  theme_bw() +
  theme(plot.title = element_text(size = 12))

Model amb `Beach` i `NAP` com a factors aleatoris, sense correlació recíproca

Per ajustar un model com l’anterior però, que ignorés la correlació entre efectes aleatoris, l’haurem d’escriure utilitzant la doble barra ||, de la següent manera:

lmer2.2 <- lme4::lmer(riquesa ~ NAP + (NAP||Beach), data = rikz)
summary(lmer2.2)

## Linear mixed model fit by REML ['lmerMod']
## Formula: riquesa ~ NAP + ((1 | Beach) + (0 + NAP | Beach))
##    Data: rikz
## 
## REML criterion at convergence: 238
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5983 -0.3995 -0.1062  0.1372  3.3093 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Beach    (Intercept) 12.295   3.506   
##  Beach.1  NAP          2.711   1.647   
##  Residual              6.793   2.606   
## Number of obs: 45, groups:  Beach, 9
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.7870     1.2488   5.435
## NAP          -2.6791     0.7034  -3.809
## 
## Correlation of Fixed Effects:
##     (Intr)
## NAP -0.075

Veiem que, com era d’esperar, desapareix el paràmetre de correlació i desapareix també el missatge d’advertiment.

Inferència amb models mixts

Després d’estimar els paràmetres del model, el pas següent és la inferència estadística, és a dir, derivar conclusions estadístiques i biològiques.

La forma clàssica fe fer inferència es amb tests específics per comprovar la significació de cada variable i la conveniència de incloure-la al model. En el cas dels models mixts la cosa es complica perquè hi tenim dos tipus d’efectes, els fixos i els aleatoris, que han de ser avaluats separadament.

L’avaluació dels efectes aleatoris presenta el problema addicional i és que l’estimació per màxima versemblança ML, subestima la variància aleatòria i residual, cosa que pot afectar-ne la significació en la inferència. L’opció més habitual per a solucionar aquest problema és utilitzar la màxima versemblança restringida REML. Això fa que en l’avaluació dels models mixts es procedeixi en dues fases: utilitzant l’algoritme REML per estimar els efectes aleatoris i l’algoritme ML per als efectes fixos. Per optimitzar tant l’estructura d’efectes aleatoris com la dels efectes fixos en un models es recomana seguir aquest procediment:

Començar amb un model on la component fixa contingui totes les variables fixes i interaccions possibles (model més allà de l’òptim).
Sobre aquest model buscar de la component aleatòria òptima, provant diferent models alternatius, amb diverses components aleatòries, sense canviar la component fixa. En en aquesta fase utilitzar REML.
Un cop obtinguda l’estructura òptima de la component aleatòria, buscarem l’estructura òptima de la component fixa. En aquesta fase utilitzar ML.
Quan s’ha seleccionat l’estructura de la component fixa, es presenta el model final utilitzant un estimador REML.

Aquí detallarem diferents eines a utilitzar per seleccionar els efectes fixos i aleatoris en un model mixt.

Eines d’inferència: Test de la F

Aquests test s’empren únicament per a la selecció d’efectes fixos. En el cas del models mixts, calcular aquest test té la dificultat d’estimar adequadament els graus de llibertat corresponents als paràmetres fixos en presència d’efectes aleatoris. El que resulta particularment delicat en el cas de disseny no equilibrat.

En el cas de models ajustats amb lmer, la funció anova() proporciona el valor de l’F-estadístic però, intencionalment, omet els valors p. Per exemple, tornant al nostre primer model

anova(lmer1)

## Analysis of Variance Table
##     npar Sum Sq Mean Sq F value
## NAP    1 252.33  252.33  26.952

Aquest comportament és subsanable més directament amb la funció Anova() que estima els graus de llibertat amb l’aproximació Kenward-Roger, realitzant el test de la F (amb el paràmetre test.statistic="F"):

car::Anova(lmer1, test.statistic="F")

## Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
## 
## Response: riquesa
##          F Df Df.res    Pr(>F)    
## NAP 26.481  1 37.028 8.961e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Eines d’inferència: Test d’hipòtesi de models niats

Comparació de models ambdós mixts

El test d’hipòtesi per a models niats, també conegut com a test de quocient de versemblança (LRT). És una mesura relativa de bondat d’ajustament perduda o guanyada inserint un nou terme al model. No es recomana l’LRT per fer inferència sobre els efectes fixos del model (és molt anticonservatiu, propens a l’error de tipus I).

anova(lmer1,lmer2)

## Data: rikz
## Models:
## lmer1: riquesa ~ NAP + (1 | Beach)
## lmer2: riquesa ~ NAP + (NAP | Beach)
##       npar    AIC    BIC  logLik -2*log(L) Chisq Df Pr(>Chisq)  
## lmer1    4 249.83 257.06 -120.92    241.83                      
## lmer2    6 246.66 257.50 -117.33    234.66 7.173  2    0.02769 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Com funciona el test de quocient de versemblança (LRT)

Per realitzar un LRT en R, es comparen dos models en funció de la seva versemblança (likelihood), és a dir, la probabilitat que el model hagi generat les dades observades:

Model nul (H₀): el model més simple i niat.
Model alternatiu (H₁): El model més complex que conté tots els termes del model nul, a més de termes addicionals.

La prova calcula el quocient d’aquestes probabilitats. Com que les probabilitats logarítmiques són les que es s’utilitzen habitualment en models mixts, la raó es converteix en una simple resta (la diferència de desviació):

L.ratio = -2 × (Log-Likelihood of Null Model - Log-Likelihood of Alternative

Per tant, cal fixar l’atenció en el valor del quocient que serà Chisq=241.83-234.66=7.17 que és el que tenim a la taula. Finalment, s’obtindrà el valor p comparant amb la distribució de \(\chi^2\).

pchisq(7.173, , df=2, lower.tail = F)

## [1] 0.02769509

Comparació d’un model amb efectes aleatoris i un que no en té

També podem voler comparar un model mixt i un model que només tingui efectes fixos. En aquest cas, el model de només efectes fixos hauria d’ajustar-se amb la funció gls() de la llibreria nlme, i el model amb efectes fixos amb la funció lme() de nlme.

lme2 <- nlme::lme(riquesa ~NAP, data = rikz, random = ~ 1|Beach)
gls.1 <- nlme::gls(riquesa ~ NAP, data = rikz)
anova(lme2, gls.1)

##       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## lme2      1  4 247.4802 254.5250 -119.7401                        
## gls.1     2  3 258.2010 263.4846 -126.1005 1 vs 2 12.72075   4e-04

Eines d’inferència: Comparació de models amb AIC

El criteris d’informació es poden fer servir per comparar dos o més models de forma simultània. Entre els indicadors basats en criteris d’informació hi ha el Criteri d’informació d’Akaiké (AIC) que sabem que es deineix amb \[ \small AIC=2k-2L(\theta) \] on \(\small L(\theta)\) és la log-versemblança del model i \(\small k\) el nombre de paràmetres ajustats.

En general es defineix que com més petit l’AIC, millor és l’ajustament, sempre en terme relatius, és clar.

Aquest mètode es pot emprar per comparar bateries de models alternatius amb diferents efectes i combinacions d’efectes. En principi, es pot comparar qualsevol conjunt de models, complint les següents condicions:

La variable de resposta és la mateixa en tots els models.
Els models s’han construït i ajustat utilitzant el mateix nombre d’observacions.
L’estructura de distribució dels errors és la mateixa en tots els models.
Els paràmetres dels models han d’haver-se estimat amb el mateix mètode (ML vs. REML).

En general una diferència \(\small \Delta>2\) es considera suficient per a considerar un model millor o pitjor d’un altre.

Comparació automàtica de models

Per comparar un gran nombre de models podem fes servir una funció que ho fa de manera automàtica (per a que coneix la temàtica de sèries temporals, seria equivalent a auto.arima()): la funció dredge() de MuMln. Aquesta funció no permet generar models amb diferents estructures d’efectes aleatoris. Genera un resultat en forma de taula on tots els models estan ordenats de menor (el millor) a major AIC. Els valors a sota de cada variable són coeficients estimats. La funció dredg() utilitza pr defecte l’AIC corregit (AICc), per a mostres de petites dimensions (com en el nostre exemple).

Per exemple, podem pensar ajustar ajustar una sèrie de models, començant pel més complex entre ells, en el qual la riquesa del bentos pot estar en funció de NAP, salinitat i temperatura:

lmer5 <- lme4::lmer(riquesa ~ NAP + temperature + salinity + (1|Beach), 
                    data = rikz, 
                    REML = F,
                    na.action = "na.fail")
dd <-MuMIn::dredge(lmer5)
dd

## Global model call: lme4::lmer(formula = riquesa ~ NAP + temperature + salinity + 
##     (1 | Beach), data = rikz, REML = F, na.action = "na.fail")
## ---
## Model selection table 
##   (Intrc)    NAP slnty  tmprt df   logLik  AICc delta weight
## 8 -60.910 -2.818 1.820 0.8753  6 -114.320 242.8  0.00  0.774
## 4 -48.720 -2.523 1.968         5 -117.032 245.6  2.75  0.196
## 6 -12.080 -2.700       0.9963  5 -119.567 250.7  7.82  0.016
## 2   6.584 -2.576               4 -120.915 250.8  7.98  0.014
## 3 -57.420        2.246         4 -127.694 264.4 21.54  0.000
## 7 -64.970        2.181 0.4998  5 -127.111 265.8 22.91  0.000
## 1   5.689                      3 -131.652 269.9 27.04  0.000
## 5  -6.746              0.6623  4 -131.231 271.5 28.61  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | Beach

Es poden seleccionar els millors models que podem especificar. Am la funció base subset. En el nostre cas si volem només \(\small \Delta<2\), és a dir només els models que tinguin una diferència d’AIC amb el millor model menor que 2:

subset(dd, delta <2)

## Global model call: lme4::lmer(formula = riquesa ~ NAP + temperature + salinity + 
##     (1 | Beach), data = rikz, REML = F, na.action = "na.fail")
## ---
## Model selection table 
##   (Intrc)    NAP slnty  tmprt df  logLik  AICc delta weight
## 8  -60.91 -2.818  1.82 0.8753  6 -114.32 242.8     0      1
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | Beach

En aquest cas, dredge() ens ha fet tota la feina, deixant només un model, clarament millor que els altres. Però no necessàriament ha de sempre així. Aquí doncs poden entrar consideracions de diversa índole per considerar un model respecte a un altre. En qualsevol cas, sigui quina suigui la decisió es recomana justificar de manera clara, aportant el detall que faci falta.

Avaluació del model

Un cop identificats els millors models el pas final necessita l’avaluació del model. Al contrari d’altres casos, en aquest no disposem d’una funció que generi gràfics de valors residuals de manera automàtica. És a dir que els haurem de crear manualment. No saltres ho farem amb ggplot2:

lmer2.aval <- data.frame(res = residuals(lmer2, type = "response"),
           fit = fitted(lmer2))

p1 <- lmer2.aval %>% 
  ggplot(aes(x=fit,y=res)) + 
  geom_point() +
  geom_hline(yintercept = 0) +
  labs(x = "Fitted values",
       y="Residuals") +
  theme_bw()

p2 <- rikz %>% 
  ggplot(aes(x=NAP, y=lmer2.aval$res)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  labs(x = "NAP",
       y="Residuals") +
  theme_bw()

p3 <- lmer2.aval %>% 
  ggplot(aes(x=res)) + 
  geom_histogram(bins = 9, fill= "grey", colour="black") +
  labs(x = "Residuals",
       y="Frequency") +
  theme_bw()

p4 <- qplot(data=lmer2.aval,sample=res)+
  stat_qq_line()+
  labs(x = "Theoretical quantiles",
       y="Sample quantiles") +
  theme_bw()

library(patchwork)

(p1+p2) / (p3+p4)

Estimació dels components de la variància en models mixts

Un cop ajustat el model, feta la inferència i controlats els supòsits, pot ser interessant valorar la bondat de l’ajustament del model a les dades.

En el cas dels models mixts, trobar una mesura de la bondat d’ajustament similar al coeficient de determinació (\(\small R^2\)) del models lineals o a la deviància explicada \(\small D^2\) del models GLM és més complicat. Això es deu a què hi ha múltiples prediccions dels efectes fixos, atès que els coeficients se’n veuen afectats per un error aleatori determinat per les variàncies estimades per a cada un dels factors aleatoris.

Sortosament, es pot estimar un anàleg de l’\(\small R^2\), gràcies a les equacions proposades Nakagawa i Schielzeth (2013) i Nakagawa (2017), que es poden impementar ara amb la funció r.sqaredGLMM() de la llibreria MuMIn. L’esmentada funció calcula dos \(\small R^2\): la condicional, \(\small R_c^2\) i la marginal, \(\small R_m^2\).

\(\small R^2\) marginal (\(\small R_m^2\)): es refereix només a la part de variància explicada només pels efectes fixos del model.
\(\small R^2\) condicional (\(\small R_c^2\)): Representa la proporció de variància explicada pel model sencer, efectes fixos inclosos.

Per exemple, utilitzant un model ja ajustat abans,

MuMIn::r.squaredGLMM(lmer2)

##            R2m       R2c
## [1,] 0.2951295 0.7276647

es pot veure que els efectes fixos (en aquest cas el NAP) expliquen un 30% de la variància, però que més d’un 43% (0.728-0.295) depenent del factor aleatori.

Models lineals generalitzats mixts (GLMM)

Com ben sabem, amb dades de recomptes (com “abundancia”, “riquesa”, etc.), la distribució de Poisson o la binomial negativa solen ser més adequades, utilitzant per tant mètodes lineals generalitzats (GLM). Quan a aquests aquests mètodes afegim factors aleatoris, tenim mètodes lineals generalitzats mixts (GLMM).

Ajustament del model

En R, ajustar GLMM és relativament fàcil utilitzant la funció glmer() del paquet lme4. La funció necessita que s’especifiqui la família de distribució d’errors (com amb glm()), però la resta d’arguments resten iguals als que ja hem vist per als ajustaments amb lmer(). Ja veurem que a més a més que, a diferència de lmer(), glmer() executa models només de màxima versemblança. Per tant, caldrà tenir present que l’estructura òptima, tant del efectes fixos com dels aleatoris, s’estarà analitzant amb el log-lik.

Considerem un model GLM equivalent al model lmer2 que ja hem executat

glmer2 <- lme4::glmer(riquesa ~ NAP + (NAP | Beach), family = poisson, 
            data = rikz)
summary(glmer2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: riquesa ~ NAP + (NAP | Beach)
##    Data: rikz
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     218.7     227.8    -104.4     208.7        40 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.35846 -0.51129 -0.21846  0.09802  2.45384 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr
##  Beach  (Intercept) 0.2630   0.5128       
##         NAP         0.0891   0.2985   0.18
## Number of obs: 45, groups:  Beach, 9
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.6942     0.1868   9.071  < 2e-16 ***
## NAP          -0.6074     0.1374  -4.421 9.81e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## NAP 0.121

Naturalment els coeficients dels efectes fixos dels models glmer2 i glmer2 no coincideixen, en bona part perquè l’elecció de la distribució d’errors de tipus Poisson selecciona per defecte una funció de vincle logarítmica que linearitza la relació entre variable de resposta i les explicatives.

Avaluació del model

Comprovarem ara els gràfics de valors residuals de manera anàloga a la que ja hem vist a dalt

glmer2.aval <- data.frame(
  res = residuals(glmer2, type = "response"),
  fit = fitted(glmer2)
  )

p1g <- glmer2.aval %>% 
  ggplot(aes(x=fit,y=res)) + 
  geom_point() +
  geom_hline(yintercept = 0) +
  labs(x = "Fitted values",
       y="Residuals") +
  theme_bw()

p2g <- rikz %>% 
  ggplot(aes(x=NAP, y=glmer2.aval$res)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  labs(x = "NAP",
       y="Residuals") +
  theme_bw()

p3g <- glmer2.aval %>% 
  ggplot(aes(x=res)) + 
  geom_histogram(bins = 9, fill= "grey", colour="black") +
  labs(x = "Residuals",
       y="Frequency") +
  theme_bw()

p4g <- qplot(data=glmer2.aval,sample=res)+
  stat_qq_line()+
  labs(x = "Theoretical quantiles",
       y="Sample quantiles") +
  theme_bw()

library(patchwork)

(p1g+p2g) / (p3g+p4g)

En general veiem hi ha un lleuger millorament dels gràfics de valors residuals, una lleugera millora de la normalitat, se deixa de veure el patró de con que s’albirava en el model lineal anterior.

Altres distribucions d’errors disponibles amb `glmer()`

La funció `glmer() permet ajustar GLMM amb distribucions diverses (Poisson, binomial, gamma, quasi-Poisson, etc.). Fins fa poc no disposava de cap funció per ajustar models amb distribució d’errors binomial negativa. Recentment s’ha resolt inserint la funció glmer.nb(), específica per aq aquesta distribució que funciona de la següent manera:

glmer2.nb <- lme4::glmer.nb(riquesa ~ NAP + (NAP | Beach), data = rikz)
summary(glmer2.nb)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(680.3001)  ( log )
## Formula: riquesa ~ NAP + (NAP | Beach)
##    Data: rikz
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     220.7     231.6    -104.4     208.7        39 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.36077 -0.50851 -0.21766  0.09643  2.44649 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr
##  Beach  (Intercept) 0.26221  0.5121       
##         NAP         0.08803  0.2967   0.18
## Number of obs: 45, groups:  Beach, 9
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.6937     0.1871   9.053  < 2e-16 ***
## NAP          -0.6076     0.1384  -4.391 1.13e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## NAP 0.123

En tot cas, els autors del paque lmer4 destaquen que algunes parts de glmer.nb() encara són experimentals i els mètodes encara falten o són subòptims. En particular, encara no hi ha cap inferència disponible per al paràmetre de dispersió θ.

Per ajustar un model binomial negatiu amb un paràmetre de sobredispersió conegut (per exemple, com a part d’un exercici de comparació de models), recomanen utilitzar glmer() amb la família negative.binomial del paquet MASS, per exemple, glmer(...,family=MASS::negative.binomial(theta=1.75)).

Recordem que és +1 perquè incloem també el vector corresponent a l’ordenada a l’origen.↩︎
Bentos marí: Organismes i comunitats que viuen als fons marins.↩︎

Models mixts (LMM, GLMM, GAMM) en ecologia

Notes d’estadística aplicada a les ciències ambientals: Nota Núm. 9

Stefano Sanfilippo

Juny 2026

Introducció

Definició de factor fix i factor aleatori

Formulació teòrica del model d’efectes mixts

Referències

Models mixts amb R: nlme i nlme4

La funció lmer()

Model amb un factor aleatori

Model amb un factor aleatori amb efecte sobre un factor fix o covariable

Model amb més d’un factor aleatori

Model amb factors anidats

Entendre els models lineals mixts (LMM) amb un exemple

El dataset RIKZ

Preparació de les dades

Aproximació amb model lineal

Cas 0: independència entre mostres

Cas 1: pendent constant

Cas 2: intercept constant

Cas 3: intercept i pendent variables

Conclusions

Aproximació amb model lineal mixt

Model amb Beach com a factor aleatori

Descripció del model

Ajustament del model amb lmer()

Model amb Beach i NAP com a factors aleatoris

Descripció del model

Ajustament del model amb lmer()

Model amb Beach i NAP com a factors aleatoris, sense correlació recíproca

Inferència amb models mixts

Eines d’inferència: Test de la F

Eines d’inferència: Test d’hipòtesi de models niats

Comparació de models ambdós mixts

Comparació d’un model amb efectes aleatoris i un que no en té

Eines d’inferència: Comparació de models amb AIC

Comparació automàtica de models

Avaluació del model

Estimació dels components de la variància en models mixts

Models lineals generalitzats mixts (GLMM)

Ajustament del model

Avaluació del model

Altres distribucions d’errors disponibles amb glmer()

Models mixts amb R: `nlme` i `nlme4`

La funció `lmer()`

El dataset `RIKZ`

Model amb `Beach` com a factor aleatori

Ajustament del model amb `lmer()`

Model amb `Beach` i `NAP` com a factors aleatoris

Ajustament del model amb `lmer()`

Model amb `Beach` i `NAP` com a factors aleatoris, sense correlació recíproca

Altres distribucions d’errors disponibles amb `glmer()`