Modelos Estadísticos. Grado Biotecnología

Introducción

Son frecuentes los datos discretos que provienen de conteos de sucesos que se producen por azar con cierta frecuencia y son modelizables en términos de tasas de incidencia que dependen de ciertas variables explicativas. Para modelizar este tipo de datos se utiliza la distribución de Poisson, \(X \sim Po(\mu)\), donde \(\mu\) representa el número medio de ocurrencias, de forma que: \[E(X) = \mu \text{ y } V(X) = \mu.\]

El parámetro \(\mu\) requiere una definición cuidadosa. A menudo es necesario describirlo como una tasa; por ejemplo, el número promedio de clientes que compran un particular producto de cada 100 clientes que ingresan a la tienda. De manera más general, la tasa se especifica en términos de unidades de “exposición” o en términos de años-persona ‘en riesgo’, por ejemplo, el número de personas que sufren cierta enfermedad sobre el total de personas para un instante de tiempo determinado. El efecto de las variables predictoras sobre la respuesta se modeliza valorando su efecto sobre \(\mu\).

Sea \(Y_1, Y_2,..., Y_n\) un conjunto de variables aleatorias Poisson, donde \(Y_i\) representa el número de eventos observados a partir de la exposición \(n_i\) para el i-esimo patrón de covariables. En esta situación el valor esperado de \(Y_i\) se puede escribir como: \[E(Y_i) = \mu_i = n_i \lambda_i,\] donde \(\lambda_i\) representa la tasa de incidencia para el patrón i-ésimo.

Veamos dos situaciones prácticas:

Telas: Se recogen los resultados de un experimento para determinar el efecto del tipo de lana (A o B) y la tensión (baja, media o alta) en el número de roturas de deformación en la fabricación de telas. Se recopilaron datos de nueve telares para cada combinación de configuraciones. Las variables consideradas son breaks (número de roturas), wool (tipo de lana), y tension (tensión de la lana). La pregunta de interés es si ¿resulta posible predecir el número de roturas en función del tipo de lana y de la tensión?

En este caso la exposición a riesgo es constante (\(n_i =1\)), y por tanto se considera que las posibles predictoras influyen directamente sobre \(\mu_i = breaks_i\)

Barcos: En este banco de datos se analizan lo datos proporcionados por la aseguradora Lloyd·s con respecto al número de incidentes causados por la olas en cierto componente de los buques cargueros. El banco de datos contiene información sobre: type (el tipo de barco, codificado de A hasta E), year (el año de construcción del barco en los periodos 1960–64, 65–69, 70–74, 75–79 codificados como 60, 65, 70, 75), period (período de operación del barco 1960-74, 75-79 codificados como 60, 75), service (meses de servicio), incidents (número de incidentes)

En este caso la exposición a riesgo no es constante, ya que caunto mayor sea el tiempo de servicio de un barco parece más probable que sufra un incidente. En este caso \(n_i = service_i\) y modelizamos la tasa de riesgo \(\lambda_i\) con la información de las variables predictoras. En la práctica \(\lambda_i = incidents_i/service_i\) y represetamos dicha variable con respecto a las predictoras para el planteamiento del modelo.

GLM Poisson

El modelo lineal generalizado para variables aleatorias Poisson o más conocido como regresión de Poisson se puede especificar siguiendo la definición establecida en el punto anterior como:

La distribución de cada variable aleatoria \(Y_i\) viene dada por \[Y_i \sim Po(\mu_i)\]
El valor medio de cada variable viene dado por: \[E(Y_i) = \mu_i = n_i\lambda_i\]
Si consideramos un conjunto \(x_{i1}, x_{i2},...,x_{ip}\) de \(p\) variables predictoras que pueden afectar el comportamiento de \(Y_i\) asumimos que: \[\lambda_i = exp(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + ...+ \beta_p x_{ip}),\] donde cada \(\beta\) representa el efecto de la correspondiente variable predictora.

En esta situación y tomando como función de enlace el logaritmo tenemos que: \[log(E(Y_i)) = log(\mu_i) = log(n_i\lambda_i) = log(n_i) + log(\lambda_i)= log(n_i) + \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + ...+ \beta_p x_{ip}\]

donde tenemos la expresión del modelo lineal generalizado. La única diferencia de este modelo con los vistos hasta ahora es la inclusión de un término fijo en el modelo (\(log(n_i)\)) que actúa como “offset” y hace referencia al logaritmo de sujetos en riesgo para la combinación de predictoras correspondiente al índice \(i\). Este términos se debe incorporar en la estimación del modelo como una constante para cada \(i\). Para el ejemplo de las telas cada \(n_i = 1\) y cada \(log(n_i) = 0\) desapareciendo de la estimación del modelo. En el caso de los barcos tenemos que \(log(n_i) = log(service_i)\) remarcando el hecho de que el riesgo aumenta cuando lo hace el tiempo de servicio. El offset no se puede estimar ya que es fijo, y la estimación del modelo se centra únicamente en los parámetros asociados con cada una de las predictoras.

Respecto de la interpretación de dichos coeficientes se puede demostrar que \(exp(\beta_j)\) representa el riesgo relativo (\(RR\)) sobre la tasa de incidencia de los sucesos asociado a un incremento de una unidad en la covariable \(x_j\). Es similar a la interpretación que hacíamos en los modelos lineales, salvo el cambio de escala debido a la transformación logaritmo sobre la media de la respuesta.

Para una variable de tipo numérico se define el riesgo relativo al incrementar una unidad el valor de la variable predictora \(X\) como: \[RR(X) = \frac{\lambda(X = x+1)}{\lambda(X = x+1)} = exp(\beta_x)\] mientras que riesgo relativo para valorar el cambio entre dos niveles de una variable categórica (A y B) se define como:

\[RR(X_{AB}) = \frac{\lambda(X = B)}{\lambda(X = A)} = exp(\text{coeficiente asociado a B} - \text{coeficiente asociado a A})\]

El riesgo relativo aumenta cuando \(RR > 1\) y disminuye cuando \(RR < 1\). Tenemos el mismo riesgo cuando \(RR = 1\)

Análisis preliminar

En este punto realizamos un análisis preliminar de los bancos de datos de los ejemplos para estudiar el comportamiento de la tasa de incidencia con respecto a cada predictora y determinar los posible modelos a considerar en la fase de estimación. En primer lugar cargamos las librerías de trabajo

Telas

Comenzamos con la lectura del banco de datos y realizamos un estudio descriptivo inicial

require(datasets)
glm_poi_01 <- warpbreaks
# Descriptivo
summary(glm_poi_01)

##      breaks      wool   tension
##  Min.   :10.00   A:27   L:18   
##  1st Qu.:18.25   B:27   M:18   
##  Median :26.00          H:18   
##  Mean   :28.15                 
##  3rd Qu.:34.00                 
##  Max.   :70.00

En este tipo de modelos es muy habitual trabajar con variables categóricas de tipo ordinal. Dichas varaible spueden ser incorporadas en el modelo como factores pero en ocasiones interesa introducirlas como numéricas para captar las posibles tendencias de la respuesta conforme aumneta la relevancia del factor. Recordemos que una variable categórica ordinal se puede obtener siempre a partir de una variable de tipo numérico. Para poder hacer esto introducimos un código numérico artificial asociado con cada categoría. Disponemos de dos alternativas:

1. Asignar un código continuo: 1, 2, 3,… donde los valores más bajos se asocian con los niveles más bajos de la variable categórica.

2. Cuando la variable categórica viene dada en términos de un intervalo se puee usar el punto medio de dicho intervalo para establecer el código numerico

En este caso optamos por asignar un código continuo:

# Creamos varuables par acada categoria con su código numérico. Construimos la variable conjunta y 
# nos quedamos sólo con las variables de interés
glm_poi_01 <- glm_poi_01 %>% 
  mutate(t1 = ifelse(tension == "L",1,0)) %>%
  mutate(t2 = ifelse(tension == "M",2,0)) %>%
  mutate(t3 = ifelse(tension == "H",3,0)) %>%
  mutate(tension.num = t1+t2+t3) 
# Selecionamos el conjunto de variables final
glm_poi_01 <- glm_poi_01[,c("breaks","wool","tension","tension.num")]

Representamos la incidencia (en este caso \(n_i = 1\)) que corresponde con la variable breaks, es decir, representamos el logaritmo de la incidencia con respecto a las predictoras, para valorar la linealidad del efecto.

# Utilizamos las dos varible categóricas
ggplot(glm_poi_01,aes(x = tension, y = log(breaks), color = wool)) + geom_boxplot() + theme_bw()

En el gráfico se aprecia cierta tendencia con la tensión, pero no tanto con el tipo de lana. Realizamos el gráfico de interacción con las medias para las combinaciones de ambos factores.

ggplot(glm_poi_01, aes(x = tension, y = log(breaks), group = wool, color = wool)) +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.y = mean, geom = "line") +
  theme_bw()

Se puede ver un efecto conjunto entre tension y wool. Parece bastante claro que deberemos considerar un efecto de interacción para este modelo. De hecho, se aprecian curvas descendientes con tension (indicando que las telas sufren menos roturas a mayor tensión) no apreciándose efecto lineal lo que puede provocar que no se verifiquen las hipótesis del modelo, y que sea necesario introducir ciertos tendencias en las variables predictoras. En el tipo de lana resulta imposible hacer esto porque es una variable categórica nominal, pero si podríamos hacerlo con tension si la consideramos como numérica a partir de los códigos introducidos para cada categoría. Ajustaremos ambos modelos y decidiremos cual de ellos produce resultados más interesantes.

Barcos

Comenzamos con la lectura del banco de datos y realizamos un estudio descriptivo inicial

require(MASS)
glm_poi_07 <- ships
# generamos factores a partir de los códigos numéricos
glm_poi_07$year.f <- as.factor(glm_poi_07$year)
glm_poi_07$period.f <- as.factor(glm_poi_07$period)
summary(glm_poi_07)

##  type       year           period        service          incidents   
##  A:8   Min.   :60.00   Min.   :60.0   Min.   :    0.0   Min.   : 0.0  
##  B:8   1st Qu.:63.75   1st Qu.:60.0   1st Qu.:  175.8   1st Qu.: 0.0  
##  C:8   Median :67.50   Median :67.5   Median :  782.0   Median : 2.0  
##  D:8   Mean   :67.50   Mean   :67.5   Mean   : 4089.3   Mean   : 8.9  
##  E:8   3rd Qu.:71.25   3rd Qu.:75.0   3rd Qu.: 2078.5   3rd Qu.:11.0  
##        Max.   :75.00   Max.   :75.0   Max.   :44882.0   Max.   :58.0  
##  year.f  period.f
##  60:10   60:20   
##  65:10   75:20   
##  70:10           
##  75:10           
##                  
##

Como se puede ver en el resumen hay barcos que tiene 0 meses de servicio, y por tanto deben ser eliminados del análisis, ya que al tratar de ajustar la tasa dividiríamos por un valor de 0, lo que provocaría la imposibilidad de un ajuste adecuado. Seleccionamos todos los barcos con tiempo de servicio mayor que cero y obtenemos la tasa de incidencia a partir de la relación entre el número de incidentes y las horas de servicio de cada barco.

glm_poi_07 <- glm_poi_07 %>% 
  filter(service>0) %>%
  mutate(tasa = incidents/service)

Realizamos el gráfico de la media del logaritmo de la tasa para type, year y period. En realidad trabajamos con los factores para representar adecuadamente el gráfico de interacción.

ggplot(glm_poi_07, aes(x = year.f, y = log(tasa), group = period.f, color = period.f)) +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.y = mean, geom = "line") +
  facet_grid(~type) +
  theme_bw()

Se observa que los tipos de barco se comportan de forma diferente en cuanto a la tasa de incidentes. Los barcos del tipo “E” son los que tienen tasas más altas. En cuanto al año de fabricación y el periodo de funcionamiento varían con el tipo de barco. Mientras que para los tipos “A” y “B” se observan comportamientos muy similares no ocurre los mismo para el tipo “C”. De nuevo habrá que vigilar el efecto para asegurar la linealidad del GLM.

En este caso nos podemos plantear un modelo con tres predictoras de tipo factor y considerar sus posible interacciones (lo que nos da un modelo bastante complejo), u optar por un modelo más sencillo donde consideramos year y periodcomo numéricas. Obtenemos la tabla cruzada de ambos factores para ver si hay alguna combinación que no se da, y por tanto desechar el modelo a que darían lugar.

table(glm_poi_07$year.f,glm_poi_07$period.f)

##     
##      60 75
##   60  5  4
##   65  5  5
##   70  5  5
##   75  0  5

Se observa que la combinación year_f 75 y period_f 60 no tienen ningún dato. Si consideramos la interacción de ambos factores en el modelo obtendremos un error ya que no disponemos de datos para esa combinación. Procedemos por tanto a considerarlas como numéricas para evitar esos problemas con la estimación del modelo.

Estimación

Como hemos hecho hasta ahora, en un primer momento estimamos el modelo más complejo posible para proceder posteriormente con la selección de aquellos efectos que de verdad resultan significativos. Para la estimación utilizamos la función glm con las especificaciones siguientes:

glm(modelo,family = poisson, offset = value, data_set)

La diferencia principal con los modelos de regresión logística es la introducción de los offsets en el ajuste del modelo. Si no disponemos de offsets (como en el ejemplo de la lana) no aparecerán en el ajuste del modelo.

Telas

Para este banco de datos ajustamos dos modelos completos considerando la tensión como factor o como variable numérica. Comenzamos con el modelo considerando todas las predictoras como factores. Consideramos tanto los efectos asociados a cada factor como la posible interacción entre ellos. En este caso no hay offsets. El modelo propuesto es: \[breaks \sim wool*tension\]

telas.cat <-glm(breaks ~ wool*tension,family = poisson, glm_poi_01)
# Resumen del modelo
summary(telas.cat)

Los coeficientes del modelo resultan significativos salvo en al combinación wool B y tension H. Ajustamos ahora el otro modelo asumiendo la variable tensión con el código numérico: \[breaks \sim wool*tension.num\]

telas.num <-glm(breaks ~ wool*tension.num,family = poisson, glm_poi_01)
# Resumen del modelo
summary(telas.num)

Los coeficientes del modelo resultan significativos pero el valor del AIC es superior al del modelo con los dos factores. Por el momento, nos quedamos con el modelo telas.cat y escribimos las expresiones de dicho modelo:

\[\left\{\begin{array}{lll} log(breaks)_{woolA;tensionL} & = 3.79674 & \to exp(3.79674) = 44.6\\ log(breaks)_{woolA;tensionM} & = 3.79674 - 0.61868 = 3.17806 & \to exp(3.17806) = 24.0\\ log(breaks)_{woolA;tensionH} & = 3.79674 - 0.59580 = 3.20094 & \to exp(3.20094) = 24.6\\ log(breaks)_{woolB;tensionL} & = 3.79674 - 0.45663 = 3.34011 & \to exp(3.34011) = 28.2\\ log(breaks)_{woolB;tensionM} & = 3.79674 - 0.45663 - 0.61868 + 0.63818 = 3.35961 & \to exp(3.35961) = 28.8\\ log(breaks)_{woolB;tensionH} & = 3.79674 - 0.45663 - 0.59580 + 0.18836 = 2.93267 & \to exp(2.93267) = 18.8\\ \end{array}\right.\]

El proceso de estimación nos proporciona que el mayor número de roturas se produce para la combinación wool = A y tensión = L, mientras que el menor número se estima para la combinación wool = B y tensión = H.

Barcos

Ajustamos el modelo utilizando las variables con los códigos numéricos en lugar de los factores. En este caso si disponemos de offset y el modelo planteado viene dado por: \[incidents \sim offset(service) + type *(year + period)\]

barcos.num <-glm(incidents ~ type *(year + period),family = poisson, offset = log(service), glm_poi_07)
# Resumen del modelo
summary(barcos.num)

El modelo presenta muchos coeficientes que parecen no significativos, y por tanto dejamos la interpretación hasta la selección del modelo definitivo.

Selección del modelo

Para la construcción del mejor modelo utilizaremos los mismos procedimientos que presentamos en el tema de regresión logística. Estudiamos a continuación cada uno de los ejemplos.

Telas

Utilizamos el modelo con los dos factores. Dado que existe un efecto de interacción, el proceso de selección debe establecer si es posible prescindir de dicho efecto, o si por otro debemos considerarlo en el modelo.

telas.cat <-glm(breaks ~ wool*tension,family = poisson, glm_poi_01)
# Resumen del modelo
step(telas.cat, test = "Chisq")
telas.fit <- telas.cat

Se puede ver que tanto el AIC con la significatividad el test indican que el efecto de interacción no puede ser eliminado del modelo. La conclusión es que para predecir el número de roturas debemos tener en cuenta las posibles combinaciones que se producen al mezclar los tipos de lana y su tensión (comportamientos distintos en cada combinación). El resumen del modelo es el obtenido en el punto anterior.

Sin embargo, con las ecuaciones estimadas anteriormente podemos construir las comparaciones de riesgo de rotura entre los tipos de tela o entre las diferentes tensiones que nos proporcionan una medida más realista del riesgo relativo (\(RR\)) de cada combinación. Por ejemplo:

El riesgo relativo asociado con los tipos de tela A y B para una tensión fija (L) es de 0.54 lo que implica un riesgo relativo de rotura para la tela B de casi la mitad que para la tela A para una tensión L. El cálculo del \(RR\) viene dado por:

\[RR_{AB;tensionL} = \frac{\lambda_{woolB;tensionL}}{\lambda_{woolA;tensionL}} = exp(3.17806 - 3.79674) = 0.54\]

El riesgo relativo asociado con las tensiones L y H para e tipo de tela A es de 0.55 lo que implica un riesgo relativo de rotura para la tensión H de casi la mitad que para la tensión L para el tipo de tela A. El cálculo del \(RR\) viene dado por:

\[RR_{LHC;woolA} = \frac{\lambda_{woolA;tensionH}}{\lambda_{woolA;tensionL}} = exp(3.20094 - 3.79674) = 0.55\]

El riesgo relativo asociado con las tensiones L y H para e tipo de tela B es de 0.66 lo que implica un riesgo relativo de rotura para la tensión H un 33% inferior que para la tensión L con el tipo de tela B. El cálculo del \(RR\) viene dado por:

\[RR_{LHC;woolB} = \frac{\lambda_{woolB;tensionH}}{\lambda_{woolB;tensionL}} = exp(2.93267 - 3.34011) = 0.66\]

En este caso los riesgos relativos estimados se pueden obtener a partir de la predicción del predictor lineal una vez hemos desecho la transformación logaritmo.

Barcos

Tomamos como modelo de partida el que considera las variables de clasificación con sus códigos numéricos (las tratamos como numéricas). De nuevo debemos valorar la presencia o no de las interacciones entre el tipo de barco y las variables codificadas.

barcos.num <-glm(incidents ~ type *(year + period),family = poisson, offset = log(service), glm_poi_07)
step(barcos.num, test = "Chisq")

## Start:  AIC=167
## incidents ~ type * (year + period)
## 
##               Df Deviance    AIC     LRT Pr(>Chi)   
## - type:period  4   45.138 165.00  6.0015 0.199039   
## <none>             39.136 167.00                    
## - type:year    4   54.046 173.91 14.9095 0.004893 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step:  AIC=165
## incidents ~ type + year + period + type:year
## 
##             Df Deviance    AIC     LRT Pr(>Chi)   
## <none>           45.138 165.00                    
## - type:year  4   59.375 171.24 14.2369 0.006576 **
## - period     1   54.589 172.46  9.4517 0.002110 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Call:  glm(formula = incidents ~ type + year + period + type:year, family = poisson, 
##     data = glm_poi_07, offset = log(service))
## 
## Coefficients:
## (Intercept)        typeB        typeC        typeD        typeE  
##    -9.40739     -1.99170     -4.80174      2.57742     13.61537  
##        year       period   typeB:year   typeC:year   typeD:year  
##     0.03166      0.02474      0.02134      0.06103     -0.03907  
##  typeE:year  
##    -0.19193  
## 
## Degrees of Freedom: 33 Total (i.e. Null);  23 Residual
## Null Deviance:       146.3 
## Residual Deviance: 45.14     AIC: 165

barcos.fit <- glm(incidents ~ type + year +period + type:year,family = poisson, offset = log(service), glm_poi_07)

El proceso de selección de efectos indica que la interacción type:period no es relevante para el modelo considerado ya que su eliminación disminuye la deviance y además el test resulta no significativo (pvalor = 0.1990). En principio es el único efecto que puede ser eliminado del modelo. El resumen del modelo final obtenido es:

summary(barcos.fit)

## 
## Call:
## glm(formula = incidents ~ type + year + period + type:year, family = poisson, 
##     data = glm_poi_07, offset = log(service))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0704  -1.1705  -0.2069   0.4373   2.7276  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -9.407387   3.016649  -3.118  0.00182 **
## typeB       -1.991695   3.124006  -0.638  0.52377   
## typeC       -4.801741   5.647285  -0.850  0.39517   
## typeD        2.577417   4.878123   0.528  0.59725   
## typeE       13.615368   4.911409   2.772  0.00557 **
## year         0.031664   0.042913   0.738  0.46059   
## period       0.024744   0.008074   3.065  0.00218 **
## typeB:year   0.021342   0.044519   0.479  0.63166   
## typeC:year   0.061030   0.081800   0.746  0.45561   
## typeD:year  -0.039075   0.068676  -0.569  0.56938   
## typeE:year  -0.191928   0.071220  -2.695  0.00704 **
## ---
## 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:  45.138  on 23  degrees of freedom
## AIC: 165
## 
## Number of Fisher Scoring iterations: 6

Las conclusiones de este modelo son:

Con respecto al tipo de barco (sumando los coeficientes de interceptación y tipo de barco) podemos concluir que el riesgo de sufrir un incidente aumenta desde el tipo A (menor riesgo) al E (mayor riesgo).
Con respecto al año de fabricación (sumando los coeficientes de año y la interacción con tipo) podemos concluir que el riesgo aumenta con el año de fabricación, pero dicho efecto se atenúa para los distintos tipos de barco. Los barcos del tipo E disminuyen el riesgo cuando el año de fabricación es más alto. Por otro lado los barcos de tipo C muestran una mayor riesgo cuando aumenta el año de fabricación.
Con respecto al periodo de operación podemos ver que el riesgo aumenta cuando aumenta el periodo de operación.

Veamos las ecuaciones estimadas (reducimos los decimales para una mejor visualización) para cada tipo de barco \[\left\{\begin{array}{ll} log(\lambda_{typeA}) & = -9.41 + 0.03*year + 0.02*period\\ log(\lambda_{typeB}) & = (-9.41 - 1.99) + (0.03 + 0.02)*year + 0.02*period\\ log(\lambda_{typeC}) & = (-9.41 - 4.80) + (0.03 + 0.06)*year + 0.02*period\\ log(\lambda_{typeD}) & = (-9.41 + 2.58) + (0.03 - 0.04)*year + 0.02*period\\ log(\lambda_{typeE}) & = (-9.41 + 13.62) + (0.03 - 0.19)*year + 0.02*period\\ \end{array}\right.\]

Comparamos ahora los riesgos de diferentes situaciones:

El riesgo relativo de los tipos de barco A y E para un mismo año y período de operación se puede obtener a partir de las ecuaciones anteriores como

\[RR_{AE} = \frac{\lambda_{typeE}}{\lambda_{typeA}} = exp(\{4.21 - 0.16*year +0.02*period\} - \{-9.41 + 0.03*year + 0.02*period\})\] \[\text{de forma que } RR_{AE} = exp(13.62 - 0.19*year)\] Vemos el comportamiento del riesgo representándolo gráficamente con respecto al año de fabricación:

# Representamos dicho riesgo relativo
anyo <- seq(60,75,0.5)
riesgo <- exp(13.62 - 0.19*anyo)
datos <- data.frame(anyo,riesgo)
ggplot(datos,aes(x = anyo, y = riesgo)) + 
  geom_line() + 
  labs(x = "Año fabricación", y = "Riesgo Relativo", title ="Tipo E versus Tipo A") + 
  theme_bw()

El riesgo relativo de tener un incidente disminuye conforme aumenta el año de fabricación, pero para los instantes iniciales (año 60) el riesgo de los barcos de tipo E es más de siete veces superior al de los barcos de tipo A. EL hecho más destacable es que el periodo no tiene efecto sobre este riesgo debido a que la interacción entre tipo y periodo no está presente en el modelo.

El riesgo se aumentar en un año el periodo de operación para un mismo tipo de barco (A) y el mismo año de fabricación se puede obtener como:

\[RR_{period} = \frac{\lambda_{period+1,typeA}}{\lambda_{period,typeA}} = exp(\{-9.41 + 0.03*year + 0.02*(period+1)\} - \{-9.41 + 0.03*year + 0.02*period\})\] \[\text{de forma que } RR_{period} = exp(0.02) = 1.02\] El riesgo se mantiene constante cuando aumentamos en un año el periodo de operaciones. El riesgo asociado con aumentar dos años vendría dado por \(exp(2*0.02)\). Veamos el gráfico:

# Representamos dicho riesgo relativo
anyo <- seq(1,15,0.5) #incrementos de periodo desde 60 a 75
riesgo <- exp(anyo*0.024744)
datos <- data.frame(anyo,riesgo)
ggplot(datos,aes(x = anyo, y = riesgo)) + 
  geom_line() + 
  labs(x = "Periodo", y = "Riesgo Relativo", title ="Period para tipo A") + 
  theme_bw()

El riesgo aumenta cuando aumenta el periodo de servicio. Cuando este periodo se sitúa en una diferencia de 12 años el riesgo de sufrir una incidente es 1.35, es decir, un 33% más de riesgo en el barco más antiguo que en el más moderno.

Diagnóstico

Las herramientas de diagnóstico para este tipo d modelos son las mismas que el caso de los GLM con respuesta binaria. Los gráficos que vamos a utilizar son:

Residuos vs Ajustados (predictor lineal)
Residuos vs variables en el modelo (predictor lineal)
Valores influyentes

Telas

Usamos el modelo telas.fit para el diagnóstico. Obtenemos todas las cantidades de interés del modelo. Para acelerar el proceso de diagnóstico utilizamos el comando plot(modelo) que nos proporciona el gráfico de residuos versus ajustados, el gráfico qq de los residuos, el gráfico de linealidad (raíz cuadrada de la deviance residual vs valores predichos), y residuos vs variables predictoras. Obtenemos los estadísticos asociados con el modelo para valorar la existencia de valores influyentes.

telas.fit <-glm(breaks ~ wool*tension,family = poisson, glm_poi_01)
# Obtención de valores para el diagnóstico
fortify(telas.fit)

# Gráficos
par(mfrow=c(2,2))
plot(telas.fit)

No se aprecian ningún comportamiento digno de interés que pueda hacer dudar del modelo obtenido. Dado que tenemos dos factores sólo estimamos las medias correspondientes a las posibles combinaciones de ellos dos.

Barcos

Utilizamos el modelo ajustado en la selección del modelo

barcos.fit <- glm(incidents ~ type + year +period + type:year,family = poisson, offset = log(service), glm_poi_07)# Obtención de valores para el diagnóstico
fortify(barcos.fit)

# Gráficos
par(mfrow=c(2,2))
plot(barcos.fit)

Se aprecia cierta tendencia en los residuos versus predichos que debería analizarse en más detalle. Realizamos el gráfico de las residuos con respecto a las variables numéricas para valorar la posible inclusión de otro tipo de efectos.

diagnostico <-fortify(barcos.fit)

# Residuos vs year
ggplot(diagnostico) +
geom_jitter(aes(year,.stdresid, color = type),) +
  geom_smooth(aes(year,.stdresid, color = type), method=loess, se=FALSE) +
  theme_bw()

# Residuos vs period
ggplot(diagnostico) +
geom_jitter(aes(period,.stdresid, color = type),) +
  geom_smooth(aes(period,.stdresid, color = type), method=loess, se=FALSE) +
  theme_bw()

No se aprecian tendencias destacables para period pero si para year. Proponemos un modelo de suavizado con respecto a esta variable y ajustamos un nuevo modelo. Para el ajuste de ese modelo utilizamos la función gam(). Dado que tenemos pocos datos para tipo de barco, rebajamos el valor de \(k\) a 5 en la función de suavizado.

# Modelo con interacción type:period
M1 <- gam(incidents ~ type + period + 
                    s(year, k = 5, m = 2, bs = "ps", by = type) +
                    type:period
                  ,family = poisson, offset = log(service), glm_poi_07)
# Modelo sin interacción type:period
M2 <- gam(incidents ~ type + period + 
                    s(year, k = 5, m = 2, bs = "ps", by = type) 
                  ,family = poisson, offset = log(service), glm_poi_07)
AIC(M1,M2)

El modelo sin interacción de type:period tiene un menor valor de AIC y lo utilizamos como modelo final.

barcos.fit <- gam(incidents ~ type + period + 
                    s(year, k = 5, m = 2, bs = "ps", by = type) 
                  ,family = poisson, offset = log(service), glm_poi_07)
# Resumen del modelo
summary(barcos.fit)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## incidents ~ type + period + s(year, k = 5, m = 2, bs = "ps", 
##     by = type)
## 
## Parametric coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -7.393678   0.606365 -12.193  < 2e-16 ***
## typeB       -0.553323   0.255014  -2.170  0.03002 *  
## typeC       -0.696173   0.397851  -1.750  0.08015 .  
## typeD       -5.246512   7.967352  -0.659  0.51022    
## typeE        0.870856   0.313020   2.782  0.00540 ** 
## period       0.025545   0.007951   3.213  0.00131 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                 edf Ref.df Chi.sq  p-value    
## s(year):typeA 1.166  1.306  0.771  0.52117    
## s(year):typeB 1.973  2.226 22.167 2.43e-05 ***
## s(year):typeC 1.000  1.000  1.737  0.18758    
## s(year):typeD 1.871  1.983  7.753  0.02457 *  
## s(year):typeE 1.000  1.000  7.847  0.00509 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.971   Deviance explained = 86.7%
## UBRE = 0.33848  Scale est. = 1         n = 34

El modelo obtenido tiene buena capacidad explicativa. En cuanto a la interpretación de los coeficientes podemos ver que el riesgo de incidentes aumenta con el periodo y es mayor para los barcos de tipo E que para el resto. El tipo de barco con menor riesgo es el D. En cuanto a las funciones de suavizado, representamos sus efectos para observar el comportamiento del año para cada tipo de barco:

par(mfrow = c(2,3))
plot.gam(barcos.fit)

Se observan efectos casi lineales para todos los tipos salvo para los de tipo D. En este caso el riesgo de incidentes aumenta cuando lo hace año.

En cuanto al diagnóstico del modelo no se encuentra ningún comportamiento relevante (pvalores mayores que 0.05), salvo la existencia de un par de observaciones cuya predicción del predictor lineal está muy alejada del resto. Se podría intentar eliminar esas dos observaciones y ajustar un nuevo modelo, pero dado que el número de observaciones es bajo eso podría ocasionar que el modelo no se podría ajustar.

gam.check(barcos.fit)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 10 iterations.
## Gradient range [-2.348054e-07,1.792604e-12]
## (score 0.3384798 & scale 1).
## Hessian positive definite, eigenvalue range [1.827615e-07,0.02021138].
## Model rank =  26 / 26 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                 k'  edf k-index p-value
## s(year):typeA 4.00 1.17    1.17    0.82
## s(year):typeB 4.00 1.97    1.17    0.85
## s(year):typeC 4.00 1.00    1.17    0.85
## s(year):typeD 4.00 1.87    1.17    0.80
## s(year):typeE 4.00 1.00    1.17    0.81

El problema con las funciones de suavizado es que resulta más complicado la evaluación de los riesgos relativos. Al no disponer de una forma específica para las predictoras que forman parte del suavizado, no podemos obtener directamente dicho riesgo.

Para evaluar el riesgo del tipo de barco A con respecto al tipo E para un mismo valor de año y periodo utilizamos el código siguiente

# Para predecir
n <-20
year <- c(seq(60,75,length = n),seq(60,75,length = n))
period <- c(seq(60,75,length = n),seq(60,75,length = n))
type <- c(rep("A",n),rep("E",n))
newdata <- data.frame(year,period,type)
# Predicción
prediccion <- predict.gam(barcos.fit,newdata,type = "response",se.fit = FALSE)  
# Riesgo
RR <- prediccion[(n+1):(2*n)]/prediccion[1:n]
datos <- data.frame(anyo = year[1:n], riesgo = RR)
ggplot(datos,aes(x = anyo, y = riesgo)) + 
  geom_line() + 
  labs(x = "Año", y = "Riesgo", title ="Tipo E vs tipo A") + 
  theme_bw()

El riesgo disminuye cuando aumenta el año de fabricación.

Predicción

Vamos a construir la predicción para cada uno de los ejemplos que hemos venido trabajando a lo largo del tema.

Telas

Puesto que las dos predictoras son de tipo categórico solo se muestra la tabla de predicción con intervalo de confianza aproximado del 95%.

telas.fit <-glm(breaks ~ wool*tension,family = poisson, glm_poi_01)
grid <- glm_poi_01 %>% 
  data_grid(wool,tension) 
# Obtenemos la predicción para el modelo ajusatdo
prediccion <- data.frame(grid, predict(telas.fit, grid, type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza 
prediccion <- prediccion %>% mutate(
  lwr = fit - (1.96 * se.fit),
  upr = fit + (1.96 * se.fit)
)
prediccion

El gráfico con las predicciones viene dado por:

# Creamos código de las combinaciones
newdata <- prediccion %>% mutate(grupo = paste(wool,tension)) 
# Representamos ordenando por varaible de interés
ggplot(newdata, aes(x = fct_reorder(grupo,fit), y= fit, color = wool)) + 
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1) +
  geom_line() +
  geom_point() +
  labs(x = "Wool - Tension", y = "Breaks") +
  theme_bw()

Barcos

Antes de proceder con la predicción para este ejemplo es necesario realizar la siguiente anotación. La inclusión de offset altera en cierta forma el proceso de predicción. En el caso de ajustar un glm clásico la función predict(modelo.glm, type="response") nos proporciona la predicción de la respuesta para ese modelo, pero no ocurre así cuando trabajamos con un gam. En este caso la función predict(modelo.gam, type="response") nos proporciona la predicción de la tasa de incidencia \(\lambda\) y no de la respuesta, por lo que es necesario multiplicar por el valor del offset para obtener la predicción de la respuesta. Si queremos predecir un valor en particular tendremos que asumir un valor para el offset con le objeto de calcular la predicción final de la respuesta.

En este caso tenemos un offset en la variable service. Calculamos las predicciones para todo el banco de datos y verificamos cuanto se parecen a los valores reales observados.

# Modelo
barcos.fit <- gam(incidents ~ type + period + 
                    s(year, k = 5, m = 2, bs = "ps", by = type) 
                  ,family = poisson, offset = log(service), glm_poi_07)
# Predicción sin offset
prediccion <- glm_poi_07 %>% mutate(fit = predict(barcos.fit,type = "response", se.fit=FALSE))
#prediccion
# Gráfico comparativo tasa real con tasa de predicción
ggplot(prediccion,aes(x = tasa,y = fit,color = type)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  facet_grid(period ~ year) +
  labs(x = "Tasas observadas" , y = "Tasas predichas", title = "Tasas en función de year (columnas) y period (filas)") +
  theme_bw()

La predicción será buena cuando los datos se distribuyan a lo largo de la diagonal. Para el año de construcción 60 podemos ver que los datos no se colocan a lo largo de la diagonal. Esto es de esperar, dado que un supuesto natural es que cuando el barco es nuevo el riesgo de sufrir incidentes será bajo dado que sus horas de servicio no serán muy altas. Cuando relacionamos el año de fabricación y el de operación ya podemos ver que la predicción de la tasa obtenida se parece bastante a los datos reales observados, indicando un buen ajuste.

GLM para respuesta Poisson

Introducción

GLM Poisson

Análisis preliminar

Telas

Barcos

Estimación

Telas

Barcos

Selección del modelo

Telas

Barcos

Diagnóstico

Telas

Barcos

Predicción

Telas

Barcos

Bibliografía