Regresión logística: Supuestos y precisión del modelo
David Alexander Tejada Peña
2024-09-28
Dr. David Tejada
Dr. Cesar Gavidia
Especialistas en Epidemiología e Investigación en Salud
Preparación del entorno de trabajo
Cargar los paquetes necesarios
## Warning: package 'readxl' was built under R version 4.3.3## Warning: package 'tseries' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Warning: package 'ggplot2' was built under R version 4.3.3## Warning: package 'corrplot' was built under R version 4.3.3## corrplot 0.94 loaded## Warning: package 'broom' was built under R version 4.3.3## Warning: package 'car' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeConocer y establecer el directorio de trabajo
## [1] "C:/Users/Dr. David Tejada/OneDrive/Escritorio/Clases a R2"Cargar la base de datos
Conocer la estructura de la base
## tibble [174 × 59] (S3: tbl_df/tbl/data.frame)
## $ cor : num [1:174] 6 9 10 11 12 15 20 21 25 26 ...
## $ inic : chr [1:174] "DAVM" "MDC" "MAD" "LAS" ...
## $ erc : chr [1:174] "Si" "No" "No" "No" ...
## $ estad_erc : chr [1:174] "est5" "NA" "NA" "NA" ...
## $ ano_diag : chr [1:174] "2021" "NA" "NA" "NA" ...
## $ munic : chr [1:174] "Tejutla" "Dulce Nombre de María" "Dulce Nombre de María" "Santa Rita" ...
## $ area_ur : chr [1:174] "rur" "urb" "urb" "rur" ...
## $ edad : num [1:174] 52 57 80 44 63 42 60 79 67 59 ...
## $ sex : chr [1:174] "mas" "fem" "fem" "mas" ...
## $ educ : chr [1:174] "Si" "Si" "Si" "Si" ...
## $ niv_educ : chr [1:174] "bachi" "basi1_6" "basi1_6" "basi1_6" ...
## $ est_civ : chr [1:174] "casad" "viud" "viud" "vivi_parej" ...
## $ act_lab : chr [1:174] "emplea" "amacas" "desemp" "agric" ...
## $ act_lab_otr : chr [1:174] "NA" "NA" "NA" "NA" ...
## $ agric : chr [1:174] "No" "No" "No" "Si" ...
## $ tabac : chr [1:174] "Si" "No" "No" "No" ...
## $ alcoh : chr [1:174] "No" "No" "No" "No" ...
## $ sal : chr [1:174] "No" "No" "Si" "No" ...
## $ agua_adec : chr [1:174] "No" "No" "No" "Si" ...
## $ No_agua_adec : chr [1:174] "Si" "Si" "Si" "No" ...
## $ frut_verd : chr [1:174] "Si" "Si" "Si" "No" ...
## $ No_frut_verd : chr [1:174] "No" "No" "No" "Si" ...
## $ ejerc : chr [1:174] "Si" "Si" "Si" "Si" ...
## $ No_ejerc : chr [1:174] "No" "No" "No" "No" ...
## $ act_temper : chr [1:174] "Si" "No" "No" "Si" ...
## $ bajpes_prem : chr [1:174] "No" "No" "No" "No" ...
## $ dm : chr [1:174] "Si" "No" "Si" "No" ...
## $ dm_tx : chr [1:174] "Si" "NA" "Si" "NA" ...
## $ chol : chr [1:174] "Si" "Si" "No" "No" ...
## $ chol_tx : chr [1:174] "Si" "Si" "NA" "NA" ...
## $ hta : chr [1:174] "Si" "Si" "Si" "Si" ...
## $ hta_tx : chr [1:174] "Si" "Si" "Si" "Si" ...
## $ enf_autoin : chr [1:174] "No" "No" "No" "No" ...
## $ obesi : chr [1:174] "Si" "Si" "No" "No" ...
## $ ivu_rec : chr [1:174] "Si" "No" "No" "No" ...
## $ lit_ren : chr [1:174] "No" "No" "No" "No" ...
## $ af_dm : chr [1:174] "Si" "Si" "Si" "No" ...
## $ af_hta : chr [1:174] "Si" "Si" "No" "Si" ...
## $ af_acv : chr [1:174] "No" "Si" "Si" "No" ...
## $ af_erc : chr [1:174] "No" "No" "No" "No" ...
## $ af_iam : chr [1:174] "No" "No" "No" "No" ...
## $ af_dislip : chr [1:174] "Si" "Si" "Si" "No" ...
## $ af_otra : chr [1:174] "Si" "No" "No" "No" ...
## $ af_otra_especif : chr [1:174] "Cáncer cérvix" "NA" "NA" "NA" ...
## $ uso_plagui : chr [1:174] "No" "No" "Si" "Si" ...
## $ tip_activiplagui: chr [1:174] "NA" "NA" "otr_act" "mas_act" ...
## $ fum_nfum : chr [1:174] "NA" "NA" "0" "1" ...
## $ epp_plagui : chr [1:174] "NA" "NA" "No" "Si" ...
## $ aines : chr [1:174] "No" "No" "No" "No" ...
## $ per_abdo : num [1:174] 103 95 95 82 95 94 88 75 95 115 ...
## $ per_abdories : chr [1:174] "Si" "Si" "Si" "No" ...
## $ estatur_cm : num [1:174] 162 153 148 167 151 163 168 141 165 174 ...
## $ estat_mt : num [1:174] 1.62 1.53 1.48 1.67 1.51 1.63 1.68 1.41 1.65 1.74 ...
## $ peso : num [1:174] 87 70 62 75 70.5 74 79.8 57 63.5 98 ...
## $ IMC : num [1:174] 33.2 29.9 28.3 26.9 30.9 ...
## $ est_nutr : chr [1:174] "obesidad" "obesidad" "sobrepeso" "sobrepeso" ...
## $ obesib : chr [1:174] "Si" "Si" "No" "No" ...
## $ sobrepes : chr [1:174] "No" "No" "Si" "Si" ...
## $ sobpes_obes : chr [1:174] "Si" "Si" "Si" "Si" ...Transformación de las variables a factor
Verificar la estructura de la base a factor
## tibble [174 × 59] (S3: tbl_df/tbl/data.frame)
## $ cor : num [1:174] 6 9 10 11 12 15 20 21 25 26 ...
## $ inic : Factor w/ 174 levels "AAM","AAMV","ACFZ",..: 22 98 86 77 117 11 51 149 41 87 ...
## $ erc : Factor w/ 2 levels "No","Si": 2 1 1 1 2 1 2 2 2 1 ...
## $ estad_erc : Factor w/ 7 levels "est1","est2",..: 6 7 7 7 2 7 2 4 2 7 ...
## $ ano_diag : Factor w/ 12 levels "2010","2012",..: 9 12 12 12 9 12 10 10 10 12 ...
## $ munic : Factor w/ 16 levels "Arcatao","Azacualpa",..: 16 5 5 15 6 14 8 8 8 8 ...
## $ area_ur : Factor w/ 2 levels "rur","urb": 1 2 2 1 1 1 2 1 2 2 ...
## $ edad : num [1:174] 52 57 80 44 63 42 60 79 67 59 ...
## $ sex : Factor w/ 2 levels "fem","mas": 2 1 1 2 1 1 2 1 2 2 ...
## $ educ : Factor w/ 2 levels "No","Si": 2 2 2 2 2 1 2 1 2 2 ...
## $ niv_educ : Factor w/ 6 levels "bachi","basi1_6",..: 1 2 2 2 2 4 2 4 2 2 ...
## $ est_civ : Factor w/ 5 levels "casad","sep_divor",..: 1 4 4 5 3 3 2 1 1 2 ...
## $ act_lab : Factor w/ 7 levels "agric","amacas",..: 5 2 4 1 2 2 1 2 7 4 ...
## $ act_lab_otr : Factor w/ 10 levels "abogada","agente_seguridad",..: 8 8 8 8 8 8 8 8 7 8 ...
## $ agric : Factor w/ 2 levels "No","Si": 1 1 1 2 1 1 2 1 1 1 ...
## $ tabac : Factor w/ 2 levels "No","Si": 2 1 1 1 1 1 1 1 2 1 ...
## $ alcoh : Factor w/ 2 levels "No","Si": 1 1 1 1 1 1 1 1 2 1 ...
## $ sal : Factor w/ 2 levels "No","Si": 1 1 2 1 1 1 1 1 1 2 ...
## $ agua_adec : Factor w/ 2 levels "No","Si": 1 1 1 2 2 2 2 1 1 2 ...
## $ No_agua_adec : Factor w/ 2 levels "No","Si": 2 2 2 1 1 1 1 2 2 1 ...
## $ frut_verd : Factor w/ 2 levels "No","Si": 2 2 2 1 1 2 1 1 1 1 ...
## $ No_frut_verd : Factor w/ 2 levels "No","Si": 1 1 1 2 2 1 2 2 2 2 ...
## $ ejerc : Factor w/ 2 levels "No","Si": 2 2 2 2 1 2 2 2 1 2 ...
## $ No_ejerc : Factor w/ 2 levels "No","Si": 1 1 1 1 2 1 1 1 2 1 ...
## $ act_temper : Factor w/ 2 levels "No","Si": 2 1 1 2 1 1 2 2 2 2 ...
## $ bajpes_prem : Factor w/ 1 level "No": 1 1 1 1 1 1 1 1 1 1 ...
## $ dm : Factor w/ 2 levels "No","Si": 2 1 2 1 2 1 2 1 2 2 ...
## $ dm_tx : Factor w/ 3 levels "NA","No","Si": 3 1 3 1 3 1 3 1 3 3 ...
## $ chol : Factor w/ 2 levels "No","Si": 2 2 1 1 2 1 2 2 2 2 ...
## $ chol_tx : Factor w/ 3 levels "NA","No","Si": 3 3 1 1 3 1 3 3 3 3 ...
## $ hta : Factor w/ 2 levels "No","Si": 2 2 2 2 2 1 2 2 2 2 ...
## $ hta_tx : Factor w/ 3 levels "NA","No","Si": 3 3 3 3 3 1 3 3 3 3 ...
## $ enf_autoin : Factor w/ 2 levels "No","Si": 1 1 1 1 1 1 1 1 1 1 ...
## $ obesi : Factor w/ 2 levels "No","Si": 2 2 1 1 2 1 1 1 1 2 ...
## $ ivu_rec : Factor w/ 2 levels "No","Si": 2 1 1 1 1 1 1 1 1 1 ...
## $ lit_ren : Factor w/ 2 levels "No","Si": 1 1 1 1 1 1 1 1 1 1 ...
## $ af_dm : Factor w/ 2 levels "No","Si": 2 2 2 1 2 1 2 2 2 2 ...
## $ af_hta : Factor w/ 2 levels "No","Si": 2 2 1 2 2 1 2 1 1 2 ...
## $ af_acv : Factor w/ 2 levels "No","Si": 1 2 2 1 1 2 1 1 2 1 ...
## $ af_erc : Factor w/ 2 levels "No","Si": 1 1 1 1 1 1 2 1 1 2 ...
## $ af_iam : Factor w/ 2 levels "No","Si": 1 1 1 1 1 1 1 1 1 2 ...
## $ af_dislip : Factor w/ 2 levels "No","Si": 2 2 2 1 2 1 2 1 1 2 ...
## $ af_otra : Factor w/ 2 levels "No","Si": 2 1 1 1 1 1 1 1 1 1 ...
## $ af_otra_especif : Factor w/ 10 levels "Asma","Cáncer",..: 3 7 7 7 7 7 7 7 7 7 ...
## $ uso_plagui : Factor w/ 2 levels "No","Si": 1 1 2 2 1 1 2 2 1 2 ...
## $ tip_activiplagui: Factor w/ 4 levels "fum","mas_act",..: 3 3 4 2 3 3 2 4 3 2 ...
## $ fum_nfum : Factor w/ 3 levels "0","1","NA": 3 3 1 2 3 3 2 1 3 2 ...
## $ epp_plagui : Factor w/ 3 levels "NA","No","Si": 1 1 2 3 1 1 2 2 1 2 ...
## $ aines : Factor w/ 2 levels "No","Si": 1 1 1 1 1 1 1 1 1 1 ...
## $ per_abdo : num [1:174] 103 95 95 82 95 94 88 75 95 115 ...
## $ per_abdories : Factor w/ 2 levels "No","Si": 2 2 2 1 2 2 1 1 1 2 ...
## $ estatur_cm : num [1:174] 162 153 148 167 151 163 168 141 165 174 ...
## $ estat_mt : num [1:174] 1.62 1.53 1.48 1.67 1.51 1.63 1.68 1.41 1.65 1.74 ...
## $ peso : num [1:174] 87 70 62 75 70.5 74 79.8 57 63.5 98 ...
## $ IMC : num [1:174] 33.2 29.9 28.3 26.9 30.9 ...
## $ est_nutr : Factor w/ 4 levels "bajo_peso","normal",..: 3 3 4 4 3 4 4 4 2 3 ...
## $ obesib : Factor w/ 2 levels "No","Si": 2 2 1 1 2 1 1 1 1 2 ...
## $ sobrepes : Factor w/ 2 levels "No","Si": 1 1 2 2 1 2 2 2 1 1 ...
## $ sobpes_obes : Factor w/ 2 levels "No","Si": 2 2 2 2 2 2 2 2 1 2 ...Transformación de las variables de factor a numéricas
# Convertir variables categóricas a numéricas
baseerc_numeric <- baseerc %>%
mutate(across(where(is.factor), as.numeric))
# Verificar la estructura después de la conversión
str(baseerc_numeric)## tibble [174 × 59] (S3: tbl_df/tbl/data.frame)
## $ cor : num [1:174] 6 9 10 11 12 15 20 21 25 26 ...
## $ inic : num [1:174] 22 98 86 77 117 11 51 149 41 87 ...
## $ erc : num [1:174] 2 1 1 1 2 1 2 2 2 1 ...
## $ estad_erc : num [1:174] 6 7 7 7 2 7 2 4 2 7 ...
## $ ano_diag : num [1:174] 9 12 12 12 9 12 10 10 10 12 ...
## $ munic : num [1:174] 16 5 5 15 6 14 8 8 8 8 ...
## $ area_ur : num [1:174] 1 2 2 1 1 1 2 1 2 2 ...
## $ edad : num [1:174] 52 57 80 44 63 42 60 79 67 59 ...
## $ sex : num [1:174] 2 1 1 2 1 1 2 1 2 2 ...
## $ educ : num [1:174] 2 2 2 2 2 1 2 1 2 2 ...
## $ niv_educ : num [1:174] 1 2 2 2 2 4 2 4 2 2 ...
## $ est_civ : num [1:174] 1 4 4 5 3 3 2 1 1 2 ...
## $ act_lab : num [1:174] 5 2 4 1 2 2 1 2 7 4 ...
## $ act_lab_otr : num [1:174] 8 8 8 8 8 8 8 8 7 8 ...
## $ agric : num [1:174] 1 1 1 2 1 1 2 1 1 1 ...
## $ tabac : num [1:174] 2 1 1 1 1 1 1 1 2 1 ...
## $ alcoh : num [1:174] 1 1 1 1 1 1 1 1 2 1 ...
## $ sal : num [1:174] 1 1 2 1 1 1 1 1 1 2 ...
## $ agua_adec : num [1:174] 1 1 1 2 2 2 2 1 1 2 ...
## $ No_agua_adec : num [1:174] 2 2 2 1 1 1 1 2 2 1 ...
## $ frut_verd : num [1:174] 2 2 2 1 1 2 1 1 1 1 ...
## $ No_frut_verd : num [1:174] 1 1 1 2 2 1 2 2 2 2 ...
## $ ejerc : num [1:174] 2 2 2 2 1 2 2 2 1 2 ...
## $ No_ejerc : num [1:174] 1 1 1 1 2 1 1 1 2 1 ...
## $ act_temper : num [1:174] 2 1 1 2 1 1 2 2 2 2 ...
## $ bajpes_prem : num [1:174] 1 1 1 1 1 1 1 1 1 1 ...
## $ dm : num [1:174] 2 1 2 1 2 1 2 1 2 2 ...
## $ dm_tx : num [1:174] 3 1 3 1 3 1 3 1 3 3 ...
## $ chol : num [1:174] 2 2 1 1 2 1 2 2 2 2 ...
## $ chol_tx : num [1:174] 3 3 1 1 3 1 3 3 3 3 ...
## $ hta : num [1:174] 2 2 2 2 2 1 2 2 2 2 ...
## $ hta_tx : num [1:174] 3 3 3 3 3 1 3 3 3 3 ...
## $ enf_autoin : num [1:174] 1 1 1 1 1 1 1 1 1 1 ...
## $ obesi : num [1:174] 2 2 1 1 2 1 1 1 1 2 ...
## $ ivu_rec : num [1:174] 2 1 1 1 1 1 1 1 1 1 ...
## $ lit_ren : num [1:174] 1 1 1 1 1 1 1 1 1 1 ...
## $ af_dm : num [1:174] 2 2 2 1 2 1 2 2 2 2 ...
## $ af_hta : num [1:174] 2 2 1 2 2 1 2 1 1 2 ...
## $ af_acv : num [1:174] 1 2 2 1 1 2 1 1 2 1 ...
## $ af_erc : num [1:174] 1 1 1 1 1 1 2 1 1 2 ...
## $ af_iam : num [1:174] 1 1 1 1 1 1 1 1 1 2 ...
## $ af_dislip : num [1:174] 2 2 2 1 2 1 2 1 1 2 ...
## $ af_otra : num [1:174] 2 1 1 1 1 1 1 1 1 1 ...
## $ af_otra_especif : num [1:174] 3 7 7 7 7 7 7 7 7 7 ...
## $ uso_plagui : num [1:174] 1 1 2 2 1 1 2 2 1 2 ...
## $ tip_activiplagui: num [1:174] 3 3 4 2 3 3 2 4 3 2 ...
## $ fum_nfum : num [1:174] 3 3 1 2 3 3 2 1 3 2 ...
## $ epp_plagui : num [1:174] 1 1 2 3 1 1 2 2 1 2 ...
## $ aines : num [1:174] 1 1 1 1 1 1 1 1 1 1 ...
## $ per_abdo : num [1:174] 103 95 95 82 95 94 88 75 95 115 ...
## $ per_abdories : num [1:174] 2 2 2 1 2 2 1 1 1 2 ...
## $ estatur_cm : num [1:174] 162 153 148 167 151 163 168 141 165 174 ...
## $ estat_mt : num [1:174] 1.62 1.53 1.48 1.67 1.51 1.63 1.68 1.41 1.65 1.74 ...
## $ peso : num [1:174] 87 70 62 75 70.5 74 79.8 57 63.5 98 ...
## $ IMC : num [1:174] 33.2 29.9 28.3 26.9 30.9 ...
## $ est_nutr : num [1:174] 3 3 4 4 3 4 4 4 2 3 ...
## $ obesib : num [1:174] 2 2 1 1 2 1 1 1 1 2 ...
## $ sobrepes : num [1:174] 1 1 2 2 1 2 2 2 1 1 ...
## $ sobpes_obes : num [1:174] 2 2 2 2 2 2 2 2 1 2 ...Supuestos de la regresión logística
1. Linealidad
1.1. Gráfico de dispersión con línea de regresión
# Paso 1: Ajustar el modelo logístico
modelo_log <- glm(erc ~ edad, data = baseerc, family = binomial)
# Paso 2: Calcular los residuos
residuos <- residuals(modelo_log, type = "deviance")
# Paso 3: Crear el gráfico de residuos vs predicciones
predicciones <- predict(modelo_log, type = "response")
plot(predicciones, residuos, xlab = "Predicciones", ylab = "Residuos Deviance",
main = "Residuos vs Predicciones", pch = 20, col = "blue")
abline(h = 0, lty = 2, col = "red")
Interpretación: Si los residuos se distribuyen de manera aleatoria alrededor de la línea de regresión (que debería estar cerca de 0), es probable que la suposición de linealidad se cumpla. Si se observa un patrón DISTINTO Y ALEJADOS DE 0, puede ser una señal de no linealidad.
1.2. Gráfico del logit estimado vs edad
# Obtener los valores ajustados del logit
logit_ajustado <- predict(modelo_log, type = "link")
# Graficar el logit ajustado contra 'edad'
plot(baseerc$edad, logit_ajustado, xlab = "Edad", ylab = "Logit Ajustado",
main = "Logit Ajustado vs Edad", pch = 20, col = "darkgreen")
Interpretación:
Este gráfico debe mostrar una relación lineal entre el logit ajustado y edad si la suposición de linealidad se cumple. Si se observa una curva u otro patrón, puede que haya una relación no lineal.
1.3. Test de Box Tidwell
El test se plantea tomando como hipótesis nula el cumplimineto del supuesto de linealidad.
H0: Existe linealidad entre la edad y la ERC
library(car)
# Convertir la variable de respuesta 'erc' en numérica (0 y 1)
baseerc <- baseerc %>%
mutate(erc_num = as.numeric(erc) - 1) # Restar 1 para que quede como 0 y 1
# Ejecutar el Test de Box-Tidwell
resultado_box_tidwell <- boxTidwell(erc_num ~ edad, data = baseerc)
print(resultado_box_tidwell)## MLE of lambda Score Statistic (t) Pr(>|t|)
## 0.71941 -0.3466 0.7293
##
## iterations = 6Interpretación:
El valor p de 0.7293 indica que no hay evidencia significativa para rechazar la suposición de linealidad entre la variable edad y la log-odds de la variable dependiente erc. Esto sugiere que la relación puede ser considerada lineal, y no es necesario aplicar una transformación no lineal a la variable edad en este caso.
1.4. Linealidad entre IMC y circunferencia abdominal
Este ejercicio en el caso que la variable de resultado sera distinta a ERC
# Gráfico de dispersión para verificar la linealidad sin acentos
plot(baseerc$IMC, baseerc$per_abdo,
main = "Relacion entre IMC y Circunferencia abdominal",
xlab = "IMC", ylab = "Circunferencia abdominal")
Linea de tendencia
# Convertir el título a UTF-8
titulo <- iconv("Relación entre IMC y Circunferencia abdominal", from = "UTF-8", to = "UTF-8")
# Ajustar el modelo de regresión lineal
modelo <- lm(per_abdo ~ IMC, data = baseerc)
# Gráfico de dispersión
plot(baseerc$IMC, baseerc$per_abdo,
main = titulo,
xlab = "IMC", ylab = "Circunferencia abdominal")
# Agregar la línea de tendencia
abline(modelo, col = "red", lwd = 2) # Línea roja con mayor grosor
1.5 Evaluación de la linealidad por grupo
ggplot(baseerc, aes(x = IMC, y = per_abdo, color = erc)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
facet_wrap(~ erc) +
labs(title = "Relacion entre IMC y Circunferencia abdominal por ERC",
x = "IMC",
y = "Circunferencia abdominal") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
2. Ausencia de multicolinealidad
2.1. Correlación entre variables numéricas / Matriz de correlación
# Seleccionar solo las variables numéricas
numeric_vars <- baseerc_numeric %>%
select(where(is.numeric))
# Calcular la matriz de correlación
cor_matrix <- cor(numeric_vars, use = "complete.obs")## Warning in cor(numeric_vars, use = "complete.obs"): the standard deviation is
## zero# Visualizar la matriz de correlación
corrplot(cor_matrix, method = "circle", type = "lower", tl.cex = 0.7, diag = FALSE)
2.2 Filtrar las Correlaciones en el rango deseado
# Función para extraer las correlaciones en el rango deseado
filter_correlations <- function(cor_matrix, min_cor, max_cor) {
# Extraer las correlaciones en el rango deseado
cor_matrix[cor_matrix < min_cor | cor_matrix > max_cor] <- NA
return(cor_matrix)
}
# Aplicar el filtro a la matriz de correlación
filtered_cor_matrix <- filter_correlations(cor_matrix, -0.7, 0.7)
# Visualizar la matriz de correlación filtrada
library(corrplot)
corrplot(filtered_cor_matrix, method = "circle", type = "lower", tl.cex = 0.7, diag = FALSE, na.label = "NA")
2.3 Identificar las Correlaciones Fuera del Rango
# Crear una tabla con las correlaciones
cor_df <- as.data.frame(as.table(cor_matrix))
# Filtrar correlaciones fuera del rango -0.7 a 0.7
cor_out_of_range <- cor_df %>%
filter(Freq < -0.7 | Freq > 0.7)
# Mostrar las correlaciones fuera del rango
print(cor_out_of_range)## Var1 Var2 Freq
## 1 cor cor 1.0000000
## 2 inic inic 1.0000000
## 3 erc erc 1.0000000
## 4 estad_erc erc -0.9313039
## 5 ano_diag erc -0.7639696
## 6 erc estad_erc -0.9313039
## 7 estad_erc estad_erc 1.0000000
## 8 ano_diag estad_erc 0.7249331
## 9 erc ano_diag -0.7639696
## 10 estad_erc ano_diag 0.7249331
## 11 ano_diag ano_diag 1.0000000
## 12 munic munic 1.0000000
## 13 area_ur area_ur 1.0000000
## 14 edad edad 1.0000000
## 15 sex sex 1.0000000
## 16 act_temper sex 0.7208758
## 17 educ educ 1.0000000
## 18 niv_educ niv_educ 1.0000000
## 19 est_civ est_civ 1.0000000
## 20 act_lab act_lab 1.0000000
## 21 act_lab_otr act_lab_otr 1.0000000
## 22 agric agric 1.0000000
## 23 act_temper agric 0.7120117
## 24 uso_plagui agric 0.7381646
## 25 epp_plagui agric 0.7376856
## 26 tabac tabac 1.0000000
## 27 alcoh alcoh 1.0000000
## 28 sal sal 1.0000000
## 29 agua_adec agua_adec 1.0000000
## 30 No_agua_adec agua_adec -1.0000000
## 31 agua_adec No_agua_adec -1.0000000
## 32 No_agua_adec No_agua_adec 1.0000000
## 33 frut_verd frut_verd 1.0000000
## 34 No_frut_verd frut_verd -1.0000000
## 35 frut_verd No_frut_verd -1.0000000
## 36 No_frut_verd No_frut_verd 1.0000000
## 37 ejerc ejerc 1.0000000
## 38 No_ejerc ejerc -1.0000000
## 39 ejerc No_ejerc -1.0000000
## 40 No_ejerc No_ejerc 1.0000000
## 41 sex act_temper 0.7208758
## 42 agric act_temper 0.7120117
## 43 act_temper act_temper 1.0000000
## 44 bajpes_prem bajpes_prem 1.0000000
## 45 dm dm 1.0000000
## 46 dm_tx dm 0.9927944
## 47 dm dm_tx 0.9927944
## 48 dm_tx dm_tx 1.0000000
## 49 chol chol 1.0000000
## 50 chol_tx chol 0.9804473
## 51 chol chol_tx 0.9804473
## 52 chol_tx chol_tx 1.0000000
## 53 hta hta 1.0000000
## 54 hta_tx hta 0.9970836
## 55 hta hta_tx 0.9970836
## 56 hta_tx hta_tx 1.0000000
## 57 enf_autoin enf_autoin 1.0000000
## 58 obesi obesi 1.0000000
## 59 IMC obesi 0.7455104
## 60 obesib obesi 1.0000000
## 61 ivu_rec ivu_rec 1.0000000
## 62 lit_ren lit_ren 1.0000000
## 63 af_dm af_dm 1.0000000
## 64 af_hta af_hta 1.0000000
## 65 af_acv af_acv 1.0000000
## 66 af_erc af_erc 1.0000000
## 67 af_iam af_iam 1.0000000
## 68 af_dislip af_dislip 1.0000000
## 69 af_otra af_otra 1.0000000
## 70 af_otra_especif af_otra_especif 1.0000000
## 71 agric uso_plagui 0.7381646
## 72 uso_plagui uso_plagui 1.0000000
## 73 fum_nfum uso_plagui -0.9434048
## 74 epp_plagui uso_plagui 0.9639650
## 75 tip_activiplagui tip_activiplagui 1.0000000
## 76 uso_plagui fum_nfum -0.9434048
## 77 fum_nfum fum_nfum 1.0000000
## 78 epp_plagui fum_nfum -0.8994195
## 79 agric epp_plagui 0.7376856
## 80 uso_plagui epp_plagui 0.9639650
## 81 fum_nfum epp_plagui -0.8994195
## 82 epp_plagui epp_plagui 1.0000000
## 83 aines aines 1.0000000
## 84 per_abdo per_abdo 1.0000000
## 85 peso per_abdo 0.7409217
## 86 IMC per_abdo 0.7510180
## 87 per_abdories per_abdories 1.0000000
## 88 estatur_cm estatur_cm 1.0000000
## 89 estat_mt estatur_cm 1.0000000
## 90 estatur_cm estat_mt 1.0000000
## 91 estat_mt estat_mt 1.0000000
## 92 per_abdo peso 0.7409217
## 93 peso peso 1.0000000
## 94 IMC peso 0.8101746
## 95 obesi IMC 0.7455104
## 96 per_abdo IMC 0.7510180
## 97 peso IMC 0.8101746
## 98 IMC IMC 1.0000000
## 99 obesib IMC 0.7455104
## 100 est_nutr est_nutr 1.0000000
## 101 sobrepes est_nutr 0.8599721
## 102 sobpes_obes est_nutr 0.8264970
## 103 obesi obesib 1.0000000
## 104 IMC obesib 0.7455104
## 105 obesib obesib 1.0000000
## 106 est_nutr sobrepes 0.8599721
## 107 sobrepes sobrepes 1.0000000
## 108 est_nutr sobpes_obes 0.8264970
## 109 sobpes_obes sobpes_obes 1.0000000Interpretación
Se excluirán aquellas variables con una alta correlación entre ellas (<0.7 ó >0.7)
3. Independencia entre las observaciones
La independencia entre observaciones es un supuesto clave. En estudios retrospectivos, esto se verifica asegurando que no hay dependencia entre individuos. Si las observaciones provienen de individuos únicos, puedes asumir independencia. Se puede revisar el diseño del estudio para verificar esta independencia.
3.1. Gráfico de Residuos vs. Predicciones
## Warning: package 'lmtest' was built under R version 4.3.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# Ajustar el modelo logístico
modelo_log <- glm(erc ~ area_ur + edad + sex + educ + agric +
tabac + alcoh + sal + No_agua_adec + No_frut_verd + No_ejerc +
act_temper + dm + hta + obesi + ivu_rec + uso_plagui + aines, data = baseerc, family = binomial)
# Residuos estandarizados
residuos <- residuals(modelo_log, type = "pearson")
# Predicciones
predicciones <- predict(modelo_log, type = "response")
# Gráfico de Residuos vs Predicciones
plot(predicciones, residuos,
main = "Residuos vs Predicciones",
xlab = "Predicciones", ylab = "Residuos")
abline(h = 0, col = "red")
INTERPRETACIÓN
El gráfico de Residuos vs Predicciones proporciona información para evaluar la calidad del ajuste de un modelo de regresión logística.
Interpretación del gráfico:
- Residuos dispersos aleatoriamente:
Si los residuos se distribuyen aleatoriamente alrededor de la línea horizontal (en 0), indica que el modelo se ajusta bien. La ausencia de patrones visibles sugiere que los residuos son independientes y que el modelo describe adecuadamente la relación entre las variables.
- Patrones sistemáticos en los residuos:
Si se observan patrones como tendencias, curvaturas o agrupamientos en los residuos, esto podría sugerir problemas con el modelo. Algunos de estos problemas incluyen:
Curvatura: Podría indicar que la relación entre las predicciones y la variable dependiente es no lineal y no está siendo bien capturada por el modelo.
Agrupamientos: La presencia de grupos de residuos en ciertas áreas del gráfico puede ser un signo de heterocedasticidad, lo que significa que la variabilidad de los residuos no es constante.
Pendiente inclinada: Si los residuos aumentan o disminuyen sistemáticamente con las predicciones, podría haber una especificación incorrecta del modelo o variables importantes que no fueron incluidas.
- Outliers: Residuos que están alejados significativamente del resto (outliers) pueden indicar la presencia de puntos influyentes que están afectando el ajuste del modelo y que podrían necesitar ser investigados más a fondo.
Si no se identifican patrones claros en el gráfico, se puede concluir que el modelo es adecuado. Sin embargo, si aparecen tendencias o agrupamientos, es posible que se requieran ajustes adicionales en el modelo, como la inclusión de más variables o la transformación de las existentes.
Para una evaluación estadística más rigurosa, se pueden emplear pruebas adicionales, como el test de Durbin-Watson para la autocorrelación de los residuos, o análisis de heterocedasticidad.
Prueba de Durbin-Watson
H0: Existe independencia entre las variables
##
## Durbin-Watson test
##
## data: modelo_log
## DW = 2.0585, p-value = 0.6081
## alternative hypothesis: true autocorrelation is greater than 0Interpretación: Dado que el p-valor es 0.6081, mayor que el nivel de significancia común (0.05), no se rechaza la hipótesis nula. Esto indica que no hay evidencia de autocorrelación en los residuos del modelo logístico.
4. Regresión logística
# Ajustar el modelo
modelo <- glm(erc ~ (area_ur + edad + sex + educ + agric +
tabac + alcoh + sal + No_agua_adec + No_frut_verd + No_ejerc +
act_temper + dm + hta + obesi + ivu_rec + uso_plagui + aines),
data = baseerc, family = binomial)
# Resumen del modelo
summary(modelo)##
## Call:
## glm(formula = erc ~ (area_ur + edad + sex + educ + agric + tabac +
## alcoh + sal + No_agua_adec + No_frut_verd + No_ejerc + act_temper +
## dm + hta + obesi + ivu_rec + uso_plagui + aines), family = binomial,
## data = baseerc)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.82270 2.25658 -4.353 1.34e-05 ***
## area_ururb 0.29073 0.62749 0.463 0.643138
## edad 0.04162 0.02545 1.636 0.101925
## sexmas 1.71962 1.25910 1.366 0.172017
## educSi -0.46912 0.82524 -0.568 0.569715
## agricSi 2.49150 1.20727 2.064 0.039042 *
## tabacSi 1.25704 1.38694 0.906 0.364755
## alcohSi 3.06674 1.51294 2.027 0.042662 *
## salSi 1.00411 0.67220 1.494 0.135236
## No_agua_adecSi 1.07123 0.68251 1.570 0.116520
## No_frut_verdSi 2.32924 0.69676 3.343 0.000829 ***
## No_ejercSi 2.71824 0.90516 3.003 0.002673 **
## act_temperSi -1.02789 1.48016 -0.694 0.487404
## dmSi 1.85056 0.71920 2.573 0.010079 *
## htaSi 1.85723 0.70146 2.648 0.008105 **
## obesiSi -0.40857 0.72044 -0.567 0.570635
## ivu_recSi 2.66114 0.99261 2.681 0.007341 **
## uso_plaguiSi 1.84699 0.89475 2.064 0.038993 *
## ainesSi 0.07584 2.03429 0.037 0.970260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 221.507 on 173 degrees of freedom
## Residual deviance: 88.016 on 155 degrees of freedom
## AIC: 126.02
##
## Number of Fisher Scoring iterations: 7Convertir los coeficientes en OR, agregar valor de p e intervalos de confianza
# Resumen del modelo
modelo_resumen <- summary(modelo)
# Calcular los Odds Ratios
odds_ratios <- exp(coef(modelo))
# Calcular los intervalos de confianza
ci <- exp(confint(modelo))## Waiting for profiling to be done...# Extraer los valores p
valores_p <- modelo_resumen$coefficients[, 4]
# Crear un data frame con los resultados, redondeando a tres decimales
resultados <- data.frame(
Coeficiente = round(coef(modelo), 3),
OR = round(odds_ratios, 3),
IC_inf = round(ci[, 1], 3),
IC_sup = round(ci[, 2], 3),
p_value = round(valores_p, 3)
)
# Imprimir los resultados
print(resultados)## Coeficiente OR IC_inf IC_sup p_value
## (Intercept) -9.823 0.000 0.000 0.003 0.000
## area_ururb 0.291 1.337 0.392 4.753 0.643
## edad 0.042 1.043 0.994 1.100 0.102
## sexmas 1.720 5.582 0.524 79.534 0.172
## educSi -0.469 0.626 0.120 3.194 0.570
## agricSi 2.491 12.079 1.160 143.354 0.039
## tabacSi 1.257 3.515 0.219 58.137 0.365
## alcohSi 3.067 21.472 1.456 715.970 0.043
## salSi 1.004 2.729 0.745 10.791 0.135
## No_agua_adecSi 1.071 2.919 0.791 11.964 0.117
## No_frut_verdSi 2.329 10.270 2.850 45.294 0.001
## No_ejercSi 2.718 15.154 2.999 109.582 0.003
## act_temperSi -1.028 0.358 0.016 6.061 0.487
## dmSi 1.851 6.363 1.662 29.389 0.010
## htaSi 1.857 6.406 1.730 28.241 0.008
## obesiSi -0.409 0.665 0.156 2.734 0.571
## ivu_recSi 2.661 14.313 2.247 116.352 0.007
## uso_plaguiSi 1.847 6.341 1.168 41.200 0.039
## ainesSi 0.076 1.079 0.018 80.454 0.9704.1 Evaluación del desempeño del modelo
Matriz de confusión y métricas de evaluación
# Predecir probabilidades de éxito
predicciones_prob <- predict(modelo, type = "response")
# Convertir probabilidades en clases (0 o 1) usando un umbral de 0.5
umbral <- 0.5
predicciones_clases <- ifelse(predicciones_prob > umbral, 1, 0)
# Crear la matriz de confusión
matriz_confusion <- table(Prediccion = predicciones_clases, Realidad = baseerc$erc)
# Imprimir la matriz de confusión
print(matriz_confusion)## Realidad
## Prediccion No Si
## 0 105 11
## 1 11 47# Extraer valores de la matriz de confusión
TP <- matriz_confusion[2, 2] # Verdaderos positivos
TN <- matriz_confusion[1, 1] # Verdaderos negativos
FP <- matriz_confusion[2, 1] # Falsos positivos
FN <- matriz_confusion[1, 2] # Falsos negativos
# Calcular métricas
exactitud <- (TP + TN) / sum(matriz_confusion) # Exactitud
sensibilidad <- TP / (TP + FN) # Tasa de verdaderos positivos (Sensibilidad)
especificidad <- TN / (TN + FP) # Tasa de verdaderos negativos (Especificidad)
vpp <- TP / (TP + FP) # Valor predictivo positivo
vpn <- TN / (TN + FN) # Valor predictivo negativo
f1 <- 2 * (vpp * sensibilidad) / (vpp + sensibilidad) # F1 Score
# Mostrar métricas
cat("Matriz de Confusión:\n")## Matriz de Confusión:## Realidad
## Prediccion No Si
## 0 105 11
## 1 11 47##
## Exactitud: 0.874## Sensibilidad: 0.81## Especificidad: 0.905## Valor Predictivo Positivo: 0.81## Valor Predictivo Negativo: 0.905## F1 Score: 0.814.2 Pruebas de hipótesis del modelo
Test de Wald
# Cargar la biblioteca necesaria
library(lmtest)
# Ajustar el modelo completo
# Realizar el Test de Wald
wald_test <- waldtest(modelo)
# Mostrar los resultados del Test de Wald
print("Resultados del Test de Wald:")## [1] "Resultados del Test de Wald:"## Wald test
##
## Model 1: erc ~ (area_ur + edad + sex + educ + agric + tabac + alcoh +
## sal + No_agua_adec + No_frut_verd + No_ejerc + act_temper +
## dm + hta + obesi + ivu_rec + uso_plagui + aines)
## Model 2: erc ~ 1
## Res.Df Df F Pr(>F)
## 1 155
## 2 173 -18 1.9048 0.01901 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretación:
El valor de la estadística F es 1.9048 y el p-valor es 0.01901, lo que sugiere que al menos una de las variables tiene un efecto significativo en la variable dependiente erc. Por lo tanto, se concluye que el modelo completo mejora la predicción en comparación con un modelo nulo.
Test de razón de verosimilitudes (Likelihood Ratio Test, LRT)
## [1] "Resultados del Likelihood Ratio Test:"## Likelihood ratio test
##
## Model 1: erc ~ (area_ur + edad + sex + educ + agric + tabac + alcoh +
## sal + No_agua_adec + No_frut_verd + No_ejerc + act_temper +
## dm + hta + obesi + ivu_rec + uso_plagui + aines)
## Model 2: erc ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 19 -44.008
## 2 1 -110.753 -18 133.49 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Intrepretación
El test mostró un valor de chi-cuadrado de 133.49 con un p-valor menor a 2.2e-16. Esto significa que el modelo con los predictores es significativamente mejor que el modelo nulo. En otras palabras, al menos uno de los predictores es importante para predecir el estado de la enfermedad renal crónica (ERC).