Taller Modelos Pronósticos

Aclaración: Este tarller contiene análisis discutidos en grupo, sin embargo las conclusiones son individuales

Introducción

El objetivo general es desarrollar un modelo de regresión para la mortalidad a los 30 días después de un infarto agudo de miocardio (IAM), incorporando la experiencia de los datos en un pequeño conjunto de datos: La base de datos sample4, contiene 785 pacientes, de los cuales 52 fallecieron.

A continuación se realiza una descripción general de los conjuntos de datos con una descripción de las variables incluidas.

Descripción general de los predictores candidatos:

Se consideran un total de 20 predictores. 18 fueron descritos en los artículos de Lee et al y Califf et al. Además, se incluyen 2 variables que se consideraron en otros estudios: tiempo hasta el alivio del dolor torácico > 1 hora y número de derivaciones con elevación del ST en el ECG. Se crea una variable adicional llamada indice de masa corporal.

La base de datos no contiene valores perdidos.

# Leyendo la base de datos
library(haven)
mydata <- read_dta("datasets/sample4.dta")

La Tabla 1 muestra las variables de la base de datos, su descripción y codificación.

Variable	Nombre	Descripción	Códigos/Valores
1	DAY30	Mortalidad a los 30 días después del IAM	0=No, 1=Sí
2	Sex	Género femenino	0=No, 1=Sí
3	AGE	Edad	Número 19:110
4	A65	Edad > 65 años	0=No, 1=Sí
5	KILLIP	Clase Killip de función ventricular	Cat 1:4
6	SHO	Shock: Killip 3 o 4	0=No, 1=Sí
7	DIA	Diabetes	0=No, 1=Sí
8	HYP	Hipotensión (PAS<100 mmHg)	0=No, 1=Sí
9	HRT	Taquicardia (Pulso> 80 pul/min)	0=No, 1=Sí
10	ANT	Infarto de localización anterior	0=No, 1=Sí
11	PMI	Infarto de miocardio previo	0=No, 1=Sí
12	HEI	Altura en cm	Número 140:212
13	WEI	Peso en Kg	Número 36:213
14	HTN	Historia de hipertensión	0=No, 1=Sí
15	SMK	Tabaquismo	1=nunca, 2=exfumador, 3=activo
16	LIP	Hipercolesterolemia (lípidos)	0=No, 1=Sí
17	PAN	Angina de pecho previa	0=No, 1=Sí
18	FAM	Historia familiar de IAM	0=No, 1=Sí
19	STE	Elevación del segmento ST (Derivaciones)	Número 0:11
20	ST4	Elevación del segmento ST > 4 derivaciones	0=No, 1=Sí
21	TTR	Tiempo de alivio del dolor > 1h	0=No, 1=Sí
22	IMC	Indice de masa corporal	Número

Preguntas

Pregunta 1: Análisis descriptivo univariado

a) Dé algunas descripciones de las covariables.

En general se toma como codificación para los análisis 0=No (Ausencia) o 1=Sí (Presencia).

La variable DAY30, es la variable desenlace de interes e indica mortalidad a 30 dias luego de un IAM.

La clasificacion de Killip, hace referencia un índice de gravedad de la insuficiencia cardíaca en pacientes con IAM y se creó para evaluar el riesgo de muerte intrahospitalaria. Esta variable es categórica y es menos grave entre menor sea su valor. Killip 1 y 2 implica menos gravedad de Killip 3 y 4.

Hipotensión se definió como presión sistólica <100mmHg.

Taquicardia se definió como pulso > 80pulsos/min.

Se desea incluir la variable IMC que se calcula dividiendo los kilogramos del peso por el cuadrado de la estatura en metros.

Análisis exploratorio de datos

A continuación se presenta en la grafica 1 un análisis exploratorio de las variables incluidas en el estudio.

library(summarytools)

# dfSummary(mydata,valid.col = FALSE, style = "rmarkdown") # Permite hacer un análisis
# Para presentación en html

print(dfSummary(mydata,valid.col = TRUE),method = "render")

Data Frame Summary

mydata

Dimensions: 785 x 26
Duplicates: 0

No	Variable	Label	Stats / Values	Freqs (% of Valid)	Graph	Valid	Missing
1	day30 [factor]	Mortalidad a los 30 días después del IAM ('day30', n(%), 0=No, 1=Si)	1. No 2. Si	733 ( 93.4% ) 52 ( 6.6% )		785 (100.0%)	0 (0.0%)
2	sex [factor]	Género ('sex', n (%), 0=Masculino, 1=Femenino)	1. Masculino 2. Femenino	582 ( 74.1% ) 203 ( 25.9% )		785 (100.0%)	0 (0.0%)
3	age [numeric]	Edad ('age', en años)	Mean (sd) : 61.8 (11.1) min ≤ med ≤ max: 28.4 ≤ 62.8 ≤ 86 IQR (CV) : 17.2 (0.2)	710 distinct values		785 (100.0%)	0 (0.0%)
4	a65 [factor]	Edad > 65 años ('a65', n (%), 0=No, 1=Si)	1. No 2. Si	458 ( 58.3% ) 327 ( 41.7% )		785 (100.0%)	0 (0.0%)
5	killip [factor]	Clase Killip de función ventricular izq. ('killip', n (%), Clases de 1 a 4)	1. Clase 1 2. Clase 2 3. Clase 3 4. Clase 4	616 ( 78.5% ) 148 ( 18.9% ) 20 ( 2.5% ) 1 ( 0.1% )		785 (100.0%)	0 (0.0%)
6	sho [factor]	Shock: Clase Killip 3 o 4 ('sho', n (%), 0=No, 1=Si)	1. No 2. Si	764 ( 97.3% ) 21 ( 2.7% )		785 (100.0%)	0 (0.0%)
7	dia [factor]	Diabetes ('dia', n (%), 0=No, 1=Si)	1. No 2. Si	700 ( 89.2% ) 85 ( 10.8% )		785 (100.0%)	0 (0.0%)
8	hyp [factor]	Hipotensión (PAS<100 mmHg) ('hyp', n (%), 0=No, 1=Si)	1. No 2. Si	745 ( 94.9% ) 40 ( 5.1% )		785 (100.0%)	0 (0.0%)
9	hrt [factor]	Frec. cardíaca (Taquicardia: pulso>80) ('hrt', n (%), 0=No, 1=Si)	1. No 2. Si	576 ( 73.4% ) 209 ( 26.6% )		785 (100.0%)	0 (0.0%)
10	ant [factor]	Infarto de localización anterior ('ant', n (%), 0=No, 1=Si)	1. No 2. Si	506 ( 64.5% ) 279 ( 35.5% )		785 (100.0%)	0 (0.0%)
11	pmi [factor]	Infarto del miocardio previo ('pmi', n (%), 0=No, 1=Si)	1. No 2. Si	645 ( 82.2% ) 140 ( 17.8% )		785 (100.0%)	0 (0.0%)
12	hei [numeric]	Altura ('hei', en cm)	Mean (sd) : 169.2 (9) min ≤ med ≤ max: 141.3 ≤ 170 ≤ 197 IQR (CV) : 12 (0.1)	128 distinct values		785 (100.0%)	0 (0.0%)
13	wei [numeric]	Peso ('wei', en kg)	Mean (sd) : 75.2 (14) min ≤ med ≤ max: 40 ≤ 75 ≤ 128 IQR (CV) : 18 (0.2)	111 distinct values		785 (100.0%)	0 (0.0%)
14	htn [factor]	Historia de hipertensión ('htn', n (%), 0=No, 1=Si)	1. No 2. Si	490 ( 62.4% ) 295 ( 37.6% )		785 (100.0%)	0 (0.0%)
15	smk [factor]	Tabaquismo ('smk', n (%), 1=nunca, 2=exfumador, 3=activo)	1. nunca 2. exfumador 3. activo	262 ( 33.4% ) 306 ( 39.0% ) 217 ( 27.6% )		785 (100.0%)	0 (0.0%)
16	lip [factor]	Lípidos: Hipercolesterolemia ('lip', n (%), 0=No, 1=Si)	1. No 2. Si	480 ( 61.1% ) 305 ( 38.9% )		785 (100.0%)	0 (0.0%)
17	pan [factor]	Angina de pecho previa ('pan', n (%), 0=No, 1=Si)	1. No 2. Si	490 ( 62.4% ) 295 ( 37.6% )		785 (100.0%)	0 (0.0%)
18	fam [factor]	Historia familiar de Infarto del miocardio ('fam', n (%), 0=No, 1=Si)	1. No 2. Si	460 ( 58.6% ) 325 ( 41.4% )		785 (100.0%)	0 (0.0%)
19	ste [numeric]	Elevación del segmento ST (Derivaciones)	Mean (sd) : 4.3 (1.8) min ≤ med ≤ max: 0 ≤ 4 ≤ 10 IQR (CV) : 3 (0.4)	11 distinct values		785 (100.0%)	0 (0.0%)
20	st4 [factor]	Elevación ST en ECG > 4 derivaciones ('st4', n (%), 0=No, 1=Si)	1. No 2. Si	460 ( 58.6% ) 325 ( 41.4% )		785 (100.0%)	0 (0.0%)
21	ttr [factor]	Tiempo al alivio del dolor de pecho > 1 hora ('ttr', n (%), 0=No, 1=Si)	1. No 2. Si	390 ( 49.7% ) 395 ( 50.3% )		785 (100.0%)	0 (0.0%)
22	esamp [numeric]		1 distinct value	1 : 785 ( 100.0% )		785 (100.0%)	0 (0.0%)
23	grpl [numeric]		1 distinct value	4 : 785 ( 100.0% )		785 (100.0%)	0 (0.0%)
24	grps [numeric]		Mean (sd) : 9.9 (0.9) min ≤ med ≤ max: 9 ≤ 10 ≤ 11 IQR (CV) : 2 (0.1)	9 : 327 ( 41.7% ) 10 : 195 ( 24.8% ) 11 : 263 ( 33.5% )		785 (100.0%)	0 (0.0%)
25	regl [numeric]		1 distinct value	2 : 785 ( 100.0% )		785 (100.0%)	0 (0.0%)
26	hig [numeric]		Min : 0 Mean : 0.5 Max : 1	0 : 421 ( 53.6% ) 1 : 364 ( 46.4% )		785 (100.0%)	0 (0.0%)

Generated by summarytools 1.0.0 (R version 4.1.1)
2022-04-06

Del total de los pacientes, se presentó el evento de muerte a los 30 días luego de IAM en 6.6%. El 25.9% de los pacientes eran de sexo femenino. La media de edad de los pacientes fue de 61.8 años, miínimo de 28.4 años y un máximo de 86 años. El 41.7% del total de los pacientes tiene mas de 65 años. En cuanto a la clasificacion Killip, el 78.5% estan en la categoria de Killip 1, el 18.9% en killip 2, el 2.5% killip 3 y el 0.1% killip 4, por tanto la mayoria (97.3%) tiene killip 1 y 2. El 10.8% de los paciente tienen DM, el 5.1% hipotension, el 27.6% taquicardia, el 37.6% hipertensioión, el 38.9% dislipidemia. El 35.5% tienen infarto de localización anterior, el 17.8% tenían antecedente de infarto del miocardio. La altura media de los pacientes fue de 169 cm, el peso medio fue de 75.2Kg. En cuanto al consumo de tabaco, el 33.4% nunca habían fumado, el 39% eran exfumadores y el 27.6% fumadores activos. El 37.6% habían presentado angina previamente y el 41.4% tienen antecedente de historia familiar de infarto. En el EKG el 41.4% tenían más de 4 derivaciones comprometidas y más de la mitad de los pacientes (50.3%) tenían más de una hora de dolor precordial.

b) Estudiar algunas tablas cruzadas: resultado * predictor.

La Tabla 2 contiene una descripción de las características demográficas y clínicas de las variables de la muestra, discriminadas por los niveles de la variable de desenlace (Mortalidad a los 30 días después del IAM) (0=No, 1=Sí).

Tabla 2. Distribuciones bivariadas y univariadas de las variables.

Covariables	No (N=733)	Si (N=52)	Total (N=785)
Género ('sex', n (%), 0=Masculino, 1=Femenino)
Masculino	548 (74.8%)	34 (65.4%)	582 (74.1%)
Femenino	185 (25.2%)	18 (34.6%)	203 (25.9%)
Edad ('age', en años)
Mean (SD)	61.1 (11.0)	71.5 (8.72)	61.8 (11.1)
Median [Min, Max]	62.0 [28.4, 86.0]	73.8 [47.7, 85.8]	62.8 [28.4, 86.0]
Edad > 65 años ('a65', n (%), 0=No, 1=Si)
No	446 (60.8%)	12 (23.1%)	458 (58.3%)
Si	287 (39.2%)	40 (76.9%)	327 (41.7%)
Tabaquismo ('smk', n (%), 1=nunca, 2=exfumador, 3=activo)
nunca	253 (34.5%)	9 (17.3%)	262 (33.4%)
exfumador	282 (38.5%)	24 (46.2%)	306 (39.0%)
activo	198 (27.0%)	19 (36.5%)	217 (27.6%)
Altura ('hei', en cm)
Mean (SD)	170 (8.81)	165 (10.1)	169 (8.95)
Median [Min, Max]	170 [145, 197]	168 [141, 183]	170 [141, 197]
Peso ('wei', en kg)
Mean (SD)	75.7 (13.9)	68.2 (14.1)	75.2 (14.0)
Median [Min, Max]	75.0 [40.0, 128]	68.4 [40.0, 110]	75.0 [40.0, 128]
Diabetes ('dia', n (%), 0=No, 1=Si)
No	657 (89.6%)	43 (82.7%)	700 (89.2%)
Si	76 (10.4%)	9 (17.3%)	85 (10.8%)
Hipotensión (PAS<100 mmHg) ('hyp', n (%), 0=No, 1=Si)
No	698 (95.2%)	47 (90.4%)	745 (94.9%)
Si	35 (4.8%)	5 (9.6%)	40 (5.1%)
Historia de hipertensión ('htn', n (%), 0=No, 1=Si)
No	463 (63.2%)	27 (51.9%)	490 (62.4%)
Si	270 (36.8%)	25 (48.1%)	295 (37.6%)
Lípidos: Hipercolesterolemia ('lip', n (%), 0=No, 1=Si)
No	446 (60.8%)	34 (65.4%)	480 (61.1%)
Si	287 (39.2%)	18 (34.6%)	305 (38.9%)
Clase Killip de función ventricular izq. ('killip', n (%), Clases de 1 a 4)
Clase 1	589 (80.4%)	27 (51.9%)	616 (78.5%)
Clase 2	130 (17.7%)	18 (34.6%)	148 (18.9%)
Clase 3	14 (1.9%)	6 (11.5%)	20 (2.5%)
Clase 4	0 (0%)	1 (1.9%)	1 (0.1%)
Shock: Clase Killip 3 o 4 ('sho', n (%), 0=No, 1=Si)
No	719 (98.1%)	45 (86.5%)	764 (97.3%)
Si	14 (1.9%)	7 (13.5%)	21 (2.7%)
Frec. cardíaca (Taquicardia: pulso>80) ('hrt', n (%), 0=No, 1=Si)
No	551 (75.2%)	25 (48.1%)	576 (73.4%)
Si	182 (24.8%)	27 (51.9%)	209 (26.6%)
Infarto de localización anterior ('ant', n (%), 0=No, 1=Si)
No	480 (65.5%)	26 (50.0%)	506 (64.5%)
Si	253 (34.5%)	26 (50.0%)	279 (35.5%)
Infarto del miocardio previo ('pmi', n (%), 0=No, 1=Si)
No	612 (83.5%)	33 (63.5%)	645 (82.2%)
Si	121 (16.5%)	19 (36.5%)	140 (17.8%)
Angina de pecho previa ('pan', n (%), 0=No, 1=Si)
No	463 (63.2%)	27 (51.9%)	490 (62.4%)
Si	270 (36.8%)	25 (48.1%)	295 (37.6%)
Historia familiar de Infarto del miocardio ('fam', n (%), 0=No, 1=Si)
No	426 (58.1%)	34 (65.4%)	460 (58.6%)
Si	307 (41.9%)	18 (34.6%)	325 (41.4%)
Elevación del segmento ST (Derivaciones)
Mean (SD)	4.25 (1.80)	4.44 (1.83)	4.27 (1.80)
Median [Min, Max]	4.00 [0, 10.0]	4.50 [1.00, 8.00]	4.00 [0, 10.0]
Elevación ST en ECG > 4 derivaciones ('st4', n (%), 0=No, 1=Si)
No	434 (59.2%)	26 (50.0%)	460 (58.6%)
Si	299 (40.8%)	26 (50.0%)	325 (41.4%)
Tiempo al alivio del dolor de pecho > 1 hora ('ttr', n (%), 0=No, 1=Si)
No	372 (50.8%)	18 (34.6%)	390 (49.7%)
Si	361 (49.2%)	34 (65.4%)	395 (50.3%)

En cuanto a las características de la población, discriminadas por los niveles de la variable de desenlace. Específicamente para aquellos que desarrollaron el evento de muerte a los 30 días luego del IAM: El 65.4% de fueron se sexo masculino, la Edad Media fue de 71.5 años vs 61.1 del grupo que no tuvieron en desenlace de interes. El 76.9% tenian mas de 65 años. En cuanto a la clasificación Killip en los pacientes con el desenlace, el 51.9% están en la categoría de Killip 1, el 34.6% en killip 2, el 11.5% killip 3 y el 1.9% killip 4, el 86.5% tenían killip 1 y 2. El 82.7 de los paciente con el desenlace no tenían DM, el 9.6% no tenían hipotensión, el 51.9% taquicardia, el 48.1% hipertensión, el 34.6% dislipidemia. El 50% tienen infarto de localización anterior, el 34.6% tenían antecedente de infarto del miocardio. La altura media de los pacientes fue de 165 cm, el peso medio fue de 68.5Kg. En cuanto al consumo de tabaco, el 17.3% nunca habían fumado, el 46.2% eran exfumadores y el 36.5% eran fumadores activos. El 48.1% habían presentado angina previamente y el 34.6% tienen antecedente de historia familiar de infarto. En el EKG el 50% tenían mas de 4 derivaciones comprometidas y el 65.4% de los pacientes tenían más de una hora de dolor precordial.

c) ¿Cuáles son las relaciones univariadas entre los predictores y el resultado (DÍA30)?

	OR	2.5 %	97.5 %
(Intercept)	0.06	0.04	0.09
sexFemenino	1.57	0.85	2.81

	OR	2.5 %	97.5 %
(Intercept)	0.00	0.00	0.00
age	1.11	1.08	1.15

	OR	2.5 %	97.5 %
(Intercept)	0.03	0.01	0.05
a65Si	5.18	2.75	10.47

	OR	2.5 %	97.5 %
(Intercept)	0.04	0.02	0.07
smkexfumador	2.39	1.13	5.53
smkactivo	2.70	1.23	6.38

	OR	2.5 %	97.5 %
(Intercept)	355.39	2.20	54888.30
hei	0.95	0.92	0.98

	OR	2.5 %	97.5 %
(Intercept)	1.42	0.30	6.63
wei	0.96	0.94	0.98

	OR	2.5 %	97.5 %
(Intercept)	0.07	0.05	0.09
diaSi	1.81	0.80	3.70

	OR	2.5 %	97.5 %
(Intercept)	0.07	0.05	0.09
hypSi	2.12	0.70	5.23

	OR	2.5 %	97.5 %
(Intercept)	0.06	0.04	0.08
htnSi	1.59	0.90	2.80

	OR	2.5 %	97.5 %
(Intercept)	0.08	0.05	0.11
lipSi	0.82	0.45	1.47

	OR	2.5 %	97.5 %
(Intercept)	0.05	0.03	0.07
killipClase 2	3.02	1.59	5.61
killipClase 3	9.35	3.11	25.35
killipClase 4	46207703.95	0.00

	OR	2.5 %	97.5 %
(Intercept)	0.06	0.05	0.08
shoSi	7.99	2.90	20.24

	OR	2.5 %	97.5 %
(Intercept)	0.05	0.03	0.07
hrtSi	3.27	1.85	5.81

	OR	2.5 %	97.5 %
(Intercept)	0.05	0.04	0.08
antSi	1.90	1.08	3.35

	OR	2.5 %	97.5 %
(Intercept)	0.05	0.04	0.08
pmiSi	2.91	1.58	5.24

	OR	2.5 %	97.5 %
(Intercept)	0.06	0.04	0.08
panSi	1.59	0.90	2.80

	OR	2.5 %	97.5 %
(Intercept)	0.08	0.06	0.11
famSi	0.73	0.40	1.31

	OR	2.5 %	97.5 %
(Intercept)	0.06	0.03	0.11
ste	1.06	0.91	1.23

	OR	2.5 %	97.5 %
(Intercept)	0.06	0.04	0.09
st4Si	1.45	0.82	2.56

	OR	2.5 %	97.5 %
(Intercept)	0.05	0.03	0.08
ttrSi	1.95	1.09	3.58

Para indicar las relaciones univaridas se realiza una análisis de los modelos logísticos para cada predictor vs el desenlace. Se presentan las tablas de contigencia 2 x 2 para evaluar la relación entre cada una de las variables tenidas en cuenta y el desarrollo de la muerte dentro del mes luego del IAM

Tablas de contingencia

ctable(mydata$sex, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
sex * day30  
Data Frame: mydata  


----------- ------- ------------- ----------- --------------
              day30            No          Si          Total
        sex                                                 
  Masculino           548 (94.2%)   34 (5.8%)   582 (100.0%)
   Femenino           185 (91.1%)   18 (8.9%)   203 (100.0%)
      Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
----------- ------- ------------- ----------- --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    1.57        0.86       2.84   
----------------------------------

# ctable(mydata$age, mydata$day30, OR=TRUE)

ctable(mydata$a65, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
a65 * day30  
Data Frame: mydata  


------- ------- ------------- ------------ --------------
          day30            No           Si          Total
    a65                                                  
     No           446 (97.4%)   12 ( 2.6%)   458 (100.0%)
     Si           287 (87.8%)   40 (12.2%)   327 (100.0%)
  Total           733 (93.4%)   52 ( 6.6%)   785 (100.0%)
------- ------- ------------- ------------ --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    5.18        2.67      10.04   
----------------------------------

ctable(mydata$smk, mydata$day30, OR=TRUE)

OR and RR can only be used with 2 x 2 tables; parameter(s) ignored

Cross-Tabulation, Row Proportions  
smk * day30  
Data Frame: mydata  

----------- ------- ------------- ----------- --------------
              day30            No          Si          Total
        smk                                                 
      nunca           253 (96.6%)    9 (3.4%)   262 (100.0%)
  exfumador           282 (92.2%)   24 (7.8%)   306 (100.0%)
     activo           198 (91.2%)   19 (8.8%)   217 (100.0%)
      Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
----------- ------- ------------- ----------- --------------

# ctable(mydata$hei, mydata$day30, OR=TRUE)

# ctable(mydata$wei, mydata$day30, OR=TRUE)

ctable(mydata$dia, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
dia * day30  
Data Frame: mydata  


------- ------- ------------- ------------ --------------
          day30            No           Si          Total
    dia                                                  
     No           657 (93.9%)   43 ( 6.1%)   700 (100.0%)
     Si            76 (89.4%)    9 (10.6%)    85 (100.0%)
  Total           733 (93.4%)   52 ( 6.6%)   785 (100.0%)
------- ------- ------------- ------------ --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    1.81        0.85       3.86   
----------------------------------

ctable(mydata$hyp, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
hyp * day30  
Data Frame: mydata  


------- ------- ------------- ------------ --------------
          day30            No           Si          Total
    hyp                                                  
     No           698 (93.7%)   47 ( 6.3%)   745 (100.0%)
     Si            35 (87.5%)    5 (12.5%)    40 (100.0%)
  Total           733 (93.4%)   52 ( 6.6%)   785 (100.0%)
------- ------- ------------- ------------ --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    2.12        0.79       5.67   
----------------------------------

ctable(mydata$htn, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
htn * day30  
Data Frame: mydata  


------- ------- ------------- ----------- --------------
          day30            No          Si          Total
    htn                                                 
     No           463 (94.5%)   27 (5.5%)   490 (100.0%)
     Si           270 (91.5%)   25 (8.5%)   295 (100.0%)
  Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
------- ------- ------------- ----------- --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    1.59        0.90       2.79   
----------------------------------

ctable(mydata$lip, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
lip * day30  
Data Frame: mydata  


------- ------- ------------- ----------- --------------
          day30            No          Si          Total
    lip                                                 
     No           446 (92.9%)   34 (7.1%)   480 (100.0%)
     Si           287 (94.1%)   18 (5.9%)   305 (100.0%)
  Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
------- ------- ------------- ----------- --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    0.82        0.46       1.48   
----------------------------------

ctable(mydata$hrt, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
hrt * day30  
Data Frame: mydata  


------- ------- ------------- ------------ --------------
          day30            No           Si          Total
    hrt                                                  
     No           551 (95.7%)   25 ( 4.3%)   576 (100.0%)
     Si           182 (87.1%)   27 (12.9%)   209 (100.0%)
  Total           733 (93.4%)   52 ( 6.6%)   785 (100.0%)
------- ------- ------------- ------------ --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    3.27        1.85       5.78   
----------------------------------

ctable(mydata$killip, mydata$day30, OR=TRUE)

OR and RR can only be used with 2 x 2 tables; parameter(s) ignored

Cross-Tabulation, Row Proportions  
killip * day30  
Data Frame: mydata  

--------- ------- ------------- ------------- --------------
            day30            No            Si          Total
   killip                                                   
  Clase 1           589 (95.6%)   27 (  4.4%)   616 (100.0%)
  Clase 2           130 (87.8%)   18 ( 12.2%)   148 (100.0%)
  Clase 3            14 (70.0%)    6 ( 30.0%)    20 (100.0%)
  Clase 4             0 ( 0.0%)    1 (100.0%)     1 (100.0%)
    Total           733 (93.4%)   52 (  6.6%)   785 (100.0%)
--------- ------- ------------- ------------- --------------

ctable(mydata$sho, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
sho * day30  
Data Frame: mydata  


------- ------- ------------- ------------ --------------
          day30            No           Si          Total
    sho                                                  
     No           719 (94.1%)   45 ( 5.9%)   764 (100.0%)
     Si            14 (66.7%)    7 (33.3%)    21 (100.0%)
  Total           733 (93.4%)   52 ( 6.6%)   785 (100.0%)
------- ------- ------------- ------------ --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    7.99        3.07      20.78   
----------------------------------

ctable(mydata$ant, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
ant * day30  
Data Frame: mydata  


------- ------- ------------- ----------- --------------
          day30            No          Si          Total
    ant                                                 
     No           480 (94.9%)   26 (5.1%)   506 (100.0%)
     Si           253 (90.7%)   26 (9.3%)   279 (100.0%)
  Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
------- ------- ------------- ----------- --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    1.90        1.08       3.34   
----------------------------------

ctable(mydata$pmi, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
pmi * day30  
Data Frame: mydata  


------- ------- ------------- ------------ --------------
          day30            No           Si          Total
    pmi                                                  
     No           612 (94.9%)   33 ( 5.1%)   645 (100.0%)
     Si           121 (86.4%)   19 (13.6%)   140 (100.0%)
  Total           733 (93.4%)   52 ( 6.6%)   785 (100.0%)
------- ------- ------------- ------------ --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    2.91        1.60       5.29   
----------------------------------

ctable(mydata$pan, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
pan * day30  
Data Frame: mydata  


------- ------- ------------- ----------- --------------
          day30            No          Si          Total
    pan                                                 
     No           463 (94.5%)   27 (5.5%)   490 (100.0%)
     Si           270 (91.5%)   25 (8.5%)   295 (100.0%)
  Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
------- ------- ------------- ----------- --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    1.59        0.90       2.79   
----------------------------------

ctable(mydata$fam, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
fam * day30  
Data Frame: mydata  


------- ------- ------------- ----------- --------------
          day30            No          Si          Total
    fam                                                 
     No           426 (92.6%)   34 (7.4%)   460 (100.0%)
     Si           307 (94.5%)   18 (5.5%)   325 (100.0%)
  Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
------- ------- ------------- ----------- --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    0.73        0.41       1.33   
----------------------------------

ctable(mydata$st4, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
st4 * day30  
Data Frame: mydata  


------- ------- ------------- ----------- --------------
          day30            No          Si          Total
    st4                                                 
     No           434 (94.3%)   26 (5.7%)   460 (100.0%)
     Si           299 (92.0%)   26 (8.0%)   325 (100.0%)
  Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
------- ------- ------------- ----------- --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    1.45        0.83       2.55   
----------------------------------

ctable(mydata$ttr, mydata$day30, OR=TRUE)

Cross-Tabulation, Row Proportions  
ttr * day30  
Data Frame: mydata  


------- ------- ------------- ----------- --------------
          day30            No          Si          Total
    ttr                                                 
     No           372 (95.4%)   18 (4.6%)   390 (100.0%)
     Si           361 (91.4%)   34 (8.6%)   395 (100.0%)
  Total           733 (93.4%)   52 (6.6%)   785 (100.0%)
------- ------- ------------- ----------- --------------

----------------------------------
 Odds Ratio   Lo - 95%   Hi - 95% 
------------ ---------- ----------
    1.95        1.08       3.51   
----------------------------------

La Tabla 3 presenta los resultados obtenidos. Se presenta la pendiente del modelo individual y los OR de cada relación.

Tabla 3. Resumen de regresión logística de covariables vs desenlace.

No	Covariable	Pendiente Modelo Log	¿Pendiente significativa?	OR	IC 95% del OR
1.	sexFemenino	0.450	NO	1.57	0.85-2.81
2.	age	0.1062	SI	1.11	1.08-1.15
3.	a65Si	1.645	SI	5.18	2.75-10.47
4.	smkex	0.872	SI	2.39	1.33-5.53
	smkactual	0.992	SI	2.7	1.23-6.38
5.	hei	-0.0509	SI	0.95	0.92-0.98
6.	wei	-0.0417	SI	0.96	0.94-0.98
7.	diaSi	0.593	NO	1.81	0.8-3.70
8.	hypSi	0.752	NO	1.57	0.7-5.23
9.	htnSi	0.462	NO	1.59	0.90-2.8
10.	lipSi	-0.195	NO	0.82	0.45-1.47
11.	killipClase 2	1.105	SI	3.02	1.59-5.61
	killipClase 3	2.235	SI	9.35	3.11-25.35
	killipClase4	17.649	NO	46207703.95	NA
12.	shoSi	2.078	SI	7.99	2.9-20.24
13.	hrtSi	1.185	SI	3.27	1.85-5.81
14.	antSi	0.640	SI	1.9	1.08-3.35
15.	pmiSi	1.069	SI	2.9	1.58-5.24
16.	panSi	0.462	NO	1.59	0.9-2.8
17.	famSi	-0.308	NO	0.73	0.4-1.31
18.	ste	0.0574	NO	NA	NA
19.	st4Si	0.373	NO	1.45	0.82-2.56
20.	ttr	0.666	SI	1.59	1.09-3.58

Se observa que 10 de las 20 variables tuvieron regresión logística significativa y además esas mismas variables tuvieron OR importantes con IC que se correlacionan con la significia y la pendiente. Esas 10 variables por tanto son las primeras candidatas a variables pronósticas para el desenlace.

Pregunta 2: Modelado de decisiones y modelo preliminar

a) Decida si utilizará versiones continuas o dicotómicas de: edad (age/A65), clase Killip (Killip/SHO) y elevación del ST (STE/ST4).

Entre la variable age y A65, considero que se debe aceptar la recomendación de no dicotomizar variables continuas. Además ambas variables dentro del análisis univariado fueron un factor de riesgo para desarrollar el desenlace, por tanto se continúa con la versión continua de la edad.

En cuanto a las variables Killip y SHO, dentro de la caracterización se encuentra que Killip 4 solo tiene sólo un dato por lo que es difícil de interpretar y de manejar, por tanto se continúa con la versión dicotomizada.

Finalmente, cuando se comparan las variables STE y ST4 se prefiere trabajar con la variable dicotomizada por la dificultad que implica el manejo de las 12 derivaciones y sus resultados.

b) Decida su estrategia de modelado para la definición del modelo de regresión multivariable, por ejemplo, un modelo completo que incluye todos los predictores potenciales, paso a paso hacia atrás con un valor p especificado por el usuario (p. ej., 0,15), etc.

Variables para modelar e incluir

Se decidió incluir las 10 variables descritas que en la Tabla 3 que presenta el análisis univariado y resultaron significativas en sus respectivos modelos logísticos frente a la variable desenlace, además se tomó en cuenta sus OR y respectivos intervalos de confianza y por ultimo y muy importante la plausibilidad clínica. La variable fam se incluyó por ser un claro factor de riesgo para muerte por causas cardiovasculares. Adicionamos la varialbe st4 porque como se dijo, entre mas derivaciones en el EKG haya alteradas, indica mayor tejido comprometido. Se incluyó la variable indice de masa corporal (IMC) porque relaciona las variables peso y altura, además se sabe que la obesidad es un gran factor de riesgos cardivascular.

Desenlace: day30

Variables pronósticas candidatas: age, smk, hei, wei, sho, hrt, fam, ant, pmi, ttr, st4, IMC.

Insertar la variable IMC

mydata$IMC <- mydata$wei / ( (mydata$hei/100)^2 ) 
summary(mydata$IMC)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.6    23.6    25.8    26.2    28.4    44.2 

Supuesto de linealidad (para los predictores continuos significativos)

A continuación se evalúo el supuesto de linealidad para los predictores continuos significativos encontrados: age, hei, wei, IMC.

# Generemos un modelo completo para determinar el riesgo condicional de tener el evento en cada participante
modelocompleto <- glm(day30 ~ age + hei + wei + IMC, 
                      data = mydata, family = binomial)

# Calculemos la probabilidad y el log (Odds) para cada participante
mydata$logodds = predict(modelocompleto, mydata)
mydata$prob = predict(modelocompleto, mydata, type="response")

# Evaluemos el supuesto de linealidad
library(ggplot2)

par(mfrow=c(2,2))  # Función para presentar los gráficos c(nrows, ncols)

ggplot(mydata, aes(x=age, y=logodds)) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_smooth(method = "loess") + 
  theme_bw() 

`geom_smooth()` using formula 'y ~ x'

ggplot(mydata, aes(x=hei, y=logodds)) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_smooth(method = "loess") + 
  theme_bw() 

`geom_smooth()` using formula 'y ~ x'

ggplot(mydata, aes(x=wei, y=logodds)) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_smooth(method = "loess") + 
  theme_bw() 

`geom_smooth()` using formula 'y ~ x'

ggplot(mydata, aes(x=IMC, y=logodds)) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_smooth(method = "loess") + 
  theme_bw() 

`geom_smooth()` using formula 'y ~ x'

par(mfrow=c(1,1))  # Función para presentar los gráficos c(nrows, ncols)

Al evaluar el supuesto de linealidad de las variables y predictores continuos candidatos, se observa que sólo la variable edad (age) lo cumple.

Métodos de selección por pasos

Método Backward

Comenzamos con un modelo Full que incluya todas las variables que elegimos para modelar, y posteriormente se eliminan variables para llegar al modelo final basándonos en el criterio de selección de modelos AIC, y usando un nivel de significancia de 0.15.

#library(Hmisc)
#library(rms)   # Regression Modeling Strategies 

# Modelo completo
#full <- lrm(day30 ~ age + smk + hei + wei + sho + hrt + fam
#            + ant + pmi + ttr + st4 + IMC,
#            data=mydata, 
#            x=TRUE, y=TRUE,
#            linear.predictors=FALSE, tol=1e-10)

#full

## Criterio valor p < 0.15 para mantener en el modelo
#validate (full, B=100, bw=TRUE ,rule = "p" , sls=.10 , type = "individual")

## Criterio AIC:
#validate(full, B=200, bw=TRUE, rule = "aic", sls=.15, type = "individual")
# Nota: sls is the significance level for a factor to be kept in a model, or for judging the residual chi-square.

library(stats)  #Librería "base" cargada por defecto al entrar al R

## BACKWARD
modelofull <- glm(day30 ~ age + smk + hei + wei + sho + hrt + fam
            + ant + pmi + ttr + st4 + IMC,
            data=mydata, 
            family = binomial())

summary(modelofull)

Call:
glm(formula = day30 ~ age + smk + hei + wei + sho + hrt + fam + 
    ant + pmi + ttr + st4 + IMC, family = binomial(), data = mydata)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.550  -0.367  -0.215  -0.119   3.148  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -38.7688    16.3214   -2.38    0.018 *  
age            0.0871     0.0209    4.18    3e-05 ***
smkexfumador   0.1224     0.4568    0.27    0.789    
smkactivo      0.1742     0.4728    0.37    0.712    
hei            0.1828     0.0984    1.86    0.063 .  
wei           -0.2523     0.1188   -2.12    0.034 *  
shoSi          1.0696     0.5676    1.88    0.060 .  
hrtSi          0.9787     0.3297    2.97    0.003 ** 
famSi         -0.2817     0.3357   -0.84    0.401    
antSi          0.4973     0.3568    1.39    0.163    
pmiSi          0.5834     0.3482    1.68    0.094 .  
ttrSi          0.7951     0.3363    2.36    0.018 *  
st4Si         -0.0935     0.3619   -0.26    0.796    
IMC            0.6489     0.3147    2.06    0.039 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 382.78  on 784  degrees of freedom
Residual deviance: 299.37  on 771  degrees of freedom
AIC: 327.4

Number of Fisher Scoring iterations: 7

step(modelofull, direction = "backward")

Start:  AIC=327.37
day30 ~ age + smk + hei + wei + sho + hrt + fam + ant + pmi + 
    ttr + st4 + IMC

       Df Deviance   AIC
- smk   2    299.5 323.5
- st4   1    299.4 325.4
- fam   1    300.1 326.1
- ant   1    301.3 327.3
<none>       299.4 327.4
- pmi   1    302.1 328.1
- sho   1    302.7 328.7
- hei   1    302.8 328.8
- IMC   1    303.4 329.4
- wei   1    303.8 329.8
- ttr   1    305.2 331.2
- hrt   1    308.0 334.0
- age   1    319.8 345.8

Step:  AIC=323.51
day30 ~ age + hei + wei + sho + hrt + fam + ant + pmi + ttr + 
    st4 + IMC

       Df Deviance   AIC
- st4   1    299.6 321.6
- fam   1    300.2 322.2
- ant   1    301.4 323.4
<none>       299.5 323.5
- pmi   1    302.5 324.5
- sho   1    302.8 324.8
- hei   1    303.0 325.0
- IMC   1    303.7 325.7
- wei   1    304.0 326.0
- ttr   1    305.5 327.5
- hrt   1    308.1 330.1
- age   1    323.3 345.3

Step:  AIC=321.58
day30 ~ age + hei + wei + sho + hrt + fam + ant + pmi + ttr + 
    IMC

       Df Deviance   AIC
- fam   1    300.3 320.3
- ant   1    301.6 321.6
<none>       299.6 321.6
- pmi   1    302.6 322.6
- sho   1    302.9 322.9
- hei   1    303.0 323.0
- IMC   1    303.7 323.7
- wei   1    304.1 324.1
- ttr   1    305.5 325.5
- hrt   1    308.1 328.1
- age   1    323.4 343.4

Step:  AIC=320.3
day30 ~ age + hei + wei + sho + hrt + ant + pmi + ttr + IMC

       Df Deviance   AIC
- ant   1    302.2 320.2
<none>       300.3 320.3
- pmi   1    303.0 321.0
- sho   1    303.6 321.6
- hei   1    303.7 321.7
- IMC   1    304.4 322.4
- wei   1    304.7 322.7
- ttr   1    306.2 324.2
- hrt   1    308.6 326.6
- age   1    326.2 344.2

Step:  AIC=320.24
day30 ~ age + hei + wei + sho + hrt + pmi + ttr + IMC

       Df Deviance   AIC
<none>       302.2 320.2
- pmi   1    305.3 321.3
- sho   1    305.5 321.5
- hei   1    305.6 321.6
- IMC   1    306.3 322.3
- wei   1    306.5 322.5
- ttr   1    308.9 324.9
- hrt   1    312.2 328.2
- age   1    328.7 344.7

Call:  glm(formula = day30 ~ age + hei + wei + sho + hrt + pmi + ttr + 
    IMC, family = binomial(), data = mydata)

Coefficients:
(Intercept)          age          hei          wei        shoSi  
   -38.9856       0.0932       0.1805      -0.2475       1.0247  
      hrtSi        pmiSi        ttrSi          IMC  
     1.0134       0.5996       0.8311       0.6477  

Degrees of Freedom: 784 Total (i.e. Null);  776 Residual
Null Deviance:      383 
Residual Deviance: 302  AIC: 320

El método Backward arrojó el siguiente modelo como el elegido:

Variables candidatas del modelo elegido por metodología “backward” [1] age hei wei sho hrt pmi ttr IMC

\[ \hat{day30} = \beta_0 + \beta_1 age + \beta_2 hei + \beta_3 wei + \beta_4 sho + \beta_5 hrt + \beta_ 6 pmi + \beta_7 ttr + \beta_8 IMC \]

Método Forward

Comenzamos con un modelo Null o vacío que no incluye las variables que elegimos para modelar, y el método “hacia adelante” nos lleva a un modelo final, basándose en el criterio de selección de modelos AIC.

## FORWARD
modeloEmpty <- glm(day30 ~ 1, data=mydata, family=binomial())
summary(modeloEmpty)

Call:
glm(formula = day30 ~ 1, family = binomial(), data = mydata)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
 -0.37   -0.37   -0.37   -0.37    2.33  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -2.646      0.144   -18.4   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 382.78  on 784  degrees of freedom
Residual deviance: 382.78  on 784  degrees of freedom
AIC: 384.8

Number of Fisher Scoring iterations: 5

forward <- step(modeloEmpty, direction = "forward", trace = 1, 
                scope =  ~ age + smk + hei + wei + sho + hrt + fam
            + ant + pmi + ttr + st4 + IMC)

Start:  AIC=384.78
day30 ~ 1

       Df Deviance   AIC
+ age   1    334.5 338.5
+ hrt   1    366.6 370.6
+ wei   1    368.0 372.0
+ sho   1    368.9 372.9
+ pmi   1    371.7 375.7
+ hei   1    372.1 376.1
+ IMC   1    377.3 381.3
+ smk   2    375.5 381.5
+ ttr   1    377.6 381.6
+ ant   1    377.9 381.9
<none>       382.8 384.8
+ st4   1    381.1 385.1
+ fam   1    381.7 385.7

Step:  AIC=338.55
day30 ~ age

       Df Deviance   AIC
+ hrt   1    321.7 327.7
+ ttr   1    328.5 334.5
+ pmi   1    328.6 334.6
+ sho   1    328.9 334.9
+ ant   1    329.8 335.8
<none>       334.5 338.5
+ hei   1    332.8 338.8
+ st4   1    333.1 339.1
+ wei   1    333.2 339.2
+ fam   1    334.4 340.4
+ IMC   1    334.5 340.5
+ smk   2    334.3 342.3

Step:  AIC=327.72
day30 ~ age + hrt

       Df Deviance   AIC
+ ttr   1    315.0 323.0
+ pmi   1    317.0 325.0
+ sho   1    318.2 326.2
+ ant   1    319.0 327.0
<none>       321.7 327.7
+ wei   1    320.2 328.2
+ hei   1    320.5 328.5
+ st4   1    321.3 329.3
+ IMC   1    321.5 329.5
+ fam   1    321.5 329.5
+ smk   2    321.1 331.1

Step:  AIC=323.01
day30 ~ age + hrt + ttr

       Df Deviance   AIC
+ pmi   1    310.9 320.9
+ sho   1    312.2 322.2
<none>       315.0 323.0
+ ant   1    313.2 323.2
+ hei   1    313.6 323.6
+ wei   1    313.6 323.7
+ fam   1    314.9 324.9
+ IMC   1    314.9 324.9
+ st4   1    314.9 324.9
+ smk   2    314.5 326.5

Step:  AIC=320.9
day30 ~ age + hrt + ttr + pmi

       Df Deviance   AIC
+ sho   1    308.2 320.2
<none>       310.9 320.9
+ ant   1    309.3 321.3
+ hei   1    309.5 321.5
+ wei   1    309.8 321.8
+ fam   1    310.5 322.5
+ st4   1    310.8 322.8
+ IMC   1    310.9 322.9
+ smk   2    310.7 324.7

Step:  AIC=320.21
day30 ~ age + hrt + ttr + pmi + sho

       Df Deviance   AIC
<none>       308.2 320.2
+ hei   1    306.6 320.6
+ ant   1    306.7 320.7
+ wei   1    306.8 320.7
+ fam   1    307.8 321.8
+ IMC   1    308.1 322.1
+ st4   1    308.1 322.1
+ smk   2    307.8 323.8

El método Forward arrojó el siguiente modelo como el elegido:

Variables candidatas del modelo elegido por metodología “forward” Step: AIC=320.21 day30 ~ age + hrt + ttr + pmi + sho

\[ \hat{day30} = \beta_0 + \beta_1 age + \beta_2 hrt + \beta_3 ttr + \beta_4 pmi + \beta_5 sho \]

Stepwise

Comenzamos con un modelo completo que incluya todas las variables que elegimos para modelar inicialmente, y el método “dos pasos” o “hacia adelante y hacia atrás” nos lleva a un modelo elegido, basándose en el criterio de selección de modelos AIC.

library(stats)  #Librería "base" cargada por defecto al entrar al R

## STEPWISE
modelofull <- glm(day30 ~ age + smk + hei + wei + sho + hrt + fam
            + ant + pmi + ttr + st4 + IMC,
            data=mydata, 
            family = binomial())

summary(modelofull)

Call:
glm(formula = day30 ~ age + smk + hei + wei + sho + hrt + fam + 
    ant + pmi + ttr + st4 + IMC, family = binomial(), data = mydata)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.550  -0.367  -0.215  -0.119   3.148  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -38.7688    16.3214   -2.38    0.018 *  
age            0.0871     0.0209    4.18    3e-05 ***
smkexfumador   0.1224     0.4568    0.27    0.789    
smkactivo      0.1742     0.4728    0.37    0.712    
hei            0.1828     0.0984    1.86    0.063 .  
wei           -0.2523     0.1188   -2.12    0.034 *  
shoSi          1.0696     0.5676    1.88    0.060 .  
hrtSi          0.9787     0.3297    2.97    0.003 ** 
famSi         -0.2817     0.3357   -0.84    0.401    
antSi          0.4973     0.3568    1.39    0.163    
pmiSi          0.5834     0.3482    1.68    0.094 .  
ttrSi          0.7951     0.3363    2.36    0.018 *  
st4Si         -0.0935     0.3619   -0.26    0.796    
IMC            0.6489     0.3147    2.06    0.039 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 382.78  on 784  degrees of freedom
Residual deviance: 299.37  on 771  degrees of freedom
AIC: 327.4

Number of Fisher Scoring iterations: 7

step(modelofull, direction = "both")

Start:  AIC=327.37
day30 ~ age + smk + hei + wei + sho + hrt + fam + ant + pmi + 
    ttr + st4 + IMC

       Df Deviance   AIC
- smk   2    299.5 323.5
- st4   1    299.4 325.4
- fam   1    300.1 326.1
- ant   1    301.3 327.3
<none>       299.4 327.4
- pmi   1    302.1 328.1
- sho   1    302.7 328.7
- hei   1    302.8 328.8
- IMC   1    303.4 329.4
- wei   1    303.8 329.8
- ttr   1    305.2 331.2
- hrt   1    308.0 334.0
- age   1    319.8 345.8

Step:  AIC=323.51
day30 ~ age + hei + wei + sho + hrt + fam + ant + pmi + ttr + 
    st4 + IMC

       Df Deviance   AIC
- st4   1    299.6 321.6
- fam   1    300.2 322.2
- ant   1    301.4 323.4
<none>       299.5 323.5
- pmi   1    302.5 324.5
- sho   1    302.8 324.8
- hei   1    303.0 325.0
- IMC   1    303.7 325.7
- wei   1    304.0 326.0
+ smk   2    299.4 327.4
- ttr   1    305.5 327.5
- hrt   1    308.1 330.1
- age   1    323.3 345.3

Step:  AIC=321.58
day30 ~ age + hei + wei + sho + hrt + fam + ant + pmi + ttr + 
    IMC

       Df Deviance   AIC
- fam   1    300.3 320.3
- ant   1    301.6 321.6
<none>       299.6 321.6
- pmi   1    302.6 322.6
- sho   1    302.9 322.9
- hei   1    303.0 323.0
+ st4   1    299.5 323.5
- IMC   1    303.7 323.7
- wei   1    304.1 324.1
+ smk   2    299.4 325.4
- ttr   1    305.5 325.5
- hrt   1    308.1 328.1
- age   1    323.4 343.4

Step:  AIC=320.3
day30 ~ age + hei + wei + sho + hrt + ant + pmi + ttr + IMC

       Df Deviance   AIC
- ant   1    302.2 320.2
<none>       300.3 320.3
- pmi   1    303.0 321.0
+ fam   1    299.6 321.6
- sho   1    303.6 321.6
- hei   1    303.7 321.7
+ st4   1    300.2 322.2
- IMC   1    304.4 322.4
- wei   1    304.7 322.7
+ smk   2    300.2 324.2
- ttr   1    306.2 324.2
- hrt   1    308.6 326.6
- age   1    326.2 344.2

Step:  AIC=320.24
day30 ~ age + hei + wei + sho + hrt + pmi + ttr + IMC

       Df Deviance   AIC
<none>       302.2 320.2
+ ant   1    300.3 320.3
- pmi   1    305.3 321.3
- sho   1    305.5 321.5
+ fam   1    301.6 321.6
- hei   1    305.6 321.6
+ st4   1    302.1 322.1
- IMC   1    306.3 322.3
- wei   1    306.5 322.5
+ smk   2    302.1 324.1
- ttr   1    308.9 324.9
- hrt   1    312.2 328.2
- age   1    328.7 344.7

Call:  glm(formula = day30 ~ age + hei + wei + sho + hrt + pmi + ttr + 
    IMC, family = binomial(), data = mydata)

Coefficients:
(Intercept)          age          hei          wei        shoSi  
   -38.9856       0.0932       0.1805      -0.2475       1.0247  
      hrtSi        pmiSi        ttrSi          IMC  
     1.0134       0.5996       0.8311       0.6477  

Degrees of Freedom: 784 Total (i.e. Null);  776 Residual
Null Deviance:      383 
Residual Deviance: 302  AIC: 320

El método stepwise arrojó el siguiente modelo como el elegido:

Variables candidatas del modelo elegido por metodología “stepwise”

Step: AIC=320.24 day30 ~ age + hei + wei + sho + hrt + pmi + ttr + IMC

\[ \hat{day30} = \beta_0 + \beta_1 age + \beta_2 hei + \beta_3 wei + \beta_4 sho + \beta_5 hrt + \beta_6 pmi + \beta_7 ttr + \beta_8 IMC \]

Lasso (least absolute shrinkage and selection operator)

Puede llevar a coeficientes iguales a cero, lo cual puede usarse para selección de variables en el modelo. Penaliza la suma de los valores absolutos de los coeficientes de regresión

Util con grandes conjutos de datos o número de predictores.

library(glmnet)
library(tidyverse)
#https://cran.r-project.org/web/packages/glmnet/vignettes/glmnet.pdf

myvars <- c("day30", "age", "smk", "hei", "wei", "sho", "hrt",
            "fam", "ant", "pmi", "ttr", "st4", "IMC")

mydata1 <- mydata[myvars]  #Subbase

x <- model.matrix(day30~., mydata1)[,-1]
y <- mydata1$day30

# When alpha=0, Ridge Model is fit and if alpha=1, a lasso model is fit.
cv.lasso <- cv.glmnet(x, y, alpha = 1, family = "binomial")

#Minimo valor de lambda
cv.lasso$lambda.min

[1] 0.00620995

coef(cv.lasso, cv.lasso$lambda.min)

14 x 1 sparse Matrix of class "dgCMatrix"
                      s1
(Intercept)  -6.63592749
age           0.07607772
smkexfumador  .         
smkactivo     .         
hei          -0.00783367
wei          -0.00711939
shoSi         0.89676532
hrtSi         0.76295143
famSi         .         
antSi         0.26152668
pmiSi         0.49298009
ttrSi         0.44188118
st4Si         .         
IMC           .         

El método Lasso arrojó el siguiente modelo como el elegido:

Variables candidatas del modelo elegido por metodología “Lasso”

day30 ~ age + sho + hrt + fam + ant + pmi + ttr

\[ \hat{day30} = \beta_0 + \beta_1 age + \beta_2 sho + \beta_3 hrt + \beta_4 fam + \beta_5 ant + \beta_6 pmi + \beta_7 ttr \]

Imagen 1: Comparación de modelos candidatos, generados por distintos metodos.

Coeficientes Backward

c) Muestre los coeficientes de regresión de su modelo preliminar + motivación para la inclusión de la covariable en el modelo.

Luego de calcular y comparar los modelos candidatos, se decide continuar en los análisis venideros con el modelo obtenido mediante la metología forward, el cual nos pareció parsimonioso, con variables comunes en todos los modelos estudiados.

modPreliminar <- glm(day30 ~ age + hrt + ttr + pmi + sho,
                     data=mydata, family = binomial())
summary(modPreliminar)

Call:
glm(formula = day30 ~ age + hrt + ttr + pmi + sho, family = binomial(), 
    data = mydata)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.205  -0.378  -0.223  -0.137   3.139  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -10.3489     1.3697   -7.56  4.2e-14 ***
age           0.0997     0.0186    5.36  8.4e-08 ***
hrtSi         1.0128     0.3133    3.23   0.0012 ** 
ttrSi         0.7375     0.3241    2.28   0.0229 *  
pmiSi         0.6689     0.3280    2.04   0.0415 *  
shoSi         0.9155     0.5383    1.70   0.0890 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 382.78  on 784  degrees of freedom
Residual deviance: 308.21  on 779  degrees of freedom
AIC: 320.2

Number of Fisher Scoring iterations: 7

#coef(modPreliminar)

Pregunta 3: Términos de interacción

Decida si se incluirán términos de interacción; utilizar pruebas de LR.

Evaluación de interacciones

A continuación se evaluará el efecto de interacciones de orden 2 en el modelo, con las variables que se incluyeron en el modelo preliminar: age, hrt, ttr, pmi y sho.

#Test LR para comparación de dos modelos cualquiera, en este caso modelo sin y con una interacción
library(lmtest)
modPreliminar <- glm(day30 ~ age + hrt + ttr + pmi + sho,
                     data=mydata, family = binomial())
#complex es el modelo con la interacción
complex <- glm(day30 ~ age + hrt + ttr + pmi + sho + sho*hrt,
                     data=mydata, family = binomial())

lrtest(modPreliminar, complex)

Model 1: day30 ~ age + hrt + ttr + pmi + sho
Model 2: day30 ~ age + hrt + ttr + pmi + sho + sho * hrt

L.R. Chisq       d.f.          P 
 0.0433204  1.0000000  0.8351231 

Pregunta 4: Exactitud del modelo

a) Determinar la bondad de ajuste del modelo con la prueba de Hosmer-Lemeshow.

Las hipótesis a testear serían:

H0: Las proporciones de tasas de muerte a 30 días observadas y esperadas son iguales H1: Las proporciones de tasas de muerte a 30 días observadas y esperadas no son iguales

Realicemos el test HL:

R2.hl <- modelChi/modPreliminar$null.deviance
R2.hl

[1] 0.0652088

library(ResourceSelection)
set.seed(43657)
n <- 100
x <- rnorm(n)
xb <- x
pr <- exp(xb)/(1+exp(xb))
y <- 1*(runif(n) < pr)
mod <- glm(y~x, family=binomial)
hl <- hoslem.test(mydata$day30, fitted(modPreliminar), g=10)
hl

    Hosmer and Lemeshow goodness of fit (GOF) test

data:  mydata$day30, fitted(modPreliminar)
X-squared = 785, df = 8, p-value <2e-16

El valor-p mayor a 0.05 no nos permite rechazar la hipótesis nula, por lo que el test HL concluye que el modelo logístico planteado ajusta bien a los datos.

b) Determinar la capacidad de discriminación calculando el área bajo la curva ROC.

#### Capacidad de discriminación
library(Hmisc)
library(rms)
library(reportROC)

modfinal <- lrm(day30 ~ age + hrt + ttr + pmi + sho,
                        data=mydata, x=T, y=T, linear.predictors=F, tol=1e-10)

#Predicciones
mydata$p.modfinal <- predict(modfinal, type="fitted.ind")

#Curva ROC
resultROC <- reportROC(gold=mydata$day30,
                       predictor=mydata$p.modfinal,
                       important="se", plot=TRUE,
                       positive = "l")

Setting levels: control = No, case = Si

resultROC

Este modelo tiene un area bajo la curva (AUC) de 0.81, lo que indica que el modelo pronóstico propuesto es un modelo con una muy buena capacidad de discriminación entre individuos con o sin la muerte a los 30 días después del IAM.

Pregunta 5: Presentación

a) Presentar el modelo final como una fórmula de regresión

El modelo logístico final elegido tiene la siguiente forma como fórmula de regresión:

\[ \log\left(\frac{p(day30)}{1 - p(day30)}\right) = -10.3489 + 0.0997\times age + 1.0128\times I(hrt=1) + 0.7375\times I(ttr=1) + 0.6689\times I(pmi=1) + 0.9155\times I(sho=1) \]

En términos de probabilidad, sería \[ P(Y) = \frac{1}{1+e^{- (-10.3489 + 0.0997\times age + 1.0128\times I(hrt=1) + 0.7375\times I(ttr=1) + 0.6689\times I(pmi=1) + 0.9155\times I(sho=1) )}} \] otra forma \[ P(Y) = \frac{e^{ -10.3489 + 0.0997\times age + 1.0128\times I(hrt=1) + 0.7375\times I(ttr=1) + 0.6689\times I(pmi=1) + 0.9155\times I(sho=1)} } {1+e^{-10.3489 + 0.0997\times age + 1.0128\times I(hrt=1) + 0.7375\times I(ttr=1) + 0.6689\times I(pmi=1) + 0.9155\times I(sho=1) }} \]

b) Presente el modelo final como una regla de predicción clínica en un formato que crea que los médicos apreciarán. Sugerencias: un gráfico de puntuación con puntuaciones basadas en 10 veces los coeficientes de regresión redondeados, una tabla, etc.

Se presenta una calculadora en excel (ver formato adjunto), en el que se calcula con la fórmula matemática descrita la probabilidad de presentar el evento. El clínico podra ingresar las variables dicotomicas como “si” o “no”, dependiendo de su paciente y debe además incluir la edad del paciente en años. Finalmente se le arrojará en el resultado final la probabilidad que su paciente desarrolle el desenlace muerte a los 30 días después del IAM.