Regresión Logística sobre Tabla de Contingencia

Autor: Sergio García Prado

Fecha: Diciembre de 2018

Descripción:

La respuesta binaria a cierto estimulo está relacionada con unas variables x (numerica), t (cualitativa) y z (cualitativa). Se han realizado 100 experimentos aleatorios en cada una de las combinaciones correspondientes a diferentes niveles de las variables anteriores. La tabla datos contiene las variables explicativas junto con la variable y (número de ocurrencias de una de las respuestas en los 100 experimentos realizados). Identifica el modelo más interesante para ajustar. Consigue intervalos de confianza para las odds ratio (Wald y verosimilitud perfil), para la dosis letal 40% (Wald) y para las predicciones en los niveles de las variables explicativas estudiados (Wald). Representa las curvas correspondientes a la predicción.

Configuración de Entorno

In [1]:
rm(list = ls())
In [2]:
library(magrittr, warn.conflicts = FALSE)
library(dplyr, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)
In [3]:
options(repr.plot.width = 8, repr.plot.height = 4)

Lectura de Datos

In [4]:
tries <- 100
In [5]:
datos <- read.table('datos.txt', header = TRUE)
In [6]:
colnames(datos)[2] <- "x"
datos$z <- as.factor(datos$z)
datos$t <- as.factor(datos$t)

Análisis Exploratorio de los Datos

In [7]:
summary(datos)
       y               x        z      t     
 Min.   : 0.00   Min.   : 1.0   0:50   1:20  
 1st Qu.: 8.00   1st Qu.: 3.0   1:50   2:20  
 Median :36.00   Median : 5.5          3:20  
 Mean   :39.61   Mean   : 5.5          4:20  
 3rd Qu.:65.50   3rd Qu.: 8.0          5:20  
 Max.   :98.00   Max.   :10.0                
In [8]:
sum(datos$y)
3961
In [9]:
xtabs(y ~ x + t + z, datos)
, , z = 0

    t
x     1  2  3  4  5
  1   2  1  2  0  0
  2   3  2  3  2  0
  3   1  4  1  1  1
  4   7  3  6  4  1
  5   7 14  8  4  9
  6   8  8  6  8 17
  7  14 20 17 17 18
  8  27 28 26 28 36
  9  41 32 38 38 39
  10 48 49 50 50 56

, , z = 1

    t
x     1  2  3  4  5
  1  18 17 16 22 18
  2  25 29 37 27 22
  3  38 41 38 36 43
  4  54 47 54 53 47
  5  61 64 60 56 64
  6  70 76 71 76 78
  7  83 81 73 82 80
  8  87 88 90 90 87
  9  94 98 92 98 94
  10 98 96 94 96 97
In [10]:
pairs(y ~ x + t + z, datos, pch = 20)

Ajuste del Modelo

In [11]:
datos.tunned <- datos %>% 
    rename(sucesses = y) %>%
    mutate(failures = tries - sucesses)

Modelo Logístico Completo (full)

In [12]:
model.logit.full <- glm(cbind(sucesses, failures) ~ x + t + z, data = datos.tunned, 
                     family = binomial(link = logit))
In [13]:
summary(model.logit.full)
Call:
glm(formula = cbind(sucesses, failures) ~ x + t + z, family = binomial(link = logit), 
    data = datos.tunned)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.28796  -0.60134   0.04575   0.49385   2.41219  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.163053   0.117652 -43.884   <2e-16 ***
x            0.515782   0.012217  42.218   <2e-16 ***
t2           0.044675   0.086293   0.518    0.605    
t3          -0.014940   0.086429  -0.173    0.863    
t4           0.007461   0.086376   0.086    0.931    
t5           0.078046   0.086224   0.905    0.365    
z1           3.101047   0.067702  45.804   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5187.820  on 99  degrees of freedom
Residual deviance:   87.865  on 93  degrees of freedom
AIC: 504.62

Number of Fisher Scoring iterations: 4

Como vemos, en el modelo con todas las variables, la influencia de las categorías formadas por la variable cualitativa t no aporta información significativa para el ajuste de la variable $y$. Por lo tanto, se eliminará del modelo para así tratar de simplificar el mismo.

Modelo Logístico Restringido

In [14]:
model.logit <- glm(cbind(sucesses, failures) ~ x + z, data = datos.tunned, 
             family = binomial(link = logit))
In [15]:
summary(model.logit)
Call:
glm(formula = cbind(sucesses, failures) ~ x + z, family = binomial(link = logit), 
    data = datos.tunned)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.43976  -0.67094   0.02646   0.51833   2.36379  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.13896    0.10388  -49.47   <2e-16 ***
x            0.51568    0.01221   42.22   <2e-16 ***
z1           3.10041    0.06769   45.80   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5187.820  on 99  degrees of freedom
Residual deviance:   89.401  on 97  degrees of freedom
AIC: 498.15

Number of Fisher Scoring iterations: 4

Con el modelo ajustado únicamente por las variables $X$ y $Z$ se consigue que ahora todos los parámetros sean significativos. A pesar de que en este caso la deviance residual del modelo es mayor ($89.401$ frente a $87.865$ del modelo con todas las variables predictoras), el AIC es menor. Esto es porque la reducción del número de parámetro a estimar "compensa" el peor ajuste del modelo.

Modelo Normal Inverso (probit)

Para comprobar el ajuste sobre un modelo lineal generalizado con respuesta de tipo binomial diferente al logístico, se va a realizar un ajuste utilizando la función de enlace normal estándar inversa.

In [16]:
model.probit <- glm(cbind(sucesses, failures) ~ x + z, data = datos.tunned,
                    family = binomial(link = probit))
In [17]:
summary(model.probit)
Call:
glm(formula = cbind(sucesses, failures) ~ x + z, family = binomial(link = probit), 
    data = datos.tunned)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0931  -0.6003  -0.1092   0.6668   2.3895  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.937276   0.053728  -54.67   <2e-16 ***
x            0.293968   0.006519   45.09   <2e-16 ***
z1           1.777507   0.035824   49.62   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5187.820  on 99  degrees of freedom
Residual deviance:   92.839  on 97  degrees of freedom
AIC: 501.59

Number of Fisher Scoring iterations: 4

En este caso, el ajuste sigue siendo similar al logístico, obteniendo una deviance residual similiar, junto con un AIC menor que el obtenido por el modelo logístico completo (por el uso de un menor número de variables). Sin embargo, en este caso el modelo no obtiene mejores resultados que el modelo logístico restringido al uso de las variables $X$ y $Z$.

A continuación se incluyen las tablas resumen que recogen los valores del Criterio de Akaike y del Criterio Bayesiano, donde ambos indican el mejor ajuste del modelo logístico restringido.

In [18]:
AIC(model.logit.full, model.logit, model.probit)
dfAIC
model.logit.full7 504.6179
model.logit3 498.1535
model.probit3 501.5914
In [19]:
BIC(model.logit.full, model.logit, model.probit)
dfBIC
model.logit.full7 522.8541
model.logit3 505.9690
model.probit3 509.4069

Por lo tanto, el modelo que se ha elegido es el siguiente (donde $y$ en este caso representaría la observacion sin agrupar y $logit(p) = log(\frac{p}{1-p})$):

$$E[logit(y)] = X + Z$$

Intervalos de Confianza

Los intervalos de confianza se van a calcular utilizando un nivel de confianza $\alpha = 0.01$.

In [20]:
alpha <- 0.01
In [21]:
logit <- function(p) {
    log(p / (1 - p))
}
In [22]:
inv.logit <- function(x) {
    exp(x) / (1 + exp(x))
}

Odds-Ratio

Los intervalos de confianza para las odd ratios se pueden conseguir a partir de las funciones confint (verosimilitud perfil) y confint.default (Wald) de manera sencilla.

Verosimilitud Perfil

In [23]:
confint(model.logit)
Waiting for profiling to be done...
2.5 %97.5 %
(Intercept)-5.3450602-4.9378172
x 0.4919734 0.5398608
z1 2.9690507 3.2344230

Wald

In [24]:
confint.default(model.logit)
2.5 %97.5 %
(Intercept)-5.3425549-4.9353591
x 0.4917373 0.5396189
z1 2.9677440 3.2330787

Los intervalos basados en ambos métodos son muy similares, por lo que interpretaremos de manera indistinta ambos. La interpretación del término indica que cuando la variable $X$ toma el valor $0$ y la variable categórica $Z$ está en la categoría $\{Z = 0\}$, la probabilidad de $y$ es prácticamente $0$. En cuanto al oddratio para la variable $X$, podemos decir que por cada unidad que esta aumenta, la ventaja aumenta un $50\%$. En cuanto a la variable categórica $Z$, diremos que cuando esta se fija en la categoría $\{Z = 1\}$ la ventaja aumenta un $300\%$.

Dosis Letal (40%)

A continuación se muestra el intervalo de confianza de wald de la variable regresora $X$ cuando la variable $Y$ alcanza la dosis letal del $40\%$. Para ello nos podemos apoyar en la función dose.p del paquete MASS. Sin embargo, está pensada para casos sencillos en que únicamente hay término independiente y una variable continua (en este caso además tenemos una variable binaria).

In [25]:
coefs <- as.vector(model.logit$coefficients)
coefs.V <- vcov(model.logit)

Una vez tenemos los valores de los parámetros, así como su matriz de varianzas convarianzas, podemos aplicar el método delta para obtener la función de Dosis Letal del $40\%$. En este caso se ha hecho para los dos valores que puede tomar la variable $Z$.

In [26]:
g <- function(x) {
    c((logit(0.4) - x[1]) / x[2],
      (logit(0.4) - x[1] - x[3]) / x[2])
}

G <- function(x) {
    matrix(c(1 / x[2], - (logit(0.4) - x[1]) / x[2] ^ 2, 0,
             1 / x[2], - (logit(0.4) - x[1] - x[3]) / x[2] ^ 2, 1 / x[3]), 
           2, 3, byrow = TRUE)
}

Una vez definidas las transformaciones $g$ (función de transformación) y $G$ (matriz de derivadas), podemos calcular fácilmente los intervalos de confianza de Wald.

In [27]:
g(coefs) + rbind(c(-1, 1), c(-1, 1)) * qnorm(1 - alpha / 2) * (G(coefs) %*% coefs.V %*% t(G(coefs)))
8.7473339.448086
2.8979353.335498

La primera fila se refiere al intervalo de confianza para la dosis letal (al $40\%$) de la variable $X$, cuando la variable $Z$ toma el valor $\{Z = 0\}$, mientras que la segunda fila se refiere a la misma situación, cuando la variable $Z$ toma el valor $\{Z = 1\}$. Los intervalos de confianza se han calculado con un nivel de significacia del $99\%$.

Niveles de Variables Explicativas

En este caso se calcula el nivel de la variable respuesta $Y$ respecto de $X$ y el valor que toma la variable categórica $Z$. Además, se han incluido las bandas de confianza de nivel $99\%$.

In [28]:
datos.support <- data.frame(x = rep(seq(1, 10), 2), z = as.factor(c(rep(0, 10), rep(1, 10))))

predictions <- predict(model.logit, datos.support, type = "response", se.fit = TRUE)
datos.support$predictions <- predictions$fit

datos.support$lower <- predictions$fit - qnorm(1 - alpha / 2) * predictions$se.fit
datos.support$upper <- predictions$fit + qnorm(1 - alpha / 2) * predictions$se.fit
In [29]:
ggplot(datos.support, aes(y = predictions,ymin = lower, ymax = upper, x = x, colour = z)) + 
    geom_line() + 
    geom_ribbon(alpha = 0.2, colour = NA) + 
    facet_grid(~ z) +
    coord_cartesian(ylim = c(0, 1)) +
    theme(legend.position = "top")