Referencias

Las notas se han elaborado a partir de vario material:

A.García Perez, Cap. 10 de Estadística aplicada avanzada con R (EAAR), Madrid, 2022.

Base teórica

Ya conocemos los modelos de regresión lineal mediante los cuales se explica una variable teórica $Y$ en función de $k$ variables $X_1,...,X_n$ independientes o covariables, suponiendo una relación lineal del tipo \[ E[Y]=\beta_0+\beta_1x_1+...+\beta_k x_k \]

siendo $Y$ una variable continua con distribución normal.

En el capítulo sobre Regresión logística hemos visto que la variable dependiente es de tipo dicotómico con distribución binomial.

Un tipo intermedio entre la situación de datos binalrios y la de datos contínuos es el Análisis de Regresión de Poisson. Éste supone una variable independiente $Y$ que sigue la distribución de Poisson.

La relación que liga $Y$ con las variables dependientes

\[ \text{log}\;r=\beta_0+\beta_1x_1+...+\beta_k x_k \] donde $r$ es una tasa sobre la variable dependiente $Y$. Este modelo recibe el nombre de Modelo de Regresión de Poisson.

Los objetivos son los mismos:

Estimar los parámetros $\beta_j$
Contrastar cuáles, entre las variables $X_j$, son significativas para explicar $Y$ (o el logaritmo de $r$).

El Modelo de Regresión de Poisson es adecuado en situaciones en que la variable $Y$ es un recuento de datos dentro de grupos en los que se haya dividida la población. También, con otras palabras, el propósito expresar la tasa $r$, esto es el cociente entre indivíduos/unidades (a menudo, personas) a los que acontece un fenómeno (la muerte, por ejemplo) y los individuos/unidades expuestos (por ejemplo, que hayan contraído una enfermedad).

Para ello, se dividirá la población en función de las covariables que se haya seleccionado, como en el ejemplo que sigue.

Regresión de Poisson con R

Ejemplo 1

Se realizó un estudio sobre el número de fallecidos en una determinada población clasificada por $SEXO$ y $EDAD$. Los resultados fueron los siguientes:

Tasas de mortalidad por Sexo y Edad
	Hombres			Mujeres
.	Poblacion	Muertes	Tasa	Poblacion	Muertes	Tasa
30-50	300	33	0.11	250	28	0.112
50-70	200	20	0.10	200	18	0.090
70-90	100	8	0.08	150	12	0.080

En nuestro ejemplo, el modelo a ejecutar será

\[ \text{log}\;r=\beta_0+\beta_1SEXO_i+\beta_2 EDAD_j \]

La clave para poder ejecutar el modelo será que

Poder obtener la relación arrriba descrita
La variable $Y$ siga realmente una distribución de Poisson.

Para la inclusión/exclusión de variables se recomienda un proceso stepwise utilizando el $\text{p-value}$ asociado a cada variable.

Variables regresoras

El set de variables regresoras se construye como en la regresión logística, creando tantas variables cuantas clases tenga cada variable menos 1. Así pues, $EDAD$ (3 clases) tendrá 2 variables y $SEXO$ (2 clases) 1 variable.

library(dplyr)
m <- matrix(c(300,200,100,33,20,8,0.11,0.10,0.08,250,200,150,28,18,12,0.112,0.09,0.08), ncol = 6)

colnames(m) <- c("Poblacion", "Muertes","Tasa","Poblacion","Muertes", "Tasa")
rownames(m) <- c("30-50","50-70","70-90")

df <- as.data.frame(m)[,1:3] %>% 
  tibble::rownames_to_column() %>% 
  bind_rows(tibble::rownames_to_column(as.data.frame(m)[,4:6])) %>% 
  rename(Edad=rowname) %>% 
  mutate(Sexo= ifelse(rownames(.)<4, "Hombre", "Mujer")) %>% 
  select(1,5,4,2,3)

df

##    Edad   Sexo  Tasa Poblacion Muertes
## 1 30-50 Hombre 0.110       300      33
## 2 50-70 Hombre 0.100       200      20
## 3 70-90 Hombre 0.080       100       8
## 4 30-50  Mujer 0.112       250      28
## 5 50-70  Mujer 0.090       200      18
## 6 70-90  Mujer 0.080       150      12

Conversión en variables indicadoras (dummy variables)

Para convertir el dataset automáticamente, podemos usar $\tt model.matrix()$ de $\tt stats$ que ya hemos visto funcionar en otros capítulos. Solo añadiremos la variable cuantitativa $Tasa$ que es la que nos interesa. Notar que la función asigna por defecto los nombres a las variables.

dummy_df <- as.data.frame(model.matrix(~ Edad+Sexo,data=df[,1:2])[,-1]) %>% 
  bind_cols(select(df,3))
dummy_df

##   Edad50-70 Edad70-90 SexoMujer  Tasa
## 1         0         0         0 0.110
## 2         1         0         0 0.100
## 3         0         1         0 0.080
## 4         0         0         1 0.112
## 5         1         0         1 0.090
## 6         0         1         1 0.080

Ejecución del modelo de Poisson

modelPois1 <- glm(Tasa ~ `Edad50-70`+`Edad70-90`+SexoMujer, family = poisson, data = dummy_df)
modelPois1

## 
## Call:  glm(formula = Tasa ~ `Edad50-70` + `Edad70-90` + SexoMujer, family = poisson, 
##     data = dummy_df)
## 
## Coefficients:
## (Intercept)  `Edad50-70`  `Edad70-90`    SexoMujer  
##    -2.18434     -0.15565     -0.32750     -0.02797  
## 
## Degrees of Freedom: 5 Total (i.e. Null);  2 Residual
## Null Deviance:       0.01066 
## Residual Deviance: 0.0004327     AIC: Inf

Lo que permite decir que el modelo de Poisson estimado es

\[ \text{log}\; r_{ij}=-2.18434-0.02797\text{·SexoMujer}-0.15565\text{·Edad50-70}+-0.32750\text{·Edad70-90} \]

Por lo tanto, podemos predecir para cada celda de la tabla los valores correspondientes. Por ejemplo, para hombres de 50-70 años sería: \[ \text{log}\; r_{02}=-2.18434-0.02797\text{·0}-0.15565\text{·1}+-0.32750\text{·0}=-2.18434-0.15565=-2.33999 \] o, que es lo mismo \[ r=\text{exp}(-2.33999)=0.09633 \] por lo cual, si volvemos a la población original, corresponde a una frecuencia estimada de $200\times 0.09633=19.266$, contra $20$ de frecuencia observada.

I así podemos ir especificando, para todas la celdas de la tabla original.

Para eleborar la tabla de valores observados y estimados, podemos utilizar la opción $\tt \$fitted.values$ de $\tt glm()$

Tasa_estimada <- modelPois1$fitted.values

tab <- cbind(df, Tasa_estimada) %>% 
  mutate(Frec_estimada=Tasa_estimada*Poblacion)
tab

##    Edad   Sexo  Tasa Poblacion Muertes Tasa_estimada Frec_estimada
## 1 30-50 Hombre 0.110       300      33    0.11255245     33.765734
## 2 50-70 Hombre 0.100       200      20    0.09632867     19.265734
## 3 70-90 Hombre 0.080       100       8    0.08111888      8.111888
## 4 30-50  Mujer 0.112       250      28    0.10944755     27.361888
## 5 50-70  Mujer 0.090       200      18    0.09367133     18.734266
## 6 70-90  Mujer 0.080       150      12    0.07888112     11.832168

Bondad del ajuste

$\lambda=\sum \frac{(n_{ij}-n_{ij}^t)^2}{n{ij}^t}=$

tab %>% 
 mutate(lamb = (Muertes-Frec_estimada)^2/Frec_estimada) %>% 
  summarise(lambda=sum(lamb))

##       lambda
## 1 0.09293395

Con 6 celdas y 2 parámetros independientes los grados de libertad serán $GL=(3-1)(2-1)=2$.

Para 2 grados de libertad, la probabilidad del estimador será¹ \[P\{\chi_2^2\}=0.955\]

Aceptándose la hiótesis nula de adecuación del modelo regresión de Poisson estimado.

pchisq(modelPois1$deviance, modelPois1$df.residual,lower.tail = F)

## [1] 0.9997837

Ejemplo 2

El número de premios obtenidos por los estudiantes en una escuela secundaria, $num_premios$ trátase de predecirse en base a dos variables independientes: el Programa en el que estaba inscrito el estudiante y la Puntuación sobre 100 puntos obtenida en su examen final de matemáticas.

Los datos se encuentran en el dataset $\tt academic\_awards$ de la librería $\tt gorica$ creado por la Universidad UCLA.

Contiene tres variables:

$\tt num\_awards$: integer Outcome variable; indicates the number of awards earned by students at a high school in a year.
$\tt math$ integer Continuous predictor variable; represents students’ scores on their math final exam
$\tt prog$ factor Categorical predictor variable with three levels, indicating the type of program in which the students were enrolled: “General”, “Academic”, or “Vocational”.

library(gorica)

data(academic_awards)

head(academic_awards)

##    id num_awards       prog math
## 1  45          0 Vocational   41
## 2 108          0    General   41
## 3  15          0 Vocational   44
## 4  67          0 Vocational   42
## 5 153          0 Vocational   40
## 6  51          0    General   42

Análisis visuales previos

Distribución de Poisson

La variable $academic\_awards$ parece seguir aproximadamente dicha distribución:

hist(academic_awards$num_awards)

Se verifica además que media y varianza sean similares

mean(mean(academic_awards$num_awards))

## [1] 0.63

var(academic_awards$num_awards)

## [1] 1.108643

Pero, no solo los valores generales, sino también en las clases de la covariable

academic_awards %>% 
  group_by(prog) %>% 
  summarise(mean= mean(num_awards), var= round(var(num_awards),2))

## # A tibble: 3 × 3
##   prog        mean   var
##   <fct>      <dbl> <dbl>
## 1 Academic    1     1.63
## 2 General     0.2   0.16
## 3 Vocational  0.24  0.27

with(academic_awards, tapply(num_awards, prog, function(x){
  round(c("mean"=mean(x), "var"=var(x)),2)
}))

## $Academic
## mean  var 
## 1.00 1.63 
## 
## $General
## mean  var 
## 0.20 0.16 
## 
## $Vocational
## mean  var 
## 0.24 0.27

Dentro de cada programa

library(ggplot2)
ggplot(academic_awards, aes(num_awards, fill = prog)) +
  geom_histogram(bins = 5, colour= "grey40") +
  facet_wrap(~ prog)+
  scale_x_continuous(breaks = c(0,1,2,3,4,5,6))

Como se puede ver, cada programa presenta distribución de Poisson

Execución del modelo

mod_awards <- glm(num_awards ~ prog + math, 
                  family = "poisson", 
                  data = academic_awards)

summary(mod_awards)

## 
## Call:
## glm(formula = num_awards ~ prog + math, family = "poisson", data = academic_awards)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -4.16327    0.66288  -6.281 3.37e-10 ***
## progGeneral    -1.08386    0.35825  -3.025  0.00248 ** 
## progVocational -0.71405    0.32001  -2.231  0.02566 *  
## math            0.07015    0.01060   6.619 3.63e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 189.45  on 196  degrees of freedom
## AIC: 373.5
## 
## Number of Fisher Scoring iterations: 6

Análisis de los resíduos

Vualización gráfica

par(mfrow=c(1,2))

hist(resid(mod_awards), main = "Histograma de Residuals",xlab = "Residuals")

plot(resid(mod_awards), type = "h", main = "Residuals según índice",ylab = "Residuals")
abline(h = 0, col = "red",lwd=2)

Test de evaluación

Vamos a evaluar la normalidad de los resíduoa con el test Kolgomorov-Smirnov comparando con la normal correspondiente (notar que se pueden especificar la media y la desviación estándar).

ks.test(x=resid(mod_awards), y= "pnorm", mean(resid(mod_awards),sd(resid(mod_awards))))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resid(mod_awards)
## D = 0.16797, p-value = 2.509e-05
## alternative hypothesis: two-sided

El test rechaza claramente la hipótesis $H_0$ de normalidad de los resíduos. Con lo cual habría que buscar un modelo robusto de regresión de Poisson.

En cualquier caso, el test de $\chi^2$ presenta un $\text{p-value}$ que confirmaría la bondad del ajutes a un Modelo de Regresión de Poisson.

pchisq(mod_awards$deviance, mod_awards$df.residual, lower.tail = F)

## [1] 0.6182274

Sistema alternativo

La variable $\tt prog$ en el data frame original está está factorizada como en el data frame del libro EAAR. Sin embargo los $\tt levels$ aplicados en la función $\tt factor()$ son dierentes. Por eso igualamos este aspecto.

academic_awards2 <- academic_awards %>% 
  mutate(prog= factor(prog, 
                      levels = c("General", "Academic","Vocational")
                      ))

Execución del modelo

mod_awards2 <- glm(num_awards ~ prog + math, 
                  family = "poisson", 
                  data = academic_awards2)

summary(mod_awards2)

## 
## Call:
## glm(formula = num_awards ~ prog + math, family = "poisson", data = academic_awards2)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -5.24712    0.65845  -7.969 1.60e-15 ***
## progAcademic    1.08386    0.35825   3.025  0.00248 ** 
## progVocational  0.36981    0.44107   0.838  0.40179    
## math            0.07015    0.01060   6.619 3.63e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 189.45  on 196  degrees of freedom
## AIC: 373.5
## 
## Number of Fisher Scoring iterations: 6

La diferente distribución de niveles parece tener un pequeño impacto en el resultado de la ejecución del modelo que ahora iguala el libro.

Test de evaluación

Vamos a evaluar la normalidad de los resíduoa con el test Kolgomorov-Smirnov comparando con la normal correspondiente (notar que se pueden especificar la media y la desviación estándar).

ks.test(x=resid(mod_awards2), y= "pnorm", mean(resid(mod_awards2),sd(resid(mod_awards2))))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  resid(mod_awards2)
## D = 0.16797, p-value = 2.509e-05
## alternative hypothesis: two-sided

En este caso sigue la diferencia, y no sabemos por qué…

Finalmente, el test de $\chi^2$ da el mismo resultado

pchisq(mod_awards2$deviance, mod_awards2$df.residual, lower.tail = F)

## [1] 0.6182274

pchisq(0.09293395,df=2,lower.tail = F)↩︎

Regresión de Poisson

Notes d’estadística aplicada - Nota Núm. 10

Stefano Sanfilippo

Març 2026

Referencias

Base teórica

Regresión de Poisson con R

Ejemplo 1

Variables regresoras

Conversión en variables indicadoras (dummy variables)

Ejecución del modelo de Poisson

Bondad del ajuste

Ejemplo 2

Análisis visuales previos

Distribución de Poisson

Execución del modelo

Análisis de los resíduos

Vualización gráfica

Test de evaluación

Sistema alternativo

Execución del modelo

Test de evaluación