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:
Busca las variables con valores omitidos y trata de imputar un valor usando alguna de las técnicas aprendidas en clase.
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.
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.
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.
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.
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
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