rm(list = ls()) #Borra todas las variables
ls()## character(0)
library(readr)
library(GGally)
library(dplyr)
library(tidyverse)
#library(tidymodels)
library(gridExtra)
#library(rsample)
library(pROC)
library(ggplot2)
library(caTools)
library(caret)#Referencias para este trabajo
# https://rpubs.com/Edenner/854339 (Heart Disease Prediction)
# https://halimatusyak.medium.com/heart-disease-uci-logistic-regression-in-r-b95b821088e6
#setwd("D:/Dropbox/Doctorado/Cursos_2022_2do_Cuat/Aprendizaje Estadistico")
#heart<-read_csv("D:/Dropbox/Doctorado/Cursos_2022_2do_Cuat/Aprendizaje Estadistico/heart.csv")
setwd("C:/Users/jpere/Dropbox/Doctorado/Cursos_2022_2do_Cuat/Aprendizaje Estadistico")
heart<-read_csv("heart.csv")## Rows: 303 Columns: 14
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## dbl (14): age, sex, cp, trestbps, chol, fbs, restecg, thalach, exang, oldpea...
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
names = c("age", "sex", "cp", "trestbps", "chol", "fbs", "restecg", "thalach", "exang", "oldpeak", "slope", "ca", "thal", "heart_disease")
colnames(heart) <- names
heart <- heart %>%
mutate(sex = case_when(sex == 0 ~ "female",
sex == 1 ~ "male")) %>%
mutate(cp = case_when(cp == 0 ~ "typical angina",
cp == 1 ~ "atypical angina",
cp == 2 ~ "non-anginal pain",
cp == 3 ~ "asymptomatic")) %>%
mutate(fbs = case_when(fbs == 1 ~ "high",
fbs == 0 ~ "low")) %>%
mutate(exang = case_when(exang == 0 ~ "no",
exang == 1 ~ "yes")) %>%
mutate(heart_disease = case_when(heart_disease == 1 ~ "absence",
TRUE ~ "presence"))
heart <- heart %>%
mutate(sex = as.factor(sex)) %>%
mutate(cp = as.factor(cp)) %>%
mutate(fbs = as.factor(fbs)) %>%
mutate(exang = as.factor(exang)) %>%
mutate(heart_disease = as.factor(heart_disease))
heart <- heart %>%
select(age, sex, cp, trestbps, chol, fbs, thalach, exang, heart_disease) %>%
rename("max_hr" = "thalach",
"exercise_angina" = "exang") %>%
drop_na()
ggpairs(heart,
columns = 1:9,
lower=list(combo=wrap("facethist", binwidth=0.5)))summary(heart)## age sex cp trestbps
## Min. :29.00 female: 96 asymptomatic : 23 Min. : 94.0
## 1st Qu.:47.50 male :207 atypical angina : 50 1st Qu.:120.0
## Median :55.00 non-anginal pain: 87 Median :130.0
## Mean :54.37 typical angina :143 Mean :131.6
## 3rd Qu.:61.00 3rd Qu.:140.0
## Max. :77.00 Max. :200.0
## chol fbs max_hr exercise_angina heart_disease
## Min. :126.0 high: 45 Min. : 71.0 no :204 absence :165
## 1st Qu.:211.0 low :258 1st Qu.:133.5 yes: 99 presence:138
## Median :240.0 Median :153.0
## Mean :246.3 Mean :149.6
## 3rd Qu.:274.5 3rd Qu.:166.0
## Max. :564.0 Max. :202.0
En el siguiente grafico se puede ver como afectan la presencia de tener o no tener enfermedad del corazon en funcion de “trestbps” (presion de sangre en reposo) y “chol” (colesterol). Podemos ver una tendencia que para altos valores de presion de sangre en reposo aumenta la presencia de enfermedad cardiaca. Por otro lado el colesterol no parece tener influencia en la presencia de enfermedad cardiaca.
gg.plot <- ggplot(heart, aes(x=chol, y=trestbps)) +
geom_point(aes(col=heart_disease))
gg.plotEn el primer grafico abajo vemos que el grafico de ausencia de enfermedad tiene una distribucion centrada mientras que el de presencia de enfermadad, grafico de la derecha, sugiere que mayor frecuencia de presencia de enfermedad para personas de edad mas avanzada
En el segundo grafico podemos ver que la angina tipica tiene una gran influencia en la presencia de la enfermedad.
En el tercer grafico vemos que para los hombres se presenta con mayor frecuencia la enfermedad como asi tambien para la angina inducida por el ejercicio. Mientras que el Fasting Blood Sugar o glucemia en ayunas no presenta tendencia.
Finalmente en el último grafico tipo boxplot, vemos que la variable que marca una tendencia con la aparición de la enfermedad es la de máxima pulsaciones por minuto alcanzadas durante el ejercicio, para aquellas personas que alcanzan pulsaciones mas elevadas son de tendencia sana mientras que los que alcanzan menos pulsaciones tienden a tener la enfermedad cardíaca.
age.plot <- ggplot(heart, mapping = aes(x = age, fill = heart_disease)) +
geom_histogram() +
facet_wrap(vars(heart_disease)) + stat_bin(bins = 50) +
labs(title = "Presencia de la enfermedad por edad", x = "Age (years)", y = "Count", fill = "Heart Disease")
age.plotcp.plot <- ggplot(heart, mapping = aes(x=heart_disease, fill = cp)) +
geom_bar(position = "dodge") +
labs(title = "Presencia enfermedad para distintos tipos de dolores de pecho - Chest Pain", x = "Heart Disease", y = "Count", fill = "Chest Pain Type")
cp.plotsex.plot <- ggplot(heart, mapping = aes(x = sex, fill = heart_disease)) +
geom_bar(position = "fill") +
labs(x = "Sex", y = "Proportion", fill = "Heart Disease") +
theme(axis.text.x = element_text(size = 12), axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12), axis.text.y = element_text(size = 12))
fbs.plot <- ggplot(heart, mapping = aes(x=fbs, fill=heart_disease)) +
geom_bar(position = "fill") +
labs(x = "Fasting Blood Sugar", y = "Proportion", fill = "Heart Disease") +
scale_x_discrete(labels = c("low", "high"))+
theme(axis.text.x = element_text(size = 12), axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12), axis.text.y = element_text(size = 12))
exang.plot <- ggplot(heart, mapping = aes(x = exercise_angina, fill = heart_disease)) +
geom_bar(position = "fill") +
labs(x = "Exercise induced angina", y = "Proportion", fill = "Heart Disease") +
theme(axis.text.x = element_text(size = 12), axis.title.x = element_text(size = 12))
grid.arrange(sex.plot, fbs.plot, exang.plot, nrow=2)trestbps.plot <- ggplot(heart, mapping = aes(x=trestbps, y=heart_disease)) +
geom_boxplot() +
labs(x = "Resting Blood Pressure (mm Hg)", y = "Heart Disease") +
theme(axis.text.x = element_text(size = 12), axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12), axis.text.y = element_text(size = 12))
chol.plot <- ggplot(heart, mapping = aes(x=chol, y=heart_disease)) +
geom_boxplot() +
labs(x = "Serum Cholestoral (mg/dl)", y = "Heart Disease") +
theme(axis.text.x = element_text(size = 12), axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12), axis.text.y = element_text(size = 12))
maxhr.plot <- ggplot(heart, mapping = aes(x = max_hr, y = heart_disease)) +
geom_boxplot() +
labs(x = "Maximum Heart Rate (bpm)", y = "Heart Disease") +
theme(axis.text.x = element_text(size = 12), axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12), axis.text.y = element_text(size = 12))
grid.arrange(trestbps.plot, chol.plot, maxhr.plot, nrow=2)Una regresión logística se utiliza típicamente cuando hay una variable de resultado dicotómica (como ganar o perder), y una variable predictiva continua que está relacionada con la probabilidad de la variable de resultado. También se puede utilizar con predictores categóricos y con múltiples predictores.
La funcion logıstica es una funcion acotada entre 0 y 1 dada por:
\(\Large {f(t)} = \frac{1}{1+e^{-t}}\)
Para la construccion de modelo y obtencion de los parametros usaremos la fucion “glm” mediante el argumento family = “binomial”
El modelo logistico sera:
\(\Large {p(x)} = \frac{1}{1+e^{-(\beta_0+\beta_1x+...+\beta_{p-1}x_{p-1})}}\)
Finalmente los valores pronosticados del modelo pueden convertirse a probabilidades de evento de la siguiente manera:
\(\Large {p_i} = 1-\frac{1}{1+e^{-t}}\)
La matriz de confusión permite evaluar la calidad del modelo. Para esto se usaran dos tipos de conjuntos de datos, uno de entrenamiento que su función es la de input para encontrar los parametros del modelo, y un conjunto de test que es el que se usara para obtener los valores de predicción y asi poder comparalos con los valores datos.
Matriz de Confusión.
Para poder hacer los comparativos de los modelos vamos a usar la matriz de confusión y asi poder seleccionar el modelo que mejor predice la presencia y/o ausencia de enfermedad cariaca.
A partir de la tabla podemos observar que el error de clasificación empirico se calcula como L=(FN+FP)/n, con n el tamaño de la muestra (n=VP+VN+FP+FN)
Desde el punto de vista de sanos y enfermos:
La sensibilidad nos indica la capacidad de nuestro estimador para dar como casos positivos los casos realmente enfermos; proporci´on de enfermos correctamente identificados. Es decir, la sensibilidad caracteriza la capacidad de la prueba para detectar la enfermedad en sujetos enfermos.
La especificidad nos indica la capacidad de nuestro estimador para dar como casos negativos los casos realmente sanos; proporci´on de sanos correctamente identificados. Es decir, la especificidad caracteriza la capacidad de la prueba para detectar la ausencia de la enfermedad en sujetos sanos.
En este modelos se toman las 8 variables independientes para la construcción del modelo. Las variables significativas encontradas para esta regresión logística son:
Las variables que no tienen significancia en la predicción de enfermedad del corazón son:
## 75% of the sample size
smp_size <- floor(0.75 * nrow(heart))
## set the seed to make your partition reproducible
set.seed(123)
train_ind <- sample(seq_len(nrow(heart)), size = smp_size)
heart.train <- heart[train_ind, ]
heart.test <- heart[-train_ind, ]
heart.full <- glm(heart_disease~., data = heart.train, family = binomial)
summary(heart.full)##
## Call:
## glm(formula = heart_disease ~ ., family = binomial, data = heart.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7106 -0.6727 -0.2480 0.6309 2.2184
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.579922 2.646483 -1.353 0.176149
## age 0.025563 0.023551 1.085 0.277729
## sexmale 1.776222 0.446005 3.983 6.82e-05 ***
## cpatypical angina -0.222318 0.757590 -0.293 0.769175
## cpnon-anginal pain -0.054859 0.655676 -0.084 0.933320
## cptypical angina 1.724031 0.626603 2.751 0.005934 **
## trestbps 0.025979 0.010497 2.475 0.013324 *
## chol 0.005381 0.003663 1.469 0.141792
## fbslow 0.105204 0.493648 0.213 0.831237
## max_hr -0.034879 0.010481 -3.328 0.000875 ***
## exercise_anginayes 0.908019 0.413920 2.194 0.028257 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 312.74 on 226 degrees of freedom
## Residual deviance: 196.28 on 216 degrees of freedom
## AIC: 218.28
##
## Number of Fisher Scoring iterations: 5
#
# Matriz de Confusion para el Train
pdata <- predict(heart.full, heart.train, type="response")
pfn = as.factor(as.numeric(pdata>0.5))
levels(pfn) <- c("absence", "presence")
dfn = as.factor(heart.train$heart_disease)
#confusionMatrix(data = pfn, reference = dfn)
#
#Matriz de Confusion para el Test
pdata <- predict(heart.full, heart.test, type="response")
pfn = as.factor(as.numeric(pdata>0.5))
levels(pfn) <- c("absence", "presence")
dfn = as.factor(heart.test$heart_disease)
confusionMatrix(data = pfn, reference = dfn)## Confusion Matrix and Statistics
##
## Reference
## Prediction absence presence
## absence 37 9
## presence 4 26
##
## Accuracy : 0.8289
## 95% CI : (0.7253, 0.9057)
## No Information Rate : 0.5395
## P-Value [Acc > NIR] : 1.083e-07
##
## Kappa : 0.6521
##
## Mcnemar's Test P-Value : 0.2673
##
## Sensitivity : 0.9024
## Specificity : 0.7429
## Pos Pred Value : 0.8043
## Neg Pred Value : 0.8667
## Prevalence : 0.5395
## Detection Rate : 0.4868
## Detection Prevalence : 0.6053
## Balanced Accuracy : 0.8226
##
## 'Positive' Class : absence
##
#Vamos ahora a realizar un modelo de regresion logistica removiendo las variables encontradas que no tienen significancia: age, chol y fbs. Vamos a llamar a este modelo heart.fit Podemos ver ahora que todas las variables del modelo son significativas.
heart_fit <- glm(heart_disease~. - age - chol -fbs, data = heart.train, family = binomial)
summary(heart_fit)##
## Call:
## glm(formula = heart_disease ~ . - age - chol - fbs, family = binomial,
## data = heart.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5784 -0.6757 -0.2684 0.6084 2.2654
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.733402 2.023038 -0.363 0.716960
## sexmale 1.506312 0.412218 3.654 0.000258 ***
## cpatypical angina -0.201655 0.745065 -0.271 0.786656
## cpnon-anginal pain -0.099607 0.641681 -0.155 0.876641
## cptypical angina 1.710198 0.615522 2.778 0.005462 **
## trestbps 0.029821 0.009903 3.011 0.002601 **
## max_hr -0.036887 0.009628 -3.831 0.000127 ***
## exercise_anginayes 0.860198 0.403556 2.132 0.033044 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 312.74 on 226 degrees of freedom
## Residual deviance: 200.48 on 219 degrees of freedom
## AIC: 216.48
##
## Number of Fisher Scoring iterations: 5
#
# Matriz de Confusion para el Train
pdata <- predict(heart_fit, heart.train, type="response")
pfn = as.factor(as.numeric(pdata>0.5))
levels(pfn) <- c("absence", "presence")
dfn = as.factor(heart.train$heart_disease)
#confusionMatrix(data = pfn, reference = dfn)
#
#Matriz de Confusion para el Test
pdata <- predict(heart_fit, heart.test, type="response")
pfn = as.factor(as.numeric(pdata>0.5))
levels(pfn) <- c("absence", "presence")
dfn = as.factor(heart.test$heart_disease)
confusionMatrix(data = pfn, reference = dfn)## Confusion Matrix and Statistics
##
## Reference
## Prediction absence presence
## absence 37 8
## presence 4 27
##
## Accuracy : 0.8421
## 95% CI : (0.7404, 0.9157)
## No Information Rate : 0.5395
## P-Value [Acc > NIR] : 2.51e-08
##
## Kappa : 0.6796
##
## Mcnemar's Test P-Value : 0.3865
##
## Sensitivity : 0.9024
## Specificity : 0.7714
## Pos Pred Value : 0.8222
## Neg Pred Value : 0.8710
## Prevalence : 0.5395
## Detection Rate : 0.4868
## Detection Prevalence : 0.5921
## Balanced Accuracy : 0.8369
##
## 'Positive' Class : absence
##
#Vamos ahora a realizar un 3er modelo sacando la variable “cp” que tambien parece tener una significacia relativamente baja. Vamos a llamar a este modelo hert.red
heart.red <- glm(heart_disease~. - age - chol -fbs - cp, data = heart.train, family = binomial)
summary(heart.red)##
## Call:
## glm(formula = heart_disease ~ . - age - chol - fbs - cp, family = binomial,
## data = heart.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4673 -0.7687 -0.3957 0.7484 2.2525
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.435230 1.737103 0.826 0.408679
## sexmale 1.338711 0.372280 3.596 0.000323 ***
## trestbps 0.024924 0.009161 2.721 0.006516 **
## max_hr -0.042189 0.008809 -4.789 1.67e-06 ***
## exercise_anginayes 1.460648 0.359264 4.066 4.79e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 312.74 on 226 degrees of freedom
## Residual deviance: 224.90 on 222 degrees of freedom
## AIC: 234.9
##
## Number of Fisher Scoring iterations: 4
#
# Matriz de Confusion para el Train
pdata <- predict(heart.red, heart.train, type="response")
pfn = as.factor(as.numeric(pdata>0.5))
levels(pfn) <- c("absence", "presence")
dfn = as.factor(heart.train$heart_disease)
#confusionMatrix(data = pfn, reference = dfn)
#
#Matriz de Confusion para el Test
pdata <- predict(heart.red, heart.test, type="response")
pfn = as.factor(as.numeric(pdata>0.5))
levels(pfn) <- c("absence", "presence")
dfn = as.factor(heart.test$heart_disease)
confusionMatrix(data = pfn, reference = dfn)## Confusion Matrix and Statistics
##
## Reference
## Prediction absence presence
## absence 34 11
## presence 7 24
##
## Accuracy : 0.7632
## 95% CI : (0.6518, 0.8532)
## No Information Rate : 0.5395
## P-Value [Acc > NIR] : 4.725e-05
##
## Kappa : 0.5193
##
## Mcnemar's Test P-Value : 0.4795
##
## Sensitivity : 0.8293
## Specificity : 0.6857
## Pos Pred Value : 0.7556
## Neg Pred Value : 0.7742
## Prevalence : 0.5395
## Detection Rate : 0.4474
## Detection Prevalence : 0.5921
## Balanced Accuracy : 0.7575
##
## 'Positive' Class : absence
##
#Para medir la calidad global de la prueba mediante el área, se usa tanto la curva ROC como el área bajo ella (AUC). El área debajo de esta curva mide la capacidad o habilidad que tiene el modelo para discriminar.
Desde el punto de vista de diagnóstico de enfermedades, ll análisis de curvas ROC (receiver operating characteristic curve) constituye un método estadístico para determinar la exactitud diagnóstica de estos tests, siendo utilizadas con tres propósitos específicos: determinar el punto de corte de una escala continua en el que se alcanza la sensibilidad y especificidad más alta, evaluar la capacidad discriminativa del test diagnóstico, es decir, su capacidad de diferenciar sujetos sanos versus enfermos, y comparar la capacidad discriminativa de dos o más tests diagnósticos que expresan sus resultados como escalas continuas.
heart.test.pred = predict(heart.full, heart.test, type = "response")
testcomp <- data.frame(heart.test$heart_disease, heart.test.pred)
colnames(testcomp) <- c("test.response", "test.prediction")
testcomp$test.prediction <- testcomp$test.prediction > 0.5
testcomp <- testcomp %>%
mutate(test.response = factor(case_when(test.response == "absence" ~ 0,
test.response == "presence" ~ 1))) %>%
mutate(test.prediction = factor(case_when(test.prediction == FALSE ~ 0,
test.prediction == TRUE ~ 1)))
heart.roc <- roc(response = ordered(testcomp$test.response), predictor = ordered(testcomp$test.prediction))## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(heart.roc, print.thres = "best", main = "Receiver Operating Characteritic Technique Plot")print(auc(heart.roc))## Area under the curve: 0.8226
heart.test.pred = predict(heart_fit, heart.test, type = "response")
testcomp <- data.frame(heart.test$heart_disease, heart.test.pred)
colnames(testcomp) <- c("test.response", "test.prediction")
testcomp$test.prediction <- testcomp$test.prediction > 0.5
testcomp <- testcomp %>%
mutate(test.response = factor(case_when(test.response == "absence" ~ 0,
test.response == "presence" ~ 1))) %>%
mutate(test.prediction = factor(case_when(test.prediction == FALSE ~ 0,
test.prediction == TRUE ~ 1)))
heart.roc <- roc(response = ordered(testcomp$test.response), predictor = ordered(testcomp$test.prediction))## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(heart.roc, print.thres = "best", main = "Receiver Operating Characteritic Technique Plot")print(auc(heart.roc))## Area under the curve: 0.8369
heart.test.pred = predict(heart.red, heart.test, type = "response")
testcomp <- data.frame(heart.test$heart_disease, heart.test.pred)
colnames(testcomp) <- c("test.response", "test.prediction")
testcomp$test.prediction <- testcomp$test.prediction > 0.5
testcomp <- testcomp %>%
mutate(test.response = factor(case_when(test.response == "absence" ~ 0,
test.response == "presence" ~ 1))) %>%
mutate(test.prediction = factor(case_when(test.prediction == FALSE ~ 0,
test.prediction == TRUE ~ 1)))
heart.roc <- roc(response = ordered(testcomp$test.response), predictor = ordered(testcomp$test.prediction))## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(heart.roc, print.thres = "best", main = "Receiver Operating Characteritic Technique Plot")print(auc(heart.roc))## Area under the curve: 0.7575
Podemos ver que del anáalisis de las curvas ROC, similarmente a lo que vimos en las Tablas de Confusión, el Modelo Optimizado es el mejor ya que arroja el mayor AUC (0.837) y los valores de Especificidad y Sensibilidad que vemos coinciden con los calculados en las tablas de confusión.
Del análisis de los tres modelos logísticos desarrollados, podemos concluir que el mejor modelo es el Modelo Optimizado. Esta conclusión se basa en las siguientes observaciones a partir del análisis de las matrices de confusión y Curvas ROC para los valores de predicción del conjunto de datos usados como test.
https://rpubs.com/Edenner/854339 (Heart Disease Prediction)>.
https://halimatusyak.medium.com/heart-disease-uci-logistic-regression-in-r-b95b821088e6