Introducción

El presente ejercicio tiene como objetivo analizar si la presencia de satélites (machos residiendo cerca) en hembras de cangrejo horseshoe depende de características de estas (color, estado de espinas, ancho y peso). Para ello, se aplicará y ajustará un modelo de regresión logística múltiple para establecer la probabilidad de que una cangrejo tenga satélites.

#carga de librerías
library(ggplot2)   # gráficos estáticos
library(table1)    # tablas descriptivas
## 
## Adjuntando el paquete: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
library(readxl)    # para leer Excel
library(caret)
## Cargando paquete requerido: lattice
## Registered S3 method overwritten by 'plyr':
##   method    from  
##   [.indexed table1
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Adjuntando el paquete: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
library(MASS)
library(pscl)
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
#carga de datos 

load("C:/Users/Juan Camilo Patiño/Documents/Universidad/Semestre XIII/Bioestadistica/Parcial 2 - Bioestadistica/crab.RData")

head(crab)
##   crab satellites weight width color spine y
## 1    1          8   3.05  28.3     2     3 1
## 2    2          0   1.55  22.5     3     3 0
## 3    3          9   2.30  26.0     1     1 1
## 4    4          0   2.10  24.8     3     3 0
## 5    5          4   2.60  26.0     3     3 1
## 6    6          0   2.10  23.8     2     3 0
# Ver estructura del dataset 

str(crab)
## 'data.frame':    173 obs. of  7 variables:
##  $ crab      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ satellites: int  8 0 9 0 4 0 0 0 0 0 ...
##  $ weight    : num  3.05 1.55 2.3 2.1 2.6 2.1 2.35 1.9 1.95 2.15 ...
##  $ width     : num  28.3 22.5 26 24.8 26 23.8 26.5 24.7 23.7 25.6 ...
##  $ color     : int  2 3 1 3 3 2 1 3 2 3 ...
##  $ spine     : int  3 3 1 3 3 3 1 2 1 3 ...
##  $ y         : int  1 0 1 0 1 0 0 0 0 0 ...

Descripción de las variables

Se cargó la base de datos que posee información sobre 173 individuos de cangrejos hembra, junto a las siguientes variables:

  • Satellites: el número de satélites observados en las hembras.

  • Weight: la variable numérica del peso de los cangrejos hembra.

  • Width: la variable númerica del ancho de cada cangrejo hembra.

  • Color: el color de las hembras, el cual es una variable ordinal (1=claro, 2=medio claro, 3=medio oscuro, 4=oscuro).

  • Spine: el estado de las dos espinas que tienen las hembras, el cual es una variable ordinal (1=ambas espinas en buenas condiciones, 2=una gastada o quebrada, 3=ambas gastadas o quebradas).

  • Y: donde 1 es si el número de satélites es mayor a 0, y donde 0 si el número de satélites es 0.

summary(crab)
##       crab       satellites         weight          width          color      
##  Min.   :  1   Min.   : 0.000   Min.   :1.200   Min.   :21.0   Min.   :1.000  
##  1st Qu.: 44   1st Qu.: 0.000   1st Qu.:2.000   1st Qu.:24.9   1st Qu.:2.000  
##  Median : 87   Median : 2.000   Median :2.350   Median :26.1   Median :2.000  
##  Mean   : 87   Mean   : 2.919   Mean   :2.437   Mean   :26.3   Mean   :2.439  
##  3rd Qu.:130   3rd Qu.: 5.000   3rd Qu.:2.850   3rd Qu.:27.7   3rd Qu.:3.000  
##  Max.   :173   Max.   :15.000   Max.   :5.200   Max.   :33.5   Max.   :4.000  
##      spine             y         
##  Min.   :1.000   Min.   :0.0000  
##  1st Qu.:2.000   1st Qu.:0.0000  
##  Median :3.000   Median :1.0000  
##  Mean   :2.486   Mean   :0.6416  
##  3rd Qu.:3.000   3rd Qu.:1.0000  
##  Max.   :3.000   Max.   :1.0000
# Crear variables para el análisis
satelites = crab$satellites
peso = crab$weight
ancho = crab$width
color = crab$color
espina = crab$spine
y = crab$y

# Generar tablas descriptivas
table1(~satelites|color, data=crab)
1
(N=12)		 2
(N=95)		 3
(N=44)		 4
(N=22)		 Overall
(N=173)
satelites
Mean (SD) 4.08 (3.12) 3.29 (3.21) 2.23 (2.60) 2.05 (3.62) 2.92 (3.15)
Median [Min, Max] 4.50 [0, 9.00] 3.00 [0, 15.0] 1.00 [0, 10.0] 0 [0, 12.0] 2.00 [0, 15.0] 
table1(~satelites|espina,data=crab)
1
(N=37)		 2
(N=15)		 3
(N=121)		 Overall
(N=173)
satelites
Mean (SD) 3.65 (3.39) 2.00 (2.36) 2.81 (3.13) 2.92 (3.15)
Median [Min, Max] 4.00 [0, 14.0] 0 [0, 6.00] 2.00 [0, 15.0] 2.00 [0, 15.0]

Interpretación

En la base de datos de estudio, se observa que se registraron más individuos de cangrejos hembra con presencia de satélites (n=111 vs. n=62). Sin embargo, a nivel descriptivo, con respecto a las variables ordinales, se observa la tendencia a que las hembras de color más claro y con ambas espinas en buenas condiciones, tienden a tener un mayor número de satélites a su alrededor.

Análisis de correlaciones

cor(crab[, c("satellites","weight","width","color","spine")])
##             satellites     weight      width      color       spine
## satellites  1.00000000  0.3692474  0.3398903 -0.1907846 -0.08993242
## weight      0.36924744  1.0000000  0.8868715 -0.2507772 -0.16648173
## width       0.33989033  0.8868715  1.0000000 -0.2643863 -0.12189458
## color      -0.19078455 -0.2507772 -0.2643863  1.0000000  0.37850163
## spine      -0.08993242 -0.1664817 -0.1218946  0.3785016  1.00000000

De acuerdo con el análisis de correlaciones, La variable de satélites muestra correlaciones positivas moderadas con peso (r = 0.37) y ancho (r = 0.34), lo que indica que las hembras de mayor tamaño o masa corporal tienden a presentar un mayor número de satélites. En contraste, las correlaciones negativas que presentan la variable de satélites con color podría indicar que las hembras más oscuras presentan un menor número de satélites. No obstante, al tener un poder de correlación muy bajo, no parece tener una relación directa.

Ajuste del modelo logística múltiple

# Ajustar un modelo logístico
mod = glm(y ~ peso + ancho + color + espina, data = crab, family = "binomial")
mod
## 
## Call:  glm(formula = y ~ peso + ancho + color + espina, family = "binomial", 
##     data = crab)
## 
## Coefficients:
## (Intercept)         peso        ancho        color       espina  
##     -7.5994       0.7949       0.2733      -0.5915       0.2717  
## 
## Degrees of Freedom: 172 Total (i.e. Null);  168 Residual
## Null Deviance:       225.8 
## Residual Deviance: 186.7     AIC: 196.7
summary(mod)
## 
## Call:
## glm(formula = y ~ peso + ancho + color + espina, family = "binomial", 
##     data = crab)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -7.5994     3.7542  -2.024   0.0429 *
## peso          0.7949     0.6917   1.149   0.2505  
## ancho         0.2733     0.1893   1.443   0.1489  
## color        -0.5915     0.2417  -2.447   0.0144 *
## espina        0.2717     0.2410   1.127   0.2597  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 225.76  on 172  degrees of freedom
## Residual deviance: 186.66  on 168  degrees of freedom
## AIC: 196.66
## 
## Number of Fisher Scoring iterations: 4

Al aplicar el modelo de regresión logística múltiple, se puede identificar que la variable color es la única presenta diferencias significativas, por lo que es la mayor predictora para determinar la presencia o ausencia de satélites en los cangrejos hembra. Por lo tanto, al tener un color más claro, las hembras tienden a presentar satélites a su alrededor.

#probabilidades = predict(mod, list(x1="y"), type = "response")
probabilidades = mod$fitted.values


presencia_satelites_modelo = probabilidades > 0.5
presencia_satelites_real = crab$y

# Matriz de confusión
table(presencia_satelites_real, presencia_satelites_modelo)
##                         presencia_satelites_modelo
## presencia_satelites_real FALSE TRUE
##                        0    31   31
##                        1    16   95
## Construcción e interpretación de métricas: desempeño, sensibilidad y especificidad

desempeño = (95 + 31) / (31 + 31 + 16 + 95)

sensibilidad = 95 / (95 + 16)

especificidad = 31 / (31 + 31)

cat("Desempeño =", round(desempeño, 3), "\n")
## Desempeño = 0.728
cat("Sensibilidad =", round(sensibilidad, 3), "\n")
## Sensibilidad = 0.856
cat("Especificidad =", round(especificidad, 3), "\n")
## Especificidad = 0.5
roc_obj <- roc(crab$y, fitted(mod))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, col = "blue", lwd = 2, main = "Curva ROC del modelo logístico")

auc(roc_obj)
## Area under the curve: 0.7701

Se sacó una matriz de confusión para observar qué tan efectivo es el modelo para predecir y clasificar correctamente la presencia o ausencia de satélites, de acuerdo con las variables estudiadas. Se concluye que El modelo logístico logra una clasificación adecuada de la presencia de satélites, mostrando buena sensibilidad (detecta correctamente la mayoría de hembras con satélites), aunque con una especificidad más limitada, lo que indica cierta confusión al identificar correctamente las hembras sin satélites.

Propuestas de modelo (Punto 3)

# 1) Poisson
mod_pois <- glm(y ~ ancho + peso + color + espina,
                data = crab, family = poisson)
summary(mod_pois)
## 
## Call:
## glm(formula = y ~ ancho + peso + color + espina, family = poisson, 
##     data = crab)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.81217    1.99419  -1.410    0.158
## ancho        0.08433    0.10102   0.835    0.404
## peso         0.15565    0.35756   0.435    0.663
## color       -0.22757    0.14252  -1.597    0.110
## espina       0.10830    0.12372   0.875    0.381
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 98.515  on 172  degrees of freedom
## Residual deviance: 85.964  on 168  degrees of freedom
## AIC: 317.96
## 
## Number of Fisher Scoring iterations: 5
# 2) quasi-Poisson (no tiene AIC válido)
mod_qpois <- glm(y ~ ancho + peso + color + espina,
                 data = crab, family = quasipoisson)
summary(mod_qpois)
## 
## Call:
## glm(formula = y ~ ancho + peso + color + espina, family = quasipoisson, 
##     data = crab)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -2.81217    1.19661  -2.350  0.01993 * 
## ancho        0.08433    0.06062   1.391  0.16600   
## peso         0.15565    0.21455   0.725  0.46918   
## color       -0.22757    0.08552  -2.661  0.00855 **
## espina       0.10830    0.07424   1.459  0.14650   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.3600604)
## 
##     Null deviance: 98.515  on 172  degrees of freedom
## Residual deviance: 85.964  on 168  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
# 3) NegBin (MASS)
mod_nb <- glm.nb(y ~ ancho + peso + color + espina, data = crab)
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
summary(mod_nb)
## 
## Call:
## glm.nb(formula = y ~ ancho + peso + color + espina, data = crab, 
##     init.theta = 27466.90059, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.81217    1.99422  -1.410    0.158
## ancho        0.08433    0.10102   0.835    0.404
## peso         0.15565    0.35757   0.435    0.663
## color       -0.22757    0.14252  -1.597    0.110
## espina       0.10830    0.12373   0.875    0.381
## 
## (Dispersion parameter for Negative Binomial(27466.9) family taken to be 1)
## 
##     Null deviance: 98.514  on 172  degrees of freedom
## Residual deviance: 85.963  on 168  degrees of freedom
## AIC: 319.97
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  27467 
##           Std. Err.:  280048 
## Warning while fitting theta: iteration limit reached 
## 
##  2 x log-likelihood:  -307.967
# 4) Zero-Inflated Poisson (pscl)
mod_zip <- zeroinfl(y ~ ancho + peso + color + espina, 
                    data = crab, dist = "poisson")
## Warning in value[[3L]](cond): sistema es computacionalmente singular: número de
## condición recíproco = 3.76407e-21FALSE
summary(mod_zip)
## 
## Call:
## zeroinfl(formula = y ~ ancho + peso + color + espina, data = crab, dist = "poisson")
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9715 -0.6705  0.1266  0.4429  1.1213 
## 
## Count model coefficients (poisson with log link):
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.31609         NA      NA       NA
## ancho        0.05218         NA      NA       NA
## peso         0.25122         NA      NA       NA
## color       -0.13791         NA      NA       NA
## espina       0.08448         NA      NA       NA
## 
## Zero-inflation model coefficients (binomial with logit link):
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   121.75         NA      NA       NA
## ancho         -86.31         NA      NA       NA
## peso          313.78         NA      NA       NA
## color         483.15         NA      NA       NA
## espina       -195.28         NA      NA       NA
## 
## Number of iterations in BFGS optimization: 85 
## Log-likelihood: -151.2 on 10 Df

Finalmente, tras recebir diferentes propuestas de modelos por la IA generativa, se concluye que la mejor propuesta es la del modelo de Quasi-Poisson. Esto debido que, por un lado, permite tener los datos necesarios para analizar el modelo, en comparación a los últimos. Y, por otro lado , presenta valores de error inferiores a los demás. De hecho, se reduce notablemente el error, teniendo de referencia al primer modelo de regresión logística múltiple.

Conclusiones

En suma, al aplicar un modelo de regresión logística, se pudo observar que el color es la variable que mejor predice que los cangrejos hembra tengan satélites o no. Además, se propone ajustar análisis con un modelo de Quasi-Poisson para optimizar la clasificación del modelo.