Proyecto de minería de datos: Modelos predictivos para la diabetes
Autor: Francisco José Rodríguez Cerezo
Diciembre 2022
Introducción y planteamiento del problema.
El NIDDK (instituto nacional de diabetes y enfermedades digestivas y renales de Estados Unidos), nos pide, aprovechando nuestros conocimientos en minería de datos, elaborar una herramienta capaz de predecir si un paciente es diabético o no.
Claramente nos enfrentamos a un problema de regresión ya que intentamos predecir el valor de una variable que desconocemos en función de otras variables.
Para medir el grado de cumplimiento de nuestro objetivo, ya que nos enfrentamos a un problema de regresión analítica, utilizaremos el concepto de error cuadrático medio.
Vamos a crear un modelo predictivo que permita identificar si una persona es diabética o no a través de otros factores más generales sobre la condición física del paciente.
Para ello, disponemos de un juego de datos tanto de pacientes diabéticos como no diabéticos proporcionado por la propia entidad (NIDDK).
Este juego de datos contiene, entre otras variables relacionadas con la enfermedad, el resultado real sobre si el paciente es diabético o no. Con ello podremos entrenar y probar nuestro modelo.
Verificación de datos.
Cargamos y presentamos un breve resumen del juego de datos:
# Cargamos el juego de datos como de costumbre.
archivo = 'diabetes.CSV'
datos_diabetes <- read.csv(archivo, row.names=NULL)
# Hacemos un summary de los mismos para ver un resumen del conjunto total.
summary(datos_diabetes)## Pregnancies Glucose BloodPressure SkinThickness
## Min. : 0.000 Min. : 0.0 Min. : 0.00 Min. : 0.00
## 1st Qu.: 1.000 1st Qu.: 99.0 1st Qu.: 62.00 1st Qu.: 0.00
## Median : 3.000 Median :117.0 Median : 72.00 Median :23.00
## Mean : 3.845 Mean :120.9 Mean : 69.11 Mean :20.54
## 3rd Qu.: 6.000 3rd Qu.:140.2 3rd Qu.: 80.00 3rd Qu.:32.00
## Max. :17.000 Max. :199.0 Max. :122.00 Max. :99.00
## Insulin BMI DiabetesPedigreeFunction Age
## Min. : 0.0 Min. : 0.00 Min. :0.0780 Min. :21.00
## 1st Qu.: 0.0 1st Qu.:27.30 1st Qu.:0.2437 1st Qu.:24.00
## Median : 30.5 Median :32.00 Median :0.3725 Median :29.00
## Mean : 79.8 Mean :31.99 Mean :0.4719 Mean :33.24
## 3rd Qu.:127.2 3rd Qu.:36.60 3rd Qu.:0.6262 3rd Qu.:41.00
## Max. :846.0 Max. :67.10 Max. :2.4200 Max. :81.00
## Outcome
## Min. :0.000
## 1st Qu.:0.000
## Median :0.000
## Mean :0.349
## 3rd Qu.:1.000
## Max. :1.000# Mostramos el número de filas del juego de datos para saber el número de pacientes totales.
nrow(datos_diabetes)## [1] 768Vemos que tenemos los datos de un total de 768 pacientes. Para cada uno tenemos una serie de variables que procedemos a explicar brevemente:
- Pregnancies: Número de veces que la persona ha estado embarazada.
- Glucose: Concentración de glucosa pasadas 2 horas tras realizar un test oral (mg/dL).
- BloodPressure: Presión sanguínea (mmHg).
- SkinThickness: Espesor de la piel (mm).
- Insulin: Insulina registrada al realizar el test del serúm de insulina de 2 horas (mu U/ml).
- BMI: Índice de masa corporal (kg/m^2).
- DiabetesPedigreeFunction: Función que indica la predisposición del paciente a sufrir diabetes en base a su historial médico y familiar.
- Age: Edad del paciente (años).
- Outcome: Variable que indica si el paciente es diabético (1) o no (0).
Mostremos ahora los primeros datos:
# Mostramos los primeros valores del juego de datos para ver que pinta tienen.
head(datos_diabetes)## Pregnancies Glucose BloodPressure SkinThickness Insulin BMI
## 1 6 148 72 35 0 33.6
## 2 1 85 66 29 0 26.6
## 3 8 183 64 0 0 23.3
## 4 1 89 66 23 94 28.1
## 5 0 137 40 35 168 43.1
## 6 5 116 74 0 0 25.6
## DiabetesPedigreeFunction Age Outcome
## 1 0.627 50 1
## 2 0.351 31 0
## 3 0.672 32 1
## 4 0.167 21 0
## 5 2.288 33 1
## 6 0.201 30 0En conjunto con el summary anterior, vemos de entrada que tenemos algunos valores perdidos done las variables Insulin o SkinThickness toman el valor 0, algo que no tiene sentido.
Debido a esto, vamos a gestionar los valores nulos o perdidos de nuestro juego datos:.
print("NA")## [1] "NA"# Sumamos las muestras con valor NA en cada columna.
colSums(is.na(datos_diabetes))## Pregnancies Glucose BloodPressure
## 0 0 0
## SkinThickness Insulin BMI
## 0 0 0
## DiabetesPedigreeFunction Age Outcome
## 0 0 0No tenemos datos etiquetados como “NA”.
print("En blanco")## [1] "En blanco"# Sumamos las muestras con valor vacío en cada columna.
colSums(datos_diabetes=="")## Pregnancies Glucose BloodPressure
## 0 0 0
## SkinThickness Insulin BMI
## 0 0 0
## DiabetesPedigreeFunction Age Outcome
## 0 0 0No tenemos datos en blanco.
print("0")## [1] "0"# Sumamos las muestras con valor igual a cero en cada columna.
colSums(datos_diabetes==0)## Pregnancies Glucose BloodPressure
## 111 5 35
## SkinThickness Insulin BMI
## 227 374 11
## DiabetesPedigreeFunction Age Outcome
## 0 0 500Pero si tenemos varios datos con valor 0. Obviamente para Pregnancies, DiabetesPedigreeFunction y Outcome es de esperar que tengamos valores iguales a 0. Pero para Glucose, BloodPressure, SkinThickness, Insulin y BMI no es normal tener valor cero.
En el caso de BloodPressure, Glucose y BMI son muy pocos valores nulos por lo que directamente vamos a prescindir de ellos:
# Creamos un nuevo juego de datos cuyos valores cumplen
# la condición mayor que cero para las variables que ya
# citamos anteriormente.
datos_diabetes_v2 <- datos_diabetes[datos_diabetes$Glucose > 0 & datos_diabetes$BloodPressure > 0 & datos_diabetes$BMI > 0, ]
# Mostramos de nuevo el resultado la cantidad de valores 0 para verificar el resultado.
colSums(datos_diabetes_v2==0)## Pregnancies Glucose BloodPressure
## 99 0 0
## SkinThickness Insulin BMI
## 192 332 0
## DiabetesPedigreeFunction Age Outcome
## 0 0 475El siguiente paso sería realizar técnicas de imputación para los valores que hemos considerado perdidos en SkinThickness e Insulin (ya que no tiene sentido un valor “0” en esas variables).
Este paso para imputar los datos faltantes, junto al tratamiento de posibles outliers, corresponden la fase de limpieza de datos que trataremos a continuación en la siguiente sección.
Preprocesado y limpieza de datos.
Valores perdidos:
Vamos a realizar una regresión lineal para imputar los datos faltantes en Insulin y SkinThickness. Para ello es importante que encontremos una variable con una correlación lo más alta posible frente a la que realizar el método de regresión lineal que elijamos.
Por este motivo, primero calcularemos la matriz de correlaciones y decidiremos frente a que variable realizar el ajuste en cada caso. Crearemos un juego de datos auxiliar que llamaremos datos_diabetes_corr donde eliminamos los valores iguales a cero en Insulin y SkinThickness y evitar así valores sesgados en las correlaciones.
# Creamos un juego de datos auxiliar donde eliminamos los valores igulaes
# a cero en Insulin y SkinThickness para tener un valor no sesgado de
# la correlación.
datos_diabetes_corr <- datos_diabetes_v2[datos_diabetes_v2$SkinThickness > 0 & datos_diabetes_v2$Insulin > 0, ]
# Comprobamos que lo hemos hecho bien.
colSums(datos_diabetes_corr==0)## Pregnancies Glucose BloodPressure
## 56 0 0
## SkinThickness Insulin BMI
## 0 0 0
## DiabetesPedigreeFunction Age Outcome
## 0 0 262# Importamos la librería necesaria para lo que queremos hacer.
if (!require('corrplot')) install.packages('corrplot'); library('corrplot')## Loading required package: corrplot## corrplot 0.92 loaded# Calculamos la matriz de correlaciones.
M = cor(datos_diabetes_corr)
# Mostramos la matiz de correlaciones en un gráfico.
corrplot(M, method = "number")
Recordemos que hemos creado un juego de datos auxiliar (datos_diabetes_corr). Vemos que:
- Insulin: La variable con la que tiene mayor correlación es Glucose (0.58).
- SkinThickness: La variable con la que tiene mayor correlación es BMI (0.66).
Vamos a realizar la regresión lineal entre Insulin y Glucose y, a partir de los coeficientes calculados por el ajuste lineal, imputaremos los valores faltantes en Insulin. Para ello, nos serviremos de la función impute_lm().
Primero preparamos un juego de datos auxiliar:
# Librería necesaria para lo que queremos hacer.
if (!require("simputation")) install.packages("simputation")## Loading required package: simputation# Creamos un subset auxiliar para trabajar esta sección.
datos_0 <- datos_diabetes_v2
# Tomamos los valores 0 de insulin y los transformamos a NA. Esto es necesario
# para poder aplicar la función impute_lm() después. De otro modo tomaría 0
# como el "valor real" de la medida.
datos_0$Insulin <- replace(datos_0$Insulin, which(datos_0$Insulin==0), NA)
# Repetimos para SkinThickness ya que también nos hará falta.
datos_0$SkinThickness <- replace(datos_0$SkinThickness, which(datos_0$SkinThickness==0), NA)Procedemos a realizar el ajuste lineal y la imputación de datos para Insulin frente a Glucose:
# Aplicamos la función "impute_lm()" para predecir e imputar los valores de Insulin.
prediccion_insulin <- impute_lm(datos_0, Insulin ~ Glucose)
# Aplicamos el método "lm()" sin imputar porque nos devolverá los coeficientes necesarios
# para pintar la recta resultante del ajuste lineal. La recta "teórica" de nuestra
# regresión.
recta_insulin <- lm(Insulin ~ Glucose, data=datos_0)
# Mostramos el ajuste realizado de forma gráfica.
# Primero los datos ya imputados.
plot(prediccion_insulin$Glucose, prediccion_insulin$Insulin, xlab="Glucose (mg/dL)", ylab="Insulin (mu U/ml)", main="Ajuste lineal Insulin vs Glucose", col="darkblue")
# Añadimos la recta resultante de los coeficientes calculados en la regresion.
abline(recta_insulin, col="red")
Vemos que la recta predicha se ajusta bastante bien al conjunto de valores (incluyendo las imputaciones realizadas con impute_lm()). Tenemos algunos outliers pero los trataremos más adelante.
Ahora repetiremos el proceso para SkinThickness frente a BMI:
# Aplicamos la función "impute_lm()" para predecir e imputar los valores de SkinThickness.
# Notar que ahora en "impute_lm()" usamos "prediccion_insuline" como conjunto de datos
# para tener ya incorporados los valores previamente imputados para Insulin.
prediccion_skin <- impute_lm(prediccion_insulin, SkinThickness ~ BMI)
# Aplicamos el método "lm()" sin imputar porque nos devolverá los coeficientes necesarios
# para pintar la recta resultante del ajuste lineal. La recta "teórica" de nuestra
# regresión.
recta_skin <- lm(SkinThickness ~ BMI, data=datos_0)
# Mostramos el ajuste realizado de forma gráfica.
# Primero los datos ya imputados.
plot(prediccion_skin$BMI, prediccion_skin$SkinThickness, xlab="BMI (kg/m^2)", ylab="SkinThickness (mm)", main="Ajuste lineal SkinThickness vs BMI", col="darkgreen")
# Añadimos la recta resultante de los coeficientes calculados en la regresion.
abline(recta_skin, col="red")
Nuevamente, los datos (tanto medidos realmente como imputados) parecen encajar bien en la recta predicha por el método de regresión lineal. En general, el grueso de los datos parece más disperso en torno a la recta que en el caso de Insulin (donde el grueso de los datos parecen estar más cerca a la recta predicha en su mayoría).
Sin embargo, visualmente parece que tenemos menor cantidad de outliers, o al menos, no destacan tanto como en el tramo final del ajuste hecho para Insulin.
Vamos a renombrar este último juego de datos para consolidar un poco lo hecho hasta ahora:
datos_diabetes_v3 <- prediccion_skinConcluyendo este apartado y a modo de resumen:
- Teníamos un juego de datos donde hemos prescindido de unas pocas muestras que presentaban valores perdidos en algunas variables.
- Hemos tenido un proceso de limpieza más largo para dos variables con muchas muestras con valores perdidos. Demasiadas para poder prescindir de esas muestras.
- Decidimos imputar los valores perdidos mediante un método de regresión lineal.
- Hemos obtenido un juego de datos (datos_diabetes_v3) sin valores nulos o perdidos.
Ahora es el turno de tratar con los outliers de nuestro juego de datos.
Outliers:
Vamos a mostrar un resumen de nuestro juego de datos.
summary(datos_diabetes_v3)## Pregnancies Glucose BloodPressure SkinThickness
## Min. : 0.000 Min. : 44.00 Min. : 24.0 Min. : 7.00
## 1st Qu.: 1.000 1st Qu.: 99.75 1st Qu.: 64.0 1st Qu.:22.00
## Median : 3.000 Median :117.00 Median : 72.0 Median :28.10
## Mean : 3.866 Mean :121.88 Mean : 72.4 Mean :28.76
## 3rd Qu.: 6.000 3rd Qu.:142.00 3rd Qu.: 80.0 3rd Qu.:35.00
## Max. :17.000 Max. :199.00 Max. :122.0 Max. :99.00
## Insulin BMI DiabetesPedigreeFunction Age
## Min. :-19.93 Min. :18.20 Min. :0.0780 Min. :21.00
## 1st Qu.: 89.94 1st Qu.:27.50 1st Qu.:0.2450 1st Qu.:24.00
## Median :132.27 Median :32.40 Median :0.3790 Median :29.00
## Mean :154.39 Mean :32.47 Mean :0.4748 Mean :33.35
## 3rd Qu.:191.25 3rd Qu.:36.60 3rd Qu.:0.6275 3rd Qu.:41.00
## Max. :846.00 Max. :67.10 Max. :2.4200 Max. :81.00
## Outcome
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.3439
## 3rd Qu.:1.0000
## Max. :1.0000Aspectos a destacar que no tienen lógica y trataremos como outliers:
- Insulin tiene, al menos, un valor negativo.
- El valor máximo de Insulin es 846, un valor excesivamente alto e imposible en una persona. Además, la media de insulina es 154.39 que es aproximadamente 5 veces menos que ese valor máximo.
- El valor máximo de BMI es 67.10 (algo muy exagerado y me atrevería a decir imposible), mientras que la media se mantiene en 32.40 (menos de la mitad de ese valor máximo).
- SkinThickness también tiene un valor máximo exagerado (casi 10cm de grosor de piel), que se desvía demasiado del valor promedio.
En el resto de variables no haremos tratamiento de outliers ya que no presentan anomalías en este sumario de los datos. Resumiendo, haremos tratamiento de outliers en: Insulin, BMI, y SkinThickness
Veamos un plot rápido de estas variables:
plot(datos_diabetes_v3[c(4, 5, 6)])
Vemos claramente varios valores bastante alejados del grueso de datos en todos los gráficos, eso son los outliers que trataremos de eliminar o, por lo menos, reducir.
Veamos también una representación de los datos en cada variable como histogramas:
# Librería necesaria.
if (!require('ggplot2')) install.packages('ggplot2'); library('ggplot2')## Loading required package: ggplot2# Histograma para la variable "Insulin".
pl <- ggplot(datos_diabetes_v3, aes(x=Insulin))
pl + geom_histogram(binwidth = 18, col = "black", fill = "lightgreen") + ggtitle("Distribución de los niveles de insulina de los pacientes")
Vemos que sólo tenemos un valor negativo y el grueso de los datos está, aproximadamente, entre 75 y 300.
# Histograma para la variable "BMI".
pl <- ggplot(datos_diabetes_v3, aes(x=BMI))
pl + geom_histogram(binwidth = 1, col = "black", fill = "orange") + ggtitle("Distribución del índice de masa corporal de los pacientes")
En este caso, vemos un perfil de campana de gauss entre 20 y 45. Es decir la distribución normal, donde se supone que debemos tener la mayoría de nuestros datos.
# Histograma para la variable "SkinThickness".
pl <- ggplot(datos_diabetes_v3, aes(x=SkinThickness))
pl + geom_histogram(binwidth = 1.8, col = "black", fill = "lightblue") + ggtitle("Distribución del espesor de piel de los pacientes")
Por último para el espesor de la piel, el rango de valores típico se situa entre 14 y 38 aproximadamente.
Vamos a emplear el método del rango intercuantílico (como ya hicimos en la PEC2) para acotar los valores de estas variables y eliminar así los outliers. Recordemos brevemente como es el método:
- Tomamos la diferencia entre el tercer y primer cuantil de nuestros datos, esto es el rango intercuantílico: \[IQR=Q_{3}-Q{1}\]
- Aquellos valores que disten 1.5 veces este valor (IQR) por encima de Q3 o por debajo de Q1 se consideran outliers y quedan fuera de nuestra selección final de datos. Es decir: \[\forall x\ \epsilon[Q_{1}-1.5IQR,\ Q_{3}+1.5IQR] \rightarrow x\ no\ es\ outlier\]
Veamos ese valor en cada caso:
- Insulin: \(IQR=191.25-89.94=102.01\rightarrow x\ \epsilon[-63.21,\ 344.40]\) en este caso tomaremos \(x\ \epsilon[0,\ 344.40]\).
- BMI: \(IQR=36.60-27.50=9.10\rightarrow x\ \epsilon[13.85,\ 50.25]\)
- SkinThickness: \(IQR=35-22=13\rightarrow x\ \epsilon[2.50,\ 54.50]\)
Crearemos un nuevo juego de datos excluyendo aquellos valores que no estén dentro de estos rangos:
# Añadimos las condiciones necesarias para tener sólo aquellas muestras que con valores dentro
# de los rangos ya calculados para dichas variables
datos_diabetes_v4 <- datos_diabetes_v3[
datos_diabetes_v3$Insulin > 0 & datos_diabetes_v3$Insulin < 344.40 &
datos_diabetes_v3$BMI > 13.85 & datos_diabetes_v3$BMI < 50.25 &
datos_diabetes_v3$SkinThickness > 2.50 & datos_diabetes_v3$SkinThickness < 54.50
, ]
# Resumen de datos.
summary(datos_diabetes_v4)## Pregnancies Glucose BloodPressure SkinThickness
## Min. : 0.00 Min. : 56.0 Min. : 24.00 Min. : 7.00
## 1st Qu.: 1.00 1st Qu.: 99.0 1st Qu.: 64.00 1st Qu.:21.97
## Median : 3.00 Median :115.0 Median : 72.00 Median :28.00
## Mean : 3.91 Mean :120.4 Mean : 72.37 Mean :28.30
## 3rd Qu.: 6.00 3rd Qu.:139.0 3rd Qu.: 80.00 3rd Qu.:34.00
## Max. :17.00 Max. :199.0 Max. :122.00 Max. :54.00
## Insulin BMI DiabetesPedigreeFunction Age
## Min. : 9.167 Min. :18.20 Min. :0.0780 Min. :21.00
## 1st Qu.: 88.000 1st Qu.:27.40 1st Qu.:0.2450 1st Qu.:24.00
## Median :130.031 Median :32.00 Median :0.3725 Median :29.00
## Mean :142.120 Mean :32.12 Mean :0.4661 Mean :33.38
## 3rd Qu.:183.937 3rd Qu.:36.30 3rd Qu.:0.6138 3rd Qu.:41.00
## Max. :342.000 Max. :50.00 Max. :2.2880 Max. :81.00
## Outcome
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.3304
## 3rd Qu.:1.0000
## Max. :1.0000# Número de muestras.
nrow(datos_diabetes_v4)## [1] 690Tras el tratamiento de valores perdidos o nulos, y el tratamiento de outliers nos hemos quedado con 690 muestras de las 768 que teníamos inicialmente. Lo que supone una pérdida de aproximadamente un 11% del número de muestras.
Este valor de pérdida de información con respecto al conjunto inicial nos parece tolerable ya que estamos aún por encima del mínimo requerido de muestras que se nos indica en el enunciado de esta práctica (PRA1).
Por último volvemos a mostrar los gráficos e histogramas para comparar el resultado antes y después del tratamiento de outliers.
plot(datos_diabetes_v4[c(4, 5, 6)])
Vemos que obtenemos una mayor definición de los puntos al tener los datos más agrupados (más cerca entre sí), ya que de este modo no es necesario alejar tanto el “zoom” en las gráficas para visualizar todas las muestras.
Pasemos ahora a los histogramas:
# Histograma para la variable "Insulin".
pl <- ggplot(datos_diabetes_v4, aes(x=Insulin))
pl + geom_histogram(binwidth = 8, col = "black", fill = "lightgreen") + ggtitle("Distribución de los niveles de insulina de los pacientes sin outliers")
# Histograma para la variable "BMI".
pl <- ggplot(datos_diabetes_v4, aes(x=BMI))
pl + geom_histogram(binwidth = 0.8, col = "black", fill = "orange") + ggtitle("Distribución del índice de masa corporal de los pacientes sin outliers")
# Histograma para la variable "SkinThickness".
pl <- ggplot(datos_diabetes_v4, aes(x=SkinThickness))
pl + geom_histogram(binwidth = 1, col = "black", fill = "lightblue") + ggtitle("Distribución del espesor de piel de los pacientes sin outliers")
Comparando estos histogramas con sus homólogos sin tratar los outliers vemos que, en todos los casos, las envolventes de los histogramas se ajustan mucho mejor a una campana de gauss. Esto es buena señal ya que abarcando casi la totalidad de valores posibles en cada caso, las muestras se situan dentro de la distribución normal.
Con esto damos por finalizada la tarea de limpieza de datos, donde tenemos datos_diabetes_v4 tras haber:
- Limpiado de valores nulos y perdidos los datos originales.
- Limpiado la mayoría de outliers del juego de datos original.
Por último, en este apartado de preprocesamiento de datos nos queda codificar aquellas variables susceptibles de ello.
Codificación y categorización de variables:
En el enunciado de la práctica PRA1 se nos indica un mínimo de dos variables categóricas. Dado que la edad (Age) y el índice de masa corporal (BMI) son medidas cualitativas sobre la condición de las personas, vamos a categorizar en distintos rangos y codificar dichas variables creando columnas nuevas en cada caso.
Empezamos con la variable Age, teniendo en cuenta nuestro juego de datos (edad mínima 21) y según index mundi podemos dividir las edades como:
- Juventud: entre 21 y 24 años. Estos individuos los meteremos en la categoría Young y lo codificaremos con el número 1.
- Adultez: entre 25 y 54 años. Estos individuos los meteremos en la categoría Adult y lo codificaremos con el número 2.
- Madurez: entre 55 y 64 años. Estos individuos los meteremos en la categoría Mature y lo codificaremos con el número 3.
- Vejez: más de 65 años. Estos individuos los meteremos en la categoría Elder y lo codificaremos con el número 4.
# Creamos otro juego de datos idéntico al que ya teníamos para trabajar sin modificar el anterior,
# de este modo podemos deshacer nuestros pasos en caso de ser necesario.
datos_diabetes_v5 <- datos_diabetes_v4
# Creamos un vector AgeCategory copiando la columna Age del juego de datos.
AgeCategory <- datos_diabetes_v5$Age
# Sustituimos los valores con las categorías definidas anteriormente.
AgeCategory <- replace(AgeCategory, which(AgeCategory >= 21 & AgeCategory <= 24), "Young")
AgeCategory <- replace(AgeCategory, which(AgeCategory >= 25 & AgeCategory <= 54), "Adult")
AgeCategory <- replace(AgeCategory, which(AgeCategory >= 55 & AgeCategory <= 64), "Mature")
AgeCategory <- replace(AgeCategory, which(AgeCategory >= 65 & AgeCategory <= 99), "Elder")
# Comprobamos el cambio.
AgeCategory[0:10]## [1] "Adult" "Adult" "Adult" "Young" "Adult" "Adult" "Adult" "Adult"
## [9] "Adult" "Mature"# Ahora creamos el vector AgeCode como una copia de AgeCategory.
AgeCode <- AgeCategory
# Codificamos los valores como hemos dicho anteriormente.
# Ahora para AgeCode
AgeCode <- replace(AgeCode, which(AgeCode == "Young"), 1)
AgeCode <- replace(AgeCode, which(AgeCode == "Adult"), 2)
AgeCode <- replace(AgeCode, which(AgeCode == "Mature"), 3)
AgeCode <- replace(AgeCode, which(AgeCode == "Elder"), 4)
# Comprobamos el cambio.
AgeCode[0:10]## [1] "2" "2" "2" "1" "2" "2" "2" "2" "2" "3"Para el caso de BMI según el CDC (Centro de control y prevensión enfermedades de estados unidos), podemos categorizar el peso de las personas con los siguientes rangos:
- Bajo peso: Por debajo de 18.5. Estos individuos los meteremos en la categoría Under weight y lo codificaremos con el número 1.
- Normal: Entre 18.5 y 24.9. Estos individuos los meteremos en la categoría Normal y lo codificaremos con el número 2.
- Sobrepeso: Entre 25.0 y 29.9. Estos individuos los meteremos en la categoría Over weight y lo codificaremos con el número 3.
- Obesidad: Más de 30. Estos individuos los meteremos en la categoría Obesity y lo codificaremos con el número 4.
# Creamos un vector AgeCategory copiando la columna BMI del juego de datos.
BMICategory <- datos_diabetes_v5$BMI
# Sustituimos los valores con las categorías definidas anteriormente.
BMICategory <- replace(BMICategory, which(BMICategory >= 10 & BMICategory <= 18.4), "Under weight")
BMICategory <- replace(BMICategory, which(BMICategory >= 18.5 & BMICategory <= 24.9), "Normal")
BMICategory <- replace(BMICategory, which(BMICategory >= 25 & BMICategory <= 29.9), "Over weight")
BMICategory <- replace(BMICategory, which(BMICategory >= 30 & BMICategory <= 60), "Obesity")
# Comprobamos el cambio.
BMICategory[0:10]## [1] "Obesity" "Over weight" "Normal" "Over weight" "Obesity"
## [6] "Over weight" "Obesity" "Obesity" "Obesity" "Over weight"# Ahora creamos el vector AgeCode como una copia de AgeCategory.
BMICode <- BMICategory
# Codificamos los valores como hemos dicho anteriormente.
# Ahora para AgeCode
BMICode <- replace(BMICode, which(BMICode == "Under weight"), 1)
BMICode <- replace(BMICode, which(BMICode == "Normal"), 2)
BMICode <- replace(BMICode, which(BMICode == "Over weight"), 3)
BMICode <- replace(BMICode, which(BMICode == "Obesity"), 4)
#Comprobamos el cambio.
BMICode[0:10]## [1] "4" "3" "2" "3" "4" "3" "4" "4" "4" "3"Añadimos las nuevas columnas al juego de datos y lo renombramos como datos_diabetes_clean:
# Añadimos estos vectores como columnas a nuestro juego de datos.
datos_diabetes_v5$AgeCategory <- AgeCategory
datos_diabetes_v5$AgeCode <- AgeCode
datos_diabetes_v5$BMICategory <- BMICategory
datos_diabetes_v5$BMICode <- BMICode
# Renombramos el juego de datos.
datos_diabetes_clean <- datos_diabetes_v5
# Comprobamos el resultado.
colnames(datos_diabetes_clean)## [1] "Pregnancies" "Glucose"
## [3] "BloodPressure" "SkinThickness"
## [5] "Insulin" "BMI"
## [7] "DiabetesPedigreeFunction" "Age"
## [9] "Outcome" "AgeCategory"
## [11] "AgeCode" "BMICategory"
## [13] "BMICode"Veamos para finalizar un histograma más intuitivo en cuanto a la edad y el peso utilizando las nuevas variables que hemos creado categorizando sus homónimas:
# Histograma para la variable "BMICategory".
pl <- ggplot(datos_diabetes_clean, aes(x=BMICategory))
pl + geom_histogram(stat ="count") + stat_count(bindwidth = 1, col = "black", fill = "orange") + ggtitle("Distribución los individuos en función del peso")## Warning in geom_histogram(stat = "count"): Ignoring unknown parameters:
## `binwidth`, `bins`, and `pad`## Warning in stat_count(bindwidth = 1, col = "black", fill = "orange"): Ignoring
## unknown parameters: `bindwidth`
# Histograma para la variable "AgeCategory".
pl <- ggplot(datos_diabetes_clean, aes(x=AgeCategory))
pl + geom_histogram(stat ="count") + stat_count(bindwidth = 1, col = "black", fill = "pink") + ggtitle("Distribución los individuos en función de la edad")## Warning in geom_histogram(stat = "count"): Ignoring unknown parameters:
## `binwidth`, `bins`, and `pad`## Warning in stat_count(bindwidth = 1, col = "black", fill = "pink"): Ignoring
## unknown parameters: `bindwidth`
Vemos que la mayoría de pacientes son adultos (entre 25 y 54 años) que padecen obesidad (índice de masa corporal mayor que 30).
Estos resultados parecen bastante razonables ya que en la edad adulta se empiezan a desarrollar este tipo de enfermedades y, lógicamente, habrá mucha más población adulta que aquellos en edad madura o vejez debido a la esperanza de vida de las personas.
En cuanto a la obesidad, según el portal TFAH (Trust for America’s Health), la tasa de obesidad en adultos de EE.UU supera el 42%.
Conclusiones del preprocesamiento de datos:
Resumiendo, hemos tratado nuestro juego de datos original de la siguiente manera:
- Hemos eliminado muy pocas muestras por culpa de valores perdidos en BloodPressure, Glucose, y BMI.
- Para el grueso de valores perdidos en nuestro juego de datos, hemos realizado la imputación de datos en Insulin y SkinThickness mediante el método de regresión lineal.
- Tras inspeccionar el conjunto de datos, hemos tratado outliers en las variables SkinThickness, Insulin y BMI con el método del rango intercuantílico.
- Hemos considerado las variables Age y BMI como medidas cualitativas de los individuos y las hemos categorizado y codificado.
Tras todo este proceso, hemos creado el juego de datos datos_diabetes_clean que será el que usaremos de ahora en adelante.
# Resumen de datos.
summary(datos_diabetes_clean)## Pregnancies Glucose BloodPressure SkinThickness
## Min. : 0.00 Min. : 56.0 Min. : 24.00 Min. : 7.00
## 1st Qu.: 1.00 1st Qu.: 99.0 1st Qu.: 64.00 1st Qu.:21.97
## Median : 3.00 Median :115.0 Median : 72.00 Median :28.00
## Mean : 3.91 Mean :120.4 Mean : 72.37 Mean :28.30
## 3rd Qu.: 6.00 3rd Qu.:139.0 3rd Qu.: 80.00 3rd Qu.:34.00
## Max. :17.00 Max. :199.0 Max. :122.00 Max. :54.00
## Insulin BMI DiabetesPedigreeFunction Age
## Min. : 9.167 Min. :18.20 Min. :0.0780 Min. :21.00
## 1st Qu.: 88.000 1st Qu.:27.40 1st Qu.:0.2450 1st Qu.:24.00
## Median :130.031 Median :32.00 Median :0.3725 Median :29.00
## Mean :142.120 Mean :32.12 Mean :0.4661 Mean :33.38
## 3rd Qu.:183.937 3rd Qu.:36.30 3rd Qu.:0.6138 3rd Qu.:41.00
## Max. :342.000 Max. :50.00 Max. :2.2880 Max. :81.00
## Outcome AgeCategory AgeCode BMICategory
## Min. :0.0000 Length:690 Length:690 Length:690
## 1st Qu.:0.0000 Class :character Class :character Class :character
## Median :0.0000 Mode :character Mode :character Mode :character
## Mean :0.3304
## 3rd Qu.:1.0000
## Max. :1.0000
## BMICode
## Length:690
## Class :character
## Mode :character
##
##
## # Número de muestras.
nrow(datos_diabetes_clean)## [1] 690Vamos a exportar los datos ya preprocesados a un fichero .csv con la función write.csv():
# Coma como separador y punto como separador decimal
write.csv(datos_diabetes_clean, "diabetes_clean.csv")Estudio PCA.
Introducción:
Una forma de entender el PCA (Principal Component Analysis) es gracias al concepto de de autovectores y autovalores del álgebra lineal (Eigenvalues and eigenvectors). Matemáticamente se pueden describir las componentes principales como una rotación de los ejes del sistema de coordenadas de modo que los nuevos ejes maximizan la variabilidad de los datos para tener mayor información representativa pudiendo reducir dimensionalidad al problema (eliminando la información redundante de variables según la correlación con el objetivo).
Estos nuevos ejes rotados respecto a la base original, son los autovectores, que por definición, cuando constituyen una nueva base ortogonal, son unitarios (longitud 1). Al multiplicarse por su correspondiente autovalor, el autovector recuperan su tamaño original, este autovalor da una idea de la variabilidad retenida por cada componente principal.
Este paso de emplear las componentes principales, pese al incoveniente de interpretar el significado de estas nuevas dimensiones para el juego de datos original, muchas veces acelera los procesos de clasificación.
Estudio:
Primero vamos a crear un juego de datos auxiliar que sólo recoja aquellas variables sobre las que tenga sentido realizar PCA:
# Recogemos Glucose, BloodPressure, SkinThickness, Insulin y DiabetesPedigreeFunction.
datos_PCA <- datos_diabetes_clean[ , c(2:5, 7)]
# Comprobamos.
colnames(datos_PCA)## [1] "Glucose" "BloodPressure"
## [3] "SkinThickness" "Insulin"
## [5] "DiabetesPedigreeFunction"nrow(datos_PCA)## [1] 690Hemos seleccionado estas variables ya que son numéricas que funcionan como indicadores de la diabetes.
Vamos a instalar una nueva librería llamada FactoMiner para realizar el estudio PCA:
if (!require('FactoMineR')) install.packages('FactoMineR'); library(FactoMineR)## Loading required package: FactoMineRUtilizamos la función PCA de esta librería:
# Hacemos el análisis de componentes principales.
pca <- PCA(datos_PCA, graph = FALSE)
# Lo mostramos para ver que nos ofrece.
pca## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 690 individuals, described by 5 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"Veamos brevemente los resultados de los autovalores pca$eig:
head(pca$eig)## eigenvalue percentage of variance cumulative percentage of variance
## comp 1 2.0191698 40.383396 40.38340
## comp 2 1.0175080 20.350159 60.73356
## comp 3 0.9772242 19.544484 80.27804
## comp 4 0.7548647 15.097294 95.37533
## comp 5 0.2312333 4.624666 100.00000Que de forma visual y mucho más intuitiva:
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBafviz_screeplot(pca, addlabels = TRUE)
Esto nos viene a decir el porcentaje de la varianza de los datos que explica cada componente principal. Vemos que las dos primeras son las componentes más importantes ya que entre ambas se explica el 80% de la varianza de los datos.
Otra forma visual de analizar la situación, además de ver la correlación entre componentes, es mediante un gráfico circular con fviz_pca_var:
fviz_pca_var(pca, col.var = "cos2",
geom.var = "arrow",
labelsize = 1,
repel = TRUE)
Donde las flechas que representan cada componente principal, estarán más juntas entre sí, cuanto mayor sea la correlación entre ambas. De este modo vemos que las dos primeras componentes principales tienen una correlación muy alta y entre ambas explican un 80% de la varianza de los datos.
Por otro lado las componentes principales 3 y 4 también están fuertemente correlacionadas entre sí aunque explican un porcentaje menor de la varianza de los datos (sólo el 34%).
Por último, la quinta componente principal no guarda correlación alguna con las demás y sólo explica un 4.6% de la varianza de los datos, por lo que carece importancia y podríamos obviarla.
Conclusiones del estudio PCA:
En este breve estudio sobre PCA hemos introducido una definición matemática de lo que son las componentes principales y lo que representan en un conjunto de datos.
A continuación hemos realizado el PCA y visto los resultados bajo el enfoque de los autovalores de dichas componentes.
Hemos visto como las dos primeras componentes principales explican la mayoría de la varianza de los datos y que, además, están fuertemente correlacionadas.
Modelo no supervisado basado en el concepto de distancia: k-means.
Dada la naturalez de nuestro problema, emplearemos el método k-means de modo que agrupe las muestras con 2 clústers, los diabéticos y los que no lo son. Pero antes deberemos elegir que variables emplear en nuestro modelo. Para ello veamos brevemente la matriz de correlación del juego de datos.
# Creamos un juego de datos auxiliar en el que sólo tengamos las variables numéricas.
datos_corr <- datos_diabetes_clean[ , c(1:9)]
# Importamos la librería necesaria.
if (!require('corrplot')) install.packages('corrplot'); library('corrplot')
# Creamos la matriz de correlaciones.
M = cor(datos_corr)
# Mostramos gráficamente la matriz de correlaciones.
corrplot(M, method = "number")
Como podemos ver, las variables con mayor correlación con Outcome son: Glucose e Insulin. Así que haremos la clasificación de k-means con estas dos variables. Primero cargaremos la librería necesaria y prepararemos dos conjuntos de datos auxiliares idénticos, uno para cada métrica.
# Cargamos la librería necesaria y le decimos que la instale en caso de no tenerla.
if (!require('cluster')) install.packages('cluster'); library(cluster)## Loading required package: cluster# Creamos un juego de datos auxiliar solo con las variables que nos interesan.
datos_k <- datos_corr[ , c(2, 5)]
datos_k_v2 <- datos_corr[ , c(2, 5)]Realizamos el clustering con k-means con la distancia euclidiana.
# Realizamos el método con k=2 que es el número de clusters apropiado para nuestra tarea
# ya que el outcome sólo puede tener dos valores: 1 (es diabético) o 0 (no es diabético).
fit <- kmeans(datos_k, 2)
y_cluster <- fit$cluster
# Pintamos el resultado.
clusplot(datos_k, fit$cluster, color=TRUE, shade=TRUE, labels=2, lines=0, main = "Clasificación k-means con k=2 (distancia euclidiana)", xlab = "Glucose (mg/dL)", ylab = "Insulin (mu U/ml)")
Vemos que los dos clusters explican la variabilidad del 100% de los datos, es decir, no existen muestras sin clasificar. Lo más llamativo de la clasificación es la línea fuertemente marcada que atraviesa ambos clusters de lado a lado. Esto se debe a que la correlación entre las variables Glucose e Insulin es de 0.77,la cuale s la correlación má alta en todo el juego de datos (exceptuando la diagonal de la matriz obviamente).
Esto tiene lógica, ya que la insulina es una hormona que segrega el cuerpo humano como respuesta a la glucosa para poder procesarla de forma efectiva, por lo que existe una relación muy marcada entre ambas variables, seguramente lineal por lo que podemos ver a priori.
Aquellas muestras más alejadas de esta línea podríamos deducir que se trata de casos fuera de lo normal en nuestro conjunto de pacientes, es decir, aquellos en los que la respuesta de insulina ante sus niveles de glucosa en sangre se aleja del promedio. Esto no quiere decir necesariamente que sean persona diabéticas, ya que en ambos clusters tenemos muestras alejadas de esta marcada tendencia, por lo que podemos deducir que la respuesta en insulina frente a los niveles de glucosa en sangre pueden ser “anómalos” sin ser necesariamente diabéticos.
En cuanto a la distancia entre los clusters vemos que existe una pequeña zona donde se superponen ambos, por lo que tenemos un porcentaje de muestras que no podríamos clasificar con seguridad si pertenecen al grupo de pacientes diabéticos o no.
El método kmeans de la librería cluster utiliza por defecto la distancia euclidiana. Mediante el método daisy recalcularemos los puntos del conjunto de datos indicando el uso de la métrica correspondiente para la distancia de Manhattan.
# Empleamos el método daisy para recalcular los valores bajo el concepto de ditancia Manhattan.
distance <- daisy(datos_k_v2, metric = "manhattan")
# Repetimos el proceso de usar el método kmeans pero pasamos como datos el conjunto "distance" que acabamos de crear con el método "daisy".
fit_v2 <- kmeans(distance, 2, algorithm = "Forgy")
y_cluster_v2 <- fit_v2$cluster
clusplot(datos_k_v2, fit_v2$cluster, color=TRUE, shade=TRUE, labels=2, lines=0, main = "Clasificación k-means con k=2 (distancia de manhattan)", xlab = "Glucose (mg/dL)", ylab = "Insulin (mu U/ml)")
El resultado prácticamente no cambia de usar la métrica euclidiana a la métrica con el concepto de distancia Manhattan. Esto puede deberse, bien a que hemos hecho algo mal en nuestro código, o bien a que la correlación entre ambas variables está tan marcada que realmente no influye el concepto de distancia que usemos.
Segundo modelo no supervisado: métodos DBSCAN y OPTICS.
A continuación vamos a aplicar los métodos DBSCAN y OPTICS para realizar el clustering con nuestro juego de datos.
# Cargamos la librería necesaria.
if (!require('dbscan')) install.packages('dbscan'); library('dbscan')## Loading required package: dbscan##
## Attaching package: 'dbscan'## The following object is masked from 'package:stats':
##
## as.dendrogramPara determinar un criterio de vecindad (minPts), veamos los valores máximo y mínimo de cada variable.
datos_dbscan <- datos_corr[ , c(1:9)]
summary(datos_dbscan)## Pregnancies Glucose BloodPressure SkinThickness
## Min. : 0.00 Min. : 56.0 Min. : 24.00 Min. : 7.00
## 1st Qu.: 1.00 1st Qu.: 99.0 1st Qu.: 64.00 1st Qu.:21.97
## Median : 3.00 Median :115.0 Median : 72.00 Median :28.00
## Mean : 3.91 Mean :120.4 Mean : 72.37 Mean :28.30
## 3rd Qu.: 6.00 3rd Qu.:139.0 3rd Qu.: 80.00 3rd Qu.:34.00
## Max. :17.00 Max. :199.0 Max. :122.00 Max. :54.00
## Insulin BMI DiabetesPedigreeFunction Age
## Min. : 9.167 Min. :18.20 Min. :0.0780 Min. :21.00
## 1st Qu.: 88.000 1st Qu.:27.40 1st Qu.:0.2450 1st Qu.:24.00
## Median :130.031 Median :32.00 Median :0.3725 Median :29.00
## Mean :142.120 Mean :32.12 Mean :0.4661 Mean :33.38
## 3rd Qu.:183.937 3rd Qu.:36.30 3rd Qu.:0.6138 3rd Qu.:41.00
## Max. :342.000 Max. :50.00 Max. :2.2880 Max. :81.00
## Outcome
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.3304
## 3rd Qu.:1.0000
## Max. :1.0000La escala de las variables será la diferencia entre su valor máximo y mínimo:
- Pregnancies: 17 - 0 = 17
- Glucose: 199 - 56 = 143
- BloodPressure: 122 - 124 = 98
- SkinThickness: 54 - 7 = 47
- Insulin: 342 - 9.16 = 332.84
- BMI: 50 - 18.20 = 31.80
- DiabetesPedigreeFunction: 2.28 - 0.07 = 2.21
- Age: 81 - 21 = 60
A pesar de esto, debemos mantener cierta lógica en nuestro planteamiento además de una homogeneidad en los calculos para poder comparar los resultados. Aludiendo nuevamente a la matriz de correlaciones, realmente las variables más correlacionadas con Outcome son Glucose e Insulin. Prepararemos un juego de datos auxiliar en el que sólo tengamos estas dos variables.
datos_dbscan <- datos_corr[ , c(2, 5)]
summary(datos_dbscan)## Glucose Insulin
## Min. : 56.0 Min. : 9.167
## 1st Qu.: 99.0 1st Qu.: 88.000
## Median :115.0 Median :130.031
## Mean :120.4 Mean :142.120
## 3rd Qu.:139.0 3rd Qu.:183.937
## Max. :199.0 Max. :342.000La variable Glucose tiene un rango de variabilidad menor por lo que será la variable más sensible a pequeños cambios en el valor de vecindad (minPts). Dividiremos a la mitad (71.5) la longitud de su intervalo y nos acogeremos a este criterio para especificar el valor de minPts.
Ahora procedemos a ordenar el conjunto de datos por vecindad con este criterio (minPts = 71.5).
# Usamos el método optics para ordenar agrupar las muestras por vecindad con el criterio establecido de minPts.
res <- optics(datos_dbscan, minPts = 71.5)
res## OPTICS ordering/clustering for 690 objects.
## Parameters: minPts = 71.5, eps = 98.5525838231362, eps_cl = NA, xi = NA
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi# Ordenación de las muestas.
res$order## [1] 1 679 674 667 651 650 637 689 687 680 664 655 629 627 649 676 681 666
## [19] 658 648 620 617 599 597 582 571 570 567 559 554 551 549 538 537 528 521
## [37] 587 517 498 477 474 472 441 394 393 390 344 338 297 293 289 286 279 266
## [55] 253 245 243 229 228 223 222 198 195 184 181 179 175 170 152 141 132 129
## [73] 127 115 111 108 78 75 74 70 69 59 58 38 26 13 8 514 246 214
## [91] 168 596 565 529 493 422 413 377 362 300 247 230 136 56 663 639 526 499
## [109] 425 365 352 321 238 154 44 33 283 491 343 112 6 669 603 575 558 523
## [127] 438 434 403 363 326 239 203 147 145 72 57 17 552 465 395 53 25 688
## [145] 647 584 470 643 631 613 497 416 389 370 109 106 177 518 295 121 114 660
## [163] 623 439 373 341 250 196 21 508 276 23 270 665 690 659 653 641 608 605
## [181] 588 583 573 569 566 564 556 524 506 505 469 458 457 456 448 397 356 346
## [199] 337 333 330 324 310 306 280 233 216 215 187 140 383 123 95 85 68 62
## [217] 45 29 14 4 15 624 685 503 500 485 483 468 452 447 572 454 446 442
## [235] 437 412 406 396 386 376 375 336 334 314 309 304 236 226 197 173 162 159
## [253] 134 119 93 76 71 34 2 7 259 644 475 372 294 281 257 87 329 278
## [271] 400 349 153 531 495 678 473 414 307 586 453 420 151 19 609 591 522 436
## [289] 429 354 327 305 291 241 235 221 208 117 94 30 431 189 492 578 398 340
## [307] 298 328 155 163 113 626 331 555 541 539 480 325 287 104 632 275 79 581
## [325] 254 192 86 622 614 561 546 510 501 467 461 381 361 251 91 64 384 171
## [343] 128 657 616 590 576 563 562 530 512 488 466 451 427 374 366 353 316 568
## [361] 519 369 273 260 237 224 218 193 169 161 148 126 124 110 80 54 37 5
## [379] 244 445 303 550 272 63 149 401 82 671 560 504 484 317 296 256 240 225
## [397] 201 142 122 120 107 66 55 28 618 511 476 464 460 421 405 157 105 430
## [415] 392 619 462 101 574 494 313 182 156 81 20 449 482 423 419 357 311 73
## [433] 77 642 545 339 675 585 543 411 268 258 209 164 100 92 60 50 471 604
## [451] 359 242 176 96 662 635 98 200 180 32 207 261 387 385 217 10 410 43
## [469] 165 219 267 202 516 301 290 646 595 231 137 426 668 417 686 424 415 399
## [487] 265 24 46 628 589 540 205 513 409 342 654 183 234 206 645 615 600 525
## [505] 479 418 382 371 274 160 135 31 611 670 598 577 455 444 388 378 322 249
## [523] 190 150 133 103 52 42 22 263 534 358 347 509 61 496 264 35 102 535
## [541] 47 308 271 602 11 520 269 174 144 90 661 139 167 255 88 579 84 683
## [559] 507 97 408 67 532 285 481 12 677 553 99 312 640 547 211 601 277 143
## [577] 116 527 607 407 335 16 621 594 350 320 199 118 625 172 40 580 368 83
## [595] 673 612 592 262 252 204 186 27 380 402 89 39 459 656 478 131 634 213
## [613] 315 636 672 487 486 533 432 158 440 348 299 630 360 9 684 544 194 210
## [631] 489 433 443 65 536 232 288 652 404 318 450 610 351 302 130 227 379 292
## [649] 212 191 41 345 36 638 248 146 502 282 185 557 3 542 323 220 188 48
## [667] 490 391 138 51 682 548 435 355 284 166 367 515 463 428 332 178 606 364
## [685] 319 18 593 633 125 49Con esto haremos un diagrama de alcanzabilidad (reachability plot) en el que podremos ver de forma visual la distancia de alcanzabilidad de cada punto. En este tipo de gráficos los valles representan los clusters (cuanto más profundos, más densos y agrupados), mientras que los picos representan puntos lejanos a los clusters, los cuales son candidatos a outliers.
plot(res)
Vemos un gráfico con un crecimiento exponencial hacia el último cuarto de la totalidad de las muestras. También es destacable el pico inicial. Esto puede deberse a que sólo tenemos dos clusters (diabéticos y no diabéticos), visualmente vimos con el método k-means que habia una fuerte tendencia a lo largo de una línea bien definida de las muestras a través de ambos clusters.
Seguramente el grueso de este gráfico corresponde a las muestras agrupadas en esa línea (aproximadamente al rango entre 0 y 500). Al estar tan marcada esa agrupación de las muestras en forma de línea, cualquier muestra que no caiga dentro de ella en comparación al resto, estará “muy” lejos del grueso de datos, por lo que las distancias, relativamente hablando, se magnifican.
Esa tendencia linear que vimos en k-means puede que se deba a las técnicas de imputación que hicimos en la primera parte del proyecto con la regresión linea, ya que, tras consultar el informe de la primera parte del proyecto, inicialmente Glucose e Insulin tenian una correlación de 0.58 y tras las labores de limpieza e imputación tenemos una correlación de 0.77.
En este gráfico podríamos decir que en el rango de 0 a 500 tenemos un cluster y de 500 a 690 otro cluster distinto. Además, el pico inicial puede indicar que, a pesar del preprocesado realizado en la primera parte del trabajo, quedan algunos outliers.
Vamos ha aplicar el método DBSCAN. En cuanto al parámetro eps_cl, usaremos la longitud del intervalo de posibles valores para Glucose, que recordemos, es 144.
res <- extractDBSCAN(res, eps_cl = 144)
res## OPTICS ordering/clustering for 690 objects.
## Parameters: minPts = 71.5, eps = 98.5525838231362, eps_cl = 144, xi = NA
## The clustering contains 1 cluster(s) and 0 noise points.
##
## 1
## 690
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusterHemos probado con múltiples valores pero sólo somos capaces de conseguir un único cluster.
plot(res)
Realmente este resultado nos ha dejado un poco bloqueados. Quizá el preprocesado de datos fue muy agresivo o quizá el juego de datos no ha sido la elección más adecuada. También pueda deberse ha algún fallo en la metodología del que no somos conscientes. Tras revisar exhaustivamente el código no somos capaces de detectar nuestro error.
Probamos a utilizar todas las variables numéricas del juego de datos.
# Esta vez el intervalo de menor longitud era 17 para la variable "Pregnancies" cuya mitad es 8.5 opara establecer el valor de "minPts".
res2 <- optics(datos_corr[ , c(1:8)], minPts = 8.5)
plot(res2)
Ahora si somos capaces de distinguir algunos valles, pero igualmente el resultado es igual de malo, ya que sólo deberíamos tener dos clusters.
# Establecemos 17 como valor de "eps_cl" ya que es la longitud más corta del rango de valores de nuestro juego de datos.
res2 <- extractDBSCAN(res2, eps_cl = 17)
res2## OPTICS ordering/clustering for 690 objects.
## Parameters: minPts = 8.5, eps = 58.8305808997328, eps_cl = 17, xi = NA
## The clustering contains 2 cluster(s) and 301 noise points.
##
## 0 1 2
## 301 383 6
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusterEn este caso, el algoritmo detecta dos clusters pero con uno de ellos muchísimo más pequeño que el otro y una gran cantidad de ruido, puede ser que el método detecte como outliers todas aquellas muestras fuera de la línea que vimos con el método k-means.
plot(res2)
Esto se contradice con los resultados en k-means ya que vimos el método explicaba el 100% de la variabilidad de los datos, es decir, todas las muestras quedaban clasificadas ya sea en un cluster, en el otro, o en la superposición de ambas. Esto nos hace pensar que nuestra teoría de que la fuerte agrupación lineal que vimos este influyendo en como el algoritmo identifica, erróneamente, las muestras fuera de esa agrupación como ruido ya que esa agrupación tan peculiar hace que las distancias relativas entre puntos se vean magnificadas.
Probemos de nuevo con un valor distinto para eps_cl.
res3 <- optics(datos_corr[ , c(1:8)], minPts = 8.5)
res3 <- extractDBSCAN(res3, eps_cl = 10)
res3## OPTICS ordering/clustering for 690 objects.
## Parameters: minPts = 8.5, eps = 58.8305808997328, eps_cl = 10, xi = NA
## The clustering contains 2 cluster(s) and 657 noise points.
##
## 0 1 2
## 657 20 13
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusterTras varias pruebas, el resultado más interesante lo encontramos en torno a 10 para eps_cl. Si el valor eps_cl aumenta mucho, pasa a identificar un único cluster y si disminuye mucho más de 10, detecta todo como ruido. En este caso tenemos dos clusters aunque muy pequeños los dos, al menos las cifras parecen estar equilibradas entre ambos. De todos modos el resultado sigue siendo malo ya que casi la totalidad de las muestras las considera ruido.
plot(res3)
Seguimos igual de desconcertados que antes, tras volver a revisar el código no somos capaces de detectar nuestro error, por lo que muy posiblemente el error que hemos cometido sea o bien conceptual, o bien lo vengamos arrastrando desde el preprocesado de datos.
Otras posibilidades, aunque más remotas, sean las limitaciones del propio juego de datos o bien lo ya explicado con la peculiar agrupación linear entre Glucose e Insulin ya que entre ambas se realizó un ajuste linear para imputar los valores perdidos en el preprocesado de datos.
Dado que no entendemos estos resultados del todo, y es muy probable que estemos cometiendo algún error. En resumen, descartaríamos el uso de estos resultados.
Modelo supervisado: Árboles de decisión.
Definiendo el árbol de decisión:
Vamos a proceder con los modelos supervisados, en este caso crearemos distintos árboles de decisión. El primer paso es el de preparar los subconjuntos de entrenamiento y prueba para nuestro modelo. Usaremos una regla típica de asignar 2/3 del total de las muestras al subconjunto de entrenamiento y reservamos el tercio restanto para el subconjunto de prueba.
Pero primero necesitamos transformar nuestros datos a factor.
datos_arbol <- datos_corr[ , c(1:9)]
datos_arbol$Outcome <- as.factor(datos_arbol$Outcome)Debemos asegurarnos que la división de estos subconjuntos no sea desigual para no inducir sesgos en nuestro modelo, lo mejor en este caso es hacerlo aleatoriamente.
# Fijamos una semilla para poder reproducir los mismos resultados.
set.seed(23)
# Creamos un parámetro con el controlar la división de los subconjuntos de forma dinámica.
split_prop <- 3
# Generamos aleatoriamente los índices de las muestras para el conjunto de entrenamiento.
indices = sample(1:nrow(datos_arbol), size=floor(((split_prop-1)/split_prop)*nrow(datos_arbol)))
# Creamos los subconjuntos de datos (entrenamiento y prueba).
conjunto_entrenamiento <- datos_arbol[indices, ]
conjunto_prueba <- datos_arbol [-indices, ]Comprobemos la distribución de personas diabéticas en ambos subconjuntos.
table(conjunto_entrenamiento$Outcome)##
## 0 1
## 306 154table(conjunto_prueba$Outcome)##
## 0 1
## 156 74- Conjunto de entrenamiento: Tenemos 154 diabéticos frente a 306 personas no diabéticas, esto supone una proporción de: 154/306=0.503… O lo que es lo mismo, un porcentaje aproximado del 50% de personas diabéticas.
- Conjunto de prueba: En este caso tenemos 74 personas diabéticas frente a 156 que no lo son, esto supone una porporción de: 74/156=0.474… O lo que es lo mismo, un porcentaje aproximado del 47% de personas diabéticas.
La diferencia en los ratios de positivos frente a negativos en diabétes para ambos subconjuntos sólo difiere en un 3%, lo cual consideramos una diferencia tolerable.
Pasemos a la creación en sí de los árboles de decisión. Para ello necesitamos la librería C50.
if (!require('C50')) install.packages('C50'); library('C50')## Loading required package: C50Creamos nuestro primer árbol de decisión.
# Variable "target" (outcome).
trainy <- conjunto_entrenamiento$Outcome
# Resto de variables.
trainx <- conjunto_entrenamiento[, -c(9)]
# Modelo del primer árbol.
# Indicamos el parámetro rules = TRUE para poder ver las reglas de decisión.
arbol_1 <- C5.0(trainx, trainy, rules=TRUE)
summary(arbol_1)##
## Call:
## C5.0.default(x = trainx, y = trainy, rules = TRUE)
##
##
## C5.0 [Release 2.07 GPL Edition] Fri Feb 16 09:34:18 2024
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 460 cases (9 attributes) from undefined.data
##
## Rules:
##
## Rule 1: (149/8, lift 1.4)
## Glucose <= 127
## Insulin <= 106
## -> class 0 [0.940]
##
## Rule 2: (14, lift 1.4)
## Pregnancies <= 6
## Glucose <= 127
## BloodPressure > 88
## -> class 0 [0.938]
##
## Rule 3: (211/20, lift 1.4)
## Pregnancies <= 6
## Glucose <= 127
## Age <= 34
## -> class 0 [0.901]
##
## Rule 4: (181/31, lift 1.2)
## BMI <= 30.2
## -> class 0 [0.825]
##
## Rule 5: (16/3, lift 2.3)
## Pregnancies <= 6
## Glucose <= 127
## BloodPressure <= 88
## Insulin > 106
## Insulin <= 185
## BMI > 26.9
## Age > 34
## -> class 1 [0.778]
##
## Rule 6: (63/14, lift 2.3)
## Pregnancies > 6
## Insulin > 106
## BMI > 26.9
## -> class 1 [0.769]
##
## Rule 7: (113/28, lift 2.2)
## Glucose > 127
## BMI > 30.2
## -> class 1 [0.748]
##
## Rule 8: (104/35, lift 2.0)
## Insulin > 106
## BMI > 26.9
## Age > 34
## -> class 1 [0.660]
##
## Default class: 0
##
##
## Evaluation on training data (460 cases):
##
## Rules
## ----------------
## No Errors
##
## 8 84(18.3%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 264 42 (a): class 0
## 42 112 (b): class 1
##
##
## Attribute usage:
##
## 82.17% Glucose
## 70.65% BMI
## 68.48% Age
## 64.35% Pregnancies
## 57.17% Insulin
## 6.52% BloodPressure
##
##
## Time: 0.0 secsNuestro árbol ha clasificado mal 84 muestras, lo cual supone una tasa de error del 18.3%. También vemos que los errores de etiquetado no están sesgados, es decir, tenemos el mismo número de falsos positivos como de falsos negativos.
En cuanto al uso de las variables, la más empleada ha sido Glucose que se ha empleado en la mayoría de los casos (82.17%), seguida de BMI con un 70.65%, Age con un 68.48%, Pregnancies con 64.35%, Insulin con 57.17% y finalmente BloodPressure con un 6.52% de uso. Por lo que vemos, el uso de las variables para la clasificación de las muestras está repartido de una forma bastante uniforme a excepción del caso de BloodPressure.
Podemos prescindir de la variable DiabetesPedigreeFunction ya que no aparece en el summary de nuestro árbol de decisión ni en ninguna de las reglas del mismo.
Vamos con las reglas que ha generado nuestro árbol de decisión:
- Regla 1: Si Glucose es menor o igual a 127 e Insulin es menor o iguala 106, el paciente se etiqueta como no diabético.
- Regla 2: Si Pregnancies es menor o igual a 6, Glucose es menor o igual a 127 y BloodPressure es mayor a 88, el paciente se clasifica como no diabético.
- Regla 3: Si Pregnancies es menor o igual a 6, Glucose es menor o igual a 127 y Age es menor o igual a 34, el paciente se clasifica como no diabético.
- Regla 4: Si BMI es menor o igual a 30.2, el paciente se clasifica como no diabético.
- Regla 5: Si Pregnancies es menor o igual a 6, Glucose es menor o igual a 127 y BloodPressure es menor o igual a 88, Insulin es mayor a 106 pero menor o igual a 185, BMI es mayor a 26.9 y Age es mayor a 34, el paciente se clasifica como diabético.
- Regla 6: Si Pregnancies es mayor a 6, Insulin es mayor a 106 y BMI es mayor a 26.9, el paciente se considera diabético.
- Regla 7: Si Glucose es mayor a 127 y BMI es mayor a 30.2, el pasciente se considera diabético.
- Regla 8: Si Insulin es mayor a 106, BMI es mayor a 26.9 y Age es mayor a 34, el paciente se etiqueta como diabético.
Podriamos resumir este conjunto de reglas en que el umbral de Glucose a partir del cual un paciente es sospechoso de padecer diabetes es a partir del valor 127, de forma similar, el modelo considera que una persona no diabética debería tener valores de Insulin entre 106 y 185. El valor 26.9 BMI también es crítico así como la edad (Age) que a partir de los 34 el riesgo de diabetes aumenta. El número de embarazos también presenta una clara división en 6, a partir de ahí, la persona es más propensa a ser diabética.
Vamos a mostrar nuestro árbol de forma visual.
# Vamos a crear un nuevo trainy, en el que simplemente cambiamos la nomenclatura de los valores
# como "Si" para "Seguro" y "No" para "inseguro".
# Esto nos ayudará a visualizar mejor el gráfico.
trainy_arbol_1 <- trainy
trainy_arbol_1 <- as.factor(trainy_arbol_1)
# Recreamos el árbol, necesitamos hacerlo de nuevo con "rules=FALSE" (valor por defecto)
# para que nos dejé crear y mostrar el gráfico.
arbol_1_plot <- C5.0(trainx, trainy_arbol_1)
# Creamos y mostramos el árbol esquemáticamente.
plot(arbol_1_plot)
Ahora si, es momento de poner a prueba nuestro árbol de decisión.
Validación del árbol de decisión:
Vamos a utilizar el método predict con el conjunto de datos de prueba para ver como se comporta nuestro modelo.
# Variable "target" (outcome).
testy <- conjunto_prueba$Outcome
# Resto de variables.
testx <- conjunto_prueba[, -c(9)]
# Aplicamos el método "predict" para que el model haga sus predicciones.
prediccion_arbol_1 <- predict(arbol_1, testx, type = "class")Calculamos la matriz de confución
matriz_confusion <- table(testy, prediccion_arbol_1)
matriz_confusion## prediccion_arbol_1
## testy 0 1
## 0 123 33
## 1 20 54Veamos el porcentaje de acierto, es decir, la precisión de nuestro modelo.
precision_arbol_1 <- 100 * (sum(diag(matriz_confusion)) / sum(matriz_confusion))
print(sprintf("El %% de registros correctamente clasificados es: %.2f %%", precision_arbol_1))## [1] "El % de registros correctamente clasificados es: 76.96 %"Podemos ver que la precisión de nuestro modelo no es muy mala, ya que ha acertado, aproximadamente, en el 77% de los casos. Además, si nos fijamos en la matriz de confusión, vemos que ha dado más falsos positivos (33) que falsos negativos (20). Con lo cual, a pesar de que nuestro modelo está ligeramente sesgado, no es una mala noticia ya que en este caso, es preferible un falso positivo que pueda llevar a realizar pruebas más exhaustivas al potencial paciente diabético para determinar en que grado si es realmente es diabético, que el caso contrario, ya que si fuera al revés podría haber pacientes que si son diabéticos y ni ellos ni sus médicos lo saben y no se cuidan o tratan al respecto.
Segundo modelo supervisado: Regresión lineal.
Para crear nuestro modelo de regresión lineal, vamos a emplear el método lm(). Pero primero prepararemos un conjuto de datos auxiliar.
datos_regresion <- datos_corr[ , c(1:9)]Creamos el modelo, haremos una regresión de Outcome frente a Glucose y, en lugar de hacer otra distinta de Outcome frente a Insulin, la haremos frente a BMI. Ya que en el preprocesado de datos de la primera parte de este proyecto hemos hecho una regresión lineal entre la glucosa y la insulina, no tiene mucho sentido que hagamos ahora una regresión de Outcome frente a ambas variables que ya de por sí vienen fuertemente correlacionadas entre ellas.
Hemos elegido hacer la segunda regresión frente BMI ya que en la primera parte del trabajo hicimos una regresión de SkinThickness frente a BMI que ambas estaban correlacionadas entre sí pero también lo estaban con Outcome aunque en menor medida que en el caso de Glucose e Insulin.
Hemos elegido hacer los ajustes lineales con Glucose y BMI en lugar de hacerlos frente a Insulin y **SkinThickness* ya que creemos que la medida de la glucosa o el índice de masa corporal son medidas más cotidianas, directas y sencillas de obtener en una situación real.
regresion_glucose <- lm(datos_regresion$Outcome ~ datos_regresion$Glucose, data = datos_regresion)
regresion_bmi <- lm(datos_regresion$Outcome ~ datos_regresion$BMI, data = datos_regresion)Ahora sólo nos queda hacer predicciones las cuales las haremos con el método impute_lm().
prediccion_glucose <- impute_lm(datos_regresion, Outcome ~ Glucose)
prediccion_bmi <- impute_lm(datos_regresion, Outcome ~ BMI)Mostramos el resultado gráficamente.
plot(prediccion_glucose$Glucose, prediccion_glucose$Outcome, xlab="Glucose (mg/dL)", ylabel="Outcome", main="Regresión lineal Outcome vs Glucose", col="darkblue")
Hemos tenido que llegar hasta aquí, para darnos cuenta de que no tiene mucho sentido usar una regresión lineal para tratar de predecir una variable binaria, que en esencia no es más que una variable categórica (Outcome vale 1 o 0 nada más).
Esto se debe, a que una variable categórica no es una variable continua y esto, a su vez, significa que puede no tener una distribución normal.
En lugar de continuar por este sendero, buscaremos una alternativa de otro método supervisado para nuestro problema.
Alternativa como segundo modelo supervisado: Regresión logística.
Tras el problema con el que nos hemos topado en la sección anterior, tras una investigación en internet, decidimos probar a realizar una regresión logística.
Una regresión logística es una alternativa más adecuada para predecir una variable binaria como es nuestro caso. Este modelo estadístico se utiliza para predecir una variable binaria en función de uno o más predictores de forma similar a como lo hace la regresión lineal. La diferencia entre la regresión logística y la regresión lineal, es que la regresión logística no asume una distribución normal y se ajusta mejor a los datos categóricos.
Volvemos a crear nuestro juego de datos auxiliar.
datos_regresion <- datos_corr[ , c(1:9)]Vamos a usar el método glm() que se usa para realizar este tipo de regresión. Primero haremos el estudio de Outcome frente a Glucose. Para esta labor tamibén hemos consultado con este enlace de la UCLA.
# Empleamos el método "glm()" apoyándonos en la documentación, agregar el argumento "family" que hace alusión al tipo de distribución estadística de los datos, como "Outcome" es una variable que solo tiene dos posbiles valores, empleamos la distribución "binomial".
regresion_logistica_g <- glm(datos_regresion$Outcome ~ datos_regresion$Glucose, data = datos_regresion, family = "binomial")Ahora vamos a evaluar nuestro modelo, para ello camos a visualizar la curva ROC (Receiver Operating Characteristics) y su área bajo la curva (AUC). Esto lo haremos mediante las funciones roc() y auc() de la librería pROC.
# Instalamos e importamos la librería ya citada.
if (!require("pROC")) install.packages("pROC"); library(pROC)## Loading required package: pROC## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, varPero primero, debemos realizar las predicciones pertinentes con nuestro nuevo modelo.
# Usamos el método "predict" de forma similar a como lo hicimos con el árbol de decisión.
prediccion_log <- predict(regresion_logistica_g, testx, type = "response")Ahora si, visionamos el resultado de nuestro modelo.
# Pintamos la curva
curva_roc <- roc(datos_regresion$Outcome, prediccion_log)## Setting levels: control = 0, case = 1## Setting direction: controls < casesplot(curva_roc)
# Mostramos el área bajo la curva ROC.
auc(curva_roc)## Area under the curve: 0.7879Aludiendo nuevamente al concepto de curva ROC. El eje x representa la tasa de falsos positivos y el eje y representa la tasa de verdaderos positivos, esta imagen nos ayuda a comprender un poco mejor el resultado obtenido.
El AUC (Area Under the Curve) es una medida del rendimiento de un modelo de clasificación binaria que representa la proporción de la zona bajo la curva ROC que está por encima de la línea de referencia (la cual representa una clasificación aleatoria). Un AUC con valor 1 indicaría el modelo perfecto, mientras que un AUC de 0.5 indica un modelo que no es mejor que una clasificación aleatoria.
Insistimos en el visionado de esta imagen junto con una breve lectura al árticulo al que pertenece lo cual nos ha ayudado a entender mejor el resultado.
En base a lo explicado, vemos que tenemos un modelo bastante bueno, podemos apreciar una curva cuya tendencia es claramente la de acercarse al modelo ideal y un AUC de 0.78, bastante cercano a 1.
Ahora repitamos el proceso esta vez frente a la variable BMI.
# Empleamos el método "glm()" apoyándonos en la documentación, agregar el argumento "family" que hace alusión al tipo de distribución estadística de los datos, como "Outcome" es una variable que solo tiene dos posbiles valores, empleamos la distribución "binomial".
regresion_logistica_b <- glm(datos_regresion$Outcome ~ datos_regresion$BMI, data = datos_regresion, family = "binomial")# Usamos el método "predict" de forma similar a como lo hicimos con el árbol de decisión.
prediccion_log_b <- predict(regresion_logistica_b, testx, type = "response")# Pintamos la curva
curva_roc2 <- roc(datos_regresion$Outcome, prediccion_log_b)## Setting levels: control = 0, case = 1## Setting direction: controls < casesplot(curva_roc2)
# Mostramos el área bajo la curva ROC.
auc(curva_roc2)## Area under the curve: 0.6799En este caso, el modelo tiene menos precisión cuando tratamos de realizar la regresión con BMI. Esto era de esperar ya que la correlación de Glucose con Outcome es mayor que en este caso (0.77 frente a 0.30). Aún así el resultado con BMI también es bueno con un AUC de 0.68 aproximadamente.
Conclusiones del proyecto.
Que hemos hecho:
En una primera parte de este proyecto, hemos planteado el problema de minería de datos que hemos estudiado. Hemos seleccionado un juego de datos, lo hemos revisado, analizado, limpiado y preprocesado.
En una segunda parte del proyecto, hemos generado diversos modelos, tanto supervisados como no supervisados para enfrentar este proyecto.
En la sección de modelos no supervisados hemos generado un modelo de clustering basado en el concepto de distancia, en concreto hemos generado un modelo k-means en función de la distancia euclidiana y luego hemos repetido el proceso con el concepto de distancia de Manhattan.
Hemos empleado los DBSCAN y OPTICS como una alternativa no supervisada al modelo k-means pero en este caso hemos tenido un mal resultado.
La sección de modelos supervisados ha comenzado generando un árbol de decisión, hemos analizado las reglas que pudimos extraer de él, hemos generado un diagrama visual del mismo y lo hemos puesto a prueba frente a un subconjunto de datos de prueba diferente al utilizado para entrenar el modelo.
Luego pretendimos continuar con los modelos supervisados haciendo una regresión lineal, pero nos topamos con que no fue una buena idea ya que la variable objetivo de nuestro problema es una variable categórica y por lo tanto no sigue una distribución normal.
Como alternativa al imprevisto anterior investigamos sobre la regresión logística especialmente pensado para abordar problemas con variables categóricas (como nuestro caso) donde el problema no sigue una distribución normal.
Empleamos la librería pROC y los métodos roc y auc para crear un modelo supervisado, también lo pusimos a prueba realizando predicciones con un subconjunto de datos de prueba de forma similar a como lo hicimos con el árbol de decisión.
Estudiamos el concepto de curva ROC y área bajo la curva (AUC) como indicadores de la precisión del modelo de regresión logística.
Resultados:
- Matriz de correlaciones: Nuevamente en este trabajo hemos utilizado el concepto de correlación para discernir aquellas variables significativas frente a la variable target, en este caso hemos visto que la variables Glucose e Insulin son las que mayor correlación tienen con Outcome (0.77 y 0.72 respectivamente) y las hemos empleado en la mayoría del trabajo.
1. Modelos no supervisados:
- K-means:
Hemos realizado el método para k=2 ya que nuestra variable target sólo podia tener dos valores. Hemos empleado los datos con las variables Glucose e Insulin como ya hemos mencionado.
Nos ha llamado mucho la atención esa agrupación tan marcada con forma claramente lineal, muy densa y compacta, probablemente se deba a las labores de imputación realizadas en la PRA 1 ya que empleamos una regresión lineal entre ambas variables.
De todos modos, en general el modelo se ha comportado dentro de lo esperado, no ha habido muestras sin clasificar en ninguno de los clusters y el solapamiento de ambos es una región relativamente pequeña.
En cuanto a las métricas usadas, no hemos apreciado ninguna diferencia entre usar la distancia euclidiana o la de Manhattan.
- DBSCAN y OPTICS:
Este ha sido el punto más problemático del trabajo. En primer lugar decidimos trabajar con Glucose e Insulin nuevamente. El criterio empleado para definir el valor minPts que alude al umbral de vecindad de las muestras entre sí para su ordenación y posterior clasificación, ha sido el de emplear la mitad de la longitud del intervalo menor de valores entre las variables usadas. Es decir, tomamos aquella variable cuya diferencia entre su valor máximo y mínimo sea menor, luego al resultado de esa resta la dividimos a la mitad.
Lo hemos hecho así ya que ese intervalo de menor longitud será más susceptible frente a los cambios en el parámetro minPts. En nuestro caso el intervalo más corto fue el de Glucose con un dominio de longitud 143, quedando minPts en 71.5
Tras ordenar las muestras y presentar el diagrama de alcanzabilidad, a pesar de intentar con numerosos valores de eps_cl tenia mala pinta y no se distinguían los dos clusters que necesitabamos. Luego probamos con todas las variables, el resultado cambio ligeramente pero seguía siendo igual de malo.
Creemos que este mal resultado se debe a la distribución tan marcada que ya vimos en k-means considerando lo densa que es esa agrupación lineal de las muestras, es posible que haya sesgado el modelo tendiendo a magnificar las distancias relativas entre las muestras, ya que aunque visual o numéricamente no estén tan lejos entre sí, esa línea tan marcada puede hacer que el modelo sea más sensible a a las distancias que nosotros consideraríamos normales.
2. Modelos supervisados:
Árbol de decisión:
Creación y descripción del modelo:
Dividimos nuestro conjunto de datos en dos subconjuntos, uno de entrenamiento para generar el modelo, el cual sería 2/3 del conjunto total, y el tercio restante fue para el subconjunto de prueba para evaluar nuestro modelo más adelante.
Para la generación del modelo usamos la librería C50 y repetimos un proceso similar al que vimos en la PEC 3. Estudiamos las reglas generadas y la tasa de error del modelo. La tasa de error fue del 18.3% durante la generación del modelo, y la cantidad de falsos positivos y falsos negativos coincidió, por lo que nuestro árbol no presenta sesgo alguno.
En cuanto a las reglas, a modo de resumen, el modelo sugiere que los pacientes mayores de 34 años son más propensos a ser diabéticos. También sugiere que la glucosa de una persona sana no debería superar los 127 mg/dL, del mismo modo la insulina debería estar en un rango comprendido entre 106 y 185 unidades para una persona sana. En cuanto al número de embarazos, el riesgo de desarrollar diabetes aumenta a partir de 6. En cuanto a la complexión de los pacientes, un índice de masa corporal superior a 30 ya es un factor de riesgo para padecer diabetes.
Por último, el porcentaje de uso de las variables que hace el modelo está bastante repartido a excepción de la presión sanguínea que solo se usó en un 6.52% de las veces y de la variable DiabetesPedigreeFunction que ni si quiera aparece.
- Validación del modelo:
Evaluamos la precisión de nuestro modelo con el subconjunto de datos de prueba y el comando predict(). En base a esto calculamos la matriz de confusión y con ello la precisión, la cual fue de un 76.96% lo cual está bastante bien.
Otro dato a destacar es que según la matriz de confusión nuestro modelo tiene un ligero sesgo inclinado hacia los falsos positivos, ya que diagnosticó erróneamente 33 pacientes como diabéticos y 20 como no diabéticos. A pesar de que el sesgo (dentro de lo leve que es dadas las cantidades), pueda parecer algo negativo, en realidad consideramos que es un modelo “precabido”, ya que en este caso es mejor diagnosticar un falso positivo y realizar un seguimiento al paciente en cuestión sea o no diabético que el caso contrario.
Con un falso positivo puede no haber problemas, pero un falso negativo si es peligroso en cuestiones de salud, ya que ni el paciente ni los médicos sabe que padece una enfermedad.
- Regresión lineal:
Como en la PRA 1 ya hicimos una regresión lineal entre la glucosa y la insulina para imputar los valores perdidos del juego de datos, no tiene mucho sentido que hagamos el ajuste lineal de la variable target frente y la glucosa y luego frente a la insulina ya que sería redundante por la fuerte correlación que existe entre ambas variables a raíz de la imputación de valores ya comentada.
Por ese motivo, vamos a emplear la glucosa y el índice de masa corporal (BMI) que también fue utilzada en la PRA 1 en conjunto a SkinThickness para imputar valores perdidos. Creemos que las medidas de la glucosa y del índice de masa corporal son medidas más fáciles de obtener y cotidianas que las de la insulina y el espesor de la piel. Por este motivo haremos una regresión lineal de “Outcome vs Glucose” y de “Outcome vs BMI”.
Tratamos de realizar una regresión lineal pero nos dimos cuenta de que nuestra variable target era una variable categórica y no tiene sentido someterla a una regresión lineal ya que al no ser una variable continua, su distribución no es normal.
Planteamos una alternativa.
- Regresión logística:
Debido al imprevisto con la regresión lineal, tras investigación en internet nos topamos con el concepto de regresión logística, un método especialmente diseñado para variables categóricas que no siguen una distribución normal como ocurre en nuestro caso.
Consultamos la documentación de la librería pROC y una guía proporcionada por la UCLA sobre el uso de las regresiones logísticas en R.
Por los mismos motivos que argumentamos en nuestro fallido intento de regresión lineal, vamos a emplear aquí un ajuste logístico de la variable target frente a a la glucosa y el índice de masa corporal (BMI).
Tras generar el modelo con el método glm() para una distribución binomial, realizamos las predicciónes con predict(). Repetimos el proceso tanto para el ajuste de nuestra variable target frente a la glucosa como el ajuste realizado frente al índice de masa corporal (BMI).
En ambos casos evaluamos el resultado del modelo con los conceptos de curva ROC y AUC (Area Under the Curve). Es de especial ayuda esta imagen para comprender los resultados obtenidos. Todo ello esta explicado en su correspondiente sección en detalle.
Obtuvimos buenos resultados en ambos casos, para la glucosa obtuvimos un AUC de 0.79 aproximadamente y la curva ROC mostraba una clara tendencia a arquearse hacia el vertice superior izquierdo de la gráfica (lo cual representa un AUC de 1, el modelo perfecto).
Por otro lado, el ajuste en base al índice de masa corporal (BMI) no fue tan bueno aunque tampoco malo, la curva obtenida es más irregular y más cercana a la diagonal, esto también se ve reflejado en el AUC obtenida que fue de 0.68 (aproximadamente un modelo predictivo un 10% peor).
Reflexiones y propuestas:
El mayor incoveniente que vemos en este proyecto fue el mal resultado obtenido con el método DBSCAN y OPTICS, además de no saber exactamente el porque de esa tendencia lineal tan fuertemente marcada con el método k-means. En general los modelos supervisados tuvieron mucho mejor desempeño y en nuestra modesta opinión personal, se pueden entender mejor aunque sabemos que no siempre serán la opción ideal.
Otro punto a señalar es el no haber aprovechado la totalidad de los datos, es decir, siempre nos apoyamos en la glucosa, la insulina el índice de masa corporal (BMI) y la variable target para realizar nuestras predicciones y crear nuestros modelos (a excepción del árbol de decisión que si usa todas las variables). Quizá sería posible mediante otras técnicas aprovechar también esa información con otros criterios más exhaustivos que la correlación matemática a la hora de elegir que variables emplear.
Por último, una posible propuesta para el futuro seria el de buscar un juego de datos similar pero de otra fuente y otra población para aplicar estos métodos y comprobar si se pueden utilizar en una situación real como apoyo en el diagnóstico clínico. Pensamos que incluso puede ser una propuesta interesante como proyecto de TFM cunado tengamos más conocimientos y herramientas más potentes para mejorar este proyecto.