Una distribución Poisson podría considerarse adecuada dado que cuenta el número de eventos dentro de un determinado espacio de tiempo, en este caso se cuentan el número de fibras de asbesto atrapadas en 200 filtros. El modelo establece que nuestra variable \(Y\) sigue la siguiente distribución: \[ \mathbb{P}[Y=k] = e^{-\lambda}\frac{\lambda^k}{k!} \] \(Y\) cuenta el número de fibras de asbesto en un filtro.
Usando este modelo supondríamos que la media es igual a la varianza. El estimador de \(\lambda\) por máxima verosimilitud sería el promedio, el cual nos dicen es \(27.7\), es decir, \(\hat{\lambda} = 27.7\).
Se realiza la prueba de la ji-cuadrada para comprobar que viene de una v.a. Poisson, en esta prueba \(H_0 = \text{Y distribuye } Poisson(27.7)~\) y \(~H_a = \text{Y no se distribuye } Poisson(27.7)\) donde nuestro estadístico de prueba es \(T = 35.39\). Entonces: \[ T = 35.39 > 11.0705 = w_{0.95} \quad \text{donde } w_{0.95} \text{ representa el cuantil 0.95 de una v.a. } \chi^2_{(7-1-1)} \]
\(\therefore\) Se rechaza \(H_0\).
Trabajamos con los datos del dataset ships de la librería MASS. Primero removemos los que tienen 0 meses de servicio ya que no aportan información ya que no han sufrido exposición al oleaje. Se ajustó el siguiente modelo:
library(MASS)
library(tidyverse)
data(ships)
data <- ships
data <- data %>% filter(service != 0)
fit <- glm(incidents ~ offset(log(service)) + type + year + period, data = data)
summary(fit)##
## Call:
## glm(formula = incidents ~ offset(log(service)) + type + year +
## period, data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -21.8803 -2.5218 -0.6538 2.3275 22.2979
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.00182 21.22377 -0.189 0.852
## typeB 27.03993 4.58782 5.894 2.8e-06 ***
## typeC -4.32491 4.58782 -0.943 0.354
## typeD -2.93408 4.58782 -0.640 0.528
## typeE -0.10767 4.79552 -0.022 0.982
## year -0.08462 0.30201 -0.280 0.781
## period 0.13171 0.20665 0.637 0.529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 73.66828)
##
## Null deviance: 6782.3 on 33 degrees of freedom
## Residual deviance: 1989.0 on 27 degrees of freedom
## AIC: 250.84
##
## Number of Fisher Scoring iterations: 2
Comprobaremos si nuestro modelo es estadístcamente significativo.
## [1] 4793.23
## [1] 0
Este valor nos indica que efectivamente podemos rechazar la hipótesis nula y concluir que nuestro modelo es estadísticamente significativo.
Observando las betas de nuestro modelo.
## (Intercept) typeB typeC typeD typeE year
## -4.00181782 27.03992866 -4.32491454 -2.93408422 -0.10767325 -0.08461571
## period
## 0.13170807
Notamos que cambiando a comparación de un barco tipo A, los tipo B sufren más incidentes dado que el coeficiente es positivo, y además también a comparación de los tipo A, los barcos tipo E, D y C son sufren menos incidentes, siendo “más seguros” los del tipo C. También notamos que conforme aumentan el año de construcción, los barcos tienen menos incidentes. Finalmente vemos que durante el primer periodo de operación se sufren menos incidentes que en el segundo, es decir, si el periodo de operación es más reciente, más incidentes tendrá.
Sin embargo, algo importante que hay que notar es que el modelo no es bueno, ya que al realizar la prueba de bondad de ajuste:
## [1] 0
## [1] 73.66828
Este valor es 0, entonces no podemos rechazar la hipótesis de que el modelo no es un buen ajuste. Por lo tanto no es un buen modelo. Además de que se tiene una sobredispersión mucho más grande que 1.
Aún así, es buena idea usar la offset en la variable service ya que los datos que tenemos no son de una misma ventana de tiempo, esto es, tenemos desde barcos con 45 meses de servicio hasta uno con 44,882. Esto nos permite comparar nuestras variables explicativas adecuadamente ya que los incidentes si van en función de los meses de servicio y cada barco tiene su propio número de meses.