Problema hipotético
El Instituto de Ciencias Aplicadas y Tecnología de la UNAM trabaja en un proyecto para desarrollar estándares para la concentración de asbesto presente en el agua destinada a consumo humano.
En un experimento controlado se disolvieron asbestos en agua y posteriormente esta se hizo pasar por un filtro con el objetivo de contar el número de fibras de asbesto que quedaban atrapadas en el mismo.
El procedimiento se realizó con 200 filtros, lo que dió como resultados un promedio de 27.7 fibras de asbestos por filtro y lo siguiente:
| Número de fibras | 0-10 | 11-15 | 16-20 | 21-24 | 25-27 | 28-30 | >=31 | suma |
|---|---|---|---|---|---|---|---|---|
| Número observado de filtros |
2 | 1 | 36 | 52 | 50 | 39 | 20 | 200 |
| Número esperado de filtros |
0.12 | 2.43 | 34.62 | 57.51 | 45.36 | 32.62 | 27.34 | 200 |
| 30.25 | 0.85 | 0.06 | 0.53 | 0.48 | 1.25 | 1.97 | 35.39 |
El consultor estadístico asumió un modelo Poisson para el número observado de fibras.
Se utilizó un modelo Poisson pues queremos contar el número de asbestos por filtro, entonces esta será nuestra variable explicativa
Los supuestos son los siguientes:
La varianza de la variable respuesta no es independiente de la media.
Los errores no siguen una distribución Normal.
Los ceros que aparecen en la variable respuesta dan problemas a la hora de transformar la variables.
Estimador para
El consultor desea probar que el modelo es el adecuado usando la prueba chi-cuadrada. Con los resultados de la tabla se obtiene que la estadística de prueba es 35.39
Realizamoss la prueba de la chi-cuadrada
Con 6 grados de libertad (7-1 categorias)
esperados <-c(2,1, 36, 52, 50, 39, 20)
observados <- c(0.12, 2.43, 34.62, 57.51, 45.36, 32.62, 27.34)
chisq.test(esperados,observados)##
## Pearson's Chi-squared test
##
## data: esperados and observados
## X-squared = 42, df = 36, p-value = 0.227
Dado que el valor del p-value menor al nivel de significancia 0.05 podemos rechazar la hipotésis nula y por lo tanto el modelo no es el adecuado.
Problema hipotético
Usted está trabajando en un proyecto del IIMAS que consiste en analizar el daño causado a embarcaciones por la exposición al oleaje.
Como parte del proyecto se recolectaron los siguientes datos:
type: tipo de embarcación (codificado con las letra A a E)
year: año de construcción de la embarcación
period: periodo de operación (1960-74, 75-79)
service: acumulado de meses en servicio
incidents: número de incidentes por daño
(Hint: ajuste el modelo incidents ~ offset(log(service))+type+year+period)
¿Diría que el modelo es bueno?
El uso de offset() en el modelo permite agregar variables cuyo parámetro se fija en el valor 1 en lugar de permitirle al algoritmo estimarlo. Explique por qué es buena idea usarlo para agregar la variable service al modelo.
N. B. Use los datos ships del paquete MASS
library(dplyr)
library(MASS)
#Filtramos los barcos para tener aquellos que tengan al menos un mes en servicio
datos_ships <- ships %>%
filter(service != 0)Ajuste del modelo
##
## Call:
## glm(formula = incidents ~ ., family = poisson, data = datos_ships)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.4598 -1.6558 -0.3323 0.6261 3.5104
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.387e+00 1.151e+00 -3.811 0.000138 ***
## typeB 9.622e-01 2.058e-01 4.675 2.94e-06 ***
## typeC -1.212e+00 3.274e-01 -3.701 0.000215 ***
## typeD -8.652e-01 2.875e-01 -3.009 0.002619 **
## typeE -1.105e-01 2.350e-01 -0.470 0.638160
## year 5.280e-02 1.378e-02 3.831 0.000127 ***
## period 3.642e-02 9.245e-03 3.939 8.19e-05 ***
## service 4.785e-05 7.050e-06 6.787 1.15e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 614.54 on 33 degrees of freedom
## Residual deviance: 120.56 on 26 degrees of freedom
## AIC: 234.43
##
## Number of Fisher Scoring iterations: 5
## Start: AIC=234.43
## incidents ~ type + year + period + service
##
## Df Deviance AIC
## <none> 120.56 234.43
## - year 1 135.85 247.72
## - period 1 137.28 249.15
## - service 1 167.80 279.67
## - type 4 208.01 313.88
##
## Call: glm(formula = incidents ~ type + year + period + service, family = poisson,
## data = datos_ships)
##
## Coefficients:
## (Intercept) typeB typeC typeD typeE
## -4.387e+00 9.622e-01 -1.212e+00 -8.652e-01 -1.105e-01
## year period service
## 5.280e-02 3.642e-02 4.785e-05
##
## Degrees of Freedom: 33 Total (i.e. Null); 26 Residual
## Null Deviance: 614.5
## Residual Deviance: 120.6 AIC: 234.4
Ahora con el modelo HINT
fit_hint <- glm(incidents ~ offset(log(service))+type+year+period, data = datos_ships, family = poisson)
summary(fit_hint)##
## Call:
## glm(formula = incidents ~ offset(log(service)) + type + year +
## period, family = poisson, data = datos_ships)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5348 -0.9319 -0.3686 0.4654 2.8833
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.079076 0.876149 -11.504 < 2e-16 ***
## typeB -0.546090 0.178415 -3.061 0.002208 **
## typeC -0.632631 0.329500 -1.920 0.054862 .
## typeD -0.232257 0.287979 -0.807 0.419951
## typeE 0.405975 0.234933 1.728 0.083981 .
## year 0.042247 0.012826 3.294 0.000988 ***
## period 0.023705 0.008091 2.930 0.003392 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 146.328 on 33 degrees of freedom
## Residual deviance: 59.375 on 27 degrees of freedom
## AIC: 171.24
##
## Number of Fisher Scoring iterations: 5
Podemos notar que el modelo HINT es mejor que el anterior seleccionado por usando el criterio de Akaike, el modelo HINT tiene un valor menor, aunque la opción offset() fija en 1 el valor de la variable service, no me queda muy claro porque tiene también la función log()
Me parece buen modelo o al menos el mejor de todos los propuestos, y fijar la variable me pareció buena idea porque al fijar en 1 estamos de cierta forma haciendo homogeneo el número de meses (o al menos eso pienso) para ver que tanto afectan los demás factores