Regresión logística
En este documento se realizara un ajuste de modelado de regresión logística para un caso de estudio especifico, de esta manera podremos estimar la probabilidad de una variable binaria en función de una variable cuantitativa. Cabe destacar que la regresión lineal 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.
Caso de estudio: Camarones
EL caso de estudio que se analizara en este documento es el de una empresa sonorense que se dedica a producir camarón, en donde a los camarones se les proporciona cierta cantidad de comida con el objetivo de que estos lleguen a pesar 12 gramos, lo cual se concederá un éxito. El fin del estudio es determinar que cantidad de comida debe ser proporcionada a los camarones diariamente durante la ultima semana de crecimiento.
Imagen representativa de camarones
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)Utilizando datos reales de 12 estanques de camarones.
En la columna “AlimentoDiario” se ilustran los kilos de comida por día durante la última semana de alimentación al estanque, mientras que en la columna que se llama “éxito” se expone con 1 cuando los camarones llegan a pesar 12 gramos y pueden ser vendidos y con un 0 cuando no, el problema radica en que solamente 3 estanques llegaron a tener éxito.
Tabla de frecuencia de los datos
tabla <- table(camarones$Exito)
tabla##
## 0 1
## 9 3Con esto podemos apreciar que solo 3 estanques de 12 posibles llegaron a ser considerados como éxitos, para finalmente ser vendidos.
Ilustrando gráficamente
colores <- NULL
colores[camarones$Exito == 1] <- "green"
colores[camarones$Exito == 0] <- "red"
plot(camarones$AlimentoDiario, camarones$Exito, pch=21, bg=colores,
xlab="Cantidad de alimento diario", ylab="Probabilidades de éxito")
legend('bottomleft', c('Éxito', 'Sin éxito'), pch = 21, col = c('green', 'red'))
Hemos usado los argumentos pch y bg para mejorar la apariencia del gráfico. También hemos utilizado el comando legend para incluir una leyenda explicativa.
Visto de manera gráfica, y observando bien de los datos, se puede razonar que la cantidad de alimento diario en algunos casos llega a ser influyente, pero en otros no. Esto se puede deber a que quizás otros factores externos son los que impiden el crecimiento de los camarones y no necesariamente todo el peso radica en la cantidad de alimento que estos reciben. En esta práctica, vamos a ajustar un modelo de regresión logística para estudiar la posible relación.
Regresión logística
Para ver si el éxito verdaderamente es una variable que responde AlimentoDiario, se empleara el modelo de regresión logística, y así comprobaremos y rechazaremos nuestra hipótesis. Para ajustar el modelo se usa el comando glm (para modelos lineales generalizados) indicando que la respuesta es binomial mediante el argumento family:
regresion <- glm(Exito ~ AlimentoDiario, data = camarones, family=binomial)
summary(regresion)##
## 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 la tabla coeficiente, la temperatura es significativa con p valor= 0.185 el cual es mayor que el valor de significancia de 0.05. Lo anterior nos permite concluir que la hipótesis nula se acepta y se rechaza la hipótesis alternativa.
$ H_0:{1}=0 $ $ H_1:{1} $
Con respecto al intercepto: También es significativo dado que p valor= 0.175 > 0.05. Lo que significa que se acepta la hipótesis nula.
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 alimentos (entre 250 y 350 kilos): 1. Primero generamos una nueva secuencia de datos, el cual llamaremos datos nuevos
datos_nuevos <- data.frame(AlimentoDiario = seq(250, 350, 1))
datos_nuevos## AlimentoDiario
## 1 250
## 2 251
## 3 252
## 4 253
## 5 254
## 6 255
## 7 256
## 8 257
## 9 258
## 10 259
## 11 260
## 12 261
## 13 262
## 14 263
## 15 264
## 16 265
## 17 266
## 18 267
## 19 268
## 20 269
## 21 270
## 22 271
## 23 272
## 24 273
## 25 274
## 26 275
## 27 276
## 28 277
## 29 278
## 30 279
## 31 280
## 32 281
## 33 282
## 34 283
## 35 284
## 36 285
## 37 286
## 38 287
## 39 288
## 40 289
## 41 290
## 42 291
## 43 292
## 44 293
## 45 294
## 46 295
## 47 296
## 48 297
## 49 298
## 50 299
## 51 300
## 52 301
## 53 302
## 54 303
## 55 304
## 56 305
## 57 306
## 58 307
## 59 308
## 60 309
## 61 310
## 62 311
## 63 312
## 64 313
## 65 314
## 66 315
## 67 316
## 68 317
## 69 318
## 70 319
## 71 320
## 72 321
## 73 322
## 74 323
## 75 324
## 76 325
## 77 326
## 78 327
## 79 328
## 80 329
## 81 330
## 82 331
## 83 332
## 84 333
## 85 334
## 86 335
## 87 336
## 88 337
## 89 338
## 90 339
## 91 340
## 92 341
## 93 342
## 94 343
## 95 344
## 96 345
## 97 346
## 98 347
## 99 348
## 100 349
## 101 350calculo de las nuevas probabilidades
probabilidades_nuevas <- predict(regresion, datos_nuevos, type="response")
probabilidades_nuevas## 1 2 3 4 5 6
## 0.005057235 0.005694741 0.006412092 0.007219149 0.008126955 0.009147866
## 7 8 9 10 11 12
## 0.010295693 0.011585859 0.013035568 0.014663985 0.016492425 0.018544562
## 13 14 15 16 17 18
## 0.020846632 0.023427652 0.026319639 0.029557817 0.033180821 0.037230872
## 19 20 21 22 23 24
## 0.041753922 0.046799751 0.052421995 0.058678083 0.065629069 0.073339313
## 25 26 27 28 29 30
## 0.081876000 0.091308444 0.101707168 0.113142702 0.125684103 0.139397160
## 31 32 33 34 35 36
## 0.154342290 0.170572147 0.188128986 0.207041842 0.227323661 0.248968480
## 37 38 39 40 41 42
## 0.271948870 0.296213806 0.321687180 0.348267152 0.375826503 0.404214109
## 43 44 45 46 47 48
## 0.433257569 0.462766946 0.492539476 0.522365017 0.552031926 0.581333034
## 49 50 51 52 53 54
## 0.610071349 0.638065168 0.665152328 0.691193399 0.716073739 0.739704392
## 55 56 57 58 59 60
## 0.762021931 0.782987355 0.802584268 0.820816499 0.837705397 0.853286965
## 61 62 63 64 65 66
## 0.867609000 0.880728353 0.892708408 0.903616831 0.913523621 0.922499471
## 67 68 69 70 71 72
## 0.930614432 0.937936854 0.944532581 0.950464371 0.955791495 0.960569500
## 73 74 75 76 77 78
## 0.964850095 0.968681138 0.972106700 0.975167190 0.977899513 0.980337264
## 79 80 81 82 83 84
## 0.982510933 0.984448119 0.986173752 0.987710299 0.989077977 0.990294946
## 85 86 87 88 89 90
## 0.991377499 0.992340232 0.993196210 0.993957115 0.994633384 0.995234333
## 91 92 93 94 95 96
## 0.995768275 0.996242620 0.996663973 0.997038216 0.997370586 0.997665745
## 97 98 99 100 101
## 0.997927840 0.998160561 0.998367188 0.998550639 0.998713504Por defecto esta probabilidad se calcula como: $ log p_{i}/(1-p_{i}) $
colores[camarones$Exito == 1] <- "green"
colores[camarones$Exito == 0] <- "red"
plot(camarones$AlimentoDiario, camarones$Exito, pch=21, bg=colores,
xlab="Cantidad de alimento diario", ylab="Probabilidades de éxito")
legend('bottomleft', c('Éxito', 'Sin éxito'), pch = 21, col = c('green', 'red'))
lines(datos_nuevos$AlimentoDiario, probabilidades_nuevas, col="blue", lwd=1)
Conclusión
Tras analizar los datos nuevos insertados se puede concluir que las cantidades óptimas para proporcionarle a los camarones diariamente es de ronda entre los 315 y 350 kilos de alimento, debido a que según los datos existe una posibilidad mayor a partir a 90% a partir de los 313 kilos diarios dados en la última semana.