Regresión
La Regresión Logística Simple, desarrollada por David Cox en 1958, es un método de regresión que permite estimar la probabilidad de una variable cualitativa binaria en función de una variable cuantitativa. Una de las principales aplicaciones de la regresión logística es la de clasificación binaria, en el que las observaciones se clasifican en un grupo u otro dependiendo del valor que tome la variable empleada como predictor. Por ejemplo, clasificar a un individuo desconocido como hombre o mujer en función del tamaño de la mandíbula.
Es importante tener en cuenta que, aunque la regresión logística permite clasificar, se trata de un modelo de regresión que modela el logaritmo de la probabilidad de pertenecer a cada grupo. La asignación final se hace en función de las probabilidades predichas.
La existencia de una relación significativa entre una variable cualitativa con dos niveles y una variable continua se puede estudiar mediante otros test estadísticos tales como t-test o ANOVA (un ANOVA de dos grupos es equivalente al t-test). Sin embargo, la regresión logística permite además calcular la probabilidad de que la variable dependiente pertenezca a cada una de las dos categorías en función del valor que adquiera la variable independiente.
Caso de estudio: Camarones.
En el siguiente caso de estudio se utilizan datos reales de 12 estanques ubicados en Huatabampo, Sonora. En este, los camarones tienen que pesar por lo menos 12 gramos para poder comercializarse, sin embargo no todos los camarones logran llegar a ese peso. Por lo que se busca realizar un análisis con el fin de saber, sus éxitos y fracasos, y poder sacar inferencias a partir de ello.
Camarones en estanque
Importar datos
library(pacman)
p_load("prettydoc", "DT", "xfun", "readr")camarones <- read_csv("camarones.csv")## Rows: 12 Columns: 2## -- Column specification --------------------------------------------------------
## Delimiter: ","
## dbl (2): AlimentoDiario, Exito##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.datatable(camarones)Tabla de frecuencia de los datos
tabla <- table(camarones$Exito)
tabla##
## 0 1
## 9 3Ilustrando gráficamente este fenómeno
colores <- NULL
colores[camarones$Exito == 0] <- "green"
colores[camarones$Exito == 1] <- "red"
plot(camarones$AlimentoDiario, camarones$Exito, pch=21, bg=colores,
xlab="Alimento diario", ylab="Probabilidades de defectos")
legend('bottomleft', c('No defecto', 'Si defecto'), pch = 21, col = c('green', 'red'))
Hemos usado los argumentos pch y bg para mejorar la apariencia del gráfico. También hemos usado el comando legend para incluir una leyenda explicativa.
Parece razonable, a la vista de los datos, pensar que la temperatura puede influir en la probabilidad de que los propulsores tengan defectos. En esta práctica, vamos a ajustar un modelo de regresión logística para estudiar la posible relación.
La regresión logística transforma el valor devuelto por la regresión lineal $ y = (_0 + _1 X) $ empleando una función cuyo resultado está siempre comprendido entre 0 y 1.
Existen varias funciones que cumplen esta descripción, una de las más utilizadas es la función logística (también conocida como función sigmoide):
\[ \text{función sigmoide} = \sigma(x) = \dfrac{1}{1 + e^{-x}} \tag{1} \]
Para valores de x muy grandes positivos, el valor de e−x es aproximadamente 0 por lo que el valor de la función sigmoide es 1. Para valores de x muy grandes negativos, el valor e−x tiende a infinito por lo que el valor de la función sigmoide es 0. Sustituyendo la x de la ecuación 1 por la función lineal ( $ (_0 + _1 X)$)
Se obtiene lo siguiente:
\[ P(Y=k|X = x)= \frac{1}{1+e^{-(\beta_0+\beta_1X)}} \]
Para ajustar el modelo se usa el comando glm (para modelos lineales generalizados) indicando que la respuesta es binomial mediante el argumento family:
reg <- glm(Exito ~ AlimentoDiario, data = camarones, family=binomial)
summary(reg)##
## Call:
## glm(formula = Exito ~ AlimentoDiario, family = binomial, data = camarones)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.28965 -0.68424 -0.39705 -0.00008 2.00729
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -35.1229 25.8776 -1.357 0.175
## AlimentoDiario 0.1194 0.0901 1.325 0.185
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13.496 on 11 degrees of freedom
## Residual deviance: 11.311 on 10 degrees of freedom
## AIC: 15.311
##
## Number of Fisher Scoring iterations: 5En el modelo de regresión logística la raíz de las desviaciones representa el papel de los residuos:
\[ D_i = \mp \sqrt{-2 [Y_i\log \hat p_i + (1-Y_i)\log(1-\hat p_i)]}, \]
donde el signo coincide con el signo de Yi−p̂i . En la salida anterior estas cantidades se denominan deviance residuals. Para calcular estos pseudo-residuos, podemos ejecutar res = resid(reg)
Nuestro modelo se vería así:
\[ P(Y=1|X)=\dfrac{e^{15.0429-0.2322x}}{1+e^{15.0429-0.2322x}} \]
En la tabla coeficientes, la temperatura es significativa con pvalor= 0.0320 el cual es menor que el valor de significancia de 0.05. Lo anterior nos permite concluir que la hipoesis nula se rechaza
$ H_0:{1}=0 $ $ H_1:{1} $
Con respecto al intercepto: También es sinificativo dado que p valor= 0.0415< 0.05. Lo que significa que se rechaza la hipotesis nula, que afirma el intercepto es cero. Más aun, significa que bajísimos, lo que demuestra que ambas variables son importantes para explicar la variable dependiente (probabilidad de éxito).
Predicción para valores nuevos con el modelo ajustado
Para representar gráficamente la función logística estimada, calculamos las probabilidades de fallo estimadas (usando el comando predict) para un vector adecuado de nuevas cantidades de alimento diario (entre 250 y 310):
- Primero generamos una nueva secuencia de datos, el cula llamaremos datos nuevos.
datos_nuevos <- data.frame(AlimentoDiario = seq(260,310 , 0.5))
datos_nuevos## AlimentoDiario
## 1 260.0
## 2 260.5
## 3 261.0
## 4 261.5
## 5 262.0
## 6 262.5
## 7 263.0
## 8 263.5
## 9 264.0
## 10 264.5
## 11 265.0
## 12 265.5
## 13 266.0
## 14 266.5
## 15 267.0
## 16 267.5
## 17 268.0
## 18 268.5
## 19 269.0
## 20 269.5
## 21 270.0
## 22 270.5
## 23 271.0
## 24 271.5
## 25 272.0
## 26 272.5
## 27 273.0
## 28 273.5
## 29 274.0
## 30 274.5
## 31 275.0
## 32 275.5
## 33 276.0
## 34 276.5
## 35 277.0
## 36 277.5
## 37 278.0
## 38 278.5
## 39 279.0
## 40 279.5
## 41 280.0
## 42 280.5
## 43 281.0
## 44 281.5
## 45 282.0
## 46 282.5
## 47 283.0
## 48 283.5
## 49 284.0
## 50 284.5
## 51 285.0
## 52 285.5
## 53 286.0
## 54 286.5
## 55 287.0
## 56 287.5
## 57 288.0
## 58 288.5
## 59 289.0
## 60 289.5
## 61 290.0
## 62 290.5
## 63 291.0
## 64 291.5
## 65 292.0
## 66 292.5
## 67 293.0
## 68 293.5
## 69 294.0
## 70 294.5
## 71 295.0
## 72 295.5
## 73 296.0
## 74 296.5
## 75 297.0
## 76 297.5
## 77 298.0
## 78 298.5
## 79 299.0
## 80 299.5
## 81 300.0
## 82 300.5
## 83 301.0
## 84 301.5
## 85 302.0
## 86 302.5
## 87 303.0
## 88 303.5
## 89 304.0
## 90 304.5
## 91 305.0
## 92 305.5
## 93 306.0
## 94 306.5
## 95 307.0
## 96 307.5
## 97 308.0
## 98 308.5
## 99 309.0
## 100 309.5
## 101 310.0calculo de las nuevas probabilidades
probabilidades_nuevas <- predict(reg, datos_nuevos, type="response")
probabilidades_nuevas## 1 2 3 4 5 6 7
## 0.01649243 0.01748895 0.01854456 0.01966261 0.02084663 0.02210034 0.02342765
## 8 9 10 11 12 13 14
## 0.02483265 0.02631964 0.02789312 0.02955782 0.03131866 0.03318082 0.03514968
## 15 16 17 18 19 20 21
## 0.03723087 0.03943025 0.04175392 0.04420823 0.04679975 0.04953532 0.05242200
## 22 23 24 25 26 27 28
## 0.05546707 0.05867808 0.06206276 0.06562907 0.06938514 0.07333931 0.07750006
## 29 30 31 32 33 34 35
## 0.08187600 0.08647586 0.09130844 0.09638260 0.10170717 0.10729096 0.11314270
## 36 37 38 39 40 41 42
## 0.11927096 0.12568410 0.13239025 0.13939716 0.14671221 0.15434229 0.16229371
## 43 44 45 46 47 48 49
## 0.17057215 0.17918254 0.18812899 0.19741469 0.20704184 0.21701153 0.22732366
## 50 51 52 53 54 55 56
## 0.23797688 0.24896848 0.26029434 0.27194887 0.28392493 0.29621381 0.30880520
## 57 58 59 60 61 62 63
## 0.32168718 0.33484621 0.34826715 0.36193332 0.37582650 0.38992709 0.40421411
## 64 65 66 67 68 69 70
## 0.41866537 0.43325757 0.44796645 0.46276695 0.47763335 0.49253948 0.50745888
## 71 72 73 74 75 76 77
## 0.52236502 0.53723142 0.55203193 0.56674082 0.58133303 0.59578431 0.61007135
## 78 79 80 81 82 83 84
## 0.62417196 0.63806517 0.65173136 0.66515233 0.67831139 0.69119340 0.70378483
## 85 86 87 88 89 90 91
## 0.71607374 0.72804983 0.73970439 0.75103029 0.76202193 0.77267519 0.78298736
## 92 93 94 95 96 97 98
## 0.79295708 0.80258427 0.81187001 0.82081650 0.82942692 0.83770540 0.84565685
## 99 100 101
## 0.85328696 0.86060205 0.86760900Por defecto esta probabilidad se calcula como: $ log p_{i}/(1-p_{i}) $
colores <- NULL
colores[camarones$Exito == 0] <- "green"
colores[camarones$AlimentoDiario == 1] <- "red"
plot(camarones$AlimentoDiario, camarones$Exito, pch=21, bg=colores,
xlab="Alimento diario", ylab="Probabilidades de defectos")
legend('bottomleft', c('No defecto', 'Si defecto'), pch = 21, col = c('green', 'red'))
lines(datos_nuevos$AlimentoDiario, probabilidades_nuevas, col="blue", lwd=2)
Conclusión.
En conclusión, de acuerdo con lo analizado y a los datos que ingresamos en el modelo, podemos decir que la cantidad más adecuada para suministrar a los camarones diariamente es de 310 kilogramos de alimento ya que es la probabilidad más alta con un 86% aproximadamente. Sin embargo, esta cantidad puede ser un poco mayor, de está manera, la probabilidad lo sería igualmente.