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

0. Referències

  • L. Cayuela & M. De la Cruz, Análisis de datos ecológicos en R, 2022.

Introducció

Els models lineals clàssics (regressió, ANOVA, ANCOVA) es basen en tres supòsits:

  • La variable de resposta té distribució normal.
  • La variança de la variable de resposta és constant.
  • En el cas de variables de variables de resposta continues, aquestes es relacionen linealment amb les covariables.

En el cas de dades ecològiques, frequentement un o més d’aquests supòsits no es compleix. Tot sovint, a mesura que augmenta el valor de \(y\), augmenta la variança. O es poden predir valors negatius o poc creïbles.

Un dels remeis que es poden dur a terme és la transformació de les dades. Però sovint no resolen les problemàtiques de heteroscedasticitat, de falta de normalitat o de linealitat. En aquests casos, els models lineals generalitzats (GLM) poden aportar alternatives d’anàlisi més convenients.

Específicament, és molt recomanable utilitzar GLMs quan la variable de resposta consisteixi en:

  • Recomptes de casos.
  • Recomptes de casos rexpressats en proporcions.
  • Una resposta binària.

En el cas de variables de resposta contínues és raonable començar amb un model linial pero podem arribar a fer servir un GLM si s’incompleix algun dels supòsits i no es corregeix per transformació.

Estructura dels GLM

Els tres components importants del GLM són:

  • El predictor linial
  • La funció de vincle
  • El component aleatori

El predictor linial

És una combinació linial d’una o més variables explicatives. Es representa així

\[ \small \eta=\beta_0+\beta_1X_1+...+\beta_kX_k \] Representa la part determinista del model.

La funció del vincle

La funció de vincle comporta la transformació de la variable \(y\). Fonamentalment aquesta transformació busca linealitzar la relació entre la \(y\) y la \(x\). A més a més, permet transformar en continues variables de resposta que no ho són (com les que són recomptes o binàries 0-1).

L’ús de la funció de vincle canvia d’alguna manera la naturalesa de la variable de resposta. Per això, sol ser convenient aplicar la funció inversa per tal de generar prediccions amb la nostra variable de resposta. Com, per exemple, en un model lineal \[ \small ln(y)=\beta_0+\beta_1x \] Para obtener los valores estimados de \(y\) hay que aplicar la función inversa a la función logarítmica, en este caso la función exponencial \[ \small e^{ln(y)}=e^{\beta_0+\beta_1x}~~~=>~~~y=e^{\beta_0+\beta_1x} \] Ambdós models són equivalents, però serà més fàcil explicar la nostra història a partir del segon.

Un altre avantatge és que permet mantenir les nostres prediccions adherents a un determinat rang de valors. Per exemple, la funció logarítmica evitarà tenir valors negatius per a dades de recompte i la funció logit és ideal per a dades 0-1.

Les funcions de vincle més comuns per a dades ecològiques amb GLM; \(\mu\) és el valor esperat de la variable de resposta; \(\eta\) és el predictor linial.

Funció de vincle Ús Funció de vincle inversa
Identitat: \(\small\mu\) Dades continues amb distribucions normals (regressió, ANOVA, ANCOVA) Identitat: \(\small \eta\)
Logaritmica: \(\small\log(\mu)\) Recomptes amb distribucions Poisson o binomial negativa Exponencial: \(\small e^\eta\)
Arrel quadrada: \(\small\sqrt{\mu}\) Recomptes amb distribucions Poisson o binomial negativa Quadrat: \(\small\eta^2\)
Logit: \(\small\log\frac{\mu}{1-\mu}\) Proporcions (dades entre 0 i 1) o variables binaries amb distr. binomial Logística: \(\small\frac{e^\eta}{1+e^\eta}\)
Recíproca: \(\small\mu^{-1}\) Dades contínues amb distribucions asimètriques negatives (ex. gamma) Recíproca: \(\small\eta^{-1}\)


El component aleatori

Una de les novetats dels models GLM és la possibilitat d’especificar és la possibilitat d’especificar erorrs no normals per al component aleatori, cosa molt útil amb dades econògiques que, si no, implicaria acceptar assumpcions poc raonables.

Les distribucions més comuns són

  • Poisson: La variança aumenta amb el valor observat. Molt útil per a recomptes (ex. nombre animals atropellats, nombre de dies de gelades duran l’any, nombre de colònies bacterianes en una placa d’agar, etc.).
  • Binomial negativa: Se puede entender como una distr. de Poisson en la que la relació entre la variança i la mitjana no sigui lineal. S’utilitza amb dades de recomptes que presentin sobredispersió.
  • Binomial: molt útil para proporcions (taxes i ratios) i dades de presència/absència.
  • Gamma: Per a dades contínues positives amb coef. de variació constant.

Tot això es pot resumir en una taula:

Família de distribució Funció de vincle usual Rang de \(y\) \(\small Var(y/x)\)
Gaussiana Identitat \(-\infty,\infty\) \(\phi\)
Poisson Logarítmica \(0,1,2,...\) \(\mu\)
Binomial negativa Logarítmica \(0,1,2,...\) \(\small\mu+\frac{\mu^2}{\phi}\)
Binomial Logit \(\small\frac{0,1,2,...,N}{N}\) \(\small\frac{\mu(1-\mu)}{N}\)

Relació entre model linial i GLM

El model linial es pot representar així:

\[ \small y=\underbrace{\beta_0+\beta_1x_1+...+\beta_kx_k}_\text{Comp. determinista} + \underbrace{\epsilon~~~ \epsilon\sim N(0,\sigma^2)}_\text{Comp. aleatori} \]

El GLM pot agafar aquesta forma general: \[ \small f(y)=\underbrace{\beta_0+\beta_1x_1+...+\beta_kx_k}_\text{Comp. determinista} + \underbrace{\epsilon~~~ Var(y|x)\sim f(\mu)}_\text{Comp. aleatori} \]

Les diferències més evidents són doncs en:

  • \(f(y)\): La funció de vincle que transforma la variable de resposta.
  • \(\epsilon\): El component aleatori que ara pot agafar valors de la Varianza \(Var(y|x)\), no constants y que varien en funció del valor esperat \(\mu\).


La funció glm() en R

Per ajustar un model GLM en R es pot utilitzar la funció glm(). Els arguments principals són formula, que funciona igual que en lm(), i family que és una funció que indica la distribució d’errors i la funció de vincle, que pot agafar els valors següents:

  • binomial(link="logit");
  • gaussian(link="identity");
  • Gamma(link="inverse");
  • inverse.gaussian(link="1/mu^2");
  • poisson(link="log");
  • quasi(link="identity",variance="constant");
  • quasibinomial(link="logit");
  • quasipoisson(link="log").

Si no s’especifica la funció de vincle, el model agafa per defecte la més usual (veure la taula a dalt).

Algunes distribucions d’errors no estan implantades en la funció glm(). Una d’aquestes és la binomial negativa que es pot especificar amb la funció glm.nb() de la liberia MASS.

Estimació de paràmetres per màxima versemblança

L’ajustament del model implica l’estimació de paràmetres. En els models linials aquesta es fa el métode de mínims quadrats1, mentre que en GLM l’estimació de paràmetres es duu a terme mitjanant mètode de màxima versemblança ML (Maximum Likelihood).

Aquest consisteix bàsicament en què, per a cada observació \(y_i\), se’n calcula la probabilitat de pertànyer a una funció de distribució normal amb mitjana \(\mu=\beta_0+\beta_1x\) y varianza \(\sigma^2\). El producte de les probabilitats per a tots els punts ens dona la versemblança \(L\) del model, tot i que sovint es fa servir la \(\log\)-versemblança (\(\log L\)), el sumatori dels logaritmes de les probabilitats, més senzill de calcular.

Funcionament de l’algoritme de càlcul

Aquest estimador, a diferèncie del de mínims quadrats, treballa de manera iterativa. Amb l’exemple seguent mirarem com funciona.

Simulació de les dades

Simularem unes dades d’un estudi que intentar predir cóm cambia la riqueza d’espècies herbàcies des del marge cap a l’interior d’un bosc. Les dades son de tipus recompte amb la qual cosa s’assumeix una distribució de Poisson.

Generarem doncs (\(y\)) unes dades hipotitzant que aquestes segueixin un model exponencial respecte a la distància (\(x\)) del tipus \(\small y=e^{\beta_0+\beta_1x}\), donde \(\small\beta_0+\beta_1x\) és el predictor linial recordem que en el context de GLM aquest model és equivalent a \(\small\log (y)=\beta_0+\beta_1x\).

Establint els valors dels coeficients en 4 i -0.02, l’equació quedaria \[ \small y=e^{4-0.02x}~~~~<=>~~~~ \small\log (y)=4-0.02x \]

Establir valors de la \(x\)

Establim un número d’observacions n i els valors de la variables explicativa (dist), simulant una obervació cada 20m.

dist <- seq(0,200, by= 20)

n <- length(dist)

Predictor linial

Es crea una matriu del model on, com sabem, la primera columna són 1, i la segona són els valors de la \(x\) (dist).

Xmat <- model.matrix(~dist)

Establim el vector dels valors dels coeficients 4 i 0.02

vector.coeficients <- c(4,-0.02)

Obtenim el predictor lineal (\(\beta_0+\beta_1x\)), fila per fila

pred.lineal <- Xmat %*% vector.coeficients

Variable de resposta

Apliquem la inversa de la funció del vincle per obtenir els valors de la variable de resposta

mu <- exp(pred.lineal)

Valors esperats de la y

set.seed(1)
y <- rpois(n, mu)

plot(y ~ dist,
     ylab = "Riquesa d'especies", 
     xlab= "Distància del marge (m)")

Estimació dels paràmetres

Ara que tenim les dades, per estimar els coeficients dels models a partir d’aquestes, provarem diferents combinacions dels valors dels paràmetres, y compararem els models resultants des del punt de vista de la seva versemblança.

Mètode per força bruta

En aquest exemple, para entendre la lògica del procés, mirarem com treballa el mètode per força bruta, encara que no s’utilitzarà mai, pel seu alt cost computacional i poca eficiència.

En tot cas, sempre és necessari definir un rang de valors per a cada paràmetre, de manera de delimitar l’espai multidimensional on es realitzarà la cerca. En el nostre cas, els paràmetres a determinar serien 3: l’ordenada a l’origen \(\beta_0\), la pendent de la funció de regressió \(\beta_1\), la variança residual \(Var(y|x)\). Com que en la distribució de Poisson la variança coincideix amb el valor esperat (\(Var(y|x)=\mu\)), només caldrà estimar els primers dos paràmetres.

Definició dels rangs dels paràmetres

Primer hem de definir (potser amb un anàlisi exploratori previ de les dades), els rangs dels dos paràmetres \(\beta_0\) i \(\beta_1\):

range.b0 <- seq(0,5,0.1)

range.b1 <- seq(-3,0,0.01)

Iteració per trobar la millor log-likelihod

Crearem una taula (ML) on guardarem els valors del logaritme de versemblança per a diverses combinacions de valors dels paràmetres

MLE <- NULL

for (i in 1:length(range.b0)) {
  for (j in 1:length(range.b1)) {
    pred.i <- exp(range.b0[i]+range.b1[j] * dist)
    logL.i <- sum(dpois(y, pred.i, log = T))
    res.i <- c(range.b0[i], range.b1[j], logL.i)
    MLE <- rbind(MLE, res.i)
  }
}

MLE <- as.data.frame(MLE)
names(MLE) <- c("b0","b1","logL")

head(MLE)
##         b0    b1      logL
## res.i    0 -3.00 -24317.88
## res.i.1  0 -2.99 -24238.28
## res.i.2  0 -2.98 -24158.68
## res.i.3  0 -2.97 -24079.08
## res.i.4  0 -2.96 -23999.48
## res.i.5  0 -2.95 -23919.88
MLE[which.max(MLE$logL),]
##              b0    b1      logL
## res.i.12639 4.1 -0.02 -26.55552

Que és bastant precís en determinar els coeficients.

Deviança residual i deviança explicada

Un indicador important per analitzar l’ajustament del model a les dades és la deviança residual (deviance) \(D\): \[ \small D=2(\log(L_{m_s})-\log(L_{m_0})) \] on \(\log(L_{m_s})\) és la log-versemblança d’un model saturat (és a dir que s’ajusta perfectament a les dades) i \(\log(L_{m_0})\) la log-versemblança del model (\(m_0\)) observat.

Per tenir un estimador similar a \(R^2\) utilitzem \(D^2\), la proporció de “variació” explicada pel model: \[ \small D^2 =1-\frac{D}{D_{nulo}} \] on \(D\) és la deviança residual del model i \(D_{nulo}\) és la deviança d’un model on tan sol s’ajusta l’ordenada d’origen.

Per a comparar models no anidats, se sol utilitzar el Criteri d’informació d’Akaiké - AIC, es defineix com \[ \small AIC= 2k-2\log(L_\theta) \]

on \(\log(L_\theta)\) és la log-versemblnça del model i \(k\) el nombre de paràmetres ajustats.




  1. Amb el mètode dels mínims quadrats per obtenir la pendent \(\beta_1\) cal calcular la covariança, \(\sigma_{xy}\), de la mostra y la variança de la variable explicativa: \[ \small\hat\beta_1=\frac{\sigma_{xy}}{\sigma_y^2} \] i \(\beta_0\) de les mitjanes de les variables \[ \small\hat\beta=\bar y-\hat\beta_1\bar x. \]↩︎