Práctica 1 Fundamentos de Minería de Datos y Big Data

En esta práctica nos preguntamos qué nivel de colesterol aumenta la probabilidad de que suframos una enfermedad coronaria. Para ello, tenemos una base de datos obtenida del repositorio de “machine learning” de UCLA: https://archive.ics.uci.edu/dataset/45/heart+disease. Puedes consultar qué variables contiene y los posibles valores que pueden tomar. Como verás, tiene “missing values”. Perio vayamos por partes:

Ejercicio 1

Busca las variables con valores omitidos y trata de imputar un valor usando alguna de las técnicas aprendidas en clase.

Ejercicio 2

1.- Discute si los valores de las variables parecen razonables (lo que desconozcas lo tendrás que buscar), y cómo de representativa es la muestra disponible.

library(readr)
heart<-read_csv("C:/Users/Hortera de la Cutre/Desktop/DwD2/heart_FINAL.csv")
## Rows: 1025 Columns: 14
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (14): age, sex, cp, trestbps, chol, fbs, restecg, thalach, exang, oldpea...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
library(pastecs)
## Warning: package 'pastecs' was built under R version 4.2.3
options(scipen=100)
options(digits=3)
stat.desc(heart)
##                    age       sex        cp   trestbps      chol       fbs
## nbr.val       1025.000 1025.0000 1025.0000   1025.000   1025.00 1025.0000
## nbr.null         0.000  312.0000  497.0000      0.000      0.00  872.0000
## nbr.na           0.000    0.0000    0.0000      0.000      0.00    0.0000
## min             29.000    0.0000    0.0000     94.000    126.00    0.0000
## max             77.000    1.0000    3.0000    200.000    564.00    1.0000
## range           48.000    1.0000    3.0000    106.000    438.00    1.0000
## sum          55795.000  713.0000  966.0000 134902.000 252150.00  153.0000
## median          56.000    1.0000    1.0000    130.000    240.00    0.0000
## mean            54.434    0.6956    0.9424    131.612    246.00    0.1493
## SE.mean          0.283    0.0144    0.0322      0.547      1.61    0.0111
## CI.mean.0.95     0.556    0.0282    0.0631      1.074      3.16    0.0219
## var             82.306    0.2119    1.0602    306.835   2661.79    0.1271
## std.dev          9.072    0.4604    1.0296     17.517     51.59    0.3565
## coef.var         0.167    0.6618    1.0925      0.133      0.21    2.3885
##                restecg    thalach     exang   oldpeak     slope        ca
## nbr.val      1025.0000   1025.000 1025.0000 1025.0000 1025.0000 1025.0000
## nbr.null      497.0000      0.000  680.0000  329.0000   74.0000  578.0000
## nbr.na          0.0000      0.000    0.0000    0.0000    0.0000    0.0000
## min             0.0000     71.000    0.0000    0.0000    0.0000    0.0000
## max             2.0000    202.000    1.0000    6.2000    2.0000    4.0000
## range           2.0000    131.000    1.0000    6.2000    2.0000    4.0000
## sum           543.0000 152842.000  345.0000 1098.3000 1420.0000  773.0000
## median          1.0000    152.000    0.0000    0.8000    1.0000    0.0000
## mean            0.5298    149.114    0.3366    1.0715    1.3854    0.7541
## SE.mean         0.0165      0.719    0.0148    0.0367    0.0193    0.0322
## CI.mean.0.95    0.0324      1.410    0.0290    0.0720    0.0379    0.0632
## var             0.2787    529.263    0.2235    1.3808    0.3816    1.0625
## std.dev         0.5279     23.006    0.4728    1.1751    0.6178    1.0308
## coef.var        0.9965      0.154    1.4046    1.0966    0.4459    1.3668
##                   thal    target
## nbr.val      1025.0000 1025.0000
## nbr.null        7.0000  499.0000
## nbr.na          0.0000    0.0000
## min             0.0000    0.0000
## max             3.0000    1.0000
## range           3.0000    1.0000
## sum          2382.0000  526.0000
## median          2.0000    1.0000
## mean            2.3239    0.5132
## SE.mean         0.0194    0.0156
## CI.mean.0.95    0.0380    0.0307
## var             0.3852    0.2501
## std.dev         0.6207    0.5001
## coef.var        0.2671    0.9745
#Es importante que sean capaces de clasificar las variables. ¿Para cuáles tienen sentido ciertos estadísticos descriptivos? Que vean que hay variables cualitativas como thal y cp que deberán identificar y tener cuidado con cómo usan los estadísticos descriptivos.

#queremos predecir la variable "target". ¿Qué significa la media? Una  buena señal es que la muestra está equilibrada: tenemos aproximadamente el mismo número de casos con problema coronario y sin  él, por lo que nos evitaremos problemas a la hora de entrenar modelos.


#Llama la atención el bajo coef.var de la variable "chol": nivel de colesterol. Esto implica que hay poca fluctuación de esta variable en la muestra. De hecho, mirando el histograma vemos que por debajo de 150 y por encima de 350 la muestra apenas si tiene casos
hist(heart$chol)

# si buscamos información (repetir lo importante que es conocer bien el contexto de las variables que manejamos), vemos que un paciente sano tiene un colesterol por debajo de 200, por lo que la muestra tiene bastante sesgo hacia individuos con colesterol alto.



# Por otro lado, llama la atención el coeficiente de variación -por pequeño- de la variable EDAD. Si vemos su histograma

hist(heart$age)

# encontramos que hay pocos casos de individuos jóvenes (por debajo de 40 años) y de mayores (por encima de 70 años) Habrá que tener en cuenta, si usamos esta variable, la poca precisión con la que contaremos al tratar de analizar esos casos.

#algo parecido encontramos en la variable "presión sanguínea en reposo" (trestbps), la cual posee valores normales entre 90 y 120. En general la muestra está compuesta por individuos hipertensos.


#remarcad la importancia del análisis descriptivo para poder hacernos una idea de a quién representa nuestra muestra:

#podremos analizar a individuos entre 40 y 70 años, con niveles de colesterol entre 150 y 350 y con tendencia a la hipertensión. ¿Está bien representado un individuo joven sano?

#NOTA: se pueden discutir más variables tratando de pronunciarse acerca del tipo de muestra de que se dispone y, por tanto, de lo que posteriormente se podrá afirmar con mayor/menor precisión.

2.-Analiza de manera descriptiva/visual qué dicen los datos sobre la relación entre tener una enfermedad coronaria y el nivel de colesterol.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
dat <- data.frame(cholesterol = heart$chol,disease = heart$target,age=heart$age)

ggplot(dat, aes(x = cholesterol, fill = as.factor(disease))) +                       # Draw overlaying histogram
  geom_histogram(position = "identity", alpha = 0.2, bins = 50)

#aquí se puede ver que no parece haber una relación clara entre el nivel de colesterol y la enfermedad cardiaca, lo cual no parece muy sensato.

Ejercicio 3

Como habrás visto, dicha relación no parece muy clara . Aunque es bien sabido el impacto del colesterol en la salud cardiaca. ¿Ocultan algo los datos? ¿Qué crees que puede estar pasando?

#se puede comentar el efecto de las posibles variables "confusoras" y recordar la paradoja de Simpson-Yule.

Ejercicio 4

Nuestra sospecha es que la edad es un factor “confusor”. Para ello, escribe un script en R que obtenga el coeficiente de la regresión logística que predice probabilidad de enfermedad cardiaca frente a nivel de colesterol para distintos grupos de edad. Como la muestra no es grande, recomendamos hacer una ventana de edades (empezamos con toda la muestra y vamos eliminando el más joven a cada paso). Haz un gráfico y analiza a partir de qué edad el modelo detecta una relación que no contradice la evidencia médica. Piensa, entonces, si este modelo se puede generalizar a la población o debería hacerse específico para una subpoblación

Pista: se muestra un posible script del que se han omitido algunas indicaciones

coef<-vector() #almacena los distintos coeficientes
age_v<-seq(40,65,1) #genera una secuencia de edades

for (i in 1:length(age_v))
{
  
  heart_sub<-subset(heart,age>age_v[i])  #elige la submuestra correspondiente
  model <- glm( target ~chol, data = heart_sub, family = binomial) #estima el logit
  coef[i]<-model$coefficients[2] #almacena el coeficiente
  
}

plot(age_v,coef) #realiza el gráfico

#podemos ver que, a partir de los 55 años el coeficiente empieza a ser positivo. Esto indica que, a priori, utilizando sólo la edad y el colesterol, nos va a ser complicado realizar una predicción para individuos con menos de 55 años.

Ejercicio 5

Vuelve a pensar tus respuestas dadas en el ejercicio 2 y en el ejercicio 3.

Para los grupos de edad donde el modelo detecta que el colesterol tiene un efecto negativo, haz un gráfico que muestre la probabilidad de enfermedad cardiaca condicionado al nivel de colesterol.

Pista: se muestra un posible script del que se han omitido algunas indicaciones

age_v2<-seq(40,65,1) #se introduce una secuencia de edades
cholesterol<-seq(150,600,1) #se introduce una secuencia de valores de coloesterol
pred<-matrix(0,length(age_v2),length(cholesterol))#se crea una matriz para almacenar las prediicciones
data_chol<-data.frame(chol=cholesterol)
for (i in 1:length(age_v2))
{
  
  heart_sub<-subset(heart,age>age_v2[i]) 
  model <- glm( target ~chol, data = heart_sub, family = binomial)
  pred[i,]<-predict(model,newdata=data_chol, type="response")
}

row.names(pred)<-age_v2
colnames(pred)<-cholesterol
col_set <- rainbow(nrow(pred))

matplot(t(pred), x = seq(150,600,1),type = "b",main="Probability of heart disease at different ages",
        ylab="Probabiilty", xlab="cholesterol",col=col_set)

legend("bottomleft", inset=0.01, legend=age_v2, col=col_set,pch=15:19,
       bg= ("white"), horiz=F)

#podemos ver una interacción clara entre el nivel de colesterol y la edad. 

Ejercicio 6

Se suele incidir en que niveles de azúcar altos y de tensión arterial son factores de riesgo para tener enfermedad cardiaca. Discute el impacto del nivel de azúcar y de tensión arterial sobre la probabilidad de enfermedad cardiaca, así como tratar de ver la interacción de estas con el colesterol y la edad de los individuos.

Pista: se muestra un posible script para visualizar las interacciones

#crea una variable "age_q" que divida a los individuos por grupos de edades como te parezca interesante


age_q<- as.factor(ifelse(heart$age<40, '40',
                         ifelse(heart$age>=40 & heart$age<50, '40-50', 
                                ifelse(heart$age>=50 & heart$age<65, '50-65', '65 and more'))))


sugar_heart<- as.factor(ifelse(heart$fbs==1&heart$target==1, 'HighSugar_Heart1',
                               ifelse(heart$fbs==1&heart$target==0,'HighSugar_Heart0',                               ifelse(heart$fbs==0&heart$target==1,'LowSugar_Heart1','LowSugar_Heart0'))))

heart<-data.frame(heart,sugar_heart,age_q)

interaction.plot(x.factor = heart$age_q,
                 trace.factor = heart$sugar_heart, 
                 response = heart$chol, fun = mean ,xlab="Age",
                 ylab="cholesterol", trace.label="Sugar and Heart Disease", 
                 col=c("green","red","blue","black"),
                 lty=4, lwd=2.5 )

#Se debe analizar con cuidado cómo el efecto conjunto de azúcar y colesterol tiene un impacto creciente con la edad en la enfermedad cardiaca (línea roja) comparado con los individuos con alto colesterol y edad que no tienen azúcar.

#se puede hacer un gráfico similar con la tensión arterial "trestbps"



press_heart<- as.factor(ifelse(heart$trestbps<122&heart$target==1, 'LowPress_Heart1',
                         ifelse(heart$trestbps<122&heart$target==0, 'LowPress_Heart0', 
                                ifelse(heart$trestbps>=122&heart$target==1, 'HighPress_Heart1', 'HighPress_Heart0'))))





heart<-data.frame(heart,press_heart)


interaction.plot(x.factor = heart$age_q,
                 trace.factor = heart$press_heart, 
                 response = heart$chol, xlab="Age",
                 ylab="cholesterol", trace.label="Blood pressure and Heart Disease", 
                 col=c("green","red","blue","black"),
                 lty=4, lwd=2.5 )

#no se ve, sin embargo, una interacción destacable entre la tensión arterial, el colesterol y la edad, pero sí entre el colesterol y la presión arterial

#respecto al nivel de azúcar

sugar_heart<- as.factor(ifelse(heart$fbs==1&heart$target==1, 'HighSugar_Heart1',
                               ifelse(heart$fbs==1&heart$target==0, 'HighSugar_Heart0', 
                                      ifelse(heart$fbs==0&heart$target==1, 'LowSugar_Heart1', 'LowSugar_Heart0'))))

heart<-data.frame(heart,sugar_heart)

interaction.plot(x.factor = heart$age_q,
                 trace.factor = heart$sugar_heart, 
                 response = heart$chol, fun = mean ,xlab="Age",
                 ylab="cholesterol", trace.label="Sugar and Heart Disease", 
                 col=c("green","red","blue","black"),
                 lty=4, lwd=2.5 )

#Sí se aprecia  interacción entre nivel de azúcar, edad y nivel de colesterol

Ejercicio 7

Estima y discute el modelo de regresión logística que crees más lógico tras lo aprendido anteriormente. ¿Cómo has incluido las interacciones? ¿Qué grupos de individuos podrían dar lugar a predicciones contraintuitivas?

#vistas las interacciones entre edad, azúcar y tensión arterial, este modelo sería insuficiente:

mod1<-glm(target~chol+age+trestbps+as.factor(fbs),data=heart,family="binomial")
summary(mod1)
## 
## Call:
## glm(formula = target ~ chol + age + trestbps + as.factor(fbs), 
##     family = "binomial", data = heart)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.690  -1.112   0.743   1.098   1.655  
## 
## Coefficients:
##                 Estimate Std. Error z value        Pr(>|z|)    
## (Intercept)      4.27379    0.61965    6.90 0.0000000000053 ***
## chol            -0.00194    0.00130   -1.49           0.135    
## age             -0.04585    0.00774   -5.92 0.0000000032216 ***
## trestbps        -0.00944    0.00392   -2.41           0.016 *  
## as.factor(fbs)1 -0.01090    0.18310   -0.06           0.953    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1420.2  on 1024  degrees of freedom
## Residual deviance: 1356.1  on 1020  degrees of freedom
## AIC: 1366
## 
## Number of Fisher Scoring iterations: 4
# el modelo más adecuado a lo analizado gráficamente, consiste en

mod2<-glm(target~chol*age+as.factor(fbs)*chol*age+trestbps*chol,data=heart,family="binomial")
summary(mod2)
## 
## Call:
## glm(formula = target ~ chol * age + as.factor(fbs) * chol * age + 
##     trestbps * chol, family = "binomial", data = heart)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
##  -1.78   -1.11    0.66    1.10    1.72  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)   
## (Intercept)               4.5547793  3.2724258    1.39   0.1640   
## chol                     -0.0033424  0.0132414   -0.25   0.8007   
## age                      -0.1134679  0.0373634   -3.04   0.0024 **
## as.factor(fbs)1          18.0584118  9.7191055    1.86   0.0632 . 
## trestbps                  0.0160340  0.0212339    0.76   0.4502   
## chol:age                  0.0002640  0.0001528    1.73   0.0840 . 
## chol:as.factor(fbs)1     -0.0770465  0.0383337   -2.01   0.0444 * 
## age:as.factor(fbs)1      -0.2849130  0.1644981   -1.73   0.0833 . 
## chol:trestbps            -0.0000975  0.0000819   -1.19   0.2341   
## chol:age:as.factor(fbs)1  0.0012209  0.0006447    1.89   0.0583 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1420.2  on 1024  degrees of freedom
## Residual deviance: 1341.1  on 1015  degrees of freedom
## AIC: 1361
## 
## Number of Fisher Scoring iterations: 4
#donde se permite que el efecto del nivel de colesterol sobre la enfermedad cardiaca dependa de los valores de la edad, de la tensión y si es diabético o no.

#se puede proponer realizar predicciones de diferentes individuos representativos (unirlo a lo visto anteriormente) y comparar lo obtenido por el modelo con la idea de criticar el uso, de nuevo, del p-valor. Por ejemplo ¿qué predice el modelo para un individuo de 70 años con 150,200,250 y 300 de colesterol, con azúcar y tensión de 100? ¿qué ocurre si el individuo tiene 40 años?


new1<-data.frame(age=c(70,70,70,70),chol=c(150,200,250,300),fbs=c(1,1,1,1),trestbps=c(100,100,100,100))

pred1<-predict(mod2, newdata=new1, type="response")
pred1
##     1     2     3     4 
## 0.168 0.287 0.446 0.616
new2<-data.frame(age=c(40,40,40,40),chol=c(150,200,250,300),fbs=c(1,1,1,1),trestbps=c(100,100,100,100))

pred2<-predict(mod2, newdata=new2, type="response")
pred2
##     1     2     3     4 
## 0.975 0.894 0.645 0.281