Curso Bioestadística y Ciencia de Datos: Ejercicio Tercera Sesión
Jonathan Patricio Baldera BSc. MSc.
07/18/2024
Contexto
El conjunto de datos “Datos_osteoporosis” ha sido fruto de prácticas desarrolladas en las secciones de la asignatura Métodos en Bioestadística II del programa de Maestría en Bioestadística impartido por el Departamento de Epidemiología Clínica y Bioestadística de la Facultad de Medicina de la Pontificia Universidad Javeriana, Bogotá, Colombia; sin embargo, la información mostrada acá es responsabilidad del tutor. Además, cabe mencionar que estos datos solo deben ser usados con fines académicos, no pueden ser utilizados para hacer inferencias ni podrán ser extrapolados para otras poblaciones u otros contextos diferentes al objetivo de la sesión.
La osteoporosis es una enfermedad en la que se produce una reducción de la masa ósea y un deterioro del tejido óseo que conlleva un aumento del riesgo de padecer una fractura. La mayoría de las fracturas debidas a la osteoporosis se producen en mujeres de edad avanzada. Los datos del archivo “Datos_osteoporosis” que se encuentra en la plataforma, buscan estudiar las características de una población a partir de una muestra de 400 mujeres entre 45 y 69 años procedentes de una determinada localidad. Esta base de datos contiene las siguientes variables:
| Variable | Descripción | Categorías |
|---|---|---|
| id | Número de identificación | |
| edad | Edad de la paciente en años | |
| imc | Índice de Masa Corporal | |
| dmo | densidad mineral ósea | |
| menopausia | Tipo de menopausia | 0. No menopausia, 1. Natural, 2. Ovariotomía/Histerectomía |
| antmaternos | Antecedentes Maternos | (0. No 1. Si) |
| nivel_ed | Nivel educativo | 1. Sin Estudios 2. Primaria_Incomp 3. Primaria_comp 4. Bachillerato 5. Superiores |
| paridad | Tuvo Hijos | 0. No 1. Si |
| numerohijos | Número de hijos nacidos vivos | |
| abortos | Número de abortos |
Instalación de las librerías requeridas
Solo si no han sido previamente instaladas, favor de instalar las siguientes librerías, las cuales serán necesarias para realizar este taller.
install.packages("readxl")
install.packages("table1")
install.packages("stargazer")
install.packages("dplyr")
install.packages("knitr")
install.packages("car")
install.packages("MASS")
install.packages("lmtest")
install.packages("cutpointr")
install.packages("caret")
install.packages("vcd")
install.packages("pROC")
install.packages("DescTools")Cargando las librerías necesarias
library(readxl)
library(table1)
library(stargazer)
library(dplyr)
library(knitr)
library(car)
library(MASS)
library(lmtest)
library(cutpointr)
library(caret)
library(vcd)
library(pROC)
library(DescTools)Categorizando las variables
datos$menopausia <- factor(datos$menopausia, labels = c("No menopausia","Natural","Ovariotomía/Histerectomía"))
datos$antmaternos <- factor(datos$antmaternos, labels = c("No","Si"))
datos$nivel_ed <- factor(datos$nivel_ed, labels = c("Sin Estudios","Primaria Incompleta","Primaria completa","Secundaria","Superior"))
datos$paridad <- factor(datos$paridad, labels = c("No","Si"))Recodificando las variables
A continuación, haremos uso de la función cut de la
librería base para recodificar la variable edad en
grupos:
- Menor de 55 años
- 55 años o más
# Recodificación de la edad en categorías
datos$edad_cat <- factor(cut(datos$edad,
breaks = c(-Inf,54.9,Inf),
labels = c("< 55 años","Mayor_o_igual_55")))
# Verificar la nueva variable
label(datos$edad_cat) <- "Grupo edad"
table1( ~ edad_cat, data = datos,
overall = "General",
rowlabelhead = "Variable",
caption = "Tabla 2: Número de mujeres según grupo edad")| Variable | General (N=400) |
|---|---|
| Grupo edad | |
| < 55 años | 248 (62.0%) |
| Mayor_o_igual_55 | 152 (38.0%) |
Recodificando la variable IMC en categorías basadas en los humbrales definidos del índice de masa corporal:
Bajo peso: IMC < 18.5 Peso normal: 18.5 <= IMC < 24.9 Sobrepeso: 25 <= IMC < 29.9 Obesidad: IMC >= 30
# Recodificación de la variable IMC en categorías
datos$imc_cat <- factor(cut(datos$imc,
breaks = c(-Inf,24.99999,29.99999,Inf),
labels = c("Normal","Sobrepeso","Obesidad")))
# Verificar la nueva variable
label(datos$imc_cat) <- "Grupo IMC"
table1( ~ imc_cat,
data = datos,
overall = "General", rowlabelhead = "Variable",
caption = "Tabla 3: Número de mujeres según grupo Índice de Masa Corporal")| Variable | General (N=400) |
|---|---|
| Grupo IMC | |
| Normal | 97 (24.3%) |
| Sobrepeso | 199 (49.8%) |
| Obesidad | 104 (26.0%) |
Recodificando la variable dmo en grupos:
- Menor o igual a 40 (Osteoporosis)
- 41 a 62 (Osteopenia)
- Más de 62 (Normal)
# Recodificación de la variable dmo en categorías
datos$dmo_cat <- factor(cut(datos$dmo,
breaks = c(-Inf,40,62,Inf),
labels = c("Osteoporosis","Osteopenia","Normal")))
# Verificar la nueva variable
label(datos$dmo_cat) <- "Grupo DMO"
table1( ~ dmo_cat,
data = datos,
overall = "General", rowlabelhead = "Variable",
caption = "Tabla 4: Número de mujeres según grupo Densidad mineral Ósea")| Variable | General (N=400) |
|---|---|
| Grupo DMO | |
| Osteoporosis | 9 (2.3%) |
| Osteopenia | 89 (22.3%) |
| Normal | 302 (75.5%) |
- Dicotomizando la variable abortos según la ocurrencia o no de este:
# Recodificación de la variable abortos en categorías
datos$abortos_cat <- factor(cut(datos$abortos,
breaks = c(-Inf,0,Inf),
labels = c("No","Si")))
# Verificar la nueva variable
label(datos$abortos_cat) <- "Aborto"
table1( ~ abortos_cat,
data = datos,
overall = "General",
rowlabelhead = "Variable",
caption = "Tabla 5: Número de mujeres según ocurrencia o no de aborto")| Variable | General (N=400) |
|---|---|
| Aborto | |
| No | 300 (75.0%) |
| Si | 100 (25.0%) |
Análisis descriptivo de los datos
En la Tabla 1 se muestran las características de las 400 pacientes incluidas en nuestra base de datos. El promedio de edad de estas fue de 53 años (d.e. 6.56), con un rango de edades desde 45 a 69 años. Con respecto al nivel educativo de estas, cerca del 11% reportaron no tener nivel educativo, alrededor del 69% indicaron tener educación primaria completa o incompleta, cerca del 16% reportaron haber tenido una educación bachiller y el 5% tiene algún tipo de educación superior.
Por otro lado, los niveles promedio respecto a la densidad mineral ósea son alrededor del 74 (d.e. 17.50), con un rango de valores desde 11 a 135. El índice de masa corporal promedio fue de 28 (d.e. 4.49), con valores que van desde 18.2 hasta 48.4.
Con respecto a los aspectos ginecológicos y maternos de las mujeres, cerca del 92% de estas indicaron no haber tenido antecedentes maternos. El 95.5% afirmó haber tenido al menos un hijo. La razón hijos fue de 2.56 hijos por mujer, con un rango de ninguno a cinco (5) hijos. Por otro lado la razón de abortos fue de 31 abortos por cada 100 mujeres.
Cabe destacar que el nivel de completitud fue de un 100% para todas las variable de la base de datos.
label(datos$edad) <- "Edad"
label(datos$nivel_ed) <- "Nivel educativo"
label(datos$imc) <- "IMC"
label(datos$dmo) <- "DMO"
label(datos$menopausia) <- "Menopausia"
label(datos$antmaternos) <- "Antecedentes maternos"
label(datos$paridad) <- "Paridad"
label(datos$numerohijos) <- "Numero de hijos"
label(datos$abortos) <- "Abortos"
table1( ~ edad + edad_cat + nivel_ed + imc + imc_cat + dmo + dmo_cat+ menopausia + antmaternos + paridad + numerohijos + abortos_cat,
data = datos,
overall = "General",
rowlabelhead = "Características",
caption = "Tabla 1: Características demográficas y clínicas de las mujeres")| Características | General (N=400) |
|---|---|
| Edad | |
| Mean (SD) | 53.1 (6.56) |
| Median [Min, Max] | 52.0 [45.0, 69.0] |
| Grupo edad | |
| < 55 años | 248 (62.0%) |
| Mayor_o_igual_55 | 152 (38.0%) |
| Nivel educativo | |
| Sin Estudios | 43 (10.8%) |
| Primaria Incompleta | 75 (18.8%) |
| Primaria completa | 199 (49.8%) |
| Secundaria | 63 (15.8%) |
| Superior | 20 (5.0%) |
| IMC | |
| Mean (SD) | 28.0 (4.49) |
| Median [Min, Max] | 27.4 [18.2, 48.4] |
| Grupo IMC | |
| Normal | 97 (24.3%) |
| Sobrepeso | 199 (49.8%) |
| Obesidad | 104 (26.0%) |
| DMO | |
| Mean (SD) | 74.1 (17.5) |
| Median [Min, Max] | 73.5 [11.0, 136] |
| Grupo DMO | |
| Osteoporosis | 9 (2.3%) |
| Osteopenia | 89 (22.3%) |
| Normal | 302 (75.5%) |
| Menopausia | |
| No menopausia | 136 (34.0%) |
| Natural | 206 (51.5%) |
| Ovariotomía/Histerectomía | 58 (14.5%) |
| Antecedentes maternos | |
| No | 367 (91.8%) |
| Si | 33 (8.3%) |
| Paridad | |
| No | 18 (4.5%) |
| Si | 382 (95.5%) |
| Numero de hijos | |
| Mean (SD) | 2.56 (1.11) |
| Median [Min, Max] | 2.00 [0, 5.00] |
| Aborto | |
| No | 300 (75.0%) |
| Si | 100 (25.0%) |
De ser necesario, reagrupe las variables que considere.
Regresión logística
Evalué si tiene sentido hacer un modelo de regresión logística para evaluar si existe asociación entre la densidad mineral ósea en categorías y el índice de masa corporal, controlando por las demás covariables.
R.¿?Si la respuesta al punto anterior es no, ¿qué debe hacer para realizar el modelo anterior?
R.¿?# Dicotomizando la variable dmo en valores normal de dmo (mayores a 62) y por
# debajo de lo normal (≤ 62)
datos$dmo_normal <- factor(cut(datos$dmo,
breaks = c(-Inf,62,Inf),
labels = c("Bajo","Normal")))
# Se cambiará la categoría de referencia dado que la variable de respuesta
# presenta "Bajo" como primera categoría, es necesario hacer un relevel para
# dejar de base "Normal".
datos$dmo_normal <- relevel(datos$dmo_normal, ref="Normal")
# Verificar la nueva variable
label(datos$dmo_normal) <- "DMO"
table1(~ dmo_normal,
data = datos,
overall = "General", rowlabelhead = "Variable",
caption = "Tabla 6: Número de mujeres según valores de densidad mineral ósea")| Variable | General (N=400) |
|---|---|
| DMO | |
| Normal | 302 (75.5%) |
| Bajo | 98 (24.5%) |
Modelo de asociación
En este punto, debemos resaltar que queremos conocer si existe asociación entre la densidad mineral ósea y y el índice de masa corporal, controlando por las demás covariables. Antes de continuar, es importante conocer el rol de cada variable antes de plantear el modelo.
Para el punto anterior, clasifique las variables según su función en el modelo. (Nota: de considerarlo necesario, incluya interacciones)*
Variables:
1. Dependiente
- Densidad mineral ósea, dicotomizada como normal (>62) y bajo (≤ 62).
2. Exposición
- Indice de masa corporar, categorizado como normal, sobrepeso y obeso.
3. Control
Edad, agrupada como menores y mayores o igual a 55 años.
Educación, dicotomizada como sin estudios, Primaria_Incomp, Primaria_comp, Bachillerato y Superiores.
Menopausia, categorizada como no menopausica, natural e ovariotomía/histerectomía
Antecedentes maternos, dicotomizada como no y si.
Paridad, dicotomizada si tuvo hijos o no.
Abortos, dicotomizada si tuvo abortos o no.
Plantee y corra el modelo de completo para el punto anterior.
Modelo:
\[logit(dmo)=\beta_{0}+\beta_{imc}+\beta_{edad}+\beta_{educación}+\beta_{menopausia}+\beta_{ant.maternos}+\beta_{paridad}+\beta_{abortos}+u\] ***
# Realizando el modelo logistico
m1 <- glm(dmo_normal ~ imc_cat + edad_cat + nivel_ed + menopausia + antmaternos + paridad + abortos_cat,
data = datos,
family = binomial("logit"))
# Obteniendo las estadísticas de resumen del modelo
summary(m1)##
## Call:
## glm(formula = dmo_normal ~ imc_cat + edad_cat + nivel_ed + menopausia +
## antmaternos + paridad + abortos_cat, family = binomial("logit"),
## data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.58965 0.72277 -0.816 0.414600
## imc_catSobrepeso -0.37704 0.33979 -1.110 0.267165
## imc_catObesidad -0.96560 0.41193 -2.344 0.019075 *
## edad_catMayor_o_igual_55 1.18801 0.32387 3.668 0.000244 ***
## nivel_edPrimaria Incompleta 0.66170 0.44086 1.501 0.133375
## nivel_edPrimaria completa -0.15630 0.41868 -0.373 0.708914
## nivel_edSecundaria -0.06999 0.53550 -0.131 0.896013
## nivel_edSuperior 0.63031 0.67242 0.937 0.348566
## menopausiaNatural 0.53088 0.38960 1.363 0.172999
## menopausiaOvariotomía/Histerectomía -0.03626 0.48672 -0.074 0.940622
## antmaternosSi -1.50257 0.60426 -2.487 0.012895 *
## paridadSi -1.12143 0.54044 -2.075 0.037983 *
## abortos_catSi 0.34297 0.28345 1.210 0.226282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 445.42 on 399 degrees of freedom
## Residual deviance: 385.55 on 387 degrees of freedom
## AIC: 411.55
##
## Number of Fisher Scoring iterations: 5##
## ===============================================================
## Dependent variable:
## ---------------------------
## dmo_normal
## ---------------------------------------------------------------
## imc_catSobrepeso -0.377
## (0.340)
##
## imc_catObesidad -0.966**
## (0.412)
##
## edad_catMayor_o_igual_55 1.188***
## (0.324)
##
## nivel_edPrimaria Incompleta 0.662
## (0.441)
##
## nivel_edPrimaria completa -0.156
## (0.419)
##
## nivel_edSecundaria -0.070
## (0.536)
##
## nivel_edSuperior 0.630
## (0.672)
##
## menopausiaNatural 0.531
## (0.390)
##
## menopausiaOvariotomía/Histerectomía -0.036
## (0.487)
##
## antmaternosSi -1.503**
## (0.604)
##
## paridadSi -1.121**
## (0.540)
##
## abortos_catSi 0.343
## (0.283)
##
## Constant -0.590
## (0.723)
##
## ---------------------------------------------------------------
## Observations 400
## Log Likelihood -192.777
## Akaike Inf. Crit. 411.553
## ===============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Otra alternativa para conocer las estadísticas de resumen del modelo
es utilizando la función stargazer de la librería
stargazer.
##
## ===============================================================
## Dependent variable:
## ---------------------------
## dmo_normal
## ---------------------------------------------------------------
## imc_catSobrepeso -0.377
## (0.340)
##
## imc_catObesidad -0.966**
## (0.412)
##
## edad_catMayor_o_igual_55 1.188***
## (0.324)
##
## nivel_edPrimaria Incompleta 0.662
## (0.441)
##
## nivel_edPrimaria completa -0.156
## (0.419)
##
## nivel_edSecundaria -0.070
## (0.536)
##
## nivel_edSuperior 0.630
## (0.672)
##
## menopausiaNatural 0.531
## (0.390)
##
## menopausiaOvariotomía/Histerectomía -0.036
## (0.487)
##
## antmaternosSi -1.503**
## (0.604)
##
## paridadSi -1.121**
## (0.540)
##
## abortos_catSi 0.343
## (0.283)
##
## Constant -0.590
## (0.723)
##
## ---------------------------------------------------------------
## Observations 400
## Log Likelihood -192.777
## Akaike Inf. Crit. 411.553
## ===============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Factores de confusión
Evaluando si `aborto` es una variable de confusión en el modelo.m2 <- glm(dmo_normal ~ imc_cat + edad_cat + nivel_ed + menopausia + antmaternos + paridad,
data = datos,
family = binomial("logit"))
stargazer(m1,m2, type = "text")##
## ================================================================
## Dependent variable:
## ----------------------------
## dmo_normal
## (1) (2)
## ----------------------------------------------------------------
## imc_catSobrepeso -0.377 -0.348
## (0.340) (0.337)
##
## imc_catObesidad -0.966** -0.912**
## (0.412) (0.408)
##
## edad_catMayor_o_igual_55 1.188*** 1.197***
## (0.324) (0.323)
##
## nivel_edPrimaria Incompleta 0.662 0.674
## (0.441) (0.441)
##
## nivel_edPrimaria completa -0.156 -0.129
## (0.419) (0.418)
##
## nivel_edSecundaria -0.070 -0.043
## (0.536) (0.533)
##
## nivel_edSuperior 0.630 0.704
## (0.672) (0.668)
##
## menopausiaNatural 0.531 0.533
## (0.390) (0.388)
##
## menopausiaOvariotomía/Histerectomía -0.036 -0.049
## (0.487) (0.486)
##
## antmaternosSi -1.503** -1.491**
## (0.604) (0.600)
##
## paridadSi -1.121** -1.110**
## (0.540) (0.539)
##
## abortos_catSi 0.343
## (0.283)
##
## Constant -0.590 -0.563
## (0.723) (0.722)
##
## ----------------------------------------------------------------
## Observations 400 400
## Log Likelihood -192.777 -193.498
## Akaike Inf. Crit. 411.553 410.996
## ================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Al incluir o no la variable aborto en el modelo podemos ver que esta no es una variable de confusión, ya que los coeficientes de la variable de exposición no cambian de manera significativa en el modelo.
Interpretación del modelo
Exponenciando los coeficientes:
## OR 2.5 % 97.5 %
## (Intercept) 0.5545195 0.13325845 2.3068412
## imc_catSobrepeso 0.6858910 0.35200735 1.3404078
## imc_catObesidad 0.3807554 0.16767158 0.8477755
## edad_catMayor_o_igual_55 3.2805411 1.76087376 6.2944646
## nivel_edPrimaria Incompleta 1.9380829 0.82675779 4.6897523
## nivel_edPrimaria completa 0.8553032 0.38046844 1.9788016
## nivel_edSecundaria 0.9324032 0.32376014 2.6696914
## nivel_edSuperior 1.8781929 0.49184556 6.9904041
## menopausiaNatural 1.7004263 0.79417596 3.6900912
## menopausiaOvariotomía/Histerectomía 0.9643941 0.36179050 2.4689526
## antmaternosSi 0.2225574 0.05852536 0.6559818
## paridadSi 0.3258126 0.11094637 0.9465153
## abortos_catSi 1.4091222 0.80258279 2.4460523Diagnóstico del modelo
Para hacer el diagnóstico del modelo, sólo evaluaremos la presencia de valores atípicos y multicolinealidad, ya que no podemos avaluar linealidad en este modelo debido a que no tenemos variables predictoras en forma continua. Con respecto a la bondad de ajuste, tampoco evaluaremos, ya que es un modelo de asociación y no de predicción, por lo que no estamos buscando si nuestro modelo se ajusta o no a lo observado, sino que nuestro interés es saber si dmo está asociada con el IMC.
Valores atípicos
Para obtener los valores atípicos utilizaremos las distancias de cook, aplicando como punto de corte el resultado de obtener la siguiente expresión:
\(\frac{4}{N-k-1}\)
\(\frac{4}{400-5-1}\) = \(\frac{4}{394}\) = 0.01015228
En la Tabla 7 podemos ver el id de las observaciones que presentan distancias de cook más altas que el punto de corte calculado anteriormente, al comparar estos valores con los arrojados por el gráfico 1 arrojado de manera automática por el software, podemos ver que este sólo nos sugiere como valores atípicos las observaciones en la posición 83, 317 y 377. Al detectar estos valores atípicos y descartar posibles influencias por estos valores, las literaturas recomiendan excluir estos valores y correr nuevamente el modelo y observar si estos valores son influyentes o no en el modelo. En esta ocasión no recurriremos a ello, sino que ejecutaremos la prueba de bonferroni para datos extremos, la cual prueba la hipótesis nula de que cada observación es un valor atípico de cambio medio. Al ejecutar dicha prueba, se rechazó la hipótesis nula. (ver debajo)
cook <- as.data.frame(cooks.distance(m1))
cook %>% filter(cooks.distance(m1) > 0.01015228) %>%
kable(caption = "Tabla 7: id de los datos con distancia de cook mayores a 0.01020408")| cooks.distance(m1) | |
|---|---|
| 26 | 0.0103105 |
| 60 | 0.0108082 |
| 77 | 0.0121062 |
| 83 | 0.0387534 |
| 91 | 0.0114629 |
| 123 | 0.0130445 |
| 128 | 0.0119662 |
| 139 | 0.0204463 |
| 157 | 0.0143609 |
| 162 | 0.0185626 |
| 163 | 0.0151593 |
| 164 | 0.0137174 |
| 169 | 0.0104983 |
| 171 | 0.0132244 |
| 173 | 0.0103409 |
| 174 | 0.0178004 |
| 181 | 0.0198339 |
| 189 | 0.0178196 |
| 200 | 0.0180880 |
| 237 | 0.0186948 |
| 239 | 0.0111831 |
| 261 | 0.0167911 |
| 298 | 0.0189736 |
| 304 | 0.0146770 |
| 306 | 0.0146252 |
| 317 | 0.0239002 |
| 318 | 0.0102728 |
| 320 | 0.0123235 |
| 350 | 0.0125171 |
| 355 | 0.0202097 |
| 360 | 0.0132244 |
| 377 | 0.0422078 |
| 378 | 0.0119662 |
| 388 | 0.0204352 |
Realizando la prueba de bonferroni para datos extremos
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 377 2.813062 0.0049072 NAMulticolinealidad
Para evaluar la presencia o no de multicolinealidad en el modelo se realizara la prueba de los valores VIF:
## GVIF Df GVIF^(1/(2*Df))
## imc_cat 1.371392 2 1.082157
## edad_cat 1.641900 1 1.281366
## nivel_ed 1.565769 4 1.057648
## menopausia 1.706238 2 1.142904
## antmaternos 1.060220 1 1.029670
## paridad 1.062906 1 1.030973
## abortos_cat 1.021628 1 1.010756Al realizar la prueba, podemos ver que los valores para los coeficientes son cercanos a uno (1), por lo que descartamos la presencia de multicolinealidad.
Modelo de predicción
A continuación, vamos a realizar un nuevo modelo, pero en este escenario partimos de un objetivo diferente. En este escenario nos interesaría conocer cuales son las mujeres más propensas a tener una densidad ósea normal o baja.
Antes de continuar, vamos a dividir los datos utilizando el 80% de estos para entrenar el modelo y el 20% restante para probar el ajustey discriminación.
# Semilla inicial para garantizar la reproducibilidad en los resultados
set.seed(16)
# Número total de observaciones
n_total <- nrow(datos)
# Tamaño del conjunto de entrenamiento
n_train <- ceiling(n_total * 0.8)
# Crear el conjunto de entrenamiento
BD_entrenamiento <- datos %>% sample_n(size = n_train, replace = FALSE)
# Crear el conjunto de prueba seleccionando las filas que no están en el conjunto de entrenamiento
BD_prueba <- datos %>% filter(!id %in% BD_entrenamiento$id)Para la construcción del modelo de predicción se tomarán en cuenta las siguientes variables, las cuales serán utilizadas en su forma continua (las que apliquen), debido a que aportan más información como variables predictoras.
\[logit(dmo)=\beta_{0}+\beta_{edad}+\beta_{imc}+\beta_{educación}+\beta_{menopausia}+\beta_{ant.maternos}+\beta_{numerohijos}+\beta_{abortos}+u\]
m3 <- glm(dmo_normal ~ edad + imc + menopausia + antmaternos + nivel_ed + numerohijos + abortos,
data = BD_entrenamiento,
family = binomial("logit"))
summary(m3)##
## Call:
## glm(formula = dmo_normal ~ edad + imc + menopausia + antmaternos +
## nivel_ed + numerohijos + abortos, family = binomial("logit"),
## data = BD_entrenamiento)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.44696 1.86230 -3.462 0.000537 ***
## edad 0.14352 0.03075 4.667 3.05e-06 ***
## imc -0.08054 0.04061 -1.983 0.047331 *
## menopausiaNatural -0.09289 0.43671 -0.213 0.831556
## menopausiaOvariotomía/Histerectomía -0.43978 0.51880 -0.848 0.396613
## antmaternosSi -1.44979 0.69097 -2.098 0.035887 *
## nivel_edPrimaria Incompleta 0.48521 0.49099 0.988 0.323040
## nivel_edPrimaria completa -0.05179 0.46707 -0.111 0.911717
## nivel_edSecundaria 0.36466 0.58063 0.628 0.529975
## nivel_edSuperior 0.91212 0.78336 1.164 0.244276
## numerohijos -0.05665 0.12194 -0.465 0.642245
## abortos 0.14714 0.23114 0.637 0.524405
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 370.47 on 319 degrees of freedom
## Residual deviance: 315.76 on 308 degrees of freedom
## AIC: 339.76
##
## Number of Fisher Scoring iterations: 5Estraegias de selección
Excluyendo del modelo las variables no significativas utilizando la metodología stepwise.
## Start: AIC=339.76
## dmo_normal ~ edad + imc + menopausia + antmaternos + nivel_ed +
## numerohijos + abortos
##
## Df Deviance AIC
## - nivel_ed 4 319.69 335.69
## - menopausia 2 316.69 336.69
## - numerohijos 1 315.98 337.98
## - abortos 1 316.16 338.16
## <none> 315.76 339.76
## - imc 1 320.08 342.08
## - antmaternos 1 321.39 343.39
## - edad 1 340.04 362.04
##
## Step: AIC=335.69
## dmo_normal ~ edad + imc + menopausia + antmaternos + numerohijos +
## abortos
##
## Df Deviance AIC
## - menopausia 2 320.72 332.72
## - abortos 1 320.07 334.07
## - numerohijos 1 320.17 334.17
## <none> 319.69 335.69
## - antmaternos 1 324.90 338.90
## - imc 1 325.64 339.64
## + nivel_ed 4 315.76 339.76
## - edad 1 351.17 365.17
##
## Step: AIC=332.72
## dmo_normal ~ edad + imc + antmaternos + numerohijos + abortos
##
## Df Deviance AIC
## - abortos 1 321.10 331.10
## - numerohijos 1 321.19 331.19
## <none> 320.72 332.72
## + menopausia 2 319.69 335.69
## - antmaternos 1 325.91 335.91
## + nivel_ed 4 316.69 336.69
## - imc 1 327.07 337.07
## - edad 1 365.40 375.40
##
## Step: AIC=331.1
## dmo_normal ~ edad + imc + antmaternos + numerohijos
##
## Df Deviance AIC
## - numerohijos 1 321.53 329.53
## <none> 321.10 331.10
## + abortos 1 320.72 332.72
## + menopausia 2 320.07 334.07
## - antmaternos 1 326.39 334.39
## + nivel_ed 4 317.10 335.10
## - imc 1 327.31 335.31
## - edad 1 365.96 373.96
##
## Step: AIC=329.53
## dmo_normal ~ edad + imc + antmaternos
##
## Df Deviance AIC
## <none> 321.53 329.53
## + numerohijos 1 321.10 331.10
## + abortos 1 321.19 331.19
## + menopausia 2 320.51 332.51
## - antmaternos 1 326.58 332.58
## + nivel_ed 4 317.27 333.27
## - imc 1 329.10 335.10
## - edad 1 366.16 372.16##
## Call:
## glm(formula = dmo_normal ~ edad + imc + antmaternos, family = binomial("logit"),
## data = BD_entrenamiento)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.83122 1.28953 -4.522 6.13e-06 ***
## edad 0.14048 0.02271 6.186 6.18e-10 ***
## imc -0.09881 0.03807 -2.596 0.00944 **
## antmaternosSi -1.32217 0.66147 -1.999 0.04563 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 370.47 on 319 degrees of freedom
## Residual deviance: 321.53 on 316 degrees of freedom
## AIC: 329.53
##
## Number of Fisher Scoring iterations: 4Comprobando la hipótesis nula de que los coeficientes de las variables excluidas por la estrategia stepwise son iguales a cero (0).
## Likelihood ratio test
##
## Model 1: dmo_normal ~ edad + imc + menopausia + antmaternos + nivel_ed +
## numerohijos + abortos
## Model 2: dmo_normal ~ edad + imc + antmaternos
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 12 -157.88
## 2 4 -160.77 -8 5.7681 0.6732Con una confianza del 95%, no podemos rechazar la hipótesis nula de que los coeficientes de las variables excluidas por la estrategia setpwise son iguales a cero (0), por lo que nos quedamos con el nuevo modelo.
Calibración y discriminación
Mejor punto de corte
BD_prueba$predicciones <- predict(step,
newdata = BD_prueba,
type=c("response"))
BD_prueba$dmo_normal_num <- as.numeric(BD_prueba$dmo_normal)
cut_point <- cutpointr(BD_prueba, # base de datos con la que trabajo
predicciones, #variable que guarda las predicciones del modelo
dmo_normal_num, # variable dependiente convertida como numérica
direction = ">=",
pos_class = 2, # valores en mi variable observada un caso positivo
neg_class = 1, # valores en mi variable observada un caso negativo
method = maximize_metric, # Metodo para maximizar (se
metric = youden) # youden = sensitivity + specificity - 1
summary(cut_point)## Method: maximize_metric
## Predictor: predicciones
## Outcome: dmo_normal_num
## Direction: >=
##
## AUC n n_pos n_neg
## 0.7463 80 13 67
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 0.4268 0.4638 0.8625 0.5385 0.9254 7 6 5 62
##
## Predictor summary:
## Data Min. 5% 1st Qu. Median Mean 3rd Qu.
## Overall 0.02278770 0.06160613 0.1212705 0.1850035 0.2354665 0.2677784
## 1 0.02278770 0.05725612 0.1109519 0.1699847 0.2040777 0.2278231
## 2 0.09693791 0.12616457 0.1636795 0.4268500 0.3972395 0.5985858
## 95% Max. SD NAs
## 0.6002206 0.7871445 0.1710144 0
## 0.4755940 0.7074194 0.1369356 0
## 0.7090066 0.7871445 0.2356923 0El punto de corte calculado en el cual se hace máximo es 0.4268, por lo que utilizaremos este punto de corte para elaborar nuestra matriz de confusión.
BD_prueba$predicciones_cat <- as.factor(ifelse(BD_prueba$predicciones > 0.4268, "Si", "No"))
BD_prueba$dmo_normal_num_2 <- factor(BD_prueba$dmo_normal_num,labels = c("No","Si") )
confusionMatrix(BD_prueba$predicciones_cat, BD_prueba$dmo_normal_num_2, positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 62 6
## Si 5 7
##
## Accuracy : 0.8625
## 95% CI : (0.7673, 0.9293)
## No Information Rate : 0.8375
## P-Value [Acc > NIR] : 0.3349
##
## Kappa : 0.4787
##
## Mcnemar's Test P-Value : 1.0000
##
## Sensitivity : 0.5385
## Specificity : 0.9254
## Pos Pred Value : 0.5833
## Neg Pred Value : 0.9118
## Prevalence : 0.1625
## Detection Rate : 0.0875
## Detection Prevalence : 0.1500
## Balanced Accuracy : 0.7319
##
## 'Positive' Class : Si
## matriz_confusion <- table(BD_prueba$predicciones_cat,BD_prueba$dmo_normal_num_2,
dnn = c("Observaciones", "Predicciones"))
# Graficando la matriz de confusión
mosaic(matriz_confusion, shade = T, colorize = T,
gp = gpar(fill = matrix(c("green", "red3", "red3", "green"), 2, 2)))
Curva roc:
## Area under the curve: 0.7463
Interprete los resultados obtenidos
El modelo construido es capaz de clasificar correctamente alrededor del 86% de las observaciones (Accuracy = 0.8625), es decir, la capacidad del modelo de predecir a una mujer con valores de dmo menor a 62 cuando realmente el valor observado es menor de 62 y de predecir una mujer con valores de dmo mayor a 62 cuando realmente es mayor a 62 en las observaciones es de 86 por cada 100 predicciones.
La sensibilidad del modelo es del 0.5385, es decir, la capacidad del modelo de identificar correctamente a una mujer que realmente presenta valores de dmo menores a 62 es del 53.85%. Con respecto a la especificidad, esta fue de 0.9254, es decir, la capacidad del modelo de identificar correctamente a una mujer que realmente tiene valores de dmo mayores a 62 es del 92.54%.
Otros aspectos relevantes es el Pseudo \(R^2\) calculado fue de 0.2067819 y el el área calclulada bajo la curva roc es de 0.7463.
En términos generales se podría decir que no estamos frente al mejor modelo para predecir las mujeres con niveles de dmo por debajo de 62, ya que alrededor del 13.75% de las predicciones no logra clasificar correctamente. Asimismo podemos observar el valor de otras pruebas como es el valor del área bajo la curva roc (0.7463), el pseudo \(R^2\) (0.2068), entre otros.
Preguntas y Discusión
Antes de concluir, tomemos unos minutos para responder cualquier pregunta que puedan tener y discutir cualquier duda o comentario sobre los temas tratados hoy.
Próximos Pasos
Aplicación de lo Aprendido
Es importante que puedan aplicar estos conceptos y técnicas en sus propios proyectos y estudios. Les animo a:
Aplicar el análisis descriptivo y exploratorio como primer paso en cualquier análisis de datos. Considerar los supuestos necesarios antes de realizar cualquier prueba estadística y elegir la prueba adecuada en función de la naturaleza de sus datos. Utilizar visualizaciones para comprender mejor las relaciones entre variables y para comunicar sus hallazgos de manera efectiva.
Referencias
Achim Zeileis, Torsten Hothorn (2002). Diagnostic Checking in Regression Relationships. R News 2(3), 7-10. URL https://CRAN.R-project.org/doc/Rnews/
Andri Signorell et mult. al. (2021). DescTools: Tools for descriptive statistics. R package version 0.99.40.
Barbara McKnight (2021). UWbe536: UW BIOST/EPI 536 Utility Functions. R package version 0.1.0.
Benjamin Rich (2021). table1: Tables of Descriptive Statistics in HTML. R package version 1.3. https://CRAN.R-project.org/package=table1
Christian Thiele (2021). cutpointr: Determine and Evaluate Optimal Cutpoints in Binary Classification Tasks. R package version 1.1.0. https://CRAN.R-project.org/package=cutpointr
David Meyer, Achim Zeileis, and Kurt Hornik (2020). vcd: Visualizing Categorical Data. R package version 1.4-8.
Díaz Socorro Cossette, Navarro Despaigne Daysi, Santana Pérez Felipe, Domínguez Alonso Emma, Bacallao Gallestey Jorge. Factores de riesgo modificables o no, relacionados con la densidad mineral ósea en mujeres de edad mediana. Rev Cubana Endocrinol [Internet]. 2012 Abr [citado 2021 Mayo 24] ; 23( 1 ): 44-55. Disponible en: http://scielo.sld.cu/scielo.php?script=sci_arttext&pid=S1561-29532012000100004&lng=es.
Guerra R José René, Urdaneta M José Ramón, Villalobos I Noren, Contreras Benítez Alfi, García I José, Baabel Z Nasser Saleh et al . Factores de riesgo para alteraciones de la densidad mineral ósea en mujeres posmenopáusicas. Rev. chil. obstet. ginecol. [Internet]. 2015 Ago [citado 2021 Mayo 24] ; 80( 5 ): 385-393. Disponible en: http://www.scielo.cl/scielo.php?script=sci_arttext&pid=S0717-75262015000500006&lng=es. http://dx.doi.org/10.4067/S0717-75262015000500006.
Hadley Wickham and Jennifer Bryan (2019). readxl: Read Excel Files. R package version 1.3.1. https://CRAN.R-project.org/package=readxl
Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2021). dplyr: A Grammar of Data Manipulation. R package version 1.0.3. https://CRAN.R-project.org/package=dplyr
Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables. R package version 5.2.1. https://CRAN.R-project.org/package=stargazer
John Fox and Sanford Weisberg (2019). An {R} Companion to Applied Regression, Third Edition. Thousand Oaks CA: Sage. URL: https://socialsciences.mcmaster.ca/jfox/Books/Companion/
Matthew Jay (2019). generalhoslem: Goodness of Fit Tests for Logistic Regression Models. R package version 1.3.4. https://CRAN.R-project.org/package=generalhoslem
Max Kuhn (2021). caret: Classification and Regression Training. R package version 6.0-88. https://CRAN.R-project.org/package=caret
Torres-Mejía Gabriela, Guzmán Pineda Rubén, Téllez-Rojo Martha María, Lazcano-Ponce Eduardo. Peak bone mass and bone mineral density correlates for 9 to 24 year-old Mexican women, using corrected BMD. Salud pública Méx [revista en la Internet]. 2009 Ene [citado 2021 Mayo 24] ; 51( Suppl 1 ): s84-s92. Disponible en: http://www.scielo.org.mx/scielo.php?script=sci_arttext&pid=S0036-36342009000700011&lng=es.
Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0
Xavier Robin, Natacha Turck, Alexandre Hainard, Natalia Tiberti, Frédérique Lisacek, Jean-Charles Sanchez and Markus Müller (2011). pROC: an open-source package for R and S+to analyze and compare ROC curves. BMC Bioinformatics, 12, p. 77. DOI: 10.1186/1471-2105-12-77 http://www.biomedcentral.com/1471-2105/12/77/
Yihui Xie (2020). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.30.