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.
rm(list = ls())
library(magrittr, warn.conflicts = FALSE)
library(dplyr, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)
options(repr.plot.width = 8, repr.plot.height = 4)
tries <- 100
datos <- read.table('datos.txt', header = TRUE)
colnames(datos)[2] <- "x"
datos$z <- as.factor(datos$z)
datos$t <- as.factor(datos$t)
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
sum(datos$y)
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
pairs(y ~ x + t + z, datos, pch = 20)
datos.tunned <- datos %>%
rename(sucesses = y) %>%
mutate(failures = tries - sucesses)
model.logit.full <- glm(cbind(sucesses, failures) ~ x + t + z, data = datos.tunned,
family = binomial(link = logit))
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.
model.logit <- glm(cbind(sucesses, failures) ~ x + z, data = datos.tunned,
family = binomial(link = logit))
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.
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.
model.probit <- glm(cbind(sucesses, failures) ~ x + z, data = datos.tunned,
family = binomial(link = probit))
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.
AIC(model.logit.full, model.logit, model.probit)
df | AIC | |
---|---|---|
model.logit.full | 7 | 504.6179 |
model.logit | 3 | 498.1535 |
model.probit | 3 | 501.5914 |
BIC(model.logit.full, model.logit, model.probit)
df | BIC | |
---|---|---|
model.logit.full | 7 | 522.8541 |
model.logit | 3 | 505.9690 |
model.probit | 3 | 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$$Los intervalos de confianza se van a calcular utilizando un nivel de confianza $\alpha = 0.01$.
alpha <- 0.01
logit <- function(p) {
log(p / (1 - p))
}
inv.logit <- function(x) {
exp(x) / (1 + exp(x))
}
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.
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 |
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\%$.
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).
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$.
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.
g(coefs) + rbind(c(-1, 1), c(-1, 1)) * qnorm(1 - alpha / 2) * (G(coefs) %*% coefs.V %*% t(G(coefs)))
8.747333 | 9.448086 |
2.897935 | 3.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\%$.
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\%$.
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
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")