Practica1 Fundamentos de Minería de Datos y Big Data
Departamento de Matemáticas y Ciencia de Datos
2023-09-28
SOLUCIONES
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.
library(readr)
heart <- read_delim("heart_MD.csv", delim = ";",
escape_double = FALSE, trim_ws = TRUE)## Rows: 1025 Columns: 14
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (1): oldpeak
## dbl (13): age, sex, cp, trestbps, chol, fbs, restecg, thalach, exang, slope,...
##
## ℹ 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.summary(is.na(heart))## age sex cp trestbps
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:1025 FALSE:1025 FALSE:1025 FALSE:1019
## TRUE :6
## chol fbs restecg thalach
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:1025 FALSE:1025 FALSE:1025 FALSE:1025
##
## exang oldpeak slope ca
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:1025 FALSE:1025 FALSE:1025 FALSE:1025
##
## thal target
## Mode :logical Mode :logical
## FALSE:1025 FALSE:1025
## #hay 6 valores NA en la variable "trestbps". Las correlaciones mayores las tiene con colesterol y edad:
continuous_var<-data.frame(heart$trestbps,heart$chol,heart$age,heart$thalach)
data_sin_NA<-na.omit(continuous_var)
mydata.cor = cor(data_sin_NA)
index<-is.na(heart$trestbps) #identifico los NA
modelo1<-lm(heart.trestbps~heart.chol+heart.age,data=data_sin_NA) #estimo el modelo
datos_predecir<-continuous_var[index,] #localizo los valores que tengo que imputar
NA_prediccion<-predict(modelo1,newdata=datos_predecir[,2:3]) #realizo la predicción
#introducimos los valores faltantes
id<-which(index==TRUE)
heart[id,4]=NA_prediccionEjercicio 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(pastecs)## Warning: package 'pastecs' was built under R version 4.2.3options(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 134890.347 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.600 246.00 0.1493
## SE.mean 0.283 0.0144 0.0322 0.546 1.61 0.0111
## CI.mean.0.95 0.556 0.0282 0.0631 1.072 3.16 0.0219
## var 82.306 0.2119 1.0602 305.762 2661.79 0.1271
## std.dev 9.072 0.4604 1.0296 17.486 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 NA 1025.0000 1025.0000
## nbr.null 497.0000 0.000 680.0000 NA 74.0000 578.0000
## nbr.na 0.0000 0.000 0.0000 NA 0.0000 0.0000
## min 0.0000 71.000 0.0000 NA 0.0000 0.0000
## max 2.0000 202.000 1.0000 NA 2.0000 4.0000
## range 2.0000 131.000 1.0000 NA 2.0000 4.0000
## sum 543.0000 152842.000 345.0000 NA 1420.0000 773.0000
## median 1.0000 152.000 0.0000 NA 1.0000 0.0000
## mean 0.5298 149.114 0.3366 NA 1.3854 0.7541
## SE.mean 0.0165 0.719 0.0148 NA 0.0193 0.0322
## CI.mean.0.95 0.0324 1.410 0.0290 NA 0.0379 0.0632
## var 0.2787 529.263 0.2235 NA 0.3816 1.0625
## std.dev 0.5279 23.006 0.4728 NA 0.6178 1.0308
## coef.var 0.9965 0.154 1.4046 NA 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 ser capaz de clasificar las variables. ¿Para cuáles tienen sentido ciertos estadísticos descriptivos?
#Por ejemplo, "fbs"- fast blood sugar- tiene una media de 0.15, aproximadamente, indicando que un 15% de los individuos tienen niveles altos de azúcar. Sin embargo, la desviación típica o el coeficiente de variación no tienen un significado claro. Lo mismo ocurre con "restecg", que es una variable cualitativa nominal. Los estadísticos descriptivos no significan nada.
barplot(prop.table(table(heart$restecg)))
#Enseguida vemos que no vamos a poder aprender mucho de los individuos con ECG tipo 2.
#queremos predecir la variable "target". ¿Qué significa la media de 0.51? Una buena señal: la muestra está equilibrada: tenemos aproximadamente el mismo número de casos con problemas coronarios y sin ellos, por lo que podremos aprender bien de ambos tipos de casos.
#Analizando el colesterol, llama la atención el bajo coef.var de esta variable. Esto implica fluctúa poco 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 (repetimos 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 que, por otro lado, pueden ser de interés.
#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.
#Aquí podemos remarcar 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, sin azúcar elevada y con electrocardiograma nivel 0 o nivel 1 . ¿Está bien representado un individuo joven y sano? ¿podrá nuestro modelo predecir qué ocurrirá con un individuo de 75 años con azúcar alta?
#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.3dat <- 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. Entonces, se podrían establecer hipótesis de variables que hagan que el impacto del colesterol sobre la probabilidad de tener enfermedad cardiaca se vean modificadas. Un ejemplo: la edad (es normal empezar por las variables más sencillas de comprender por el analista)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(50,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
Discute cuál de los dos siguientes modelos estimados representa lo que hemos aprendido previamente.
M1<-glm( target ~chol+age, data = heart, family = binomial)
summary(M1)##
## Call:
## glm(formula = target ~ chol + age, family = binomial, data = heart)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.685 -1.122 0.736 1.105 1.741
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.34069 0.47589 7.02 0.0000000000022 ***
## chol -0.00217 0.00130 -1.67 0.094 .
## age -0.05048 0.00753 -6.71 0.0000000000198 ***
## ---
## 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: 1362.2 on 1022 degrees of freedom
## AIC: 1368
##
## Number of Fisher Scoring iterations: 4M2<-glm(target~chol*age, data = heart, family = binomial)
summary(M2)##
## Call:
## glm(formula = target ~ chol * age, family = binomial, data = heart)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.818 -1.108 0.634 1.126 1.725
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.189498 2.025581 4.04 0.000053 ***
## chol -0.022579 0.008318 -2.71 0.00664 **
## age -0.136076 0.035421 -3.84 0.00012 ***
## chol:age 0.000359 0.000144 2.49 0.01267 *
## ---
## 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.0 on 1021 degrees of freedom
## AIC: 1364
##
## Number of Fisher Scoring iterations: 4#El modelo 2 es el que permite que el coeficiente asociado a la variable colesterol cambie de signo en función de la edad¿A partir de qué valor de la edad el impacto del colesterol sobre la enfermedad cardiaca se vuelve positivo?
Se puede ver que,en el modelo M2,
\[ logit(target_i)=8.18-0.022 chol_i -0.13 age_i+0.000359 chol_i \times age_i \]
El efecto del colesterol sobre la enfermedad cardiaca será positivo cuando
\[ -0.022+0.000359age_i>0 \implies age>\frac{0.022}{0.000359}=61.3 \] años.
Ejercicio 7
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 trata de ver la interacción de estas con el colesterol y la edad de los individuos.
Pista: reutiliza los gráficos del ejercicio 5, adaptándolos
age_v2<-seq(50,66,2) #se introduce una secuencia de edades
cholesterol<-seq(150,400,1) #se introduce una secuencia de valores de coloesterol
fbs1<-rep(1,length(cholesterol))
fbs0<-rep(0,length(cholesterol))
pred1<-matrix(0,length(age_v2),length(cholesterol))#se crea una matriz para almacenar las
pred2<-matrix(0,length(age_v2),length(cholesterol))#se crea una matriz para almacenar las
data_chol1<-data.frame(chol=cholesterol, fbs=fbs1)
data_chol2<-data.frame(chol=cholesterol,fbs=fbs0)
for (i in 1:length(age_v2))
{
heart_sub<-subset(heart,age>age_v2[i])
model <- glm( target ~ fbs*chol, data = heart_sub, family = binomial)
pred1[i,]<-predict(model,newdata=data_chol1, type="response") #diabetes
pred2[i,]<-predict(model,newdata=data_chol2, type="response") #no diabetes
}
row.names(pred1)<-age_v2
colnames(pred1)<-cholesterol
col_set <- rainbow(nrow(pred1))
layout(mat = matrix(c(1, 2),
nrow = 1,
ncol = 2),
heights = c(2,2), # Heights of the two rows
widths = c(1,1)) # Widths of the two columns
# Plot 1
par(mar = c(4,4,0,0))
matplot(t(pred1), x = seq(150,400,1),type = "b",
ylab="Probabiilty", xlab="cholesterol",col=col_set,ylim=c(0,1),cex.main=1)
text(200, 0.8, "Diabetes")
legend("bottomleft", inset=0.01, legend=age_v2, col=col_set,pch=15:19,
bg= ("white"), horiz=F)
# Plot 2
par(mar = c(4,4,0,0))
matplot(t(pred2), x = seq(150,400,1),type = "b",
ylab="Probabiilty", xlab="cholesterol",col=col_set,ylim=c(0,1))
text(200, 0.8, "No Diabetes")
press_d<- as.factor(ifelse(heart$trestbps<140, 0,1)) #creamos una variable que mida hipertensión (previamente se ha buscado en google cuándo se es hipertenso)
heart<-data.frame(heart,press_d) #añadimos la variable a los datos
age_v2<-seq(50,66,2) #se introduce una secuencia de edades
cholesterol<-seq(150,400,1) #se introduce una secuencia de valores de coloesterol
press1<-rep(1,length(cholesterol))
press0<-rep(0,length(cholesterol))
pred1<-matrix(0,length(age_v2),length(cholesterol))#se crea una matriz para almacenar las
pred2<-matrix(0,length(age_v2),length(cholesterol))#se crea una matriz para almacenar las
data_chol1<-data.frame(chol=cholesterol,press_d=press1)
data_chol2<-data.frame(chol=cholesterol,press_d=press0)
for (i in 1:length(age_v2))
{
heart_sub<-subset(heart,age>age_v2[i])
model <- glm( target ~ as.factor(press_d)*chol, data = heart_sub, family = binomial)
pred1[i,]<-predict(model,newdata=data_chol1, type="response") #press
pred2[i,]<-predict(model,newdata=data_chol2, type="response") #no press
}
row.names(pred1)<-age_v2
colnames(pred1)<-cholesterol
col_set <- rainbow(nrow(pred1))
layout(mat = matrix(c(1, 2),
nrow = 1,
ncol = 2),
heights = c(2,2), # Heights of the two rows
widths = c(1,1)) # Widths of the two columns
# Plot 1:
par(mar = c(4,4,0,0))
matplot(t(pred1), x = seq(150,400,1),type = "b",
ylab="Probabiilty", xlab="cholesterol",col=col_set,ylim=c(0,1),cex.main=1)
text(200, 0.8, "Hypertension")
legend("bottomleft", inset=0.01, legend=age_v2, col=col_set,pch=15:19,
bg= ("white"), horiz=F)
#plot 2:
par(mar = c(4,4,0,0))
matplot(t(pred2), x = seq(150,400,1),type = "b",
ylab="Probabiilty", xlab="cholesterol",col=col_set,ylim=c(0,1))
text(200, 0.8, "No Hypertension")
Tanto la hiptertensión como la diabetes, parecen tener interacción con el nivel de colesterol y la edad. Ambos factores inciden en la probabilidad de afección cardiaca en individuos mayores (más la hipertensión que el azúcar, aunque recordemos que no tenemos una proporción alta de individuos con azúcar alta)
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.27425 0.62026 6.89 0.0000000000055 ***
## chol -0.00193 0.00130 -1.49 0.136
## age -0.04589 0.00774 -5.93 0.0000000030471 ***
## trestbps -0.00944 0.00392 -2.41 0.016 *
## as.factor(fbs)1 -0.01073 0.18312 -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.2 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.782 -1.108 0.659 1.094 1.719
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.4826304 3.2779616 1.37 0.1715
## chol -0.0030465 0.0132609 -0.23 0.8183
## age -0.1136158 0.0373693 -3.04 0.0024 **
## as.factor(fbs)1 18.0539707 9.7166045 1.86 0.0632 .
## trestbps 0.0166460 0.0212676 0.78 0.4338
## chol:age 0.0002645 0.0001529 1.73 0.0836 .
## chol:as.factor(fbs)1 -0.0770277 0.0383234 -2.01 0.0444 *
## age:as.factor(fbs)1 -0.2849305 0.1644622 -1.73 0.0832 .
## chol:trestbps -0.0000999 0.0000821 -1.22 0.2234
## chol:age:as.factor(fbs)1 0.0012209 0.0006445 1.89 0.0582 .
## ---
## 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.0 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.167 0.286 0.445 0.617new2<-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.282