Uso de Regresión Logística para la asociación de factores y predicción del riesgo de tener una enfermedad coronaria

Universidad Nacional de Colombia - Sede Bogotá

Kevin Duván Barrera Pezca, Daniela Bohórquez, Deivis Cárdenas, Rodolfo Alejandro García, Joan Lamprea

abril 21, 2023

Introducción

La Organización Mundial de la Salud ha estimado que 12 millones de muertes se producen al año a cause de Cardiopatía. La mitad de muertes en Estados Unidos y otros países desarrollados se deben a problemas cardio vasculares. El pronóstico temprano de problemas cardiovasculares puede aportar en la toma de decisiones para generar cambios de hábitos en pacientes de alto riesgo y en vías de reducir sus complicaciones. Esta investigación busca indicar los factores más relevantes o riesgosos relacionados con la cardiopatía, asimismo predecir el riesgo general de contraer la enfermedad bajo exposición a dichos factores.

Estos datos provienen de Kaggle, son de libre acceso, corresponden a un estudio cardiovascular en habitantes de el pueblo de Framingham, Massachusetts. La variable respuesta que mediremos en este estudio es el diagnóstico de poder tener en los próximos 10 años una enfermedad coronaria. Estos datos contienen más de 4000 registros y 15 variables.

Nuestro interés es mostrar lo presentado en la metodología del artículo y las diferencias entre dos modelos: uno de asociación y uno de predicción, con el fin de mostrar las particularidades de ambas metodologías.

Nota: Los datos que analizamos no provienen del artículo que escogimos debido a que no fue posible adquirir el conjunto de datos usado para la implementación, por esta razón optamos por usar otro conjunto de datos.

Variables del estudio

Variable respuesta

TenYearCHD: Riesgo detectado a 10 años de posible enfermedad coronoria.

Variables explicativas

Male: Hombre(1) o Mujer (1)
Age: Variable Continua de edad
CurrentSmoker: Actual fumador (1), e.o.c.(0)
CigsPerDay: Cantidad de cigarrillos que se consumen por día
BPMeds: Medicación para Presión Sanguinea (1), e.o.c. (0)
PrevalentStroke: Ataque Cardiaco Previo (1), e.o.c. (0)
PrevalentHyp: Hipertensión Diagnosticada (1), e.o.c. (0)
Diabetes: Diabetes Diagnosticada (1), e.o.c. (0)
TotChol: Nivel de Colesterol, variable continua
SysBP: Presión Sistólica, variable continua
DiaBP: Presión Diastólica, variable continua
BMI: Índice de Masa Corporal, variable continua
HeartRate: Ritmo Cardiaco, variable continua
Glucose: Nivel de glucosa en la sangre, variable continua

Análisis de la información

Los datos están contenidos en el archivo framingham.csv, vamos a importar los datos y omitiremos todos los registros con datos faltantes.

library(readr)
heart_dis <- read_csv("framingham.csv", col_types = cols(male = col_factor(levels = c("0","1")), education = col_skip(), currentSmoker = col_factor(levels = c("0", "1")), BPMeds = col_factor(levels = c("0", "1")), prevalentStroke = col_factor(levels = c("0", "1")), prevalentHyp = col_factor(levels = c("0", "1")), diabetes = col_factor(levels = c("0", "1")), TenYearCHD = col_factor(levels = c("0", "1"))))
heart_dis<- na.omit(heart_dis)

Aquí pueden visualizarse los primeros 10 registros del conjunto de datos.

knitr::kable(x = head(heart_dis,10), digits = 3, align = "c")

male	age	currentSmoker	cigsPerDay	prevalentHyp	totChol	sysBP	diaBP	BMI	heartRate	glucose	TenYearCHD
1	39	0	0	0	195	106.0	70	26.97	80	77	0
0	46	0	0	0	250	121.0	81	28.73	95	76	0
1	48	1	20	0	245	127.5	80	25.34	75	70	0
0	61	1	30	1	225	150.0	95	28.58	65	103	1
0	46	1	23	0	285	130.0	84	23.10	85	85	0
0	43	0	0	1	228	180.0	110	30.30	77	99	0
0	63	0	0	0	205	138.0	71	33.11	60	85	1
0	45	1	20	0	313	100.0	71	21.68	79	78	0
1	52	0	0	1	260	141.5	89	26.36	76	79	0
1	43	1	30	1	225	162.0	107	23.61	93	88	0

Frecuencias de las variables categóricas

resum <-summary(heart_dis[,c(1,3,5,6,7,8,15)])
knitr::kable(x = resum, digits = 3, align = "c")

	male	currentSmoker	BPMeds	prevalentStroke	prevalentHyp	diabetes	TenYearCHD
	0:2080	0:1918	0:3635	0:3728	0:2580	0:3647	0:3177
	1:1669	1:1831	1: 114	1: 21	1:1169	1: 102	1: 572

Resumen de las variables continuas

resum <-summary(heart_dis[,c(2,4,9,10,11,12,13,14)])
knitr::kable(x = resum, digits = 3, align = "c")

age	cigsPerDay	totChol	sysBP	diaBP	BMI	heartRate	glucose
Min. :32.00	Min. : 0.000	Min. :113	Min. : 83.5	Min. : 48.00	Min. :15.54	Min. : 44.0	Min. : 40.00
1st Qu.:42.00	1st Qu.: 0.000	1st Qu.:206	1st Qu.:117.0	1st Qu.: 75.00	1st Qu.:23.09	1st Qu.: 68.0	1st Qu.: 71.00
Median :49.00	Median : 0.000	Median :234	Median :128.0	Median : 82.00	Median :25.41	Median : 75.0	Median : 78.00
Mean :49.58	Mean : 9.005	Mean :237	Mean :132.4	Mean : 82.93	Mean :25.81	Mean : 75.7	Mean : 81.88
3rd Qu.:56.00	3rd Qu.:20.000	3rd Qu.:264	3rd Qu.:144.0	3rd Qu.: 90.00	3rd Qu.:28.06	3rd Qu.: 82.0	3rd Qu.: 87.00
Max. :70.00	Max. :70.000	Max. :696	Max. :295.0	Max. :142.50	Max. :56.80	Max. :143.0	Max. :394.00

Modelo de asociación

Para el desarrollo de este modelo debemos seleccionar las variables que incluiremos en el estudio.

Elección de variables

Análisis de las variables originales

De manera inicial incluimos todas las variables explicativas en nuestro modelo, y vamos a ir revisando punto a punto para ir refinando el modelo de asociación que queremos usar. Excluimos la variable educación porque no se encuentra definida en el manual del conjunto de datos.

Pruebas de correlación

Eliminamos la variable Presión Sistólica, porque en la prueba de correlación de Pearson encontramos que tiene colinealidad con la variable presión diastólica, tomamos la decisión de eliminar la variable de presión sistólica debido a que su efecto es el más afectado por confusión (interacción aditiva) con otras dos variables: Diabetes y Colesterol.

En la siguiente tabla vemos resumidos los valores p de la interacción de ambas variables respecto a la variable respuesta

cor.test(e2$sysBP,e2$diaBP)

## 
##  Pearson's product-moment correlation
## 
## data:  e2$sysBP and e2$diaBP
## t = 77.801, df = 3747, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7733535 0.7978485
## sample estimates:
##       cor 
## 0.7859091

knitr::kable(x = Interacciones, digits = 3, align = "c")

Inter	age	current	cigs	BPM	prevalentStroke	prevalentHyp	diabetes	totchol	sysBP	diaBP	BMI	heartRate	glucose
male	0.13719999999999999	2.3400000000000001E-2	0.34139999999999998	0.63400000000000001	0.9415	0.121	0.79010000000000002	0.25750000000000001	0.53100000000000003	0.42499999999999999	4.3200000000000002E-2	0.182	0.95960000000000001
age	NA	0.18440000000000001	0.82499999999999996	0.30399999999999999	0.29799999999999999	0.27760000000000001	0.1007	2.7100000000000002E-3	6.1199999999999997E-2	9.2600000000000002E-2	0.14399999999999999	0.7712	0.71199999999999997
current	NA	NA	NA	0.2868	0.31440000000000001	0.15311	0.51190000000000002	0.33200000000000002	0.98699999999999999	0.93200000000000005	0.4839	0.32200000000000001	0.65800000000000003
cigs	NA	NA	NA	0.18099999999999999	0.69520000000000004	8.5199999999999998E-2	0.68420000000000003	0.54600000000000004	0.75600000000000001	0.28599999999999998	0.56200000000000006	0.122	0.747
BPM	NA	NA	NA	NA	3.8059999999999997E-2	NA	0.50800000000000001	0.92900000000000005	0.626	0.28100000000000003	1.393E-2	6.6600000000000006E-2	0.26900000000000002
prevalentStroke	NA	NA	NA	NA	NA	0.68	0.97040000000000004	0.70499999999999996	0.77300000000000002	0.89900000000000002	0.29580000000000001	7.8899999999999998E-2	0.98899999999999999
prevalentHyp	NA	NA	NA	NA	NA	NA	0.4294	2.1399999999999999E-2	0.68300000000000005	0.99890000000000001	5.3400000000000003E-2	0.71099999999999997	0.47410000000000002
diabetes	NA	NA	NA	NA	NA	NA	NA	0.96699999999999997	2.0299999999999999E-2	0.15340000000000001	0.34089999999999998	0.72099999999999997	0.91269999999999996
totchol	NA	NA	NA	NA	NA	NA	NA	NA	4.4699999999999997E-2	0.27060000000000001	0.97319999999999995	0.18179999999999999	0.24809999999999999
sysBP	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.85919999999999996	8.4099999999999994E-2	0.1137	0.31340000000000001
diaBP	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.35039999999999999	0.5242	0.5363
BMI	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.51400000000000001	0.81499999999999995
heartRate	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.48089999999999999
glucose	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA

En el siguiente apartado el código referente a las interacciones.

### Interacciones de sexo con otras variables
summary(glm(TenYearCHD~male*age,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*currentSmoker,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*cigsPerDay,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*BPMeds,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*prevalentStroke,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*prevalentHyp,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*diabetes,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*totChol,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~male*glucose,family = binomial,data=heart_dis))

### Interacciones de la edad con otras variables
summary(glm(TenYearCHD~age*currentSmoker,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*cigsPerDay,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*BPMeds,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*prevalentStroke,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*prevalentHyp,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*diabetes,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*totChol,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~age*glucose,family = binomial,data=heart_dis))

### Interacción de Fumador Actual con otras variables
summary(glm(TenYearCHD~currentSmoker*cigsPerDay,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*BPMeds,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*prevalentStroke,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*prevalentHyp,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*diabetes,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*totChol,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~currentSmoker*glucose,family = binomial,data=heart_dis))

### Interacción de Consumo de Cigarrillos con otras variables
summary(glm(TenYearCHD~cigsPerDay*BPMeds,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*prevalentStroke,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*prevalentHyp,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*diabetes,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*totChol,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~cigsPerDay*glucose,family = binomial,data=heart_dis))

### Interacción de Consumo de Medicamentos Presión con otras variables
summary(glm(TenYearCHD~BPMeds*prevalentStroke,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BPMeds*prevalentHyp,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BPMeds*diabetes,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BPMeds*totChol,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BPMeds*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BPMeds*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BPMeds*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BPMeds*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BPMeds*glucose,family = binomial,data=heart_dis))

### Interacción de Ataque Cardiaco Previo con otras variables
summary(glm(TenYearCHD~prevalentStroke*prevalentHyp,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentStroke*diabetes,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentStroke*totChol,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentStroke*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentStroke*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentStroke*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentStroke*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentStroke*glucose,family = binomial,data=heart_dis))

### Interacción de Hipertensión previa con otras variables
summary(glm(TenYearCHD~prevalentHyp*diabetes,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentHyp*totChol,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentHyp*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentHyp*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentHyp*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentHyp*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~prevalentHyp*glucose,family = binomial,data=heart_dis))

### Interacción de presencia de diabetes con otras variables
summary(glm(TenYearCHD~diabetes*totChol,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~diabetes*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~diabetes*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~diabetes*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~diabetes*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~diabetes*glucose,family = binomial,data=heart_dis))

### Interacción de Colesterol Total con otras variables
summary(glm(TenYearCHD~totChol*sysBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~totChol*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~totChol*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~totChol*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~totChol*glucose,family = binomial,data=heart_dis))

### Interacción de Presión Sistólica con otras variables
summary(glm(TenYearCHD~sysBP*diaBP,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~sysBP*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~sysBP*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~sysBP*glucose,family = binomial,data=heart_dis))

### Interacción de Presión Diastólica con otras variables
summary(glm(TenYearCHD~diaBP*BMI,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~diaBP*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~diaBP*glucose,family = binomial,data=heart_dis))

### Interacción de Índice de Masa Corporal con otras variables
summary(glm(TenYearCHD~BMI*heartRate,family = binomial,data=heart_dis))
summary(glm(TenYearCHD~BMI*glucose,family = binomial,data=heart_dis))

### Interacción de Ritmo Cardiaco con otras variables
summary(glm(TenYearCHD~heartRate*glucose,family = binomial,data=heart_dis))

Agrupación de las variables en categorias

Para realizar nuestro análisis de asociación usaremos las variables agrupadas para poder identificar si alguna de las categorías están asociadas con la presencia de Cardiopatía.

Las variables que agruparemos serán:

Edad: tomando dos grupos, mayores y menores de 50 (respecto a estudios)
Presión Diastólica, creando tres grupos: Hipotensión(<60), normal (60<diaPD<90) e Hipertensión (>90).
Índice de Masa Corporal (Obeso, Sobrepeso, Normal, Bajo).
Ritmo Cardiaco: Bradicardia (<60), Normal (60<HR<90) y Taquicardia(>90)
Colesterol: Óptimo (<200), Aceptable (200<totChol<240) y Riesgo (>240)
Edad (Mayores y menores de 50 años)
Número de cigarrillos (Poco, Moderado, Elevado y Riesgoso)

e2 <- heart_dis %>%
  mutate(Pres_dias=case_when(diaBP > 90 ~ "Alta",
                             diaBP < 60 ~ "Baja",
                             TRUE ~ "Normal"), 
         RC=case_when(heartRate>100 ~ "Taqui",
                      heartRate<60 ~ "Bradi",
                      TRUE~ "Normal"),
         Gluc=case_when(glucose>126 ~ "Hipergluc",
                        TRUE~"Normal"),
         Coles=case_when(totChol<200 ~ "Optimo",
                       totChol>240 ~ "Riesgo",
                       TRUE~ "Aceptable"),
         IMC=case_when((BMI>=30)~"Obesidad",
                       (BMI<30 & BMI>=25)~"Sobrepeso",
                       (BMI<25 & BMI>=18.5)~"Normal",
                       TRUE~"Bajo"),
         Cigcons=case_when(cigsPerDay<=1~"Poco",
                            (cigsPerDay>1 & cigsPerDay<6)~"Moderado",
                            (cigsPerDay>5 & cigsPerDay<21)~"Elevado",
                            TRUE~"Riesgoso"),
         Edad=case_when(age<50 ~ "E1",
                        TRUE ~"E2")) %>%
  as.data.frame()

En el análisis mantenemos las variables: Sexo, Consumo Medicamentos para Presión Arterial (BPMeds), Ataques Cardiacos previos e Hipertensión previa estas ya se encontraban categorizadas desde un inicio del estudio.

Pruebas Chi-Cuadrado

Tomamos las variables y analizamos con un p-valor de 0,05 para determinar si incluimos o no las variables en el análisis.

chisq.test(table(e2$male,e2$TenYearCHD))
chisq.test(table(e2$currentSmoker,e2$TenYearCHD))
chisq.test(table(e2$BPMeds,e2$TenYearCHD))
chisq.test(table(e2$prevalentStroke,e2$TenYearCHD))
chisq.test(table(e2$prevalentHyp,e2$TenYearCHD))
chisq.test(table(e2$diabetes,e2$TenYearCHD))
chisq.test(table(e2$Pres_dias,e2$TenYearCHD))
chisq.test(table(e2$RC,e2$TenYearCHD))
chisq.test(table(e2$Gluc,e2$TenYearCHD))
chisq.test(table(e2$Coles,e2$TenYearCHD))
chisq.test(table(e2$IMC,e2$TenYearCHD))
chisq.test(table(e2$Cigcons,e2$TenYearCHD))
chisq.test(table(e2$Edad,e2$TenYearCHD))

Las variables de fumador actual y Ritmo Cardiaco las eliminamos porque no fueron significativa en la prueba Chi Cuadrado respecto a nuestra variable respuesta.

Análisis estratificado

En el siguiente código encontramos ajustadas las variables dos a dos para analizar si encontramos alguna interacción aditiva entre dos variables.

### Interacción de Sexo con otras variables
summary(glm(TenYearCHD~male*BPMeds,family=binomial,data=e2))
summary(glm(TenYearCHD~male*prevalentStroke,family=binomial,data=e2))
summary(glm(TenYearCHD~male*prevalentHyp,family=binomial,data=e2))
summary(glm(TenYearCHD~male*diabetes,family=binomial,data=e2))
summary(glm(TenYearCHD~male*Pres_dias,family=binomial,data=e2))
summary(glm(TenYearCHD~male*Gluc,family=binomial,data=e2))
summary(glm(TenYearCHD~male*Coles,family=binomial,data=e2))
summary(glm(TenYearCHD~male*IMC,family=binomial,data=e2))
summary(glm(TenYearCHD~male*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~male*Edad,family=binomial,data=e2))

### Interacción de Consumo de Medicamentos Presión con otras variables
summary(glm(TenYearCHD~BPMeds*prevalentStroke,family=binomial,data=e2))
summary(glm(TenYearCHD~BPMeds*prevalentHyp,family=binomial,data=e2))
summary(glm(TenYearCHD~BPMeds*diabetes,family=binomial,data=e2))
summary(glm(TenYearCHD~BPMeds*Pres_dias,family=binomial,data=e2))
summary(glm(TenYearCHD~BPMeds*Gluc,family=binomial,data=e2))
summary(glm(TenYearCHD~BPMeds*Coles,family=binomial,data=e2))
summary(glm(TenYearCHD~BPMeds*IMC,family=binomial,data=e2))
summary(glm(TenYearCHD~BPMeds*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~BPMeds*Edad,family=binomial,data=e2))

### Interacción de Ataque Cardiaco Previo con otras variables
summary(glm(TenYearCHD~prevalentStroke*prevalentHyp,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentStroke*diabetes,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentStroke*Pres_dias,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentStroke*Gluc,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentStroke*Coles,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentStroke*IMC,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentStroke*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentStroke*Edad,family=binomial,data=e2))

### Interacción de presencia de Hipertensión con otras variables
summary(glm(TenYearCHD~prevalentHyp*diabetes,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentHyp*Pres_dias,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentHyp*Gluc,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentHyp*Coles,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentHyp*IMC,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentHyp*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~prevalentHyp*Edad,family=binomial,data=e2))

### Interacción de presencia de diabetes con otras variables
summary(glm(TenYearCHD~diabetes*Pres_dias,family=binomial,data=e2))
summary(glm(TenYearCHD~diabetes*Gluc,family=binomial,data=e2))
summary(glm(TenYearCHD~diabetes*Coles,family=binomial,data=e2))
summary(glm(TenYearCHD~diabetes*IMC,family=binomial,data=e2))
summary(glm(TenYearCHD~diabetes*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~diabetes*Edad,family=binomial,data=e2))

### Interacción de la presión diastólica con otras variables
summary(glm(TenYearCHD~Pres_dias*Gluc,family=binomial,data=e2))
summary(glm(TenYearCHD~Pres_dias*Coles,family=binomial,data=e2))
summary(glm(TenYearCHD~Pres_dias*IMC,family=binomial,data=e2))
summary(glm(TenYearCHD~Pres_dias*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~Pres_dias*Edad,family=binomial,data=e2))

### ### Interacción de Nivel Glucosa con otras variables
summary(glm(TenYearCHD~Gluc*Coles,family=binomial,data=e2))
summary(glm(TenYearCHD~Gluc*IMC,family=binomial,data=e2))
summary(glm(TenYearCHD~Gluc*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~Gluc*Edad,family=binomial,data=e2))

### Interacción de Nivel de colesterol con otras variables
summary(glm(TenYearCHD~Coles*IMC,family=binomial,data=e2))
summary(glm(TenYearCHD~Coles*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~Coles*Edad,family=binomial,data=e2))

### Interacción de Índice de Masa Corporal con otras variables
summary(glm(TenYearCHD~IMC*Cigcons,family=binomial,data=e2))
summary(glm(TenYearCHD~IMC*Edad,family=binomial,data=e2))

### Interacción de Nivel Consumo de Cigarrillos con otras variables
summary(glm(TenYearCHD~Cigcons*Edad,family=binomial,data=e2))

Vease la siguiente tabla que resume estos valores

knitr::kable(x = int2, digits = 3, align = "c")

…1	BPMeds	prev_stroke	prev_hyp	diabetes	Pres_dias	Gluc	Coles	IMC	Cigcons	Edad
male	0.63400000000000001	0.94099999999999995	1.2E-2	0.79	0.35	0.96	0.84599999999999997	0.73199999999999998	4.4999999999999998E-2	0.22
BPMeds	NA	3.7999999999999999E-2	NA	0.50800000000000001	0.55200000000000005	0.34899999999999998	0.55800000000000005	0.96899999999999997	0.48399999999999999	0.58899999999999997
prev_stroke	NA	NA	0.68	0.97	0.92700000000000005	0.52200000000000002	0.91100000000000003	0.79200000000000004	0.52900000000000003	0.96799999999999997
prev_hyp	NA	NA	NA	0.42899999999999999	0.41099999999999998	0.88100000000000001	6.8000000000000005E-2	7.0999999999999994E-2	0.25800000000000001	0.249
diabetes	NA	NA	NA	NA	0.76200000000000001	9.6000000000000002E-2	0.96199999999999997	0.54200000000000004	0.111	0.34699999999999998
Pres_dias	NA	NA	NA	NA	NA	0.93	5.6000000000000001E-2	0.50800000000000001	0.27	0.47599999999999998
Gluc	NA	NA	NA	NA	NA	NA	0.50800000000000001	0.68799999999999994	0.38600000000000001	0.51900000000000002
Coles	NA	NA	NA	NA	NA	NA	NA	7.8E-2	0.39100000000000001	7.0000000000000001E-3
IMC	NA	NA	NA	NA	NA	NA	NA	NA	0.64100000000000001	0.33500000000000002
Cigcos	NA	NA	NA	NA	NA	NA	NA	NA	NA	0.316
Edad	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA

Teniendo en cuenta estos resultados podemos ver que la variable sexo(male) puede ser un posible efecto de confusión sobre las variables de ataque cardiaco previo y consumo de cigarrillo. Asimismo pasa con la variable edad que puede generar confusión sobre la variable del nivel de colesterol. Por lo tanto, tomaremos el sexo y la edad como posibles variables de confusión del modelo.

Análisis de los resultados

Preparamos el modelo para poder calcular los OR(RRI) asociados a cada uno de los niveles de los factores y evaluar cuáles pueden estar asociados al riesgo de tener un posible diagnóstico en los próximos 10 de años de una enfermedad coronaria.

fit1<-glm(TenYearCHD~male+prevalentStroke+prevalentHyp+BPMeds+diabetes+Pres_dias+Gluc+Coles+IMC+Cigcons+Edad, family = binomial ,data=e2)
knitr::kable(x = cbind(exp(coef(fit1)),exp(confint(fit1,level=0.95))), digits = 3, align = "c")

## Waiting for profiling to be done...

		2.5 %	97.5 %
(Intercept)	0.346	0.108	1.029
male1	1.673	1.363	2.057
prevalentStroke1	1.648	0.614	4.215
prevalentHyp1	1.713	1.341	2.185
BPMeds1	1.430	0.913	2.206
diabetes1	1.019	0.517	1.923
Pres_diasBaja	2.325	0.901	5.465
Pres_diasNormal	0.753	0.586	0.969
GlucNormal	0.298	0.147	0.599
ColesOptimo	0.835	0.618	1.119
ColesRiesgo	1.164	0.946	1.434
IMCNormal	0.693	0.312	1.771
IMCObesidad	0.858	0.373	2.242
IMCSobrepeso	0.765	0.342	1.959
CigconsModerado	0.780	0.477	1.232
CigconsPoco	0.724	0.577	0.909
CigconsRiesgoso	1.357	0.996	1.841
EdadE2	2.888	2.344	3.572

Modelo de predicción

suppressMessages(library(ggplot2))
suppressMessages(library(glmtoolbox))
suppressMessages(library(dplyr))
suppressMessages(library(tidyverse)) # para manejo de datos    
suppressMessages(library(MASS))
suppressMessages(library(class))
suppressMessages(library(caret))
suppressMessages(library(knitr))
# install.packages("smotefamily")
suppressMessages(library(smotefamily))
#suppressMessages(install.packages("epiR"))

Importamos el conjunto de datos.

heart <- suppressMessages(read_delim("framingham.csv", delim = ",", escape_double = FALSE, trim_ws = TRUE))
suppressMessages(attach(heart))

Seleccionar las filas sin valores faltantes en las variables de interés.

# Seleccionar las filas sin valores faltantes en las variables de interés
heart <- na.omit(heart)

Implementamos el modelo con todas las variables.

n <- names(heart)
f <- as.formula(paste("TenYearCHD ~", paste(n[!n %in% "TenYearCHD"], collapse = " + "))) #,"-TenYearCHD"
f

## TenYearCHD ~ male + age + education + currentSmoker + cigsPerDay + 
##     BPMeds + prevalentStroke + prevalentHyp + diabetes + totChol + 
##     sysBP + diaBP + BMI + heartRate + glucose

modelo_glm <- glm(f,data = heart,family = binomial(link="logit"));summary(modelo_glm)

## 
## Call:
## glm(formula = f, family = binomial(link = "logit"), data = heart)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9581  -0.5944  -0.4267  -0.2829   2.8402  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -8.322206   0.715476 -11.632  < 2e-16 ***
## male             0.555098   0.109046   5.091 3.57e-07 ***
## age              0.063453   0.006680   9.499  < 2e-16 ***
## education       -0.047497   0.049390  -0.962  0.33621    
## currentSmoker    0.070875   0.156749   0.452  0.65115    
## cigsPerDay       0.017929   0.006238   2.874  0.00405 ** 
## BPMeds           0.162255   0.234309   0.692  0.48863    
## prevalentStroke  0.693502   0.489532   1.417  0.15658    
## prevalentHyp     0.234638   0.138037   1.700  0.08917 .  
## diabetes         0.039461   0.315483   0.125  0.90046    
## totChol          0.002324   0.001127   2.062  0.03920 *  
## sysBP            0.015398   0.003808   4.043 5.27e-05 ***
## diaBP           -0.004132   0.006438  -0.642  0.52096    
## BMI              0.006603   0.012758   0.518  0.60476    
## heartRate       -0.003250   0.004211  -0.772  0.44030    
## glucose          0.007124   0.002234   3.189  0.00143 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3120.5  on 3655  degrees of freedom
## Residual deviance: 2754.2  on 3640  degrees of freedom
## AIC: 2786.2
## 
## Number of Fisher Scoring iterations: 5

Podemos observar que las variables que resultan significativas a un nivel del 5 % son sexo, edad, cigarrillos por día, total de calorías, presión sistólica y niveles de glucosa. Con base en eso, ponemos a prueba el nivel predictivo del modelo, su precisión y, adicionalmente, verificaremos los supuestos de este.

glm.probs <- predict(modelo_glm, type = "response")
glm.pred <- rep(0, length(glm.probs))
glm.pred[glm.probs > 0.5] <- 1
confusionMatrix(data =  as.factor(as.numeric(glm.probs>0.5)), reference = as.factor(heart$TenYearCHD), positive="1")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 3080  506
##          1   19   51
##                                           
##                Accuracy : 0.8564          
##                  95% CI : (0.8446, 0.8676)
##     No Information Rate : 0.8476          
##     P-Value [Acc > NIR] : 0.07273         
##                                           
##                   Kappa : 0.1332          
##                                           
##  Mcnemar's Test P-Value : < 2e-16         
##                                           
##             Sensitivity : 0.09156         
##             Specificity : 0.99387         
##          Pos Pred Value : 0.72857         
##          Neg Pred Value : 0.85890         
##              Prevalence : 0.15235         
##          Detection Rate : 0.01395         
##    Detection Prevalence : 0.01915         
##       Balanced Accuracy : 0.54272         
##                                           
##        'Positive' Class : 1               
##

par(mfrow=c(1,2))
ROCc(modelo_glm, plot.it = T,main="Modelo propuesto")

## [1] 0 1

## 
##  Area Under ROC Curve  =  0.739 
##      Gini Coefficient  =  0.478 
##         K-S Statistic  =  0.378

set.seed(1234)
envelope(modelo_glm,main="Modelo propuesto",silent = TRUE)

## 
  |                               
  |                         |   0%
  |                               
  |+                        |   4%
  |                               
  |++                       |   8%
  |                               
  |+++                      |  12%
  |                               
  |++++                     |  16%
  |                               
  |+++++                    |  20%
  |                               
  |++++++                   |  24%
  |                               
  |+++++++                  |  28%
  |                               
  |++++++++                 |  32%
  |                               
  |+++++++++                |  36%
  |                               
  |++++++++++               |  40%
  |                               
  |+++++++++++              |  44%
  |                               
  |++++++++++++             |  48%
  |                               
  |+++++++++++++            |  52%
  |                               
  |++++++++++++++           |  56%
  |                               
  |+++++++++++++++          |  60%
  |                               
  |++++++++++++++++         |  64%
  |                               
  |+++++++++++++++++        |  68%
  |                               
  |++++++++++++++++++       |  72%
  |                               
  |+++++++++++++++++++      |  76%
  |                               
  |++++++++++++++++++++     |  80%
  |                               
  |+++++++++++++++++++++    |  84%
  |                               
  |++++++++++++++++++++++   |  88%
  |                               
  |+++++++++++++++++++++++  |  92%
  |                               
  |++++++++++++++++++++++++ |  96%
  |                               
  |+++++++++++++++++++++++++| 100%

par(mfrow=c(1,1))

car::vif(modelo_glm)

##            male             age       education   currentSmoker      cigsPerDay 
##        1.242150        1.280953        1.053098        2.578367        2.726360 
##          BPMeds prevalentStroke    prevalentHyp        diabetes         totChol 
##        1.116193        1.023355        1.985800        1.810986        1.070603 
##           sysBP           diaBP             BMI       heartRate         glucose 
##        3.616670        2.845913        1.202350        1.096558        1.814159

Podemos observar que la precisión del modelo (85.6 %) resulta bastante apropiada para los fines del mismo; sin embargo, se debe tener en cuenta que la tasa de detección de casos positivos (sensibilidad) resulta ser del 9 % y la tasa de detección de verdaderos negativos (especificidad) del 99.3 %, mostrando que el modelo puede perder la detección de muchos casos positivos.

Por otro lado, podemos ver que el área bajo la curva ROC del 73.9 %, lo cual indica que el modelo tiene una capacidad moderada para distinguir entre las dos clases: padecer una afección cardíaca y no padecerla.

No obstante, validando los supuestos del modelo, podemos evidenciar un comportamiento muy apropiado de los residuales en el q-q plot, lo que nos permite decir que se presenta normalidad. También, se analizaron las varianzas infladas descartando cualquier problema de multicolinealidad con ninguna variable (si son menores a 5 son aceptables); hay que tener en cuenta que las variables con más alto valor de varianzas infladas son “currentSmoker”, “cigsPerDay”, “sysBP” y “diaBP”, indicando una posible correlación con otras variables.

glmtoolbox::hltest(modelo_glm)

## 
##    The Hosmer-Lemeshow goodness-of-fit test
## 
##  Group Size Observed  Expected
##      1  366        7  11.00660
##      2  366       21  17.30875
##      3  366       20  23.47755
##      4  366       30  30.62616
##      5  366       40  38.52235
##      6  366       50  47.69341
##      7  366       54  59.84074
##      8  366       93  75.37047
##      9  366       98  99.24261
##     10  362      144 153.91136
## 
##          Statistic =  10.09219 
## degrees of freedom =  8 
##            p-value =  0.25861

Ahora, utilizando el test de Hosmer-Lemeshow para evaluar la bondad de ajuste de nuestro modelo de regresión logística, podemos evidenciar que, como el valor-p es mayor que el nivel de significación, no se rechaza la hipótesis nula y se concluye que el modelo ajusta bien los datos.

pasos = stepCriterion(modelo_glm, "aic", direction = "backward")

pasos$final

## [1] "~ male + age + cigsPerDay + prevalentStroke + prevalentHyp + totChol + sysBP + glucose"

step_glm <- glm(TenYearCHD ~ age + sysBP + male + cigsPerDay + glucose + totChol + prevalentHyp + prevalentStroke , data = heart,family = binomial(link="logit"));summary(step_glm)

## 
## Call:
## glm(formula = TenYearCHD ~ age + sysBP + male + cigsPerDay + 
##     glucose + totChol + prevalentHyp + prevalentStroke, family = binomial(link = "logit"), 
##     data = heart)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0061  -0.5975  -0.4294  -0.2842   2.8627  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -8.739521   0.522563 -16.724  < 2e-16 ***
## age              0.065337   0.006444  10.140  < 2e-16 ***
## sysBP            0.014219   0.002857   4.976 6.48e-07 ***
## male             0.553152   0.107037   5.168 2.37e-07 ***
## cigsPerDay       0.019574   0.004182   4.681 2.85e-06 ***
## glucose          0.007314   0.001673   4.373 1.23e-05 ***
## totChol          0.002248   0.001122   2.003   0.0452 *  
## prevalentHyp     0.226231   0.135098   1.675   0.0940 .  
## prevalentStroke  0.751412   0.483562   1.554   0.1202    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3120.5  on 3655  degrees of freedom
## Residual deviance: 2757.2  on 3647  degrees of freedom
## AIC: 2775.2
## 
## Number of Fisher Scoring iterations: 5

Ahora el modelo sin las variables “prevalentHyp” y “prevalentStroke”; es decir, el modelo con todas las variables significativas.

nuevo_glm <- glm(TenYearCHD ~ age + sysBP + male + cigsPerDay + glucose + totChol, data = heart,family = binomial(link="logit"));summary(nuevo_glm)

## 
## Call:
## glm(formula = TenYearCHD ~ age + sysBP + male + cigsPerDay + 
##     glucose + totChol, family = binomial(link = "logit"), data = heart)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0074  -0.5966  -0.4315  -0.2864   2.8811  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -9.129843   0.475530 -19.199  < 2e-16 ***
## age          0.065896   0.006426  10.254  < 2e-16 ***
## sysBP        0.017534   0.002149   8.159 3.38e-16 ***
## male         0.561446   0.106845   5.255 1.48e-07 ***
## cigsPerDay   0.019226   0.004176   4.604 4.14e-06 ***
## glucose      0.007280   0.001677   4.342 1.41e-05 ***
## totChol      0.002272   0.001123   2.024    0.043 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3120.5  on 3655  degrees of freedom
## Residual deviance: 2762.5  on 3649  degrees of freedom
## AIC: 2776.5
## 
## Number of Fisher Scoring iterations: 5

glm.probs <- predict(nuevo_glm, type = "response")
glm.pred <- rep(0, length(glm.probs))
glm.pred[glm.probs > 0.5] <- 1
cm = confusionMatrix(data =  as.factor(as.numeric(glm.probs>0.5)), reference = as.factor(heart$TenYearCHD), positive="1")
cm

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 3076  515
##          1   23   42
##                                           
##                Accuracy : 0.8528          
##                  95% CI : (0.8409, 0.8642)
##     No Information Rate : 0.8476          
##     P-Value [Acc > NIR] : 0.1977          
##                                           
##                   Kappa : 0.1066          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.07540         
##             Specificity : 0.99258         
##          Pos Pred Value : 0.64615         
##          Neg Pred Value : 0.85659         
##              Prevalence : 0.15235         
##          Detection Rate : 0.01149         
##    Detection Prevalence : 0.01778         
##       Balanced Accuracy : 0.53399         
##                                           
##        'Positive' Class : 1               
##

par(mfrow=c(1,2))
ROCc(nuevo_glm, plot.it = T,main="Nuevo modelo propuesto")

## [1] 0 1

## 
##  Area Under ROC Curve  =  0.736 
##      Gini Coefficient  =  0.472 
##         K-S Statistic  =  0.369

set.seed(1234)
envelope(nuevo_glm,main="Nuevo modelo propuesto")

## 
  |                               
  |                         |   0%
  |                               
  |+                        |   4%
  |                               
  |++                       |   8%
  |                               
  |+++                      |  12%
  |                               
  |++++                     |  16%
  |                               
  |+++++                    |  20%
  |                               
  |++++++                   |  24%
  |                               
  |+++++++                  |  28%
  |                               
  |++++++++                 |  32%
  |                               
  |+++++++++                |  36%
  |                               
  |++++++++++               |  40%
  |                               
  |+++++++++++              |  44%
  |                               
  |++++++++++++             |  48%
  |                               
  |+++++++++++++            |  52%
  |                               
  |++++++++++++++           |  56%
  |                               
  |+++++++++++++++          |  60%
  |                               
  |++++++++++++++++         |  64%
  |                               
  |+++++++++++++++++        |  68%
  |                               
  |++++++++++++++++++       |  72%
  |                               
  |+++++++++++++++++++      |  76%
  |                               
  |++++++++++++++++++++     |  80%
  |                               
  |+++++++++++++++++++++    |  84%
  |                               
  |++++++++++++++++++++++   |  88%
  |                               
  |+++++++++++++++++++++++  |  92%
  |                               
  |++++++++++++++++++++++++ |  96%
  |                               
  |+++++++++++++++++++++++++| 100%

par(mfrow=c(1,1))

car::vif(nuevo_glm)

##        age      sysBP       male cigsPerDay    glucose    totChol 
##   1.186407   1.148048   1.196158   1.228479   1.011907   1.060822

Analizando el nuevo modelo, obtenemos una precisión del 85.2 %, siendo levemente menor al del modelo general; también la sensibilidad resultó ser del 7.5 % y la especificidad del 99.2 %, mostrando que este modelo también puede perder la detección de muchos casos positivos.

No hubo un cambio significativo en el área bajo la curva de este modelo.

Ahora, validando supuestos del modelo, podemos evidenciar un comportamiento muy apropiado de los residuales en el q-q plot, lo que nos permite decir que se presenta normalidad. No hubo problemas de multicolinealidad, dado que se excluyeron variables que potencialmente podrían estar correlacionadas con otras.

glmtoolbox::hltest(nuevo_glm)

## 
##    The Hosmer-Lemeshow goodness-of-fit test
## 
##  Group Size Observed  Expected
##      1  366       10  11.07294
##      2  366       21  17.54932
##      3  366       21  23.97677
##      4  366       27  31.11584
##      5  366       45  39.23359
##      6  366       41  48.16620
##      7  366       64  59.91823
##      8  366       79  75.13304
##      9  366      109  98.78773
##     10  362      140 152.04634
## 
##          Statistic =  7.66192 
## degrees of freedom =  8 
##            p-value =  0.46717

Nuevamente estamos en presencia de un modelo que se ajusta bien a los datos. No obstante, es importante considerar que el comportamiento que evidenciamos anteriormente sea producto de datos influyentes. Así que procedemos a analizarlos.

cook_for = cooks.distance(nuevo_glm)
# par(mfrow=c(1,2))
umbral_cook <- 3*mean(cook_for)
# Graficamos la distancia de Cook
plot(cook_for, type = "h", pch = 20, main = "Distancia de Cook Modelo propuesto",
     ylab = "Distancia de Cook", xlab = "Índice de observación")
# Agregamos una línea horizontal para el umbral de corte
abline(h = umbral_cook, col = "red", lty = 2)
# Agregamos puntos rojos para las observaciones influyentes
influyentes_for <- as.numeric(which(cook_for > umbral_cook))
points(influyentes_for, cook_for[influyentes_for], pch = 20, col = "red")

modelo_glm_1 <- update(nuevo_glm,
                      subset=-influyentes_for);summary(modelo_glm_1)

## 
## Call:
## glm(formula = TenYearCHD ~ age + sysBP + male + cigsPerDay + 
##     glucose + totChol, family = binomial(link = "logit"), data = heart, 
##     subset = -influyentes_for)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3845  -0.2604  -0.1355  -0.0631   3.1834  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -19.668101   1.301061 -15.117  < 2e-16 ***
## age           0.143819   0.014610   9.844  < 2e-16 ***
## sysBP         0.038177   0.004157   9.184  < 2e-16 ***
## male          1.467940   0.228197   6.433 1.25e-10 ***
## cigsPerDay    0.036623   0.008154   4.491 7.08e-06 ***
## glucose       0.014965   0.005040   2.969  0.00298 ** 
## totChol       0.004266   0.002433   1.753  0.07955 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1254.78  on 3219  degrees of freedom
## Residual deviance:  861.14  on 3213  degrees of freedom
## AIC: 875.14
## 
## Number of Fisher Scoring iterations: 7

Se evidencian 436 datos influyentes, los cuales procedemos a excluir del modelo. Podemos observar que el cambio más importante que se presenta es con respecto a la variable “totChol” pierde significación a un nivel del 5 %. Procedemos a evaluar la predicción del modelo sin los datos influyentes, además de sus supuestos.

glm.probs <- predict(modelo_glm_1, type = "response")
glm.pred <- rep(0, length(glm.probs))
glm.pred[glm.probs > 0.5] <- 1

confusionMatrix(data =  as.factor(as.numeric(glm.probs>0.5)), reference = as.factor(heart$TenYearCHD[-influyentes_for]), positive="1")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 3059  130
##          1    4   27
##                                          
##                Accuracy : 0.9584         
##                  95% CI : (0.9509, 0.965)
##     No Information Rate : 0.9512         
##     P-Value [Acc > NIR] : 0.03054        
##                                          
##                   Kappa : 0.2756         
##                                          
##  Mcnemar's Test P-Value : < 2e-16        
##                                          
##             Sensitivity : 0.171975       
##             Specificity : 0.998694       
##          Pos Pred Value : 0.870968       
##          Neg Pred Value : 0.959235       
##              Prevalence : 0.048758       
##          Detection Rate : 0.008385       
##    Detection Prevalence : 0.009627       
##       Balanced Accuracy : 0.585334       
##                                          
##        'Positive' Class : 1              
##

par(mfrow=c(1,2))
ROCc(modelo_glm_1, plot.it = T,main="Nuevo modelo sin influyentes")

## [1] 0 1

## 
##  Area Under ROC Curve  =  0.891 
##      Gini Coefficient  =  0.781 
##         K-S Statistic  =  0.644

set.seed(1234)
envelope(modelo_glm_1,main="Nuevo modelo sin influyentes",silent = TRUE)

## 
  |                               
  |                         |   0%
  |                               
  |+                        |   4%
  |                               
  |++                       |   8%
  |                               
  |+++                      |  12%
  |                               
  |++++                     |  16%
  |                               
  |+++++                    |  20%
  |                               
  |++++++                   |  24%
  |                               
  |+++++++                  |  28%
  |                               
  |++++++++                 |  32%
  |                               
  |+++++++++                |  36%
  |                               
  |++++++++++               |  40%
  |                               
  |+++++++++++              |  44%
  |                               
  |++++++++++++             |  48%
  |                               
  |+++++++++++++            |  52%
  |                               
  |++++++++++++++           |  56%
  |                               
  |+++++++++++++++          |  60%
  |                               
  |++++++++++++++++         |  64%
  |                               
  |+++++++++++++++++        |  68%
  |                               
  |++++++++++++++++++       |  72%
  |                               
  |+++++++++++++++++++      |  76%
  |                               
  |++++++++++++++++++++     |  80%
  |                               
  |+++++++++++++++++++++    |  84%
  |                               
  |++++++++++++++++++++++   |  88%
  |                               
  |+++++++++++++++++++++++  |  92%
  |                               
  |++++++++++++++++++++++++ |  96%
  |                               
  |+++++++++++++++++++++++++| 100%

par(mfrow=c(1,1))

car::vif(modelo_glm_1)

##        age      sysBP       male cigsPerDay    glucose    totChol 
##   1.183820   1.169391   1.386489   1.331283   1.019041   1.114904

Vemos un cambio interesante en la precisión del modelo, ya que en este caso, resulta ser del 95.8 %. La especificidad no varió significativamente, como lo hizo la sensibilidad, que terminó siendo del 17.19 %; no obstante, aún sigue siendo baja.

Otro cambio significativo ocurrió en el área bajo la curva ROC que llegó a ser del 89.1 %, mostrando una discriminación entre clases bastante apropiada.

Se siguen cumpliendo el supuesto de normalidad y no hay presencia de problemas de multicolinealidad.

glmtoolbox::hltest(modelo_glm_1)

## 
##    The Hosmer-Lemeshow goodness-of-fit test
## 
##  Group Size Observed   Expected
##      1  322        0  0.1654760
##      2  322        0  0.4461232
##      3  322        0  0.9100877
##      4  322        1  1.6526806
##      5  322        2  2.9096720
##      6  322        5  4.7474009
##      7  322        7  8.0210226
##      8  322       22 14.4887643
##      9  322       40 29.8967140
##     10  322       80 93.7620592
## 
##          Statistic =  12.90892 
## degrees of freedom =  8 
##            p-value =  0.11502

Este modelo sigue siendo apropiado para los datos que tenemos.

Ahora, aunque perdemos significación en una variable, optamos por elegir el modelo sin datos influyentes, ya que presenta precisión, sensibilidad y área bajo la curva mayor, en comparación a los otros. Ahora, procedemos a analizar sus odd ratio

coef <- exp(coef(modelo_glm_1))
int <- exp(confint(modelo_glm_1,level = 0.95))

## Waiting for profiling to be done...

oddratio <- cbind(coef,int)
suppressMessages(knitr::kable(oddratio, digits = 3))

	coef	2.5 %	97.5 %
(Intercept)	0.000	0.000	0.000
age	1.155	1.123	1.189
sysBP	1.039	1.031	1.048
male	4.340	2.794	6.844
cigsPerDay	1.037	1.021	1.054
glucose	1.015	1.005	1.025
totChol	1.004	0.999	1.009

Resultados

Modelo de Asociación

Posterior a la limpieza del conjunto de datos encontramos que el riesgo de ser diagnosticado con una enfermedad coronaria puede estar asociado con el sexo hombre (OR:1.61), con el diagnóstico de hipertensión (OR:1.71), con tener presión diastólica normal (OR: 0.748), con nivel de glucosa normal (OR:0.309), con el bajo consumo de cigarrillo (OR:0.722), con el consumo riesgo de cigarrillo (OR:1.377) y con tener menos de 50 años (OR:0.353).

Modelo de Predicción

Tenemos las siguientes conclusiones: 1. Por cada año adicional, el riesgo de desarrollar una enfermedad coronaria en 10 años aumenta 0.155 veces. 2. Por cada unidad que aumente la presión sistólica, el riesgo de desarrollar una enfermedad coronaria en 10 años aumenta 0.039 veces. 3. Ser hombre aumenta 3.34 veces el riesgo de desarrollar una enfermedad coronaria en 10 años. 4. Por cada cigarrillo diario adicional que se consuma por encima del promedio, el riesgo de desarrollar una enfermedad coronaria en 10 años aumenta 0.037 veces. 5. Por cada unidad que aumente en nivel de glucosa, el riesgo de desarrollar una enfermedad coronaria en 10 años aumenta 0.015 veces. 6. Por cada unidad que aumente el colesterol, el riesgo de desarrollar una enfermedad coronaria en 10 años aumenta 0.004 veces.

Referencias

https://www.mayoclinic.org/es-es/diseases-conditions/heart-attack/symptoms-causes/syc-20373106#:~:text=Entre%20los%20factores%20de%20riesgo,y%20las%20mujeres%20m%C3%A1s%20j%C3%B3venes
https://www.cdc.gov/healthyweight/spanish/assessing/bmi/adult_bmi/index.html#:~:text=%C2%BFC%C3%B3mo%20se%20calcula%20el%20IMC%3F,-El%20IMC%20se&text=Con%20el%20sistema%20m%C3%A9trico%2C%20la,la%20estatura%20en%20metros%20cuadrados.
https://www.paho.org/es/temas/enfermedades-cardiovasculares
https://www.drmartinezcardiologia.es/cuantos-cigarrillos-riesgo-infarto/
https://www.medicalnewstoday.com/articles/es/prueba-de-glucosa-en-sangre#diabetes
https://www.kaggle.com/datasets/dileep070/heart-disease-prediction-using-logistic-regression?resource=download
https://www.mayoclinic.org/es-es/diseases-conditions/high-blood-pressure/diagnosis-treatment/drc-20373417#:~:text=Para%20medir%20la%20presi%C3%B3n%20arterial,automatizada%20de%20la%20presi%C3%B3n%20arterial.
Factores asociados al intento de suicidio en la población colombiana. C. Gómez-Restrepo, N. Rodriguez Malagón, A. Bohórquez. Revista colombiana de psiquiatría 31 (4), 283-298