Modelos Estadísticos. Grado Biotecnología

Introducción

Consideramos en este tema la modelización de variables respuesta de tipo binario, es decir, para cada individuo, la variable respuesta para cada uno de los sujetos puede tomar únicamente dos posibles valores, que denotaremos por 0 (fracaso) y 1 (éxito),

\[Pr(Y_i = 0) = 1 - \pi_i; \text{ } Pr(Y_i = 1) = \pi_i\]

Si además se han observado una serie de covariables continuas o de tipo factor, el objetivo del análisis con el modelo lineal generalizado será predecir la probabilidad asociada al éxito (o equivalentemente al fracaso), en función de dichas covariables.

Sin embargo, las situaciones experimentales nos pueden llevar a registrar esta información de dos formas diferentes:

  • Información de cada sujeto de la muestra con una varaible que indica el éxito o el fracaso. Ejemplo: En este banco de datos aparecen los datos de 81 niñoos que fueron intervenidos quirúrgicamente para corregirles problemas en la columna vertebral. La variable binaria “kifosis” indica la presencia o ausencia de una deformidad postoperatoria en la columna, denominada kifosis. Las otras tres variables son ’Age’, edad del niñoo en meses, ’Number’, número de vértebras intervenidas en la operación, e ’Start’, que define la primera vértebra involucrada en la operación. Es de interés en el análisis investigar cómo están relacionadas variables como la edad del niño, el número de vértebras dañaadas y el número de vértebra en la que empieza la deformidad a la hora de predecir la incidencia de la kifosis en el postoperatorio.
require(rpart)
attach(kyphosis)
# Convertimos la variable respuesta factor en variable numérica para el ajuste de modelos
kyphosisb = kyphosis %>% mutate(Kyphosis = 1*(Kyphosis=="present"))
kyphosisb
# Gráficos de predictoras y variable respuesta
# Utilizamos la función melt con todas las varaibles, ya que todas son numéricas
datacomp = melt(kyphosis, id.vars='Kyphosis')
ggplot(datacomp) +
  geom_boxplot(aes(Kyphosis,value, colour=variable)) + 
  facet_wrap(~variable, scales ="free_y") +
  labs(x = "", y = "Kyphosis") +
  theme_bw()

  • Información agrupda para todos los sujetos que identifican éxito o fracaso en función de diferentes varaibles de clasificación. Ejemplo: Collet (1991) presenta un experimento sobre la toxicidad de distintas dosis (en microgramos) del piretroide trnas-cipemetrín en los capullos de gusano del tabaco. Se había comenzado a detectar resistencia de esas polillas a dicho tóxico (a determinadas dosis). El experimento consitía en exponer, durante tres dias y a distintas dosis de tóxico, a series de 20 polillas de cada sexo. Se anotaron el número de polillas muertas en cada serie. Es de interés en el análisis investigar si efectivamente se demostraba tal resistencia al tóxico (en machos y hembras) y a partir de qué dosis. También se desea determinar la dosis a la cual es posible garantizar el exterminio del 50% de los insectos.
glm_bin_04 = read_csv("https://goo.gl/w23RGz", col_types = "cdii")
# Calculamos los vivos, la probabilidad de morir, y el logartimo de dosis que es la forma habitual de medir en este tipo de situaciones
glm_bin_04 = glm_bin_04 %>%  
  mutate(alive = total-dead, probabilidad = dead/total, ldosis = log(dosis))
glm_bin_04
# Representamos la probabilidad de morir en función de las covariables sex y ldosis
ggplot(glm_bin_04,aes(x = ldosis, y = probabilidad, color = sex)) + geom_point() +theme_bw()

En el primer ejemplo la variable respuesta se puede representar mediante el modelo \(Y_i \sim Bi(n = 1, \pi_i)\), dado que cada observación corresponde a un único sujeto con \(\pi_i\) la probabilidad de sufrir una deformidad postoperatoria para el sujeto \(i\), mientras que en el segundo ejemplo tenemos un modelo \(Y_i \sim Bi(n = 20, \pi_i)\), con \(\pi_i\) la probabilidad de morir para una polilla de ese grupo.

GLM respuesta binaria

Las hipótesis que debe verificar este modelo son:

  • independencia entre las observaciones
  • linealidad entre transformaciones de la proporción de éxitos y de las variables explicativas continuas (función link)
  • consistencia entre la modelización y la interpretación física.

Si \(\pi_i\) es la probabilidad de éxito asociado al i-ésimo experimento binomial con respuesta \(Y_i\), y condiciones de experimentación observadas en las variables predictoras \(X_1;X_2,...,X_p\), la formulación del GLM viene dada por:

\[g(\pi_i) = \eta_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{21} + ... + \beta_p x_{ip}\]

Estimación

Para la estimación utilizamos la función glm con las especificaciones siguientes:

# Modelo de regresión logística
glm(modelo,family = binomial(link = logit),data_set)
# Modelo de regresión probit
glm(modelo,family = binomial(link = probit),data_set)
# Modelo cloglog
glm(modelo,family = binomial(link = cloglog),data_set)

Sin embargo, la especificación del modelo varía en función del tipo de datos (individualizados o conteos). En el caso de nuestro banco de datos contenga información individualizada para cada sujeto la variable respuesta se debe codificar con un 1 para el éxito y 0 para el fracaso, y el modelo se especifica como:

respuesta ~ predictoras

En el caso de disponer del número de éxitos y fracasos el modelo se especifica e la forma siguiente:

cbind(exitos,fracasos) ~ predictoras

Estudiamos ahora le proceso de estimación del modelo completo (con todos los efectos posibles) para cada una de los ejemplos considerados.

Kyphosis

En este caso ya hemos creado un variable para identificar el éxito (tener malformación postoperatoria) o el fracaso (no tener malformación) para cada sujeto de la muestra. El banco de datos es kyphosisb con variable respuesta Kyphosis y con tres posibles variables predictoras de tipo numérico: Age, Number, Start. Para ajustar el modelo de regresión logística escribimos:

M1 <- glm(Kyphosis  ~ Age + Number + Start,family = binomial(link = logit),kyphosisb)
# Resumen del modelo ajustado
summary(M1)
## 
## Call:
## glm(formula = Kyphosis ~ Age + Number + Start, family = binomial(link = logit), 
##     data = kyphosisb)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3124  -0.5484  -0.3632  -0.1659   2.1613  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -2.036934   1.449575  -1.405  0.15996   
## Age          0.010930   0.006446   1.696  0.08996 . 
## Number       0.410601   0.224861   1.826  0.06785 . 
## Start       -0.206510   0.067699  -3.050  0.00229 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 83.234  on 80  degrees of freedom
## Residual deviance: 61.380  on 77  degrees of freedom
## AIC: 69.38
## 
## Number of Fisher Scoring iterations: 5

En el modelo ajustado podemos observar que individualmente sólo se considera significativo el efecto asociado con la variable Start. En la parte inferior aparece la deviance asociada al modelo ajustado (Null deviance) y la deviance que no es capaz de explicar el modelo ajustado (Residual deviance). Para valorar la capacidad explicativa del modelo recordemos que debemos comparar la deviance del modelo ajustado con el cuantil 0.95 de una \(\chi^2\) con grados de libertad igual a los correspondientes a la deviance residual. En lugar de comparar los estadísticos obtenemos el pvalor asociado que viene dado por

1-pchisq(M1$deviance,M1$df.residual)
## [1] 0.9033442

El valor de dicho pvalor es superior a 0.05 indicando que el modelo considerado tiene una buena capacidad explicativa. En el punto siguiente estudiaremos la selección de variables para quedarnos con el mejor modelo posible.

Collet

En este caso trabajamos con datos agrupados que contabilizan los éxitos (número de polillas muertas) y los fracasos (número de polillas vivas) para cada combinación considerada en el diseño experimental. Se consideran las variables predictoras sex (categórica) y dosis (numérica). El banco de datos esglm_bin_04` y contiene además las polillas totales, muertas y vivas (calculadas como muertas - vivas), y el logaritmo de la dosis suministrada, ya que de forma habitual en este tipo de modelos se suele trabajar en dicha escala logarítmica. Para ajustar el modelo de regresión logística consideramos los efectos asociados con cada predictora, así como la posible interacción entre ambos:

Yres <- cbind(glm_bin_04$dead,glm_bin_04$alive)
M1 <- glm(Yres  ~ sex + ldosis + sex:ldosis,family = binomial(link = logit),glm_bin_04)
# Resumen del modelo ajustado
summary(M1)
## 
## Call:
## glm(formula = Yres ~ sex + ldosis + sex:ldosis, family = binomial(link = logit), 
##     data = glm_bin_04)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.39849  -0.32094  -0.07592   0.38220   1.10375  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.9935     0.5527  -5.416 6.09e-08 ***
## sexM          0.1750     0.7783   0.225    0.822    
## ldosis        1.3071     0.2411   5.422 5.89e-08 ***
## sexM:ldosis   0.5091     0.3895   1.307    0.191    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 124.8756  on 11  degrees of freedom
## Residual deviance:   4.9937  on  8  degrees of freedom
## AIC: 43.104
## 
## Number of Fisher Scoring iterations: 4

De los resultados obtenidos parece desprenderse que no existe efecto de interacción entre sexo y dosis. Antes de tratar de explicar el resto de efectos deberemos estudiar la posibilidad de eliminar dicho efecto del modelo. Por otro lado, el pvalor para valorar la bondad del ajuste:

1-pchisq(M1$deviance,M1$df.residual)
## [1] 0.7582464

parece indicar que el ajuste obtenido es bueno, dado que es superior a 0.05.

Selección del modelo

Para la construcción del mejor modelo utilizaremos los mismos procedimientos secuenciales de selección de efectos que vimos en el tema de los modelos de regresión lineal múltiple, modelos anova y ancova. En este caso disponemos del criterio Deviance, del AIC y del test \(\chi^2\) (pvalor significativo o no) para valorar los efectos del modelo. Para realizar la selección del modelo utilizaremos el procedimiento:

modelo.final <- step(modelo, test = "Chisq")

Estudiamos a continuación cada uno de los ejemplos.

Kyphosis

M1 <- glm(Kyphosis  ~ Age + Number + Start,family = binomial(link = logit),kyphosisb)
Mfinal <- step(M1, test = "Chisq")
## Start:  AIC=69.38
## Kyphosis ~ Age + Number + Start
## 
##          Df Deviance    AIC     LRT Pr(>Chi)   
## <none>        61.380 69.380                    
## - Age     1   64.536 70.536  3.1565 0.075623 . 
## - Number  1   65.299 71.299  3.9191 0.047739 * 
## - Start   1   71.627 77.627 10.2466 0.001369 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El proceso de selección de efectos indica que no podemos eliminar ninguna variable del modelo. Los valores de Deviance siempre aumentan (deviance asociada a cada -variable) cuando eliminamos cualquiera de las variables con respecto a la del modelo que se queda con todas ellas (). Lo mismo ocurre con el estadístico AIC. En cuanto al test Chi cuadrado tenemos un pvalor superior a 0.05 para la variable Age, pero dado que estos pvalores son aproximados deben ser tomados con cautela y solo deben ser interpretados cuando claramente reflejen una elección. El análisis se realiza valorando todos los criterios y no sólo uno de ellos. El modelo resultante:

Mfinal
## 
## Call:  glm(formula = Kyphosis ~ Age + Number + Start, family = binomial(link = logit), 
##     data = kyphosisb)
## 
## Coefficients:
## (Intercept)          Age       Number        Start  
##    -2.03693      0.01093      0.41060     -0.20651  
## 
## Degrees of Freedom: 80 Total (i.e. Null);  77 Residual
## Null Deviance:       83.23 
## Residual Deviance: 61.38     AIC: 69.38

permite escribir el predictor lineal mediante la expresión siguiente: \[log\left(\frac{\pi_i}{1-\pi_i}\right) = -2.03693 + 0.01093 * Age_i + 0.41060 * Number_i - 0.20651 * Start_i\] donde \(\pi_i\) es la probabilidad de tener una malformación postoperatoria. En este caso la interpretación de los signos de los coeficientes es similar a la de los modelos lineales. Un signo positivo indica que la probabilidad de presentar la malformación aumenta con dicho efecto, pero disminuye cuando el coeficiente es negativo. En este caso la probabilidad aumenta cuando el niño tiene más edad y el número de vértebras intervenidas es mayor, pero disminuye en función de la primera vértebra involucrada en la operación.

Podemos obtener todo el proceso de inferencia (estimaciones e intervalos de confianza) del modelo con:

cbind(tidy(Mfinal),confint_tidy(Mfinal))

Collet

Yres <- cbind(glm_bin_04$dead,glm_bin_04$alive)
M1 <- glm(Yres  ~ sex + ldosis + sex:ldosis,family = binomial(link = logit),glm_bin_04)
Mfinal <- step(M1, test = "Chisq")
## Start:  AIC=43.1
## Yres ~ sex + ldosis + sex:ldosis
## 
##              Df Deviance    AIC    LRT Pr(>Chi)
## - sex:ldosis  1   6.7571 42.867 1.7633   0.1842
## <none>            4.9937 43.104                
## 
## Step:  AIC=42.87
## Yres ~ sex + ldosis
## 
##          Df Deviance     AIC     LRT  Pr(>Chi)    
## <none>         6.757  42.867                      
## - sex     1   16.984  51.094  10.227  0.001384 ** 
## - ldosis  1  118.799 152.909 112.042 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El proceso de selección de efectos indica que no debemos considerar el efecto de interacción ya que el valor de deviance es inferior cuando está presente que cuando lo eliminamos, el AIC aumenta al eliminarlo, y el pvalor es claramente no significativo. El resto de efectos del modelo no pueden ser eliminados. Esto implica que al no haber efecto de interacción la probabilidad de morir se puede obtener mediante comportamientos paralelos para machos y hembras. El modelo resultante:

Mfinal
## 
## Call:  glm(formula = Yres ~ sex + ldosis, family = binomial(link = logit), 
##     data = glm_bin_04)
## 
## Coefficients:
## (Intercept)         sexM       ldosis  
##      -3.473        1.101        1.535  
## 
## Degrees of Freedom: 11 Total (i.e. Null);  9 Residual
## Null Deviance:       124.9 
## Residual Deviance: 6.757     AIC: 42.87

permite escribir el predictor lineal mediante las expresiones siguientes (atendiendo a los diferentes valores de sex): \[log\left(\frac{\pi_i}{1-\pi_i}\right)_{Machos} = -3.473 + 1.101 + 1,535 * ldosis_{i,Machos} = -2.372 + 1,535 * ldosis_{i,machos}\] \[log\left(\frac{\pi_i}{1-\pi_i}\right)_{Hembras} = -3.473 + 1,535 * ldosis_{i,Hembras} \]

donde \(\pi_i\) es la probabilidad de muerte de la polilla. Dado que la interceptación es más pequeña para las hembras que para los machos tenemos un indicador de que la hembras son más resistentes que los machos, dado que su probabilidad de muerte para una misma dosis es inferior. La bondad de ajuste de este modelo:

1-pchisq(Mfinal$deviance,Mfinal$df.residual)
## [1] 0.6623957

muestra que el ajuste obtenido puede considerarse como bueno. Podemos obtener todo el proceso de inferencia (estimaciones e intervalos de confianza) del modelo con:

cbind(tidy(Mfinal),confint_tidy(Mfinal))

Para verificar el efecto de la dosis por sexo podemos calcular cuál es la dosis necesaria aplicar a cada sexo para conseguir una probabilidad de muerte del 50%, o como habitualmente se conoce con el nombre de dosis letal al 50% (LD50). Para ello basta con sustituir en las ecuaciones anteriores el valor de \(\pi_i\) por 0.5 y despejar el valor de ldosis: \[ldosis50_{Machos} = \frac{2.372}{1.535} = 1.54277 \to dosis50_{Machos} = exp(1.54277) = 4.68927 \] \[ldosis50_{Hembras} = \frac{3.473}{1.535} = 2.262541 \to dosis50_{Hembras} = exp(1.54277) = 9.607468 \] De los resultados obtenidos podemos ver que hay que aplicar el doble de dosis en las hembras para conseguir la misma probabilidad de muerte.

Para conseguir cualquier otra dosis letal basta con sustituir \(\pi_i\) por la correspondiente probabilidad y despejar en las ecuaciones obtenidas, o utilizar el código siguiente para obtener la dosis letal a una probabilidad “prob”:

# prob = Probabilidad buscada
predictor <- binomial(link = logit)$linkfun(prob)
machos <-exp((predictor+2.372)/1.535)
hembras <- exp((predictor+3.473)/1.535)

Utilizamos el código anterior para representar y estudiar las curvas de dosis letales para ambos sexos (machos = azul, hembras = rojo):

prob <- seq(0.05,0.95,0.01)
# Trabajamos con dosis 
predictor <- binomial(link = logit)$linkfun(prob)
machos <-exp((predictor+2.372)/1.535)
hembras <- exp((predictor+3.473)/1.535)
dosis <- data.frame(prob,machos,hembras)
ggplot(dosis) + 
  geom_line(aes(x = prob, y = machos),color="blue") + 
  geom_line(aes(x = prob, y = hembras),color="red") + 
  scale_x_continuous(breaks = seq(0,1,0.05)) +
  scale_y_continuous(breaks = seq(0,70,5)) +
  labs(x = "Probabilidad de muerte", y = "Dosis") + 
  theme_bw()

Diagnóstico

Antes de comenzar el proceso de diagnóstico hay que recordad que en este tipo de modelos tenemos dos tipos de valores ajustados y dos tipos de residuos:

  • Los correspondientes al predictor lineal obtenidos con las estimaciones de los parámetros del modelo
  • Los correspondientes a los valores originales que se obtienen al deshacer la transformación de la función de enlace.

El diagnóstico de los GLM se basa en los procedimientos gráficos, que vimos en temas anteriores, aplicado a los residuos deviance obtenidos a partir del predictor lineal considerado (\(X\widehat{\beta}\)). Tanto para obtener los valores ajustados y residuos del predictor lineal utilizaremos la función fortify. Al resultado de dicha función podemos añadir sin muchos problemas los residuos entre valores originales de la respuesta y predichos. En el caso del modelo de regresión logística el código necesario viene dado por:

# Sólo hay que sustituir "modelo" por el modelo ajustado a nuestros datos
diagnostico <- fortify(modelo)
diagnostico$fitoriginal <- predict.glm(modelo,type = "response")
diagnostico$residoriginal <- residuals.glm(modelo,type = "response")

Los gráficos que vamos a utilizar son:

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

Kyphosis

Obtenemos todas las cantidades de interés del modelo

# Modelo final obtenido
M1 <- glm(Kyphosis  ~ Age + Number + Start,family = binomial(link = logit),kyphosisb)
# Obtención de valores para el diagnóstico
diagnostico <- fortify(M1)
diagnostico$fitoriginal <- predict.glm(M1,type = "response")
diagnostico$residoriginal <- residuals.glm(M1,type = "response")
diagnostico

Realizamos los gráficos de ajuste

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- diagnostico[,c(".fitted","Age","Number","Start",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method=lm, se=FALSE) +
  geom_smooth(aes(value,.stdresid, colour=variable), method=loess, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Como se puede observar en el gráfico de ajustados versus residuos la interpretación en este tipo de modelos se hace bastante complicada, Esto es debido a que la variable respuesta sólo toma valores 0 y 1, lo que motiva que el gráfico tenga esa pinta tan extraña. Si observamos las distancias de Cook obtenidas podemos ver que no hay ninguna superior a 1, y por tanto no tenemos ninguna observación influyente.

En cuanto a los gráficos con respecto a las variables predictoras parece apreciarse ciertas tendencias con respecto a ellas. Este comportamiento parece más evidente con la variable edad donde se ve una parábola indicando la posible existencia de un efecto polinómico de grado 2 con respecto a ella. En las otras dos predictoras el efecto es menos apreciable. La opción más habitual sería ajustar un nuevo modelo indicando dicha tendencia de grado 2, pero se podría optar también por la inclusión de efectos de suavizado sobre cada una de ellas para evitar esas tendencias observadas. En principio optamos por incluir efectos de suavizado por edad y start (ambos con tendencias parabólicas) y comparamos ambos modelos. Para el ajuste de estos modelos únicamente debemos cambiar el comando glm por gam:

# Modelo sin suavizado
M1 <- glm(Kyphosis  ~ Age +  Number + Start,family = binomial(link = logit),kyphosisb)
# Modelo con suavizado en edad y start
M2 <- gam(Kyphosis  ~ s(Age, k = 10, m = 2, bs = "ps")  + Number + s(Start, k = 10, m = 2, bs = "ps"),family = binomial(link = logit),kyphosisb)
AIC(M1,M2)

El menor valor de AIC indica que el modelo con efectos de suavizado es preferible al que no los tiene. Estudiamos con un poco más detalle dicho modelo:

# Resumen del modelo ajustado
summary(M2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## Kyphosis ~ s(Age, k = 10, m = 2, bs = "ps") + Number + s(Start, 
##     k = 10, m = 2, bs = "ps")
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -3.6011     1.1482  -3.136  0.00171 **
## Number        0.3333     0.2324   1.434  0.15160   
## ---
## 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(Age)   2.151  2.671  6.345   0.071 .
## s(Start) 1.956  2.409  9.747   0.013 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.354   Deviance explained = 39.3%
## UBRE = -0.22572  Scale est. = 1         n = 81

En los efectos paramétricos se comprueba que el efecto asociado con number no resulta significativo, mientras que en los no paramétricos vemos que el suavizado con edad no resulta significativo al 95%. Construimos todos los posibles modelos y vemos cual de ellos es el mejor:

# Modelo con suavizados y number
M1 <- gam(Kyphosis  ~ s(Age, k = 10, m = 2, bs = "ps")  + Number + s(Start, k = 10, m = 2, bs = "ps"),family = binomial(link = logit),kyphosisb)
# Modelo con suavizados y sin number
M2 <- gam(Kyphosis  ~ s(Age, k = 10, m = 2, bs = "ps") + s(Start, k = 10, m = 2, bs = "ps"),family = binomial(link = logit),kyphosisb)
# Modelo con suavizados con Start unicamente
M3 <- gam(Kyphosis  ~ s(Start, k = 10, m = 2, bs = "ps"),family = binomial(link = logit),kyphosisb)
# Valores de AIC
AIC(M1,M2,M3)
# Valores de GCV
c(M1$gcv.ubre,M2$gcv.ubre,M3$gcv.ubre)
##     GCV.Cp     GCV.Cp     GCV.Cp 
## -0.2257187 -0.2223853 -0.1562664

El menor valor de AIC corresponde al modelo con los dos suavizados y la variable number. Hay que tener en cuenta que las significatividades que parecen en las tablas son aproximaciones asintóticas, y por tanto siempre se deben tomar con cautela y proceder con otro tipo de comparación. El estadístico GCV también proporciona la misma conclusión.

# Resumen del modelo ajustado
summary(M1)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## Kyphosis ~ s(Age, k = 10, m = 2, bs = "ps") + Number + s(Start, 
##     k = 10, m = 2, bs = "ps")
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -3.6011     1.1482  -3.136  0.00171 **
## Number        0.3333     0.2324   1.434  0.15160   
## ---
## 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(Age)   2.151  2.671  6.345   0.071 .
## s(Start) 1.956  2.409  9.747   0.013 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.354   Deviance explained = 39.3%
## UBRE = -0.22572  Scale est. = 1         n = 81

Por tanto el modelo final ajustado viene dado por:

\[log\left(\frac{\pi_i}{1-\pi_i}\right) = -3.6011 + 0.3333 * Number_i + s(Age_i) + s(Start_i)\] Podemos representar el efectos de los suavizados para observar el comportamiento de dichas predictoras:

par(mfrow = c(1,2))
plot.gam(M1)

Se observa la parábola que describe la variable Age y la función monótona decreciente asociada con la variable Star. En conclusión la probabilidad de malformación es más alta en valores de edad próximos a los 100 meses, y aunque se muestra constante para cuando la vértebra de inicio se encuentra entre las cinco primeras, dicha probabilidad disminuye cuando la vértebra está más alejada. En cuanto a la variable number la probabilidad de malformación aumenta cuando el número de vértebras operadas es mayor.

Realizamos el diagnóstico de este modelo:

# Suavizado
gam.check(M1)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-5.641377e-09,6.154328e-08]
## (score -0.2257187 & scale 1).
## Hessian positive definite, eigenvalue range [0.007746516,0.01021528].
## Model rank =  20 / 20 
## 
## 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(Age)   9.00 2.15    1.11    0.89
## s(Start) 9.00 1.96    1.14    0.91

Ningún pvalor resulta significativo indicando que el suavizado utilizado es adecuado.

Collet

Obtenemos las cantidades de interés del modelo

Yres <- cbind(glm_bin_04$dead,glm_bin_04$alive)
M1 <- glm(Yres  ~ sex + ldosis ,family = binomial(link = logit),glm_bin_04)
# Obtención de valores para el diagnóstico
diagnostico <- fortify(M1)

Realizamos los gráficos de residuos versus ajustados y de residuos vs ldosis

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- diagnostico[,c(".fitted","ldosis","sex",".stdresid")] 
# Residuos vs ajustados
ggplot(diagnosis) +
geom_jitter(aes(.fitted,.stdresid, color = sex),) +
  theme_bw()

# Residuos vs predictora
ggplot(diagnosis) +
geom_jitter(aes(ldosis,.stdresid,color = sex),) +
  theme_bw()

En este caso no hemos ajustado las tendencias, ya que dado el tamaño de muestra tan pequeño en cada uno de los grupos, los resultados podrían mostrar tendencias ficticias que nos podrían hacer dudar de la validez del modelo. En estas situaciones es preferible no realizar cambios, salvo un comportamiento claramente extraño, y seguir con el modelo obtenido. No se aprecian tendencias destacables con respecto a ajustados y a las variables predictoras. Dado el número tan bajo de observaciones el planteamiento de modelos con suavizado sobre dosis no resulta una opción. Además no se detecta ninguna observación influyente ya que la distancia de Cook es inferior a 1 para todos los sujetos.

Predicción

La predicción en este tipo de modelos se divide en dos fases: i) predicción del predictor lineal, y ii) predicción de la respuesta, aunque evidentemente la más interesante es la correspondiente a la respuesta.

Sin embargo, cuando la respuesta viene individualizada, la predicción obtenida es un valor comprendido entre 0 y 1 y hay que definir una regla para clasificar finalmente al sujeto con un 1 o un 0, para identificar el éxito o el fracaso. En esta situación se toma la regla:

  • si la predicción es mayor o igual a 0.5 asignamos un éxito como resultado de la predicción
  • si la predicción es menor a 0.5 asignamos un fracaso como resultado de la predicción

El resultado de la predicción es una tabla de doble entrada donde comparamos los valores observados frente a los predichos calculados con la regla anterior. En esta situación estamos interesados en el porcentaje de observaciones que son correctamente clasificadas como éxitos o fracasos.

Este problema no aparece cuando tenemos datos agrupados ya que la variable respuesta ya es un valor entre 0 y 1 al modelizar directamente la probabilidad de éxito.

Veamos la aplicación a cada uno de los ejemplos. Utilizaremos la función fortify y completaremos las predicciones de la respuesta.

Kyphosis

Dado que en este caso tenemos datos individualizamos, calcularemos la predicción de la respuesta y utilizaremos la regla de clasificación presentada anteriormente. Finalmente calcularemos el porcentaje de clasificación correcta.

# Modelo final obtenido
M1 <- gam(Kyphosis  ~ s(Age, k = 10, m = 2, bs = "ps")  + Number + s(Start, k = 10, m = 2, bs = "ps"),family = binomial(link = logit),kyphosisb)
# Obtención de predicción de la respuesta
prediccion <- predict.gam(M1,type = "response")
# Clasificacmos a cada sujeto como éxito o fracaso
clasificado <- 1*(prediccion>=0.5)
tabla <- table(kyphosis$Kyphosis,clasificado)
tabla
##          clasificado
##            0  1
##   absent  60  4
##   present  7 10
# Porcentaje de calsificación correcta
round(100*sum(diag(tabla))/sum(tabla),2)
## [1] 86.42

Tenemos un 86.42% de clasificar correctamente a un individuo como éxito o fracaso en función de las variables predictoras consideradas.

Collet

Obtenemos las cantidades de interés del modelo

Yres <- cbind(glm_bin_04$dead,glm_bin_04$alive)
M1 <- glm(Yres  ~ sex + ldosis ,family = binomial(link = logit),glm_bin_04)
# Obtención de predicción de la respuesta y el error estándar
prediccion <- predict.glm(M1,type = "response", se.fit = TRUE)

Dado que tenemos datos agrupados, y que en cada grupo hay 20 sujetos, podemos calcular el número esperado de éxitos en cada combinación como la probabilidad obtenida con la predicción multiplicado por el número de sujetos (20).

# Número esperado de éxitos por combinación
exitos <- round(prediccion$fit*glm_bin_04$total,0)
exitos
##  1  2  3  4  5  6  7  8  9 10 11 12 
##  2  4  9 14 17 19  1  2  4  9 14 17
# Comparamos exitos originales con los predichos y añadimos la comnbinación del diseño correspondiente
comparativa <- data.frame(sex = glm_bin_04$sex, dosis = glm_bin_04$dosis, original = glm_bin_04$dead, predicho = exitos)
comparativa

En la mayoría de combinaciones el número esperado de polillas muertas coincide con el número observado indicando que el modelo construido es adecuado. Por otro lado vamos a obtener las curvas de probabilidad de éxito para las variables predictoras. Representamos además la probabilidad asociada con la variable original dosis.

# Modelo
Yres <- cbind(glm_bin_04$dead,glm_bin_04$alive)
M1 <- glm(Yres  ~ sex + ldosis ,family = binomial(link = logit),glm_bin_04)
# Utlizamos la opción each para indicar el número de repeticiones (tantas como niveles del factor)
newdata <- data.frame(ldosis = rep(seq(min(glm_bin_04$ldosis), max(glm_bin_04$ldosis), .05),each=2), sex = factor(c("M", "F")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(M1, newdata, type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza y añadimos la dosis desahaciendo el cambio del logaritmo
newdata <- newdata %>% mutate(
  upr = fit + (1.96 * se.fit),
  lwr = fit - (1.96 * se.fit),
  dosis = exp(ldosis)
)

# Gráfico de la predicción con log dosis
ggplot(newdata, aes(x = ldosis, y = fit, color = sex)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr, fill = sex), alpha = 1/5) +
  scale_x_continuous(breaks = seq(0,4,0.25)) +
  labs(x = "Logaritmo Dosis", y = "Probabilidad de morir", title = "Bandas de predicción") +
  theme_bw()

# Gráfico de la predicción con dosis
ggplot(newdata, aes(x = dosis, y = fit, color = sex)) +
  geom_line() +
  scale_x_continuous(breaks = seq(0,60,2.5)) +
  labs(x = "Dosis", y = "Probabilidad de morir", title = "") +
  theme_bw()

La conclusión general es que las hembras son más resistentes que los machos.

Bibliografía

Copyright © 2018 Javier Morales. Universidad Miguel Hernández de Elche.