Comparación de la pendiente de dos modelos (normal y poisson) de riqueza de especies

En un estudio que está llevando a cabo Fabio Suzart de Alburquerque (Departamento de Ecología, Universidad de Alcalá) se quiere comparar la pendiente de la recta de regresión de la riqueza de especies nativas y exóticas en el Reino Unido frente a una variable explicativa relacionada con la entrada de energía en el sistema (el primer eje de un análisis de componente principales), correlacionado fundamentalmente con la temperatura media anual, la temperatura media del mes de enero, y la evapotranspiración potencial). El principal problema con el que se enfrenta es que la distribución de la riqueza de nativas sigue una distribución normal, mientras que la distribución de especies exóticas sigue una distribución más de tipo Poisson o log-normal (ver Figura 1).

El siguiente código de R permite la lectura de los datos (accesibles con conexión a internet), el ajuste de dos modelos "normales" para la riqueza de nativas y exóticas, y la creación de: (1) histogramas de los residuos de los modelos de riqueza de especies nativas y exóticas; y (2) gráficos de correlación de estas variables con el primer eje del análisis de componentes principales (PCA).

data <- read.table(url("http://archivos-para-subir.googlegroups.com/web/datos+regresion.txt?gda=x8KpMUUAAACkOfwqHAVd4YqgfIB09GDRRDpDrIn9MKnJm3pWH7Uz3yEZ0eBGq9HKqr7NTMdj8uts7Hm76flVp_vYI1lbK8kaGu1iLHeqhw4ZZRj3RjJ_-A"), header = T, sep = "\t")

glm.native <- glm(data$NAT_RICH ~ data$Factor1.energy) glm.exotic <- glm(data$ALIEN_RICH ~ data$Factor1.energy)
par(mfcol= c(2, 2))
hist(residuals(glm.native), main = "Histogram of native SR residuals", xlab = "Number of species")
par(new = T)
plot(density(residuals(glm.native)), xlab = "", ylab = "", main = "", axes = F)
hist(residuals(glm.exotic), main = "Histogram of exotic SR residuals", xlab = "Number of species")
par(new = T)
plot(density(residuals(glm.exotic)), xlab = "", ylab = "", main = "", axes = F)
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
plot(data$ALIEN_RICH ~ data$Factor1.energy, main = "Exotic species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")

Figura 1. Histogramas de los residuos de los modelos de especies nativas (arriba izquierda) y exóticas (abajo izquierda) y gráficos de correlación de la riqueza de especies nativas (arriba derecha) y exóticas (abajo derecha) con el primer eje del análisis de componentes principales (PCA).

Una posible solución a este problema es utilizar modelos lineales generalizados (GLM) con una distribución de errores de tipo Poisson y una función de vínculo de tipo logarítmica para linealizar la relación entre la riqueza de especies nativas y la variable independiente. Esto impide comparar la pendiente de ambas rectas de regresión de manera directa pero nos permite, en cambio, establecer hipótesis sobre los procesos que generan la respuesta de las variables dependientes (riqueza de nativas frente a riqueza de exóticas).

glm.native <- glm(data$NAT_RICH ~ data$Factor1.energy)
glm.exotic <- glm(data$ALIEN_RICH ~ data$Factor1.energy, family = poisson)

El primer modelo es un modelo normal. Hay que considerar que los GLM no son más que una extensión de los modelos lineales, por lo que los modelos lineales que comúnmente se utilizan (ANOVA, regresión) no son más que un caso particular de GLM en dónde la distribución de errores es de tipo gaussiana (normal). Con la función summary() se obtiene un resumen de los modelos ajustados.

summary(glm.native)

Call:
glm(formula = data$NAT_RICH ~ data$Factor1.energy)

Deviance Residuals:
Min 1Q Median 3Q Max
-428.470 -45.604 5.102 47.498 218.190

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 421.489 1.515 278.12 <2e-16>
data$Factor1.energy 59.039 1.516 38.95 <2e-16>
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 5172.153)

Null deviance: 19483365 on 2251 degrees of freedom
Residual deviance: 11637345 on 2250 degrees of freedom
AIC: 25652

Number of Fisher Scoring iterations: 2

summary(glm.exotic)

Call:
glm(formula = data$ALIEN_RICH ~ data$Factor1.energy, family = poisson)

Deviance Residuals:
Min 1Q Median 3Q Max
-18.7702 -4.5256 -0.3740 3.1827 23.9153

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.900846 0.001945 2519.6 <2e-16 ***
data$Factor1.energy 0.543437 0.002224 244.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 158010 on 2251 degrees of freedom
Residual deviance: 87121 on 2250 degrees of freedom
AIC: 101806

Number of Fisher Scoring iterations: 5

Como se puede observar, las pendientes de la recta de regresión de ambos modelos no son comparables, entre otras cosas, porque en el modelo Poisson la variable respuesta (riqueza de exóticas) está transformada logarítmicamente. Esto ya nos da una idea sobre el comportamiento de ambos tipos de riqueza de especies. Así, mientras que la riqueza de especies nativas varía linealmente con el primer eje del PCA (energía), la riqueza de exóticas cambia exponencialmente, es decir, que para valores bajos de la variable explicativa, el incremento de la riqueza será relativamente pequeño, pero para valores altos de la variable explicativa, el incremento de la riqueza de especies será muy alto.

Esto se puede observar en el modelo ajustado a los puntos para cada uno de los dos tipos de riqueza.

par(mfrow=c(1,2))
yv.native <- predict(glm.native, type= "response")
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.native))
yv.exotic <- predict(glm.exotic, type= "response")
plot(data$ALIEN_RICH ~ data$Factor1.energy, main = "Exotic species richness (response) vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.exotic))


Y alternativamente podemos dibujar las predicciones del modelo en escala logarítmica. Esta relación es lineal, pero hay que tener en cuenta que la variable respuesta (riqueza de especies) está en escala logarítmica, por lo que la comparación de las pendientes de ambas rectas no es posible.

par(mfrow=c(1,2))
yv.native <- predict(glm.native, type= "response")
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.native))
yv.exotic2 <- predict(glm.exotic, type= "link")
plot(log(data$ALIEN_RICH) ~ data$Factor1.energy, main = "Exotic species richness (link) vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "log(Number of species)")
lines(sort(data$Factor1.energy), sort(yv.exotic2))

Este estudio forma parte del siguiente trabajo: Suzart de Alburquerque, F., Castro, P. Cayuela, L. & Rodríguez, M.A. Assessing the influence of environmental and human-related factors on native and exotic species richness in Great Britain. In prep.