Simulación de Datos

El siguiente código fue simulado con la ayuda de ChatGPT, su objetivo es generar un conjunto de datos sintético que modela información sobre el riesgo genético y el género de una muestra ficticia de individuos. Se generan 200 observaciones con dos variables explicativas, la primera definida como el genero al cual se le asigna 0 y 1, que representa mujer y hombre respectivamente, la segunda gen_risk sera una probabilidad entre 0 y 1 del nivel de riesgo genético de desarrollar cáncer.

# Sembrar una semilla
set.seed(123)

# Número de observaciones
n <- 200

# Variables explicativas
gen_risk <- runif(n, 0, 1) # riesgo genético (continua entre 0 y 1)
genero <- sample(c("1", "0"), n, replace = TRUE) # género como factor

# Convertir género a variable binaria para el modelo (0 = F, 1 = M)
genero_bin <- ifelse(genero == "M", 1, 0)

A partir de estas dos variables se calcula la probabilidad de tener cáncer usando la función logística. A partir de esto se genera la variable respuesta cancer, que tomara valores de 0 (no) o 1 (sí), simulando si la persona desarrolla o no cáncer.

# Función logística para la probabilidad de cáncer
logit_p <- -4 + 5 * gen_risk + 1.2 * genero_bin
prob_cancer <- 1 / (1 + exp(-logit_p))

# Generar variable dependiente
cancer <- rbinom(n, 1, prob_cancer)

# Crear data frame
datosca <- data.frame(cancer = cancer, gen_risk = gen_risk, genero = as.factor(genero))

Exploración de los Datos

Se verifica que la variable genero esté definida como un factor es decir, una variable categórica, lo cual sera importante para que el modelo sea interpretado correctamente,

library(dplyr)

# Asegurar que genero es factor
datosca <- datosca %>% mutate(genero = as.factor(datosca$genero))

# Ver primeras filas
head(datosca)
##   cancer  gen_risk genero
## 1      1 0.2875775      0
## 2      0 0.7883051      0
## 3      1 0.4089769      0
## 4      1 0.8830174      1
## 5      1 0.9404673      0
## 6      0 0.0455565      0

Tablas de Frecuencia

Se crean tablas que muestran la distribución entre la variable respuesta (cancer) y el género. Posteriormente se convierte la tabla de frecuencias absolutas (número de casos) en una tabla de frecuencias relativas donde se observaran porcentajes sobre el total.

tabla <- table(datosca$cancer, datosca$genero, dnn = c("cancer","genero")) 
addmargins(tabla)
##       genero
## cancer   0   1 Sum
##    0    68  78 146
##    1    24  30  54
##    Sum  92 108 200
tabla_frecuencias <- prop.table(tabla) * 100 
addmargins(tabla_frecuencias)
##       genero
## cancer   0   1 Sum
##    0    34  39  73
##    1    12  15  27
##    Sum  46  54 100

Visualización: Boxplot

A continuación se genera un boxplot que permite visualizar cómo varía el riesgo genético (gen_risk) según el diagnóstico de cáncer (cancer) y el género.

library(ggplot2)

ggplot(data = datosca, aes(x = factor(cancer), y = gen_risk, colour = genero)) + 
  geom_boxplot(outlier.shape = NA) + 
  theme(legend.position = "bottom") + 
  theme_bw() + 
  labs(x = "Cáncer (0 = No, 1 = Sí)", y = "Riesgo Genético", title = "Boxplot de Riesgo Genético por Cáncer y Género")

La relación entre gen_risk y cáncer es fuerte y visible es decir se puede deducir que personas que desarrollaron cáncer generalmente tiene un riesgo genético más alto. No hay datos atípicos, por eso no hay puntos.

Modelo de Regresión Logística

modelo_glm <- glm(cancer ~ gen_risk + genero, data = datosca, family = "binomial")
summary(modelo_glm)
## 
## Call:
## glm(formula = cancer ~ gen_risk + genero, family = "binomial", 
##     data = datosca)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.9335     0.6027  -6.527 6.71e-11 ***
## gen_risk      4.8583     0.8090   6.006 1.91e-09 ***
## genero1       0.2149     0.3695   0.582    0.561    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 233.30  on 199  degrees of freedom
## Residual deviance: 183.13  on 197  degrees of freedom
## AIC: 189.13
## 
## Number of Fisher Scoring iterations: 5

Se observa que gen_risk es el predictor más fuerte debido a su coeficiente positivo (4.86) y su valor p muy bajo indican que a medida que aumenta el riesgo genético, la probabilidad de tener cáncer crece exponencialmente. Mientras que genero no es significativo por su valor p = 0.561, lo que indica que una vez se controla el riesgo genético, el género no tiene un efecto sobre la probabilidad de cáncer en estos datos simulados.

Intervalos de Confianza

Otra manera de verificar cual es el predictor más fuerte sera observando los intervalos de confianza:

confint(modelo_glm) 
##                  2.5 %     97.5 %
## (Intercept) -5.2116208 -2.8359830
## gen_risk     3.3596161  6.5475371
## genero1     -0.5055583  0.9494584
confint.default(modelo_glm)  # Basado en el error estándar
##                  2.5 %     97.5 %
## (Intercept) -5.1146694 -2.7522972
## gen_risk     3.2728020  6.4438852
## genero1     -0.5092868  0.9391303

Donde el intervalo de gen_risk [3.36, 6.55] no contiene el cero, por lo tanto los valores posibles de coeficiente sera positivo y mayo a 3.36, lo que confirma que el riesgo genético es un predictor fuerte y estable.

Por otra parte el intervalo de genero [-0.5, 0.95] contiene el cero y el coeficiente podria ser positivo, negativo o nulo, por lo tanto no es estadísticamente significativo.

Visualización de la Curva de Probabilidad

datosca$cancer <- as.numeric(as.character(datosca$cancer)) 

# Nuevos valores interpolados
nuevos_gen_risk <- seq(from = min(datosca$gen_risk), to = max(datosca$gen_risk), length.out = 100)

# Para genero = 1 (hombres)
genero <- as.factor(rep("1", length(nuevos_gen_risk)))
predicciones <- predict(modelo_glm, newdata = data.frame(gen_risk = nuevos_gen_risk, genero = genero), type = "response")
datos_curva_M <- data.frame(gen_risk = nuevos_gen_risk, genero = genero, cancer = predicciones)

# Para genero = 0 (mujeres)
genero <- as.factor(rep("0", length(nuevos_gen_risk)))
predicciones <- predict(modelo_glm, newdata = data.frame(gen_risk = nuevos_gen_risk, genero = genero), type = "response")
datos_curva_F <- data.frame(gen_risk = nuevos_gen_risk, genero = genero, cancer = predicciones)

# Unir los dos
datos_curva <- rbind(datos_curva_F, datos_curva_M)

Gráfico de Resultados

ggplot(data = datosca, aes(x = gen_risk, y = as.numeric(cancer), color = genero)) + 
  geom_point() + 
  geom_line(data = datos_curva, aes(y = cancer)) +  
  theme_bw() + 
  labs(title = "P. cáncer en función del riesgo genético y género", y = "P(Cáncer)")

Se observa como a medida que aumenta el riesgo genético (gen_risk), la probabilidad de tener cáncer aumenta fuertemente para ambos géneros. Eso confirma que gen_risk es un predictor fuerte.

Tambien se evidencia como la línea azul (hombres) está apenas levemente por encima de la roja (mujeres), lo que muestra el coeficiente positivo pequeño y no significativo de genero en el modelo es decir no hace una gran diferencia en la predicción.

Los puntos en 1 (casos con cáncer) se concentran en valores altos de gen_risk, y los 0 (sin cáncer) en valores bajos por lo tanto el modelo está capturando bien la tendencia real.

Evaluacion del modelo

Likelihood Ratio

# Diferencia de residuos 

dif_residuos <- modelo_glm$null.deviance - modelo_glm$deviance 

# Grados libertad 
df <- modelo_glm$df.null - modelo_glm$df.residual 

# p-value 

p_value <- pchisq(q = dif_residuos,df = df, lower.tail = FALSE) 

paste("Diferencia de residuos:", round(dif_residuos, 4))
## [1] "Diferencia de residuos: 50.1763"
paste("Grados de libertad:", df) 
## [1] "Grados de libertad: 2"
paste("p-value:", round(p_value, 4)) 
## [1] "p-value: 0"

Reforzando los p-values encontrados en los coeficientes del modelo, es posible afirmar que el modelo en su conjunto es significativo, pues el P-value del mismo es 0.

Conclusiones

predicciones <- ifelse(test = modelo_glm$fitted.values > 0.5, yes = 1, no = 0) 
matriz_confusion <- table(modelo_glm$model$cancer, predicciones, 
dnn = c("observaciones", "predicciones")) 
matriz_confusion
##              predicciones
## observaciones   0   1
##             0 133  13
##             1  26  28

El modelo es capaz de clasificar correctamente el 80% de las observaciones, dado por la fórmula: \(\frac{133+28}{133+28+26+28}=0.8\)

Sin embargo, se observa que el modelo identificó a 28 pacientes con cáncer cuando en realidad solo había 26, lo cual indica que ha sobreestimado los casos positivos. Es decir, predijo más casos de cáncer de los que efectivamente se encontraban en los datos originales.

Por tanto, se descarta ajustar el umbral (threshold) del modelo, ya que el nivel actual proporciona una buena capacidad predictiva. En consecuencia, se puede concluir que el modelo se ajusta de manera adecuada para predecir los odds de padecer cáncer.