RLM: TP #1

Carga de librerías.

library(readxl)

## Warning: package 'readxl' was built under R version 4.4.2

library(data.table)

## Warning: package 'data.table' was built under R version 4.4.3

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.4.3

## Warning: package 'tidyr' was built under R version 4.4.2

## Warning: package 'dplyr' was built under R version 4.4.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between()     masks data.table::between()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ dplyr::first()       masks data.table::first()
## ✖ lubridate::hour()    masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ dplyr::last()        masks data.table::last()
## ✖ lubridate::mday()    masks data.table::mday()
## ✖ lubridate::minute()  masks data.table::minute()
## ✖ lubridate::month()   masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second()  masks data.table::second()
## ✖ purrr::transpose()   masks data.table::transpose()
## ✖ lubridate::wday()    masks data.table::wday()
## ✖ lubridate::week()    masks data.table::week()
## ✖ lubridate::yday()    masks data.table::yday()
## ✖ lubridate::year()    masks data.table::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(broom)
library(vcd)

## Warning: package 'vcd' was built under R version 4.4.3

## Cargando paquete requerido: grid

library(sjPlot)

## Warning: package 'sjPlot' was built under R version 4.4.3

library(epiR)

## Warning: package 'epiR' was built under R version 4.4.2

## Cargando paquete requerido: survival
## Package epiR 2.0.80 is loaded
## Type help(epi.about) for summary information
## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses

library(parameters)

## Warning: package 'parameters' was built under R version 4.4.3

Carga y preparación de datos.

setwd("C:/Users/Usuario/Desktop/iecs/ano_2/regresion_logistica_multiple/actividades_asincronicas_del_18_05_al_5_06")
data <- fread("framdata_ej_1.csv")
data <- data %>% mutate(across(c("sex", "newchd", "male", "muerte"), as.factor))

Ejercicios

Ajuste un modelo logístico entre enfermedad coronaria (newchd) y edad (age).

Modelo logístico con edad como predictor: Se ajusta un modelo de regresión logística para predecir la probabilidad de enfermedad coronaria (newchd) en función de la edad (age). La familia “binomial” indica que se trata de un modelo para variable dependiente dicotómica (0 = no, 1 = sí). Este es el modelo base, sin otras variables.

modelo_1 <- glm(newchd ~ age, family = "binomial", data = data)

options(scipen = 999)

summary(modelo_1)

## 
## Call:
## glm(formula = newchd ~ age, family = "binomial", data = data)
## 
## Coefficients:
##             Estimate Std. Error z value       Pr(>|z|)    
## (Intercept) -4.88431    0.77372  -6.313 0.000000000274 ***
## age          0.06581    0.01446   4.550 0.000005374208 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1330.2  on 1361  degrees of freedom
## AIC: 1334.2
## 
## Number of Fisher Scoring iterations: 4

¿Cómo se interpreta?

0.0658 es el B de la variable age.
El encabezado Log-Odds indica que la tabla muestra la ecuación lineal en escala logit, es decir, antes de aplicar la exponencial para interpretar en términos de odds ratio.
El coeficiente B1 = 0.065810 indica que por cada año adicional de edad, el logaritmo de las odds del evento (enfermedad coronaria) aumenta en 0.0658.
Para transformarlo a una odds ratio (OR), se puede aplicar la exponencial:
- OR = e^0.0658 ≈ 1.07

Esto implica que cada año de edad incrementa las odds de enfermedad coronaria en un 7%.

model_parameters(modelo_1, ci = 0.95, digits = 6)

## Parameter   |  Log-Odds |       SE |                 95% CI |         z |      p
## --------------------------------------------------------------------------------
## (Intercept) | -4.884311 | 0.773723 | [-6.415876, -3.380611] | -6.312742 | < .001
## age         |  0.065810 | 0.014465 | [ 0.037598,  0.094343] |  4.549623 | < .001

## 
## Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
##   computed using a Wald z-distribution approximation.

## 
## The model has a log- or logit-link. Consider using `exponentiate =
##   TRUE` to interpret coefficients as ratios.

Presentación tabular del modelo 1: Esta línea genera una tabla con los resultados del modelo modelo_1 usando tab_model() (de sjPlot). Se muestra el error estándar (show.se), desviación, log-verosimilitud, pero se oculta el intercepto y otros detalles menos relevantes para simplificar la presentación. El título especifica que se trata del modelo con edad como única variable explicativa.

La línea con el código knitr permite que la tabla generada por tab_model() se imprima correctamente en el documento HTML o PDF generado por R Markdown. asis_output() le indica a knitr que la salida debe tratarse como código HTML/R ya procesado.

t1 <- tab_model(modelo_1, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Y = Enfermedad coronaria; X<sub>1</sub> = Edad") 

knitr::asis_output(t1$knitr)

Y = Enfermedad coronaria; X₁ = Edad
	newchd
Predictors	Odds Ratios	std. Error	CI	p
age	1.07	0.02	1.04 – 1.10	<0.001
Observations	1363
Deviance	1330.165
log-Likelihood	-665.082

¿Cómo se interpreta el efecto de la edad en el modelo simple?

Por cada unidad de incremento en la edad (expresada en años), las odds de presentar enfermedad coronaria aumentan en un 7%, en comparación con personas de menor edad. Este resultado es estadísticamente significativo (p < 0.001), y el intervalo de confianza del 95% para el odds ratio (OR) se encuentra entre 1.04 y 1.10, lo que indica que con un 95% de certeza el efecto verdadero está dentro de ese rango y no cruza la unidad.

Agregamos la variable sex (hombre = 1)

modelo_2 <- glm(newchd ~ age + sex, family = "binomial", data = data)

summary(modelo_2)

## 
## Call:
## glm(formula = newchd ~ age + sex, family = "binomial", data = data)
## 
## Coefficients:
##             Estimate Std. Error z value        Pr(>|z|)    
## (Intercept) -5.30821    0.78626  -6.751 0.0000000000147 ***
## age          0.06672    0.01458   4.575 0.0000047515111 ***
## sex1         0.71613    0.14052   5.096 0.0000003461213 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1303.5  on 1360  degrees of freedom
## AIC: 1309.5
## 
## Number of Fisher Scoring iterations: 4

model_parameters(modelo_2, ci = 0.95, digits = 6)

## Parameter   |  Log-Odds |       SE |                 95% CI |         z |      p
## --------------------------------------------------------------------------------
## (Intercept) | -5.308213 | 0.786263 | [-6.865548, -3.781093] | -6.751194 | < .001
## age         |  0.066718 | 0.014582 | [ 0.038282,  0.095485] |  4.575470 | < .001
## sex [1]     |  0.716135 | 0.140517 | [ 0.442399,  0.993692] |  5.096428 | < .001

## 
## Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
##   computed using a Wald z-distribution approximation.

t2 <- tab_model(modelo_2, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = F, title = "Y = Enfermedad coronaria; X<sub>1</sub> = Edad, X<sub>2</sub> = Sexo") 

knitr::asis_output(t2$knitr)

Y = Enfermedad coronaria; X₁ = Edad, X₂ = Sexo
	newchd
Predictors	Odds Ratios	std. Error	CI	p
age	1.07	0.02	1.04 – 1.10	<0.001
sex1	2.05	0.29	1.56 – 2.70	<0.001
Observations	1363
Deviance	1303.531
log-Likelihood	-651.766

¿Cómo se interpreta el efecto de la variable sex en el modelo múltiple?

El odds de presentar enfermedad coronaria entre los hombres (sex = 1) es aproximadamente el doble que entre las mujeres (categoría de referencia) ajustado o controlado o independientemente del resto de las variables. Este hallazgo es estadísticamente significativo (p < 0.001), y el intervalo de confianza del 95% para el odds ratio (OR) se encuentra entre 1.56 y 2.70, lo que indica que con un 95% de certeza el valor real del OR se encuentra dentro de ese rango y no cruza la unidad.

¿La variable sex agrega valor predictivo al modelo?, ¿cómo lo explica?

Para evaluar si la variable sex agrega valor predictivo al modelo, una manera correcta y formal es comparar el modelo que contiene solo age (modelo_1) con el modelo que incluye age + sex (modelo_2), mediante un test de verosimilitud (LRT). Esto permite contrastar estadísticamente si la inclusión de sex mejora el ajuste del modelo.

anova(modelo_1, modelo_2, test = "LRT")

## Analysis of Deviance Table
## 
## Model 1: newchd ~ age
## Model 2: newchd ~ age + sex
##   Resid. Df Resid. Dev Df Deviance    Pr(>Chi)    
## 1      1361     1330.2                            
## 2      1360     1303.5  1   26.633 0.000000246 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se comparó el modelo con solo age frente al modelo que agrega sex mediante un test de verosimilitud (LRT). El resultado fue estadísticamente significativo (p < 0.001), indicando que sex mejora el ajuste del modelo, es decir, agrega valor predictivo. Esto se condice con el hallazgo de un OR significativo para sex en el modelo multivariable.

Entonces, concluimos que la variable sex agrega valor predictivo al modelo. Esto se evidencia por:

- Un p-valor significativo en el test de Wald (p < 0.001), lo que indica que sex está asociada a la probabilidad de enfermedad coronaria.

- El intervalo de confianza del OR (1.56 a 2.70) no incluye la unidad, reforzando la significancia clínica.

- La deviance residual del modelo múltiple (1303.5) es menor que la del modelo simple (1351.2).

-Al realizar una comparación formal entre el modelo simple (age) y el modelo con age + sex mediante un test de verosimilitud (LRT), se obtuvo un p-valor < 0.001, confirmando que el modelo con sex se ajusta significativamente mejor a los datos.

¿Cómo se interpreta el efecto de la variable age en el modelo múltiple?

Ajustado por sexo, el coeficiente para age fue 0.0667, lo que implica que por cada año adicional de edad (entre 45 y 62 años, rango de la base de datos), el log odds de enfermedad coronaria aumenta 0.0667 unidades, y el odds ratio (OR) asociado es:

exp(0.0667)≈1.07

Es decir, el riesgo de enfermedad coronaria aumenta un 7% por cada año adicional.

Si se desea interpretar el efecto de la edad en intervalos de 10 años (o qué podría ocurrir con una persona con 10 años más de edad), se puede calcular:

exp(0.667)≈1.95

lo que indica que por cada 10 años adicionales, los odds de enfermedad coronaria se duplican aproximadamente.

¿Existe efecto confundidor de sex sobre la relación entre edad y enfermedad coronaria?

Si edad se ingresa primero, y el análisis se plantea como “¿qué efecto tiene la edad sobre la probabilidad de enfermedad?”, entonces es la exposición principal. Lo que se evalúa es si esa relación está confundida o modificada por sex. Al comparar el modelo simple (edad) con el modelo múltiple (edad + sexo), el OR para la variable edad no se modifica (1.07 en ambos modelos). Esto indica que la variable sexo no actúa como confundidor en la asociación entre edad y enfermedad coronaria. El efecto de la edad es estable independientemente del ajuste por sexo, por lo cual no se observa confusión.

¿Existe interacción entre edad y sexo?

modelo_3 <- glm(newchd ~ age + sex + age:sex, family = "binomial", data = data)

summary(modelo_3)

## 
## Call:
## glm(formula = newchd ~ age + sex + age:sex, family = "binomial", 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value     Pr(>|z|)    
## (Intercept) -6.99749    1.25445  -5.578 0.0000000243 ***
## age          0.09820    0.02319   4.234 0.0000229071 ***
## sex1         3.54459    1.60431   2.209       0.0271 *  
## age:sex1    -0.05297    0.02987  -1.773       0.0762 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1300.4  on 1359  degrees of freedom
## AIC: 1308.4
## 
## Number of Fisher Scoring iterations: 4

model_parameters(modelo_3, ci = 0.95, digits = 6)

## Parameter     |  Log-Odds |       SE |                 95% CI |         z |      p
## ----------------------------------------------------------------------------------
## (Intercept)   | -6.997492 | 1.254449 | [-9.509017, -4.581859] | -5.578139 | < .001
## age           |  0.098198 | 0.023190 | [ 0.053266,  0.144350] |  4.234490 | < .001
## sex [1]       |  3.544593 | 1.604311 | [ 0.419079,  6.715232] |  2.209418 | 0.027 
## age × sex [1] | -0.052967 | 0.029871 | [-0.111867,  0.005345] | -1.773211 | 0.076

## 
## Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
##   computed using a Wald z-distribution approximation.

t3<-tab_model(modelo_3, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = T, title = "Interacción") 

knitr::asis_output(t3$knitr)

Interacción
	newchd
Predictors	Odds Ratios	std. Error	p
age	1.10 (1.05 – 1.16)	0.03	<0.001
age:sex1	0.95 (0.89 – 1.01)	0.03	0.076
sex1	34.63 (1.52 – 824.87)	55.55	0.027
Observations	1363
Deviance	1300.364
log-Likelihood	-650.182

En el modelo con término de interacción edad:sexo, el valor p asociado a la interacción fue 0.076. Si bien este valor no alcanza significación estadística al umbral clásico (p < 0.05), se encuentra cercano, y el OR de interacción (0.95) sugiere una posible atenuación del efecto de la edad en varones. Desde el punto de vista epidemiológico y biológico, el sexo podría modificar la asociación entre edad y enfermedad coronaria. Sin embargo, en este análisis no se encontró evidencia estadísticamente significativa de interacción.

Modelo sin interacción

data %>% ggplot()+geom_point(aes(age,predict(modelo_2, data), color= sex))+ylab("logodds EC")

En el modelo sin interacción (modelo_2), estás viendo dos rectas paralelas (una para cada sexo), con pendientes similares (porque no hay interacción) pero diferente intercepto. El log-odds (logaritmo del odds) es una escala lineal del riesgo. Si los hombres tienen log-odds más altos (menos negativos) que las mujeres para la misma edad, eso significa que los hombres tienen mayor probabilidad de enfermedad coronaria que las mujeres a cualquier edad. Como las rectas son paralelas, el efecto de la edad es el mismo en ambos sexos (misma pendiente). Pero el nivel de riesgo basal es diferente (interceptos distintos): los hombres parten de un riesgo más alto a cualquier edad. Esto refuerza lo que se ve en el modelo_2: sex1 es un predictor independiente significativo (OR ≈ 2.05; p < 0.001) pero no modifica el efecto de la edad → no hay interacción estadística significativa.

Modelo con término de interacción

data %>% ggplot()+geom_point(aes(age,predict(modelo_3, data), color= sex))+ylab("logodds EC")

En este modelo (modelo_3), al incluir edad:sex1, las dos líneas ya no son paralelas: La pendiente para hombres y la de mujeres se calculan por separado, porque hay un término de interacción. Si bien las líneas no se cruzan, sus pendientes difieren levemente. Log-odds más altos en hombres a todas las edades → se mantiene la diferencia de riesgo basal. La pendiente es levemente menor en hombres (la línea sube más lento), reflejando que el incremento del riesgo por edad sería un poco menor en varones que en mujeres. Como no se cruzan, el efecto de la edad no se invierte entre sexos. El término de interacción (edad:sexo1) tiene un OR de 0.95, lo que indica una leve atenuación del efecto de la edad en hombres, pero no significativa (p = 0.076). Este gráfico refuerza que el sexo masculino tiene mayor riesgo absoluto, pero su riesgo crece a un ritmo ligeramente menor con la edad comparado con mujeres. Esto podría tener valor clínico o biológico, aunque el efecto no sea significativo estadísticamente. Aunque las líneas parecen casi paralelas, hay una ligera diferencia de pendiente. En mujeres: la pendiente de la recta representa el efecto completo de la variable edad. En hombres: la pendiente es edad + edad:sexo1, es decir, un poco menos pronunciada. No obstante, debido a que el valor de p está cerca de 0.1 y las curvas de riesgo según edad difieren visualmente entre hombres y mujeres, se puede justificar mantener el término de interacción en modelos posteriores o reportar análisis estratificados.

Existe evidencia epidemiológica que sugiere algunas diferencias importantes entre hombres y mujeres con respecto al riesgo de EC. Una de ellas podría ser el efecto que la hipertensión sistólica tiene sobre la relacion entre sexo y EC.

modelo_4 <- glm(newchd ~ sbp, family = "binomial", data = data)

summary(modelo_4)

## 
## Call:
## glm(formula = newchd ~ sbp, family = "binomial", data = data)
## 
## Coefficients:
##              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept) -3.769129   0.362313 -10.403 < 0.0000000000000002 ***
## sbp          0.015623   0.002309   6.767      0.0000000000131 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1305.0  on 1361  degrees of freedom
## AIC: 1309
## 
## Number of Fisher Scoring iterations: 4

model_parameters(modelo_4, ci = 0.95, digits = 6)

## Parameter   |  Log-Odds |       SE |                 95% CI |          z |      p
## ---------------------------------------------------------------------------------
## (Intercept) | -3.769129 | 0.362313 | [-4.488853, -3.066568] | -10.402961 | < .001
## sbp         |  0.015623 | 0.002309 | [ 0.011127,  0.020193] |   6.767211 | < .001

## 
## Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
##   computed using a Wald z-distribution approximation.

t4 <- tab_model(modelo_4, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = T, title = "Y = Enfermedad coronaria; X<sub>3</sub> = Hipertensión sistólica") 

knitr::asis_output(t4$knitr)

Y = Enfermedad coronaria; X₃ = Hipertensión sistólica
	newchd
Predictors	Odds Ratios	std. Error	p
sbp	1.02 (1.01 – 1.02)	0.00	<0.001
Observations	1363
Deviance	1305.048
log-Likelihood	-652.524

Por cada 1 mmHg de incremento en la presión arterial sistólica, las odds (razón de chances) de presentar enfermedad coronaria aumentan en un 2% en comparación con una persona con 1 mmHg menos. Esta asociación es estadísticamente significativa (p < 0.001) y, con un 95% de certeza, el OR verdadero se encuentra entre 1,01 y 1,02.

Transforme sbp en una variable dicotómica “sbp_cat” que tome valor 0 si sbp < 141 y 1 si sbp >140. Ajuste un modelo entre newchd y sbp_cat y otro agregando sex.

data$sbp_cat <- ifelse(data$sbp > 140, 1,0)

modelo_5 <- glm(newchd ~ sbp_cat, family = "binomial", data = data)

summary(modelo_5)

## 
## Call:
## glm(formula = newchd ~ sbp_cat, family = "binomial", data = data)
## 
## Coefficients:
##             Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)  -1.8770     0.1145 -16.398 < 0.0000000000000002 ***
## sbp_cat       0.8162     0.1435   5.689         0.0000000128 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1317.3  on 1361  degrees of freedom
## AIC: 1321.3
## 
## Number of Fisher Scoring iterations: 4

model_parameters(modelo_5, ci = 0.95, digits = 6)

## Parameter   |  Log-Odds |       SE |                 95% CI |          z |      p
## ---------------------------------------------------------------------------------
## (Intercept) | -1.877033 | 0.114467 | [-2.107840, -1.658580] | -16.398002 | < .001
## sbp cat     |  0.816161 | 0.143462 | [ 0.537786,  1.100717] |   5.689040 | < .001

## 
## Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
##   computed using a Wald z-distribution approximation.

t5 <- tab_model(modelo_5, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = T, title = "Y = Enfermedad coronaria; X<sub>4</sub> = Hipertensión sistólica [CAT]") 

knitr::asis_output(t5$knitr)

Y = Enfermedad coronaria; X₄ = Hipertensión sistólica [CAT]
	newchd
Predictors	Odds Ratios	std. Error	p
sbp_cat	2.26 (1.71 – 3.01)	0.32	<0.001
Observations	1363
Deviance	1317.253
log-Likelihood	-658.626

El odds del evento (enfermedad coronaria) entre hipertensos es 2.2 veces mayor que entre no hipertensos, siendo esta asociación estadísticamente significativa (p < 0.001) y, con 95% de certeza, el verdadero valor del OR está ubicado entre 1,71 y 3,01.

modelo_6 <- glm(newchd ~ sbp_cat + sex, family = "binomial", data = data)

summary(modelo_6)

## 
## Call:
## glm(formula = newchd ~ sbp_cat + sex, family = "binomial", data = data)
## 
## Coefficients:
##             Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)  -2.3641     0.1498 -15.779 < 0.0000000000000002 ***
## sbp_cat       0.9171     0.1466   6.254         0.0000000004 ***
## sex1          0.8176     0.1429   5.720         0.0000000106 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1351.2  on 1362  degrees of freedom
## Residual deviance: 1283.5  on 1360  degrees of freedom
## AIC: 1289.5
## 
## Number of Fisher Scoring iterations: 4

model_parameters(modelo_6, ci = 0.95, digits = 6)

## Parameter   |  Log-Odds |       SE |                 95% CI |          z |      p
## ---------------------------------------------------------------------------------
## (Intercept) | -2.364110 | 0.149823 | [-2.665479, -2.077781] | -15.779404 | < .001
## sbp cat     |  0.917058 | 0.146639 | [ 0.632682,  1.208051] |   6.253834 | < .001
## sex [1]     |  0.817583 | 0.142924 | [ 0.539394,  1.100113] |   5.720420 | < .001

## 
## Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
##   computed using a Wald z-distribution approximation.

t6 <- tab_model(modelo_5,modelo_6, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = T, title = "Y = Enfermedad coronaria; X<sub>2</sub> = Sexo, X<sub>4</sub> = Hipertensión sistólica [CAT]") 

knitr::asis_output(t6$knitr)

Y = Enfermedad coronaria; X₂ = Sexo, X₄ = Hipertensión sistólica [CAT]
	newchd			newchd
Predictors	Odds Ratios	std. Error	p	Odds Ratios	std. Error	p
sbp_cat	2.26 (1.71 – 3.01)	0.32	<0.001	2.50 (1.88 – 3.35)	0.37	<0.001
sex1				2.27 (1.71 – 3.00)	0.32	<0.001
Observations	1363			1363
Deviance	1317.253			1283.467
log-Likelihood	-658.626			-641.733

El ajuste por sexo aumentó el OR de SBP categorizado. El cambio es del 10.6%, y ambos intervalos de confianza son estadísticamente significativos. Además, el OR para sex1 en modelo_6 es 2.27, mostrando que el sexo sigue siendo un predictor independiente.

¿Mejora la predicción del modelo múltiple?

anova(modelo_6, test = "LRT")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: newchd
## 
## Terms added sequentially (first to last)
## 
## 
##         Df Deviance Resid. Df Resid. Dev       Pr(>Chi)    
## NULL                     1362     1351.2                   
## sbp_cat  1   33.993      1361     1317.2 0.000000005532 ***
## sex      1   33.786      1360     1283.5 0.000000006152 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El análisis de devianza mediante LRT mostró que tanto la hipertensión sistólica categorizada (sbp_cat) como el sexo (sex1) contribuyen significativamente al modelo (p < 0.001 para ambas variables). La inclusión secuencial de cada predictor redujo de manera significativa la devianza, lo que indica que el modelo múltiple tiene mejor capacidad predictiva en comparación con el modelo nulo. Por lo tanto, ambos predictores aportan información independiente y estadísticamente significativa sobre la probabilidad de presentar enfermedad coronaria. El modelo nulo (null model) es el modelo más simple posible, que no incluye ninguna variable predictora. Solo estima una constante (intercepto), es decir, predice la probabilidad de enfermedad coronaria (newchd) usando únicamente la proporción general del evento en la muestra. Este modelo asume que todos los individuos tienen la misma probabilidad de presentar enfermedad coronaria, sin importar sus características (edad, sexo, hipertensión, etc.). Sirve como referencia base para comparar si un modelo con predictores (como modelo_5, modelo_6, etc.) mejora la predicción. Si el modelo con predictores reduce significativamente la devianza con respecto al nulo (como este test anova(modelo_6, test = "LRT")), entonces se considera que aporta valor.

Se comparó el modelo que incluye solo sex frente al modelo que agrega también sbp_cat mediante un test de verosimilitud (LRT). El resultado fue estadísticamente significativo (p < 0.001), indicando que sbp_cat mejora el ajuste del modelo, es decir, agrega valor predictivo.

Esto se evidencia por:

- Un p-valor significativo en el test de Wald para sbp_cat (p < 0.001), lo que indica que esta variable está asociada de manera independiente a la probabilidad de enfermedad coronaria.

- El intervalo de confianza del OR para sbp_cat (1.83 a 3.35) no incluye la unidad, reforzando su significación clínica.

- La deviance residual del modelo que incluye sbp_cat (1283.5) es menor que la del modelo sin sbp_cat (1317.3).

- Al realizar la comparación formal mediante el test de verosimilitud (LRT), se obtuvo un p-valor < 0.001, confirmando que el modelo con sbp_cat se ajusta significativamente mejor a los datos.

Interpretación de columnas del output:

Columna	Descripción
Df	Grados de libertad del modelo ajustado hasta ese punto.
Deviance	Medida de falta de ajuste (menor = mejor ajuste). Se calcula como: -2 × log-verosimilitud.
Resid. Df	Grados de libertad residuales (N - número de parámetros).
Resid. Dev	Deviance residual del modelo ajustado hasta ese punto.
Pr(>Chi)	Valor p del LRT. Compara el modelo actual con el anterior. Si es < 0.05, indica que la nueva variable mejora significativamente el modelo.

Evalúe si sex es un confundidor, un modificador de efecto o ambos de la relación entre sbp_cat y newchd. Justifique su respuesta usando tablas de 2x2 y regresión logística.

Exposed: sbp_cat; Outcome: newchd; Strata: sex.

a <- epi.2by2(table(exposure = data$sbp_cat,outcome = data$newchd, strata= data$sex), method = "cohort.count", digits = 3, conf.level = 0.95, outcome = "as.columns", interpret = F)
a

##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +          575           88        663  86.727 (83.905 to 89.217)
## Exposed -          520          180        700  74.286 (70.878 to 77.487)
## Total             1095          268       1363  80.337 (78.126 to 82.417)
## 
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio (crude)                         1.167 (1.107, 1.231)
## Inc risk ratio (M-H)                           1.185 (1.125, 1.249)
## Inc risk ratio (crude:M-H)                     0.985
## Inc Odds ratio (crude)                         2.262 (1.707, 2.996)
## Inc Odds ratio (M-H)                           2.526 (1.890, 3.375)
## Inc Odds ratio (crude:M-H)                     0.895
## Attrib risk in the exposed (crude) *           12.441 (8.300, 16.583)
## Attrib risk in the exposed (M-H) *             13.651 (9.246, 18.056)
## Attrib risk (crude:M-H)                        0.911
## -------------------------------------------------------------------
##  M-H test of homogeneity of IRRs: chi2(1) = 0.103 Pr>chi2 = 0.749
##  M-H test of homogeneity of ORs: chi2(1) = 3.400 Pr>chi2 = 0.065
##  Test that M-H adjusted OR = 1:  chi2(1) = 40.471 Pr>chi2 = <0.001
##  Wald confidence limits
##  M-H: Mantel-Haenszel; CI: confidence interval
##  * Outcomes per 100 population units

b <- a$massoc.detail$OR.strata.wald
row.names(b)<-c("sex_0", "sex_1")
b

##            est    lower    upper
## sex_0 3.655952 2.208210 6.052860
## sex_1 2.019627 1.408432 2.896052

Aunque no hay evidencia estadística de modificación de efecto por sexo (p = 0.065), el cambio mayor al 10% en el OR tras el ajuste sugiere que sexo actúa como una variable confusora, ya que su inclusión en el modelo altera sustancialmente la estimación del efecto de sbp_cat sobre newchd.

modelo_7 <- glm(newchd ~ sbp_cat + sex + sbp_cat:sex, family = "binomial", data = data)

t7 <- tab_model(modelo_7, show.se = T, show.dev = T, show.loglik = T, show.intercept = F,show.stat = F, show.reflvl = T, digits = 2, digits.p = 3, show.r2 = F,collapse.ci = T, title = "Y = Enfermedad coronaria; X2 = Sexo, X4 = Hipertensión sistólica [CAT]; Interacción") 

knitr::asis_output(t7$knitr)

Y = Enfermedad coronaria; X2 = Sexo, X4 = Hipertensión sistólica [CAT]; Interacción
	newchd
Predictors	Odds Ratios	std. Error	p
sbp_cat	3.66 (2.25 – 6.19)	0.94	<0.001
sbp_cat:sex1	0.55 (0.29 – 1.02)	0.17	0.061
sex1	3.38 (2.05 – 5.80)	0.89	<0.001
Observations	1363
Deviance	1279.818
log-Likelihood	-639.909

La diferencia entre los resultados de modelo_6 y los que aparecen ahora se debe al orden en que se ingresan las variables en el modelo, en especial cuando usás anova(modelo_6, test = "LRT") o cuando hacés comparaciones de mejora secuencial.

¿Por qué cambian los OR?

En términos simples:

Cuando hacés modelo_6 <- glm(newchd ~ sbp_cat + sex1, ...), el efecto de cada variable se estima ajustando por la otra.
Pero cuando usás tab_model(modelo_5, modelo_6, ...), o bien comparás modelos con anova(), el orden de entrada importa si estás estimando efectos marginales secuenciales, como en el test LRT (Likelihood Ratio Test).

En el modelo más reciente con la interacción (modelo_7):

sbp_cat ya no representa el OR crudo ni ajustado general, sino el OR cuando sex = 0 (o sea, en mujeres).
Y el coeficiente sex representa el efecto del sexo cuando sbp_cat = 0 (sin hipertensión).
Mientras que sbp_cat:sex1 mide cuánto varía el efecto de la hipertensión en hombres respecto de mujeres.
Los OR cambian porque en modelo_6 los efectos son ajustados, mientras que en modelo_7 con la interacción los efectos son estratificados por nivel de la otra variable.

¿Qué significa cada término del modelo?

sbp_cat

Este OR refiere al efecto de la hipertensión en mujeres (la categoría de referencia de sex).
Es decir, el OR de sbp_cat cuando sex = 0.
Entonces: El odds del evento (enfermedad coronaria), entre mujeres con hipertensión, es 3.66 veces mayor comparado con mujeres sin hipertensión.

sex

Este OR refiere al efecto del sexo masculino cuando no hay hipertensión (sbp_cat = 0).
Es decir, el OR de sex cuando sbp_cat = 0.
Entonces: El odds del evento, entre hombres sin hipertensión, es 3.38 veces mayor comparado con mujeres sin hipertensión.

sbp_cat:sex

Este es el término de interacción: mide cuánto cambia el efecto de sbp_cat en hombres respecto al efecto en mujeres.
En términos prácticos: si el OR de sbp_cat en mujeres es 3.66 y el OR de la interacción es 0.55, entonces el OR de sbp_cat en hombres es:

OR sbp_cat en hombres = 3.66 × 0.55 = 2.01

El odds ratio del efecto de la hipertensión sobre la enfermedad coronaria en hombres es 0.55 veces el efecto de la hipertensión en mujeres. O, si querés decirlo en relación al evento: El efecto de la hipertensión en hombres es 45% menor (OR = 0.55) que en mujeres.

COMENTARIOS:

El odds del evento entre hipertensos y no hipertensos, ¿es igual entre hombres y mujeres? No, ¿cuánto? 0.5 (término de interacción), hombres * término de interacción y mujeres x 1.*

El efecto de sbp_cat sobre el riesgo de enfermedad coronaria no es homogéneo entre hombres y mujeres, lo que sugiere una posible interacción. El término de interacción (sbp_cat:sex1) mostró un p-valor de 0.063, cercano al umbral convencional de significación estadística.

Aunque no alcanza el criterio clásico de p < 0.05, la diferencia en los OR —3.61 en mujeres versus 2.01 en hombres— sugiere un efecto clínicamente relevante, compatible con una modificación del efecto por sexo.

Desde un enfoque epidemiológico, esta interacción puede interpretarse como un posible efecto de subgrupo, en el cual el impacto de la hipertensión sistólica categorizada (sbp_cat) sobre la enfermedad coronaria es más marcado en mujeres.

Por último, aunque el sexo actúe también como un confusor en el modelo, su efecto de confusión es limitado y queda en segundo plano ante la presencia de una posible interacción clínicamente significativa.

RLM: TP #1

Juan Pablo Vinicki

2025-05-18