# Enfermedad cardiaca y frecuencia cardiaca
# https://rpubs.com/JOJOQUIMO/Heart_Disease_and_thalach
#C:\Users\HP\Downloads\Heart_Disease_Prediction.xlsx
datos <- read.csv("Heart_Disease_Prediction.csv", head=T, sep=";")
str(datos)
## 'data.frame': 270 obs. of 14 variables:
## $ Age : int 70 67 57 64 74 65 56 59 60 63 ...
## $ Sex : int 1 0 1 1 0 1 1 1 1 0 ...
## $ cp : int 4 3 2 4 2 4 3 4 4 4 ...
## $ trestbps : int 130 115 124 128 120 120 130 110 140 150 ...
## $ Chol : int 322 564 261 263 269 177 256 239 293 407 ...
## $ fbs : int 0 0 0 0 0 0 1 0 0 0 ...
## $ EKG : int 2 2 0 0 2 0 2 2 2 2 ...
## $ thalach : int 109 160 141 105 121 140 142 142 170 154 ...
## $ exang : int 0 0 0 1 1 0 1 1 0 0 ...
## $ oldpeak : chr "2,4" "1,6" "0,3" "0,2" ...
## $ Slope : int 2 2 1 2 1 1 2 2 2 2 ...
## $ ca : int 3 0 0 1 1 0 1 1 2 3 ...
## $ Thal : int 3 7 7 7 3 7 6 7 7 7 ...
## $ Heart.Disease: int 1 0 1 0 0 0 1 1 1 1 ...
# convirtiendo a factor la variable target y etiquetando
datos$Heart.Disease <- factor(datos$Heart.Disease,levels = c("0","1"),
labels=c("no", "yes"))
str(datos)
## 'data.frame': 270 obs. of 14 variables:
## $ Age : int 70 67 57 64 74 65 56 59 60 63 ...
## $ Sex : int 1 0 1 1 0 1 1 1 1 0 ...
## $ cp : int 4 3 2 4 2 4 3 4 4 4 ...
## $ trestbps : int 130 115 124 128 120 120 130 110 140 150 ...
## $ Chol : int 322 564 261 263 269 177 256 239 293 407 ...
## $ fbs : int 0 0 0 0 0 0 1 0 0 0 ...
## $ EKG : int 2 2 0 0 2 0 2 2 2 2 ...
## $ thalach : int 109 160 141 105 121 140 142 142 170 154 ...
## $ exang : int 0 0 0 1 1 0 1 1 0 0 ...
## $ oldpeak : chr "2,4" "1,6" "0,3" "0,2" ...
## $ Slope : int 2 2 1 2 1 1 2 2 2 2 ...
## $ ca : int 3 0 0 1 1 0 1 1 2 3 ...
## $ Thal : int 3 7 7 7 3 7 6 7 7 7 ...
## $ Heart.Disease: Factor w/ 2 levels "no","yes": 2 1 2 1 1 1 2 2 2 2 ...
# Representación de las observaciones
require(ggplot2)
## Cargando paquete requerido: ggplot2
require(gridExtra)
## Cargando paquete requerido: gridExtra
require(caret)
## Cargando paquete requerido: caret
## Cargando paquete requerido: lattice
require(tidyverse)
## Cargando paquete requerido: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::lift() masks caret::lift()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
datos %>%
group_by(datos$Heart.Disease) %>%
summarise( numero_casos = n(),
porcentaje = numero_casos / nrow(datos)
)
## # A tibble: 2 × 3
## `datos$Heart.Disease` numero_casos porcentaje
## <fct> <int> <dbl>
## 1 no 150 0.556
## 2 yes 120 0.444
qplot(Heart.Disease, data = datos, col = Heart.Disease,
main = " ",geom = "bar") +
geom_text(aes(label=scales::percent(..count../sum(..count..))),
stat='count',position=position_stack(0.5))+
geom_text(aes(label=..count..),
stat="count",position=position_stack())
## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Histograma de la variable cuantitativa
ggplot(datos, aes(x=thalach))+
geom_histogram(color="steelblue", fill="lightblue")+
ggtitle("Histograma de notas de matemáticas")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# grafico de cajas
p1 <- ggplot(data = datos, aes(x = thalach, y = Heart.Disease)) +
geom_point(aes(color = Heart.Disease)) +
theme_bw() + theme(legend.position = "null")
p2 <- ggplot(data = datos, aes(x = Heart.Disease, y = thalach , color=Heart.Disease)) +
geom_boxplot(outlier.shape=NA) +
geom_jitter (width=0.1) +
theme_bw() +
theme(legend.position = "null")
grid.arrange(p1, p2, nrow = 1)
# Estadísticos descriptivos
datos %>% filter(!is.na(thalach)) %>% group_by(Heart.Disease) %>%
summarise(media = mean(thalach),
mediana = median(thalach),
min = min(thalach),
max = max(thalach))
## # A tibble: 2 × 5
## Heart.Disease media mediana min max
## <fct> <dbl> <dbl> <int> <int>
## 1 no 158. 161 96 202
## 2 yes 139. 142. 71 195
# estimación del modelo
modelo <- glm(Heart.Disease ~ thalach, data = datos, family = "binomial")
summary(modelo)
##
## Call:
## glm(formula = Heart.Disease ~ thalach, family = "binomial", data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.159089 1.020518 6.035 1.59e-09 ***
## thalach -0.042753 0.006767 -6.318 2.64e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 370.96 on 269 degrees of freedom
## Residual deviance: 320.04 on 268 degrees of freedom
## AIC: 324.04
##
## Number of Fisher Scoring iterations: 4
# exponenciando para facilitar la interpretación de los coeficientes
round(exp(coefficients(modelo)),6)
## (Intercept) thalach
## 472.996901 0.958148
# Intervalos de Confianza para los coeficientes
round(confint(modelo, level=0.95),5)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 4.23925 8.25190
## thalach -0.05662 -0.03002
# Tasa de ventajas e IC 95%
round(exp(cbind(OR = coef(modelo), confint(modelo,level=0.95))),5)
## Waiting for profiling to be done...
## OR 2.5 % 97.5 %
## (Intercept) 472.99690 69.35609 3834.91473
## thalach 0.95815 0.94496 0.97043
# El mismo calculo se puede obtener directamente con:
anova(modelo, test = "Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Heart.Disease
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 269 370.96
## thalach 1 50.919 268 320.04 9.627e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(ResourceSelection)
## ResourceSelection 0.3-6 2023-06-27
hoslem.test(datos$Heart.Disease, fitted(modelo))
## Warning in Ops.factor(1, y): '-' no es significativo para factores
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: datos$Heart.Disease, fitted(modelo)
## X-squared = NA, df = 8, p-value = NA
# medidas tipo R2
# Pseudo R2 de McFadden
library(pscl)
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
pscl :: pR2 (modelo) ["McFadden"]
## fitting null model for pseudo-r2
## McFadden
## 0.1372624
# medidas tipo R2
# Pseudo R2 de McFadden
(RsqrMcFadden <- 1 - modelo$deviance/modelo$null.deviance)
## [1] 0.1372624
# Pseudo R2 de Cox y Snell
LR <- modelo$null.deviance - modelo$deviance
N <- sum(weights(modelo))
(RsqrCN <- 1 - exp(-LR/N))
## [1] 0.1718724
# Pseudo R2de Nlkerke
L0.adj <- exp(-modelo$null.deviance/N)
(RsqrNal <- RsqrCN/(1 - L0.adj))
## [1] 0.230118
# criterio AIC (para comparar modelos)
AIC(modelo)
## [1] 324.0405
# BIC (para comparar modelos) criterio de inform, bayesiana
BIC(modelo)
## [1] 331.2373
# CURVA ROC Y AREA BAJO LA CURVA
library(stats)
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Adjuntando el paquete: 'pROC'
##
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
prob=predict(modelo,type="response")
rocA1C <- with(datos,roc(Heart.Disease, prob))
## Setting levels: control = no, case = yes
## Setting direction: controls < cases
plot(rocA1C, col="red", print.auc=TRUE)
closest <- coords(rocA1C,"b",ret=c("threshold","specificity","sensitivity","npv","ppv"),
best.method="closest.topleft")
closest
## threshold specificity sensitivity npv ppv
## threshold 0.4633215 0.7666667 0.6333333 0.7232704 0.6846847
#Calculando area bajo la curva
areaROC<-auc(roc(datos$Heart.Disease,prob))
## Setting levels: control = no, case = yes
## Setting direction: controls < cases
ROC<-plot.roc(datos$Heart.Disease,prob, xlab="1- Especificidad", ylab="Sensibilidad",
main = paste('Area Bajo la Curva =',round(areaROC,4)), col="blue")
## Setting levels: control = no, case = yes
## Setting direction: controls < cases
#library(InformationValue)
#plotROC(sjlabelled::as_numeric(datos$Heart.Disease,start.at = 0), prob)
# probabilidades y grupos estimados
prob=predict(modelo,type="response")
head(prob,20)
## 1 2 3 4 5 6 7 8
## 0.8174253 0.3359470 0.5326784 0.8415775 0.7282897 0.5433046 0.5220226 0.5220226
## 9 10 11 12 13 14 15 16
## 0.2480686 0.3953472 0.3264772 0.8043164 0.1770489 0.4899733 0.3455505 0.6931633
## 17 18 19 20
## 0.7366668 0.3851742 0.5006603 0.1898531
library(vcd)
## Cargando paquete requerido: grid
predicciones <- ifelse(test = modelo$fitted.values > 0.5, yes = 1, no = 0)
head(predicciones,20)
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 0 1 1 1 1 1 1 0 0 0 1 0 0 0 1 1 0 1 0
matriz_confusion <- table(modelo$model$Heart.Disease, predicciones,
dnn = c("observaciones", "predicciones"))
matriz_confusion
## predicciones
## observaciones 0 1
## no 119 31
## yes 52 68
# Convertir Heart.Disease a factor con niveles 0 y 1
if (!is.factor(datos$Heart.Disease)) {
datos$Heart.Disease <- as.factor(datos$Heart.Disease)
}
levels(datos$Heart.Disease) <- c(0, 1) # Asegurar los niveles correctos
# Convertir las predicciones a factor con los mismos niveles
predicciones <- factor(predicciones, levels = c(0, 1))
# Crear la matriz de confusión usando la función confusionMatrix
matriz_confusion <- confusionMatrix(predicciones, datos$Heart.Disease)
# Extraer la tabla de contingencia de la matriz de confusión
tabla_contingencia <- matriz_confusion$table
# Calcular el error 1
error1 <- sum(tabla_contingencia[1, 2], tabla_contingencia[2, 1]) / sum(tabla_contingencia)
print(error1)
## [1] 0.3074074
# Preparar la base de datos para guardar valores predichos y probabilidades
prob <- predict(modelo, type = "response")
finaldata <- cbind(datos, prob, predicciones)
head(finaldata, 20)
## Age Sex cp trestbps Chol fbs EKG thalach exang oldpeak Slope ca Thal
## 1 70 1 4 130 322 0 2 109 0 2,4 2 3 3
## 2 67 0 3 115 564 0 2 160 0 1,6 2 0 7
## 3 57 1 2 124 261 0 0 141 0 0,3 1 0 7
## 4 64 1 4 128 263 0 0 105 1 0,2 2 1 7
## 5 74 0 2 120 269 0 2 121 1 0,2 1 1 3
## 6 65 1 4 120 177 0 0 140 0 0,4 1 0 7
## 7 56 1 3 130 256 1 2 142 1 0,6 2 1 6
## 8 59 1 4 110 239 0 2 142 1 1,2 2 1 7
## 9 60 1 4 140 293 0 2 170 0 1,2 2 2 7
## 10 63 0 4 150 407 0 2 154 0 4 2 3 7
## 11 59 1 4 135 234 0 0 161 0 0,5 2 0 7
## 12 53 1 4 142 226 0 2 111 1 0 1 0 7
## 13 44 1 3 140 235 0 2 180 0 0 1 0 3
## 14 61 1 1 134 234 0 0 145 0 2,6 2 2 3
## 15 57 0 4 128 303 0 2 159 0 0 1 1 3
## 16 71 0 4 112 149 0 0 125 0 1,6 2 0 3
## 17 46 1 4 140 311 0 0 120 1 1,8 2 2 7
## 18 53 1 4 140 203 1 2 155 1 3,1 3 0 7
## 19 64 1 1 110 211 0 2 144 1 1,8 2 0 3
## 20 40 1 1 140 199 0 0 178 1 1,4 1 0 7
## Heart.Disease prob predicciones
## 1 1 0.8174253 1
## 2 0 0.3359470 0
## 3 1 0.5326784 1
## 4 0 0.8415775 1
## 5 0 0.7282897 1
## 6 0 0.5433046 1
## 7 1 0.5220226 1
## 8 1 0.5220226 1
## 9 1 0.2480686 0
## 10 1 0.3953472 0
## 11 0 0.3264772 0
## 12 0 0.8043164 1
## 13 0 0.1770489 0
## 14 1 0.4899733 0
## 15 0 0.3455505 0
## 16 0 0.6931633 1
## 17 1 0.7366668 1
## 18 1 0.3851742 0
## 19 0 0.5006603 1
## 20 0 0.1898531 0
# Predicción para nuevos individuos
nuevo1 <- data.frame(thalach = 60)
print(predict(modelo, newdata = nuevo1, type = "response"))
## 1
## 0.9732447
nuevo1 <- data.frame(thalach = 20)
print(predict(modelo, newdata = nuevo1, type = "response"))
## 1
## 0.9950531
nuevo1 <- data.frame(thalach = 63)
print(predict(modelo, newdata = nuevo1, type = "response"))
## 1
## 0.9696942
# Gráfico del modelo: SIN INTERVALOS DE CONFIANZA
# Convirtiendo en numerico (0,1) la variable respuesta
datos$Heart.Disease <- as.character(datos$Heart.Disease)
datos$Heart.Disease <- as.numeric(datos$Heart.Disease)
plot(Heart.Disease ~ thalach, datos, col = "darkblue",
main = "Modelo regresion logistica Heart.Disease ~ nota thalach",
ylab = "P(Heart.Disease = 1 | thalach)", xlab = "thalach", pch = 19)
# type = 'response' devuelve las predicciones en forma de probabilidad en
# lugar de en log_ODDs
curve(predict(modelo, data.frame(thalach = x), type = "response"),
col = "firebrick", lwd = 2.5, add = TRUE)
# MEDIANTE GGPLOT2 INCLUYENDO INTERVALOS DE CONFIANZA
nuevos_puntos <- seq(from = min(datos$thalach), to = max(datos$thalach), by = 0.5)
# Predicciones de los nuevos puntos segun el modelo
predicciones <- predict(modelo, data.frame(thalach = nuevos_puntos), se.fit = TRUE)
# Mediante la funcion logit se transforman los log_ODDs a probabilidades
predicciones_logit <- exp(predicciones$fit) / (1 + exp(predicciones$fit))
# Se calcula el limite inferior y superior del IC del 95%
limite_inferior <- predicciones$fit - 1.96 * predicciones$se.fit
limite_inferior_logit <- exp(limite_inferior) / (1 + exp(limite_inferior))
limite_superior <- predicciones$fit + 1.96 * predicciones$se.fit
limite_superior_logit <- exp(limite_superior) / (1 + exp(limite_superior))
# Se crea un data frame con los nuevos puntos y sus predicciones
datos_curva <- data.frame(thalach = nuevos_puntos,
probabilidad_Heart.Disease = predicciones_logit,
limite_inferior_logit = limite_inferior_logit,
limite_superior_logit = limite_superior_logit)
ggplot(datos, aes(x = thalach, y = Heart.Disease)) +
geom_point(aes(color = as.factor(Heart.Disease)), shape = 19, size = 3) +
geom_line(data = datos_curva, aes(y = probabilidad_Heart.Disease), color = "firebrick") +
geom_line(data = datos_curva, aes(y = limite_inferior_logit), linetype = "dashed") +
geom_line(data = datos_curva, aes(y = limite_superior_logit), linetype = "dashed") +
theme_bw() +
labs(title = "Modelo regresion logistica Heart.Disease ~ nota thalach",
y = "P(Heart.Disease = 1 | thalach)", x = "thalach") +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5))
# Modelo gausiano
modelo_gaussiano <- glm(Heart.Disease ~ thalach, data = datos, family = "gaussian")
summary(modelo_gaussiano)
##
## Call:
## glm(formula = Heart.Disease ~ thalach, family = "gaussian", data = datos)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.790614 0.180564 9.917 < 2e-16 ***
## thalach -0.008994 0.001192 -7.544 7.12e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.2051856)
##
## Null deviance: 66.667 on 269 degrees of freedom
## Residual deviance: 54.990 on 268 degrees of freedom
## AIC: 342.58
##
## Number of Fisher Scoring iterations: 2