library(tidyverse)
library(flextable)
Regresion de Poisson
Ejemplo: Ecología reproductiva del cangrejo de herradura
Info sobre esta especie, Más info
Durante la temporada de desove, las hembras de esta especie migran a la costa para realizar la puesta, con un macho enganchado en su cola. Cavan en la arena, donde depositan los huevos que son fertilizados externamente, tanto por el macho enganchado a la hembra como por otros machos que se reúnen alrededor de la pareja. Estos otros machos reciben el nombre de satélites. El archivo crabs.txt contiene datos de un estudio realizado sobre esta especie en el golfo de Mexico. Las variables contenidas en el archivo son:
- Crab: Número de identificación de la hembra.
- y: Número de machos satélite.
- Weight: Peso
- Width: Anchura del caparazón
- Color: color del caparazón, de 1 a 4. El color se relaciona con la edad, siendo 1 el color más claro, correspondiente a ejemplares más jóvenes y 4 el color más oscuro.
- Spine: Estado de las espinas laterales del cangrejo: 1: buen estado en ambos lados, 2: Un lado en buen estado, el otro dañado, 3: ambos lados dañados.
Leemos el archivo de datos y mostramos sus primeras 15 filas:
=read.table("datos/Crabs.txt",header=TRUE)
cangrejoshead(cangrejos,15)
crab y weight width color spine
1 [1,] 8 3.05 28.3 2 3
2 [2,] 0 1.55 22.5 3 3
3 [3,] 9 2.30 26.0 1 1
4 [4,] 0 2.10 24.8 3 3
5 [5,] 4 2.60 26.0 3 3
6 [6,] 0 2.10 23.8 2 3
7 [7,] 0 2.35 26.5 1 1
8 [8,] 0 1.90 24.7 3 2
9 [9,] 0 1.95 23.7 2 1
10 [10,] 0 2.15 25.6 3 3
11 [11,] 0 2.15 24.3 3 3
12 [12,] 0 2.65 25.8 2 3
13 [13,] 11 3.05 28.2 2 3
14 [14,] 0 1.85 21.0 4 2
15 [15,] 14 2.30 26.0 2 1
Nuestro objetivo en este capítulo será predecir el número esperado de machos satélite alrededor de una hembra, utilizando el resto de variables como variables predictoras.
Distribución de Poisson
Una variable \(Y\) sigue una distribución de probabilidad de Poisson cuando \[P\left(Y=k\right)=\frac{\lambda^{k}}{k!}e^{-\lambda},\,\,\,\,\,k\ge0\] Aquí el valor \(\lambda\) representa el valor esperado de la variable \(Y\).
El modelo de regresión de Poisson se aplica cuando la variable \(Y\) es una variable de recuento que depende de unas variables explicativas \(X_i\), de tal modo que el valor esperado de \(Y\) para un conjunto de valores de \(X_i\) puede aproximarse mediante una distribución de Poisson de parámetro \(\lambda(\mathbf{X})=exp(\beta_0+\beta_1 X_1 + \dots + \beta_p X_p)\).
Para interpretar los parámetros de este modelo basta tener en cuenta que de acuerdo con la descripción anterior:
\[log\left(\lambda(\mathbf{X})\right)=\beta_0+\beta_1 X_1 + \dots + \beta_p X_p\]
por lo que si los valores de todas las \(X_i\) se mantienen fijos, salvo el valor de \(X_j\), que se incrementa en una unidad, entonces:
\[ \log\left(\frac{\lambda\left(X_{1},\dots,X_{j}+1,\dots,X_{p}\right)}{\lambda\left(X_{1},\dots,X_{j},\dots,X_{p}\right)}\right) =\log\left(\lambda\left(X_{1},\dots,X_{j}+1,\dots,X_{p}\right)\right)-\log\left(\lambda\left(X_{1},\dots,X_{j},\dots,X_{p}\right)\right)= \]
\[ =\beta_{j}\left(X_{j}+1\right)-\beta_{j}X_{j}=\beta_{j} \]
o lo que es lo mismo:
\[\frac{\lambda\left(X_{1},\dots,X_{j}+1,\dots,X_{p}\right)}{\lambda\left(X_{1},\dots,X_{j},\dots,X_{p}\right)}=e^{\beta_{j}}\]
Por tanto \(e^{\beta_j}\) representa el incremento proporcional en el valor de \(Y\) por cada unidad de incremento en el valor de \(X_j\). Por ejemplo, si fuese \(\beta_j=1.388\) se tiene que \(e^{1.388}=4\). Ello significa que por cada unidad que se incrementa el valor de \(X_j\) se cuadruplica el valor de \(Y\) (siempre y cuando se mantengan fijos los valores del resto de las variables).
Estimación del modelo de Poisson
Predicción del número de satélites a partir del peso del cangrejo hembra
El modelo se estima por máxima verosimilitud. Utilizaremos R para ajustar el modelo de regresión de Poisson para predecir el número de satélites macho en función del peso del cangrejo hembra:
=glm(y~weight,data=cangrejos,family=poisson(link=log))
modPoissummary(modPois)
Call:
glm(formula = y ~ weight, family = poisson(link = log), data = cangrejos)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.9307 -1.9981 -0.5627 0.9298 4.9992
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.42841 0.17893 -2.394 0.0167 *
weight 0.58930 0.06502 9.064 <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: 632.79 on 172 degrees of freedom
Residual deviance: 560.87 on 171 degrees of freedom
AIC: 920.16
Number of Fisher Scoring iterations: 5
Calculamos los valores de \(e^{\beta_i}\) y sus intervalos de confianza:
exp(coef(modPois))
(Intercept) weight
0.6515473 1.8027335
exp(confint(modPois))
2.5 % 97.5 %
(Intercept) 0.4597023 0.9268994
weight 1.5835992 2.0431615
Por tanto, por cada kg de peso adicional del cangrejo hembra, el número medio de satélites macho se multiplica por 1.8, si bien con una confianza del 95% podemos decir que dicho multiplicador es un número entre 1.58 y 2.04.
Gráficamente:
# Generamos valores de peso de cangrejo en el rango observado para
# predecir el número de satélites a partir del modelo anterior
=data.frame(weight=seq(1,5.5,by=0.05))
weightNew$ypred=predict(modPois,newdata=weightNew,type="response")
weightNew
# Mostramos los puntos observados y las predicciones en el rango
# de valores que acabamos de calcular:
%>%
cangrejos ggplot(aes(weight,y)) +
geom_point(shape=21) +
geom_line(data=weightNew,aes(x=weight,y=ypred),color="red",lwd=1.2)
Predicción del número de satélites añadiendo otras variables
=glm(y~weight+width+color+spine,data=cangrejos,family=poisson(link=log))
modPoissummary(modPois)
Call:
glm(formula = y ~ weight + width + color + spine, family = poisson(link = log),
data = cangrejos)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.0126 -1.8846 -0.5406 0.9448 4.9602
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.52848 0.94867 -0.557 0.57748
weight 0.47246 0.16490 2.865 0.00417 **
width 0.02753 0.04794 0.574 0.56588
color -0.18493 0.06652 -2.780 0.00544 **
spine 0.03998 0.05681 0.704 0.48160
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 632.79 on 172 degrees of freedom
Residual deviance: 551.85 on 168 degrees of freedom
AIC: 917.15
Number of Fisher Scoring iterations: 6
Las variables width
y spine
no resultan significativas. Analizamos si podemos eliminarlas:
=glm(y~weight+color,data=cangrejos,family=poisson(link=log))
modPois2anova(modPois,modPois2, test="Chisq")
Analysis of Deviance Table
Model 1: y ~ weight + width + color + spine
Model 2: y ~ weight + color
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 168 551.85
2 170 552.79 -2 -0.94042 0.6249
Efectivamente, las dos variables no aportan capacidad predictiva al modelo, por lo cuál nos quedamos definitivamente con el modelo modPois2
:
summary(modPois2)
Call:
glm(formula = y ~ weight + color, family = poisson(link = log),
data = cangrejos)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.9785 -1.9159 -0.5471 0.9181 4.8338
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.08855 0.25443 0.348 0.72783
weight 0.54588 0.06749 8.088 6.05e-16 ***
color -0.17282 0.06155 -2.808 0.00499 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 632.79 on 172 degrees of freedom
Residual deviance: 552.79 on 170 degrees of freedom
AIC: 914.09
Number of Fisher Scoring iterations: 6
Los coeficientes \(e^\beta_j\) así como sus intervalos de confianza son:
exp(coef(modPois2))
(Intercept) weight color
1.0925855 1.7261244 0.8412914
exp(confint(modPois2))
2.5 % 97.5 %
(Intercept) 0.6637615 1.7997895
weight 1.5087970 1.9656124
color 0.7449264 0.9482147
Así pues, para un peso fijo, el incremento de un grado de color produce que el número de machos disminuya a un 84% del valor anterior; dicho de otra forma, el envejecimiento de la hembra produce que tenga menos machos satélite.
Asimismo, para dos hembras de la misma edad, cada kg de peso adicional hace que el número de machos satélites se multiplique aproximadamente por 1.73. Este valor no es muy diferente del estimado más arriba (1.8) lo que indica que la adición del color al modelo no modifica el efecto del peso; interpretado desde un punto de vista biológico, los efectos del peso y el color actúan independientemente sobre el número de machos satélite.
Por último, la representación gráfica de este modelo sería la siguiente:
# Genero las predicciones del modelo:
=data.frame(weight=seq(1,5.5,length=100))
weightNew<- NULL
weiCol for (k in 1:4){
<- rbind(weiCol,cbind(weightNew,color=k))
weiCol
}<- weiCol %>%
weiCol mutate(color=factor(color),
ypred=predict(modPois2,newdata=weiCol,type="response"))
# Hago la gráfica:
%>%
cangrejos mutate(color=factor(color)) %>%
ggplot(aes(weight,y,color=color)) +
geom_point() +
geom_line(data=weiCol,aes(x=weight,y=ypred,color=color),lwd=1.2)
O también:
%>%
cangrejos mutate(color=factor(color)) %>%
ggplot(aes(weight,y,color=color)) +
geom_point() +
geom_line(data=weiCol,aes(x=weight,y=ypred,color=color),lwd=1.2) +
facet_wrap(~factor(color),ncol=2)
Bondad de ajuste del modelo de Poisson
El siguiente contraste permite determinar la bondad de ajuste del modelo; un p-valor mayor que 0.05 indica que hay un buen ajuste del modelo a los datos. En caso contrario, los valores observados se apartan significativamente de los valores predichos por el modelo, que por tanto debe desecharse. Dicho apartamiento puede tener diversas causas: falta de linealidad en el logaritmo de la respuesta, exceso de variabilidad en la misma, o quizás ausencia en el modelo de variables importantes que debieran tenerse en cuenta (y que el investigador ignora)
En nuestro ejemplo del análisis del número de machos satélites de cada hembra ponedora tenemos:
# Test de bondad de ajuste para el modelo de Poisson
pchisq(modPois2$deviance, df=modPois2$df.residual, lower.tail=FALSE)
[1] 4.865544e-42
Por tanto nos encontramos ante una falta de ajuste importante del modelo.
Sobredispersión
Una de las posibles causas de la falta de ajuste de un modelo de Poisson es la sobredispersión. La distribución de Poisson se caracteriza porque su esperanza y su varianza coinciden; cuando se ajusta un modelo de Poisson a un conjunto de datos puede ocurrir que ambos valores difieran de manera significativa. Se dice entonces que el modelo presenta sobredispersión. En el paquete AER de R hay un test que permite contrastar la presencia de sobredispersión en un modelo. En nuestro ejemplo el resultado de dicho test es:
# Test de sobredispersión
library(AER)
dispersiontest(modPois2)
Overdispersion test
data: modPois2
z = 5.1735, p-value = 1.149e-07
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion
3.150114
Por tanto nos encontramos ante un caso de sobredispersión; la variabilidad de las observaciones es mayor que la esperada bajo la hipótesis del modelo de Poisson
El modelo de regresión binomial negativa
Cuando se detecta sobredispersión en un modelo de Poisson, se puede utilizar como alternativa la regresión binomial negativa. Recuérdese que una variable aleatoria \(V\) sigue una distribución binomial negativa de parámetros \(r\) y \(p\) si: \[P\left(V=k\right)=\binom{k-1}{r-1}p^{r}\left(1-p\right)^{k-r}\] La esperanza y varianza de esta distribución son:
\[\mu=E\left[X\right]=r\cdot\frac{1}{p}\] \[\sigma^{2}=r\cdot\frac{1-p}{p^{2}}\]
Lo que hace interesante a esta distribución es que toma valores discretos y, al contrario que en la distribución de Poisson, la esperanza y la varianza son distintas, lo que permite modelar situaciones más generales.
En el caso de la regresión binomial negativa se asume que la variable respuesta \(Y\), dados los valores de unas variables explicativas \(X_1, X_2,\dots,X_p\) sigue una distribución binomial negativa de esperanza \(E[Y|X_1,X_2,\dots,X_p] = \beta_0 +\beta_1X_1+\dots+\beta_p X_p\). La interpretación de los coeficientes del modelo es, pues, análoga al caso del modelo de regresión de Poisson, con la ventaja en este caso de que la posible sobredispersión queda incluida en el modelo a través de la estimación de la varianza de la distribución binomial negativa.
Para ajustar un modelo de regresión binomial negativa a nuestros datos de cangrejos, el procedimiento a seguir en R sería el siguiente:
library(MASS) # Librería que contiene el procedimiento para la regresión BN
<- glm.nb(y~weight+width+color+spine,data=cangrejos) modBN
Vemos que al considerar todas las variables, ninguna resulta significativa. Ejecutamos un procedimiento paso a paso:
<- stepAIC(modBN) modBN1
Start: AIC=756.39
y ~ weight + width + color + spine
Df AIC
- spine 1 754.41
- width 1 754.43
- color 1 756.37
<none> 756.39
- weight 1 756.93
Step: AIC=754.41
y ~ weight + width + color
Df AIC
- width 1 752.45
<none> 754.41
- color 1 754.57
- weight 1 754.93
Step: AIC=752.45
y ~ weight + color
Df AIC
<none> 752.45
- color 1 752.64
- weight 1 766.68
summary(modBN1)
Call:
glm.nb(formula = y ~ weight + color, data = cangrejos, init.theta = 0.9555421418,
link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8635 -1.3849 -0.3215 0.4446 2.1361
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.3220 0.5540 -0.581 0.561
weight 0.7072 0.1612 4.387 1.15e-05 ***
color -0.1734 0.1199 -1.445 0.148
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(0.9555) family taken to be 1)
Null deviance: 219.50 on 172 degrees of freedom
Residual deviance: 196.64 on 170 degrees of freedom
AIC: 754.45
Number of Fisher Scoring iterations: 1
Theta: 0.956
Std. Err.: 0.174
2 x log-likelihood: -746.452
El método AIC ha elegido las dos variables peso y color, aunque la variable color no resulta significativa y podríamos suprimirla:
<- glm.nb(y~weight,data=cangrejos)
modBN2 anova(modBN1,modBN2,test="Chisq")
Likelihood ratio tests of Negative Binomial Models
Response: y
Model theta Resid. df 2 x log-lik. Test df LR stat.
1 weight 0.9310592 171 -748.6437
2 weight + color 0.9555421 170 -746.4524 1 vs 2 1 2.19128
Pr(Chi)
1
2 0.1387939
summary(modBN2)
Call:
glm.nb(formula = y ~ weight, data = cangrejos, init.theta = 0.9310592338,
link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8394 -1.4122 -0.3247 0.4744 2.1279
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.8647 0.4048 -2.136 0.0327 *
weight 0.7603 0.1578 4.817 1.45e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(0.9311) family taken to be 1)
Null deviance: 216.43 on 172 degrees of freedom
Residual deviance: 196.16 on 171 degrees of freedom
AIC: 754.64
Number of Fisher Scoring iterations: 1
Theta: 0.931
Std. Err.: 0.168
2 x log-likelihood: -748.644
Bondad de ajuste:
# Test de bondad de ajuste para el modelo Binomial Negativo
pchisq(modBN2$deviance, df=modBN2$df.residual, lower.tail=FALSE)
[1] 0.09100026
La presentación e interpretación de resultados del modelo binomial negativo es idéntica que en el caso de Poisson:
cbind(exp(coef(modBN1)),exp(confint(modBN1)))
2.5 % 97.5 %
(Intercept) 0.7246848 0.2340738 2.213001
weight 2.0282519 1.4428449 2.891609
color 0.8408252 0.6697334 1.056964
cbind(exp(coef(modBN2)),exp(confint(modBN2)))
2.5 % 97.5 %
(Intercept) 0.421196 0.173036 1.007211
weight 2.138872 1.523032 3.045503
Es decir, el incremento de un kg de peso en la hembra multiplica por 2 (intervalo de confianza entre 1.5 y 3) el número de machos satélite.
Otro ejemplo de regresión de Poisson
library(GGally)
En un estudio sobre la población de un parásito se hizo un recuento de parásitos en 15 localizaciones con diversas condiciones ambientales. Los datos obtenidos son los siguientes:
=data.frame(
parasitosTemperatura=c(15,16,24,13,21,16,22,18,20,16,28,27,13,22,23),
Humedad=c(70,65,71,64,84,86,72,84,71,75,84,79,80,76,88),
Recuento=c(156,157,177,145,197,184,172,187,157,169,200,193,167,170,192)
)
Representamos gráficamente estos datos:
<- function(data, mapping, ...){
my_fn <- ggplot(data = data, mapping = mapping) +
p geom_point() +
geom_smooth(method=lm, fill="blue", color="red", ...)
p
}
ggpairs(parasitos,lower = list(continuous = my_fn))
library(plot3D)
with(parasitos,
scatter3D(Temperatura, Humedad, Recuento, bty = "g",
colkey = FALSE,theta = 30, phi = 30, pch=19,
xlab="Temp",ylab="Hum",zlab="Recuento",
ticktype="detailed",col="blue",type=""))
# x, y, z variables
<- parasitos$Temperatura
x <- parasitos$Humedad
y <- parasitos$Recuento
z # Compute the linear regression (z = ax + by + d)
<- lm(z ~ x + y)
fit # predict values on regular xy grid
= 26
grid.lines <- seq(min(x), max(x), length.out = grid.lines)
x.pred <- seq(min(y), max(y), length.out = grid.lines)
y.pred <- expand.grid( x = x.pred, y = y.pred)
xy <- matrix(predict(fit, newdata = xy),
z.pred nrow = grid.lines, ncol = grid.lines)
# fitted points for droplines to surface
<- predict(fit)
fitpoints # scatter plot with regression plane
scatter3D(x, y, z, pch = 19, cex = 1.2, col=gg.col(20),
theta = 35, phi = 25, ticktype = "detailed",
xlab = "Temp", ylab = "Hum", zlab = "Recuento",
surf = list(x = x.pred, y = y.pred, z = z.pred,
facets = NA, fit = fitpoints), clab=c("Numero de","parásitos"),
main = "Recuento de Parásitos")
Ajuste de un modelo lineal
=lm(Recuento~Temperatura+Humedad,data=parasitos)
modlinsummary(modlin)
Call:
lm(formula = Recuento ~ Temperatura + Humedad, data = parasitos)
Residuals:
Min 1Q Median 3Q Max
-9.8617 -2.0406 0.4319 2.9881 8.5047
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.7115 14.3725 1.789 0.098876 .
Temperatura 1.5818 0.3203 4.939 0.000343 ***
Humedad 1.5424 0.1995 7.731 5.32e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.351 on 12 degrees of freedom
Multiple R-squared: 0.914, Adjusted R-squared: 0.8996
F-statistic: 63.75 on 2 and 12 DF, p-value: 4.051e-07
Ajuste de un modelo de Poisson
=glm(Recuento~Temperatura+Humedad,data=parasitos,family="poisson")
modPoissummary(modPois)
Call:
glm(formula = Recuento ~ Temperatura + Humedad, family = "poisson",
data = parasitos)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.72308 -0.14742 0.04359 0.23994 0.61624
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.296214 0.207132 20.741 < 2e-16 ***
Temperatura 0.009037 0.004477 2.018 0.04357 *
Humedad 0.008964 0.002834 3.163 0.00156 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 22.9606 on 14 degrees of freedom
Residual deviance: 1.9468 on 12 degrees of freedom
AIC: 112.92
Number of Fisher Scoring iterations: 3
cbind(exp(coef(modPois)),exp(confint(modPois)))
2.5 % 97.5 %
(Intercept) 73.421327 48.853974 110.039597
Temperatura 1.009077 1.000256 1.017968
Humedad 1.009005 1.003410 1.014622
El modelo significa que, para una humedad fija, el incremento de un grado de temperatura significa que el número de parásitos se multiplica por 1.009. Asimismo, para una temperatura fija, el incremento de una unidad en el valor de humedad implica que el número de parásitos se multiplica por 1.009
Comparación del ajuste conseguido con ambos modelos
=data.frame(observado=parasitos$Recuento,
predslineal=predict(modlin),
poisson=exp(predict(modPois)))
<- pivot_longer(preds,cols=2:3,names_to = "modelo",values_to = "predicho")
preds3 ggplot(preds3,aes(x=observado,y=predicho,color=modelo))+geom_point()+
geom_abline(intercept=0,slope=1,col="lightblue",lwd=1.2)
Bondad de ajuste del modelo
# Test de bondad de ajuste para el modelo de Poisson
pchisq(modPois$deviance, df=modPois$df.residual, lower.tail=FALSE)
[1] 0.9994831
Sobredispersión
No es necesario hacer el test dado que hay un buen ajuste general. En cualquier caso, si lo llevamos a cabo obtenemos:
# Test de sobredispersión
dispersiontest(modPois)
Overdispersion test
data: modPois
z = -19.935, p-value = 1
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion
0.129322