Regresion de Poisson

library(tidyverse)
library(flextable)

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:

cangrejos=read.table("datos/Crabs.txt",header=TRUE)
head(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:

modPois=glm(y~weight,data=cangrejos,family=poisson(link=log))
summary(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
weightNew=data.frame(weight=seq(1,5.5,by=0.05))
weightNew$ypred=predict(modPois,newdata=weightNew,type="response")

# 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

modPois=glm(y~weight+width+color+spine,data=cangrejos,family=poisson(link=log))
summary(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:

modPois2=glm(y~weight+color,data=cangrejos,family=poisson(link=log))
anova(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:
weightNew=data.frame(weight=seq(1,5.5,length=100))
weiCol <- NULL
for (k in 1:4){
  weiCol <- rbind(weiCol,cbind(weightNew,color=k))
}
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
modBN <- glm.nb(y~weight+width+color+spine,data=cangrejos)

Vemos que al considerar todas las variables, ninguna resulta significativa. Ejecutamos un procedimiento paso a paso:

modBN1 <- stepAIC(modBN)

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:

modBN2 <- glm.nb(y~weight,data=cangrejos)
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:

parasitos=data.frame(
  Temperatura=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:

my_fn <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    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
x <- parasitos$Temperatura
y <- parasitos$Humedad
z <- parasitos$Recuento
# Compute the linear regression (z = ax + by + d)
fit <- lm(z ~ x + y)
# predict values on regular xy grid
grid.lines = 26
x.pred <- seq(min(x), max(x), length.out = grid.lines)
y.pred <- seq(min(y), max(y), length.out = grid.lines)
xy <- expand.grid( x = x.pred, y = y.pred)
z.pred <- matrix(predict(fit, newdata = xy), 
                 nrow = grid.lines, ncol = grid.lines)
# fitted points for droplines to surface
fitpoints <- predict(fit)
# 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

modlin=lm(Recuento~Temperatura+Humedad,data=parasitos)
summary(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

modPois=glm(Recuento~Temperatura+Humedad,data=parasitos,family="poisson")
summary(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

preds=data.frame(observado=parasitos$Recuento,
                 lineal=predict(modlin),
                 poisson=exp(predict(modPois)))
preds3 <- pivot_longer(preds,cols=2:3,names_to = "modelo",values_to = "predicho")
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