Introducción

En este documento se presenta una serie de códigos estandarizados en R para facilitar y sistematizar el proceso de manejo y análisis de datos. El flujo de trabajo incluye:

  1. Preparación de los datos: carga de paquetes, exploración inicial y limpieza.
  2. Análisis descriptivo: variables categóricas y numéricas, y detección de valores atípicos.
  3. Análisis bivariado: categórica × categórica, categórica × numérica y numérica × numérica.
  4. Análisis multivariado: regresión logística binaria y regresión lineal múltiple.

Para propósitos de demostración, se utilizará la base de datos birthwt, disponible en el paquete MASS de R.

1 Preparación de los datos

1.1 Carga de paquetes

En esta sección se cargan los paquetes necesarios para el manejo, exploración, limpieza y análisis de los datos.

# Cargar paqueterías
library(MASS)       # Example dataset
library(dplyr)      # Data manipulation
library(janitor)    # Data cleaning and standardization
library(skimr)      # Data exploration and descriptive summaries
library(GGally)     # Correlation matrix plots (ggpairs)
library(car)        # VIF (multicollinearity)
library(lmtest)     # Breusch-Pagan test

1.2 Importación de la base de datos

# Importar base de datos que utilizaremos como ejemplo
data("birthwt")
data <- birthwt

1.3 Exploración inicial

Antes de realizar modificaciones en la base de datos, se recomienda llevar a cabo una exploración inicial para conocer su estructura, dimensiones, tipos de variables, valores faltantes y distribución general de los datos.

# Dimensiones de la base de datos
dim(data)
## [1] 189  10
# Nombres de las variables
names(data)
##  [1] "low"   "age"   "lwt"   "race"  "smoke" "ptl"   "ht"    "ui"    "ftv"  
## [10] "bwt"
# Estructura y tipo de las variables
str(data)
## 'data.frame':    189 obs. of  10 variables:
##  $ low  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : int  2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: int  0 0 1 1 1 0 0 0 1 1 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ui   : int  1 0 0 1 1 0 0 0 0 0 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
# Primeras observaciones
head(data)

1.3.1 Resumen general con skim()

La función skim() ofrece un resumen general de la base de datos: tipo de variable, número de valores faltantes (missing values), media y desviación estándar para variables numéricas, percentiles, mínimo, máximo y distribución general.

skim(data)
Data summary
Name data
Number of rows 189
Number of columns 10
_______________________
Column type frequency:
numeric 10
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
low 0 1 0.31 0.46 0 0 0 1 1 ▇▁▁▁▃
age 0 1 23.24 5.30 14 19 23 26 45 ▇▇▃▁▁
lwt 0 1 129.81 30.58 80 110 121 140 250 ▆▇▂▁▁
race 0 1 1.85 0.92 1 1 1 3 3 ▇▁▂▁▆
smoke 0 1 0.39 0.49 0 0 0 1 1 ▇▁▁▁▅
ptl 0 1 0.20 0.49 0 0 0 0 3 ▇▁▁▁▁
ht 0 1 0.06 0.24 0 0 0 0 1 ▇▁▁▁▁
ui 0 1 0.15 0.36 0 0 0 0 1 ▇▁▁▁▂
ftv 0 1 0.79 1.06 0 0 0 1 6 ▇▂▁▁▁
bwt 0 1 2944.59 729.21 709 2414 2977 3487 4990 ▁▅▇▆▁

2 Limpieza y estandarización de los datos

Antes de comenzar los análisis estadísticos, se recomienda verificar la consistencia de la base de datos e identificar posibles problemas relacionados con nombres de variables, valores duplicados, valores faltantes y categorías.

2.1 Estandarización de nombres de variables

# Estandarizar los nombres de las variables
# clean_names() convierte los nombres a minúsculas y utiliza "_" para separar palabras
data <- data %>%
  clean_names()

# Verificar los nombres de las variables luego de la estandarización
names(data)
##  [1] "low"   "age"   "lwt"   "race"  "smoke" "ptl"   "ht"    "ui"    "ftv"  
## [10] "bwt"

2.2 Revisión de duplicados

Primero se crea una variable id (colocada como primera columna) y luego se identifican observaciones completamente duplicadas.

# Para crear variable id y colocarla como primera variable
data <- data %>%
  mutate(id = row_number()) %>%
  relocate(id)

# Identificar observaciones completamente duplicadas
# get_dupes() muestra las filas que aparecen más de una vez en la base de datos
data %>%
  get_dupes(-id)

Si quisiéramos eliminar los duplicados, usamos la siguiente función:

# Eliminar observaciones completamente duplicadas
# Se conserva únicamente la primera observación de cada duplicado
data <- data %>%
  distinct(across(-id), .keep_all = TRUE)

# Verificar que no permanezcan observaciones completamente duplicadas
data %>%
  get_dupes(-id)

2.3 Valores faltantes

Para identificar los valores faltantes por columna:

# Identificar la cantidad de valores faltantes por variable
colSums(is.na(data))
##    id   low   age   lwt  race smoke   ptl    ht    ui   ftv   bwt 
##     0     0     0     0     0     0     0     0     0     0     0
# Calcular el porcentaje de valores faltantes por variable
round(colMeans(is.na(data)) * 100, 2)
##    id   low   age   lwt  race smoke   ptl    ht    ui   ftv   bwt 
##     0     0     0     0     0     0     0     0     0     0     0

Nota: para este ejemplo no hay valores faltantes en ninguna variable.

Comprobación general de valores faltantes:

# Calcular el total de valores faltantes en la base de datos
sum(is.na(data))
## [1] 0
# Identificar si existe al menos un valor faltante en la base de datos
anyNA(data)
## [1] FALSE

2.4 Tipo de variables y recodificación

Luego de evaluar los valores faltantes, se recomienda verificar que cada variable esté almacenada con el tipo de dato correspondiente. Las variables categóricas codificadas numéricamente deben identificarse y, cuando sea apropiado, convertirse a factores con etiquetas descriptivas antes de realizar los análisis estadísticos.

2.4.1 Conversión a factor

# Revisar los valores únicos de las variables categóricas
unique(data$low)
## [1] 0 1
unique(data$race)
## [1] 2 3 1
unique(data$smoke)
## [1] 0 1
unique(data$ht)
## [1] 0 1
unique(data$ui)
## [1] 1 0
# Convertir las variables categóricas a tipo factor
data$low   <- as.factor(data$low)
data$race  <- as.factor(data$race)
data$smoke <- as.factor(data$smoke)
data$ht    <- as.factor(data$ht)
data$ui    <- as.factor(data$ui)

# Verificar la estructura de las variables
str(data)
## 'data.frame':    184 obs. of  11 variables:
##  $ id   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ low  : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : Factor w/ 3 levels "1","2","3": 2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: Factor w/ 2 levels "0","1": 1 1 2 2 2 1 1 1 2 2 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ui   : Factor w/ 2 levels "0","1": 2 1 1 2 2 1 1 1 1 1 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...

2.4.2 Renombrar niveles de las categorías

# Renombrar los niveles de las variables categóricas

# Bajo peso al nacer
levels(data$low) <- c("No", "Yes")

# Raza
levels(data$race) <- c("White", "Black", "Other")

# Consumo de tabaco durante el embarazo
levels(data$smoke) <- c("No", "Yes")

# Historial de hipertensión
levels(data$ht) <- c("No", "Yes")

# Presencia de irritabilidad uterina
levels(data$ui) <- c("No", "Yes")

# Verificar los nuevos niveles
levels(data$low)
## [1] "No"  "Yes"
levels(data$race)
## [1] "White" "Black" "Other"
levels(data$smoke)
## [1] "No"  "Yes"
levels(data$ht)
## [1] "No"  "Yes"
levels(data$ui)
## [1] "No"  "Yes"

2.4.3 Verificación de la recodificación

# Verificar la estructura final de la base de datos
str(data)
## 'data.frame':    184 obs. of  11 variables:
##  $ id   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ low  : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : Factor w/ 3 levels "White","Black",..: 2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: Factor w/ 2 levels "No","Yes": 1 1 2 2 2 1 1 1 2 2 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ui   : Factor w/ 2 levels "No","Yes": 2 1 1 2 2 1 1 1 1 1 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
# Obtener un resumen de la base de datos luego de la recodificación
skim(data)
Data summary
Name data
Number of rows 184
Number of columns 11
_______________________
Column type frequency:
factor 5
numeric 6
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
low 0 1 FALSE 2 No: 125, Yes: 59
race 0 1 FALSE 3 Whi: 93, Oth: 66, Bla: 25
smoke 0 1 FALSE 2 No: 113, Yes: 71
ht 0 1 FALSE 2 No: 172, Yes: 12
ui 0 1 FALSE 2 No: 157, Yes: 27

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
id 0 1 96.21 54.74 1 49.75 96.5 143.25 189 ▇▇▇▇▇
age 0 1 23.39 5.29 14 19.75 23.0 26.25 45 ▇▇▃▁▁
lwt 0 1 130.25 30.71 80 110.00 121.5 140.25 250 ▅▇▂▁▁
ptl 0 1 0.20 0.50 0 0.00 0.0 0.00 3 ▇▁▁▁▁
ftv 0 1 0.81 1.07 0 0.00 0.0 1.00 6 ▇▂▁▁▁
bwt 0 1 2939.67 736.79 709 2410.00 2977.0 3501.25 4990 ▁▅▇▆▁

3 Análisis descriptivo

3.1 Variables categóricas

Luego de convertir y recodificar las variables categóricas, se recomienda verificar la distribución de sus categorías mediante tablas de frecuencia y porcentajes. Esto permite confirmar que los niveles fueron definidos correctamente e identificar posibles inconsistencias en los datos.

# Análisis descriptivo para variables categóricas
tabyl(data, low)
tabyl(data, race)
tabyl(data, smoke)
tabyl(data, ht)
tabyl(data, ui)

3.2 Variables numéricas

Para las variables numéricas, se recomienda evaluar medidas de tendencia central y dispersión, así como los valores mínimos, máximos y percentiles. Estas estadísticas permiten conocer la distribución general de los datos e identificar valores que requieran una evaluación adicional.

3.2.1 Con summary()

# Análisis descriptivo para variables numéricas (una por una)
summary(data$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   14.00   19.75   23.00   23.39   26.25   45.00
summary(data$lwt)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    80.0   110.0   121.5   130.2   140.2   250.0
summary(data$bwt)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     709    2410    2977    2940    3501    4990
# Análisis descriptivo para variables numéricas (todas juntas)
data %>%
  select(age, lwt, ptl, ftv, bwt) %>%
  summary()
##       age             lwt             ptl              ftv        
##  Min.   :14.00   Min.   : 80.0   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:19.75   1st Qu.:110.0   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :23.00   Median :121.5   Median :0.0000   Median :0.0000  
##  Mean   :23.39   Mean   :130.2   Mean   :0.2011   Mean   :0.8098  
##  3rd Qu.:26.25   3rd Qu.:140.2   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :45.00   Max.   :250.0   Max.   :3.0000   Max.   :6.0000  
##       bwt      
##  Min.   : 709  
##  1st Qu.:2410  
##  Median :2977  
##  Mean   :2940  
##  3rd Qu.:3501  
##  Max.   :4990

3.2.2 Con skim()

La función skim() permite obtener de manera sencilla estadísticas como la media, desviación estándar, mínimo, máximo y percentiles.

# Resumen descriptivo de las variables numéricas
data %>%
  select(age, lwt, ptl, ftv, bwt) %>%
  skim()
Data summary
Name Piped data
Number of rows 184
Number of columns 5
_______________________
Column type frequency:
numeric 5
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
age 0 1 23.39 5.29 14 19.75 23.0 26.25 45 ▇▇▃▁▁
lwt 0 1 130.25 30.71 80 110.00 121.5 140.25 250 ▅▇▂▁▁
ptl 0 1 0.20 0.50 0 0.00 0.0 0.00 3 ▇▁▁▁▁
ftv 0 1 0.81 1.07 0 0.00 0.0 1.00 6 ▇▂▁▁▁
bwt 0 1 2939.67 736.79 709 2410.00 2977.0 3501.25 4990 ▁▅▇▆▁

3.3 Valores atípicos (outliers)

Luego del análisis descriptivo, se recomienda evaluar la presencia de posibles valores atípicos (outliers) en las variables numéricas. Estos pueden identificarse mediante estadísticas descriptivas y herramientas gráficas como los boxplots.

Importante: la presencia de un valor atípico no implica necesariamente que sea un error, por lo que se recomienda investigar la observación antes de modificarla o eliminarla.

3.3.1 Exploración gráfica (boxplots)

# Crear boxplots para identificar posibles valores atípicos
boxplot(data$age,
        main = "Edad materna",
        ylab = "Edad")

boxplot(data$lwt,
        main = "Peso materno",
        ylab = "Peso")

boxplot(data$bwt,
        main = "Peso al nacer",
        ylab = "Peso")

3.3.2 Método del rango intercuartílico (IQR)

Una manera común de identificar posibles valores atípicos es mediante el rango intercuartílico (IQR). Se consideran potencialmente atípicos aquellos valores que se encuentran por debajo de Q1 − 1.5 × IQR o por encima de Q3 + 1.5 × IQR.

# Calcular cuartiles e IQR para la variable edad
Q1 <- quantile(data$age, 0.25, na.rm = TRUE)
Q3 <- quantile(data$age, 0.75, na.rm = TRUE)
IQR_age <- IQR(data$age, na.rm = TRUE)

# Calcular límites para identificar posibles outliers
limite_inferior <- Q1 - 1.5 * IQR_age
limite_superior <- Q3 + 1.5 * IQR_age

limite_inferior
## 25% 
##  10
limite_superior
## 75% 
##  36

Para identificar las observaciones que se presentan como posibles outliers:

data %>%
  filter(age < limite_inferior | age > limite_superior) %>%
  select(id, age)

3.3.3 Método directo con boxplot.stats()

Otra manera más directa de hacer esto:

# Identificar posibles valores atípicos de manera sencilla
boxplot.stats(data$age)$out
## [1] 45
boxplot.stats(data$lwt)$out
##  [1] 202 215 189 250 229 190 235 241 187 200 187 190
boxplot.stats(data$bwt)$out
## [1] 709
# Guardar los valores identificados como posibles outliers
outliers_bwt <- boxplot.stats(data$bwt)$out

# Identificar las observaciones correspondientes
data %>%
  filter(bwt %in% outliers_bwt) %>%
  select(id, bwt)

Las observaciones identificadas como posibles outliers deben ser evaluadas antes de tomar una decisión sobre su manejo. Se recomienda verificar si corresponden a errores de entrada, valores imposibles o valores extremos válidos. No se recomienda eliminar observaciones únicamente por haber sido identificadas como outliers.

3.4 Verificación de rangos plausibles

Además de evaluar posibles valores atípicos, se recomienda verificar que los valores de las variables se encuentren dentro de rangos razonables según el contexto de los datos. Un valor extremo no necesariamente representa un error; sin embargo, valores imposibles o fuera de los rangos esperados deben ser investigados antes de continuar con el análisis.

3.4.1 Valores mínimos y máximos

data %>%
  summarise(
    age_min = min(age, na.rm = TRUE),
    age_max = max(age, na.rm = TRUE),

    lwt_min = min(lwt, na.rm = TRUE),
    lwt_max = max(lwt, na.rm = TRUE),

    bwt_min = min(bwt, na.rm = TRUE),
    bwt_max = max(bwt, na.rm = TRUE)
  )

3.4.2 Observaciones fuera de un rango definido

Para investigar directamente observaciones que no cumplan un rango que tú hayas definido según el contexto del estudio:

# Ejemplo: identificar observaciones fuera de un rango esperado
data %>%
  filter(age < 10 | age > 60) %>%
  select(id, age)

Los límites utilizados para identificar valores fuera de rango deben establecerse utilizando el conocimiento de las variables, el diseño del estudio y, cuando esté disponible, el diccionario de datos. No deben establecerse arbitrariamente.

3.5 Verificación final de la limpieza

Una vez completado el proceso de limpieza, se recomienda realizar una última revisión de la estructura y distribución de los datos antes de comenzar los análisis estadísticos.

# Verificación general de la base de datos luego del proceso de limpieza
dim(data)
## [1] 184  11
str(data)
## 'data.frame':    184 obs. of  11 variables:
##  $ id   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ low  : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : Factor w/ 3 levels "White","Black",..: 2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: Factor w/ 2 levels "No","Yes": 1 1 2 2 2 1 1 1 2 2 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ui   : Factor w/ 2 levels "No","Yes": 2 1 1 2 2 1 1 1 1 1 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
skim(data)
Data summary
Name data
Number of rows 184
Number of columns 11
_______________________
Column type frequency:
factor 5
numeric 6
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
low 0 1 FALSE 2 No: 125, Yes: 59
race 0 1 FALSE 3 Whi: 93, Oth: 66, Bla: 25
smoke 0 1 FALSE 2 No: 113, Yes: 71
ht 0 1 FALSE 2 No: 172, Yes: 12
ui 0 1 FALSE 2 No: 157, Yes: 27

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
id 0 1 96.21 54.74 1 49.75 96.5 143.25 189 ▇▇▇▇▇
age 0 1 23.39 5.29 14 19.75 23.0 26.25 45 ▇▇▃▁▁
lwt 0 1 130.25 30.71 80 110.00 121.5 140.25 250 ▅▇▂▁▁
ptl 0 1 0.20 0.50 0 0.00 0.0 0.00 3 ▇▁▁▁▁
ftv 0 1 0.81 1.07 0 0.00 0.0 1.00 6 ▇▂▁▁▁
bwt 0 1 2939.67 736.79 709 2410.00 2977.0 3501.25 4990 ▁▅▇▆▁

4 Análisis bivariado

Una vez completado el proceso de limpieza y análisis descriptivo de los datos, se pueden realizar análisis bivariados para evaluar la relación o asociación entre dos variables. La prueba estadística utilizada dependerá del tipo de variables que se estén analizando y de los supuestos correspondientes.

En esta sección se presentan ejemplos para evaluar la relación entre:

Tipo de variables Prueba(s)
Categórica × categórica Chi-cuadrado / Fisher
Categórica × numérica Prueba t / Wilcoxon
Numérica × numérica Correlación de Pearson / Spearman

4.1 Categórica × categórica: Chi-cuadrado / Fisher

Para evaluar la asociación entre dos variables categóricas se puede utilizar la prueba de Chi-cuadrado de independencia. Antes de realizar la prueba, se recomienda construir una tabla de contingencia para observar la distribución de las variables. Si las frecuencias esperadas son demasiado pequeñas para cumplir con los supuestos de Chi-cuadrado, se puede considerar la prueba exacta de Fisher.

4.1.1 Tabla de contingencia

# Crear tabla de contingencia entre dos variables categóricas
tabyl(data, low, smoke)
# Crear tabla de contingencia con porcentajes por fila
data %>%
  tabyl(low, smoke) %>%
  adorn_percentages("row") %>%
  adorn_pct_formatting(digits = 1)

4.1.2 Prueba de Chi-cuadrado

# Realizar prueba de Chi-cuadrado de independencia
chi_result <- chisq.test(data$low, data$smoke)

# Mostrar resultados
chi_result
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data$low and data$smoke
## X-squared = 4.7738, df = 1, p-value = 0.0289

4.1.3 Verificación de supuestos

Para verificar el supuesto de frecuencias esperadas por celda:

# Verificar las frecuencias esperadas de la prueba de Chi-cuadrado
chi_result$expected
##         data$smoke
## data$low      No     Yes
##      No  76.7663 48.2337
##      Yes 36.2337 22.7663

Importante:

  • Ninguna celda debería tener una frecuencia esperada < 1.
  • No más del 20% de las celdas deberían tener frecuencias esperadas < 5.

Si no se cumple con el supuesto, se recomienda utilizar la prueba exacta de Fisher.

4.1.4 Prueba exacta de Fisher

# Realizar prueba exacta de Fisher
fisher_result <- fisher.test(data$low, data$smoke)

# Mostrar resultados
fisher_result
## 
##  Fisher's Exact Test for Count Data
## 
## data:  data$low and data$smoke
## p-value = 0.02328
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.071142 4.184248
## sample estimates:
## odds ratio 
##   2.110414

4.2 Categórica × numérica: prueba t / Wilcoxon

Para evaluar diferencias en una variable numérica entre dos grupos independientes, se puede utilizar la prueba t de Student. Antes de realizar la prueba, se recomienda evaluar la distribución de la variable numérica y los supuestos correspondientes, incluyendo normalidad y homogeneidad de varianzas.

4.2.1 Exploración descriptiva y gráfica

# Estadísticas descriptivas de la variable numérica por grupo
data %>%
  group_by(smoke) %>%
  summarise(
    n = n(),
    media = mean(bwt, na.rm = TRUE),
    sd = sd(bwt, na.rm = TRUE),
    mediana = median(bwt, na.rm = TRUE),
    min = min(bwt, na.rm = TRUE),
    max = max(bwt, na.rm = TRUE)
  )
# Evaluar visualmente la distribución de la variable numérica por grupo
boxplot(bwt ~ smoke,
        data = data,
        main = "Peso al nacer según consumo de tabaco",
        xlab = "Consumo de tabaco",
        ylab = "Peso al nacer")

4.2.2 Normalidad

La normalidad de la variable numérica dentro de cada grupo puede evaluarse mediante la prueba de Shapiro-Wilk. Un valor p > 0.05 indica que no existe evidencia estadísticamente significativa para rechazar el supuesto de normalidad.

# Evaluar normalidad dentro de cada grupo
shapiro.test(data$bwt[data$smoke == "No"])
## 
##  Shapiro-Wilk normality test
## 
## data:  data$bwt[data$smoke == "No"]
## W = 0.98626, p-value = 0.3049
shapiro.test(data$bwt[data$smoke == "Yes"])
## 
##  Shapiro-Wilk normality test
## 
## data:  data$bwt[data$smoke == "Yes"]
## W = 0.98301, p-value = 0.4515

La prueba de normalidad no debe evaluarse de manera aislada. También se recomienda considerar la distribución de los datos, los valores atípicos y el tamaño de los grupos.

4.2.3 Homogeneidad de varianzas

Para evaluar si las varianzas de los dos grupos pueden considerarse iguales, se puede utilizar una prueba de igualdad de varianzas.

# Evaluar igualdad de varianzas entre los dos grupos
var_test_result <- var.test(bwt ~ smoke,
                            data = data)

# Mostrar resultados
var_test_result
## 
##  F test to compare two variances
## 
## data:  bwt by smoke
## F = 1.3058, num df = 112, denom df = 70, p-value = 0.2286
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.8445385 1.9768528
## sample estimates:
## ratio of variances 
##           1.305806

4.2.4 Prueba t de Student (varianzas iguales)

Si el valor p es mayor de 0.05, no se rechaza la hipótesis de igualdad de varianzas y se puede realizar la prueba t asumiendo varianzas iguales.

# Prueba t asumiendo igualdad de varianzas
t_test_equal <- t.test(bwt ~ smoke,
                       data = data,
                       var.equal = TRUE)

t_test_equal
## 
##  Two Sample t-test
## 
## data:  bwt by smoke
## t = 2.735, df = 182, p-value = 0.006855
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
##   83.54809 516.26246
## sample estimates:
##  mean in group No mean in group Yes 
##          3055.398          2755.493

4.2.5 Prueba t de Welch (varianzas desiguales)

Si el valor p es menor o igual a 0.05, existe evidencia de que las varianzas son diferentes. En este caso se puede utilizar la prueba t de Welch, que no asume igualdad de varianzas.

# Prueba t de Welch para varianzas desiguales
t_test_welch <- t.test(bwt ~ smoke,
                       data = data,
                       var.equal = FALSE)

t_test_welch
## 
##  Welch Two Sample t-test
## 
## data:  bwt by smoke
## t = 2.8196, df = 163.29, p-value = 0.005405
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
##   89.87443 509.93611
## sample estimates:
##  mean in group No mean in group Yes 
##          3055.398          2755.493

4.2.6 Alternativa no paramétrica: Wilcoxon

Cuando los supuestos de una prueba paramétrica no sean apropiados, se puede considerar la prueba de Wilcoxon como alternativa no paramétrica.

# Prueba de Wilcoxon para dos grupos independientes
wilcox_result <- wilcox.test(bwt ~ smoke,
                             data = data)

wilcox_result
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  bwt by smoke
## W = 4990.5, p-value = 0.005396
## alternative hypothesis: true location shift is not equal to 0

4.3 Numérica × numérica: correlación de Pearson / Spearman

Para evaluar la relación entre dos variables numéricas se puede utilizar un análisis de correlación. La selección del coeficiente de correlación dependerá de las características y distribución de los datos. Entre los métodos más utilizados se encuentran la correlación de Pearson y la correlación de Spearman.

Para este ejemplo se evaluará la relación entre la edad de la madre (age) y su peso (lwt).

4.3.1 Exploración gráfica

Antes de realizar el análisis de correlación, se recomienda visualizar la relación entre las variables mediante un diagrama de dispersión (scatterplot). Esto permite evaluar la dirección y forma de la relación e identificar posibles valores atípicos.

# Diagrama de dispersión entre las dos variables numéricas
plot(data$age, data$lwt,
     xlab = "Edad materna",
     ylab = "Peso materno",
     main = "Relación entre edad y peso materno")

4.3.2 Evaluación de normalidad

Para la correlación de Pearson se recomienda evaluar la distribución de las variables y la presencia de valores atípicos. La prueba de Shapiro-Wilk puede utilizarse como una herramienta adicional para evaluar normalidad.

# Evaluar normalidad de las variables
shapiro.test(data$age)
## 
##  Shapiro-Wilk normality test
## 
## data:  data$age
## W = 0.96335, p-value = 9.726e-05
shapiro.test(data$lwt)
## 
##  Shapiro-Wilk normality test
## 
## data:  data$lwt
## W = 0.89304, p-value = 3.201e-10

Un valor p > 0.05 indica que no existe evidencia estadísticamente significativa para rechazar el supuesto de normalidad. Sin embargo, esta prueba no debe utilizarse de manera aislada y se recomienda considerar también las características de la distribución y la presencia de valores atípicos.

4.3.3 Correlación de Pearson

La correlación de Pearson se utiliza para evaluar la fuerza y dirección de una relación lineal entre dos variables numéricas. El coeficiente de correlación (r) puede tomar valores entre -1 y 1.

# Correlación de Pearson
pearson_result <- cor.test(data$age,
                           data$lwt,
                           method = "pearson")

# Mostrar resultados
pearson_result
## 
##  Pearson's product-moment correlation
## 
## data:  data$age and data$lwt
## t = 2.3055, df = 182, p-value = 0.02227
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.02438334 0.30566275
## sample estimates:
##       cor 
## 0.1684502

4.3.4 Matriz de correlación

# Seleccionar variables numéricas
variables_cor <- data %>%
  select(age, lwt, ptl, ftv, bwt)

# Matriz de correlación
cor(variables_cor,
    use = "complete.obs",
    method = "pearson")
##            age        lwt         ptl         ftv         bwt
## age 1.00000000  0.1684502  0.06151932  0.20471278  0.09815339
## lwt 0.16845018  1.0000000 -0.14736982  0.13470219  0.18570067
## ptl 0.06151932 -0.1473698  1.00000000 -0.05092996 -0.15277343
## ftv 0.20471278  0.1347022 -0.05092996  1.00000000  0.06412657
## bwt 0.09815339  0.1857007 -0.15277343  0.06412657  1.00000000

4.3.5 Correlación de Spearman

Cuando los supuestos para la correlación de Pearson no sean apropiados, o cuando se desee evaluar una relación monotónica utilizando los rangos de las observaciones, se puede utilizar la correlación de Spearman.

# Correlación de Spearman
spearman_result <- cor.test(data$age,
                            data$lwt,
                            method = "spearman",
                            exact = FALSE)

# Mostrar resultados
spearman_result
## 
##  Spearman's rank correlation rho
## 
## data:  data$age and data$lwt
## S = 861914, p-value = 0.02119
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.169816

4.3.6 Gráfico de matriz de correlaciones (ggpairs)

# Seleccionar variables numéricas
variables_cor <- data %>%
  select(age, lwt, bwt)

# Crear matriz de correlaciones
ggpairs(variables_cor)

5 Análisis multivariado

El análisis multivariado permite evaluar simultáneamente la relación entre una variable dependiente y múltiples variables independientes. Estos modelos permiten estimar la asociación de cada variable independiente con el outcome mientras se controla estadísticamente por las demás variables incluidas en el modelo.

La selección del modelo dependerá principalmente del tipo de variable dependiente, el diseño del estudio, el objetivo del análisis y los supuestos correspondientes.

5.1 Regresión logística binaria

La regresión logística binaria se utiliza cuando la variable dependiente presenta dos categorías. Este modelo permite evaluar la asociación entre múltiples variables independientes y la probabilidad de que ocurra el evento de interés.

En este ejemplo se utilizará la variable low como variable dependiente, la cual indica la presencia o ausencia de bajo peso al nacer.

5.1.1 Verificación de la variable dependiente

Antes de realizar el modelo, se recomienda verificar las categorías de la variable dependiente y confirmar cuál será utilizada como categoría de referencia.

# Verificar niveles de la variable dependiente
levels(data$low)
## [1] "No"  "Yes"
# Verificar distribución de la variable dependiente
tabyl(data, low)

En este caso, la variable low contiene las categorías "No" y "Yes". La categoría "No" se encuentra como primera categoría y "Yes" representa el evento de interés: bajo peso al nacer.

5.1.2 Categoría de referencia

Cuando sea necesario, se puede establecer explícitamente la categoría de referencia utilizando relevel().

# Establecer "No" como categoría de referencia
data$low <- relevel(data$low, ref = "No")

# Verificar niveles
levels(data$low)
## [1] "No"  "Yes"

Esto permite interpretar el modelo en términos de la ocurrencia del evento "Yes" en comparación con la categoría de referencia "No".

5.1.3 Ajuste del modelo

Para crear el modelo se utiliza la función glm() con la familia binomial y el enlace logit.

# Crear modelo de regresión logística binaria
modelo_logistico <- glm(
  low ~ age + lwt + race + smoke + ht + ui,
  data = data,
  family = binomial(link = "logit")
)

# Mostrar resultados del modelo
summary(modelo_logistico)
## 
## Call:
## glm(formula = low ~ age + lwt + race + smoke + ht + ui, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.951311   1.229366   0.774  0.43904   
## age         -0.030047   0.035940  -0.836  0.40314   
## lwt         -0.017642   0.007019  -2.514  0.01195 * 
## raceBlack    1.281103   0.534062   2.399  0.01645 * 
## raceOther    0.817943   0.437370   1.870  0.06146 . 
## smokeYes     1.042921   0.396316   2.632  0.00850 **
## htYes        1.848820   0.692818   2.669  0.00762 **
## uiYes        0.943965   0.457588   2.063  0.03912 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 230.87  on 183  degrees of freedom
## Residual deviance: 198.46  on 176  degrees of freedom
## AIC: 214.46
## 
## Number of Fisher Scoring iterations: 4

El resultado de summary() presenta los coeficientes del modelo en escala logarítmica (log-odds), junto con sus errores estándar, estadísticos de prueba y valores p.

5.1.4 Odds Ratios (OR)

Para facilitar la interpretación de los coeficientes de una regresión logística, estos pueden transformarse mediante la función exponencial para obtener Odds Ratios (OR).

# Calcular Odds Ratios
exp(coef(modelo_logistico))
## (Intercept)         age         lwt   raceBlack   raceOther    smokeYes 
##   2.5891029   0.9704002   0.9825130   3.6006075   2.2658344   2.8374937 
##       htYes       uiYes 
##   6.3523223   2.5701524

Interpretación general:

  • OR = 1: no se observa cambio en los odds del outcome.
  • OR > 1: mayores odds del outcome.
  • OR < 1: menores odds del outcome.

La interpretación específica dependerá del tipo de variable independiente y de su categoría de referencia.

5.1.5 Odds Ratios e intervalos de confianza

Además de los Odds Ratios, se recomienda calcular sus intervalos de confianza del 95%.

# Calcular Odds Ratios e intervalos de confianza del 95%
OR_IC <- exp(
  cbind(
    OR = coef(modelo_logistico),
    confint(modelo_logistico)
  )
)

# Mostrar resultados
OR_IC
##                    OR     2.5 %     97.5 %
## (Intercept) 2.5891029 0.2419222 30.6159150
## age         0.9704002 0.9030621  1.0403580
## lwt         0.9825130 0.9682452  0.9954047
## raceBlack   3.6006075 1.2706104 10.4957734
## raceOther   2.2658344 0.9724586  5.4573304
## smokeYes    2.8374937 1.3237817  6.3195287
## htYes       6.3523223 1.6939674 27.0431266
## uiYes       2.5701524 1.0449218  6.3699230

De manera general, si el intervalo de confianza del 95% de un Odds Ratio incluye el valor 1, no existe evidencia estadísticamente significativa de una asociación al nivel de significancia de 0.05.

5.1.6 Probabilidades predichas

El modelo también puede utilizarse para calcular la probabilidad estimada del outcome para cada observación.

# Calcular probabilidades predichas
data$probabilidad_predicha <- predict(
  modelo_logistico,
  type = "response"
)

# Visualizar algunas probabilidades
head(
  data %>%
    select(id, low, probabilidad_predicha)
)

Las probabilidades predichas toman valores entre 0 y 1 y representan la probabilidad estimada de presentar el evento de interés según las variables incluidas en el modelo.

5.1.7 Evaluación general del modelo

Finalmente, se recomienda revisar el modelo de manera integral, considerando los coeficientes estimados, Odds Ratios, intervalos de confianza, significancia estadística, posibles problemas de multicolinealidad y la plausibilidad de los resultados.

# Resumen final del modelo
summary(modelo_logistico)
## 
## Call:
## glm(formula = low ~ age + lwt + race + smoke + ht + ui, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.951311   1.229366   0.774  0.43904   
## age         -0.030047   0.035940  -0.836  0.40314   
## lwt         -0.017642   0.007019  -2.514  0.01195 * 
## raceBlack    1.281103   0.534062   2.399  0.01645 * 
## raceOther    0.817943   0.437370   1.870  0.06146 . 
## smokeYes     1.042921   0.396316   2.632  0.00850 **
## htYes        1.848820   0.692818   2.669  0.00762 **
## uiYes        0.943965   0.457588   2.063  0.03912 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 230.87  on 183  degrees of freedom
## Residual deviance: 198.46  on 176  degrees of freedom
## AIC: 214.46
## 
## Number of Fisher Scoring iterations: 4
# Odds Ratios e intervalos de confianza
OR_IC
##                    OR     2.5 %     97.5 %
## (Intercept) 2.5891029 0.2419222 30.6159150
## age         0.9704002 0.9030621  1.0403580
## lwt         0.9825130 0.9682452  0.9954047
## raceBlack   3.6006075 1.2706104 10.4957734
## raceOther   2.2658344 0.9724586  5.4573304
## smokeYes    2.8374937 1.3237817  6.3195287
## htYes       6.3523223 1.6939674 27.0431266
## uiYes       2.5701524 1.0449218  6.3699230

La interpretación de los resultados debe realizarse considerando el diseño del estudio, las características de las variables, los posibles factores de confusión y el contexto de la investigación.

5.2 Evaluación adicional del modelo logístico

Luego de ajustar el modelo de regresión logística, se recomienda evaluar diferentes aspectos del modelo antes de realizar la interpretación final. Entre estos se encuentran la multicolinealidad, posibles observaciones influyentes y el ajuste general del modelo.

5.2.1 Multicolinealidad

La multicolinealidad ocurre cuando dos o más variables independientes presentan una relación fuerte entre sí. Esto puede aumentar la variabilidad de los coeficientes estimados y dificultar la interpretación de los efectos individuales de las variables.

El Variance Inflation Factor (VIF) puede utilizarse para evaluar la presencia de multicolinealidad.

# Evaluar multicolinealidad
vif(modelo_logistico)
##           GVIF Df GVIF^(1/(2*Df))
## age   1.030894  1        1.015330
## lwt   1.288174  1        1.134978
## race  1.489296  2        1.104702
## smoke 1.293516  1        1.137328
## ht    1.154384  1        1.074423
## ui    1.023924  1        1.011891

Valores elevados de VIF pueden indicar problemas de multicolinealidad. No se recomienda utilizar un punto de corte de manera automática; los resultados deben evaluarse considerando las características del modelo y las variables incluidas.

5.2.2 Observaciones influyentes

Algunas observaciones pueden tener una influencia considerable sobre los coeficientes estimados del modelo. Una herramienta que puede utilizarse para identificar posibles observaciones influyentes es la distancia de Cook.

# Calcular distancia de Cook
cook_logistico <- cooks.distance(modelo_logistico)

# Visualizar distancia de Cook
plot(cook_logistico,
     type = "h",
     main = "Distancia de Cook",
     xlab = "Observación",
     ylab = "Distancia de Cook")

También se pueden identificar las observaciones con los valores más altos:

# Identificar observaciones con mayor distancia de Cook
head(
  sort(cook_logistico, decreasing = TRUE),
  10
)
##         98         28         77         10        144        138         59 
## 0.05579945 0.04593139 0.03396932 0.02778999 0.02567499 0.02447404 0.02372116 
##         57        119        197 
## 0.02289425 0.02228267 0.02169551

Una observación identificada como potencialmente influyente no debe eliminarse automáticamente. Primero se recomienda verificar la observación, determinar si existe algún error en los datos y evaluar su posible influencia sobre los resultados.

5.2.3 Ajuste mediante AIC

El Akaike Information Criterion (AIC) puede utilizarse para comparar modelos ajustados utilizando los mismos datos y la misma variable dependiente.

# Obtener AIC del modelo
AIC(modelo_logistico)
## [1] 214.4631

Valores menores de AIC indican un mejor balance relativo entre ajuste y complejidad al comparar modelos. El valor de AIC no debe interpretarse de manera aislada, sino en comparación con otros modelos apropiados.

5.3 Modelos crudos y ajustados

En algunos análisis puede ser útil comenzar evaluando la asociación entre el outcome y una variable independiente mediante un modelo crudo, y posteriormente construir un modelo ajustado que incluya otras covariables relevantes.

5.3.1 Modelo crudo

Como ejemplo, se puede evaluar la asociación entre el consumo de tabaco durante el embarazo (smoke) y el bajo peso al nacer (low).

# Modelo crudo
modelo_crudo <- glm(
  low ~ smoke,
  data = data,
  family = binomial(link = "logit")
)

# Mostrar resultados
summary(modelo_crudo)
## 
## Call:
## glm(formula = low ~ smoke, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0635     0.2154  -4.938  7.9e-07 ***
## smokeYes      0.7511     0.3227   2.328   0.0199 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 230.87  on 183  degrees of freedom
## Residual deviance: 225.43  on 182  degrees of freedom
## AIC: 229.43
## 
## Number of Fisher Scoring iterations: 4

Los Odds Ratios y sus intervalos de confianza pueden obtenerse de la siguiente manera:

# Odds Ratios e intervalos de confianza del modelo crudo
exp(
  cbind(
    OR = coef(modelo_crudo),
    confint(modelo_crudo)
  )
)
##                    OR     2.5 %    97.5 %
## (Intercept) 0.3452381 0.2227119 0.5198449
## smokeYes    2.1194281 1.1277800 4.0090138

5.3.2 Modelo ajustado

Posteriormente se pueden incluir otras variables relevantes para estimar la asociación entre smoke y low controlando estadísticamente por las demás covariables incluidas en el modelo.

# Modelo ajustado
modelo_ajustado <- glm(
  low ~ smoke + age + lwt + race + ht + ui,
  data = data,
  family = binomial(link = "logit")
)

# Mostrar resultados
summary(modelo_ajustado)
## 
## Call:
## glm(formula = low ~ smoke + age + lwt + race + ht + ui, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.951311   1.229366   0.774  0.43904   
## smokeYes     1.042921   0.396316   2.632  0.00850 **
## age         -0.030047   0.035940  -0.836  0.40314   
## lwt         -0.017642   0.007019  -2.514  0.01195 * 
## raceBlack    1.281103   0.534062   2.399  0.01645 * 
## raceOther    0.817943   0.437370   1.870  0.06146 . 
## htYes        1.848820   0.692818   2.669  0.00762 **
## uiYes        0.943965   0.457588   2.063  0.03912 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 230.87  on 183  degrees of freedom
## Residual deviance: 198.46  on 176  degrees of freedom
## AIC: 214.46
## 
## Number of Fisher Scoring iterations: 4
# Odds Ratios ajustados e intervalos de confianza
OR_ajustado <- exp(
  cbind(
    OR = coef(modelo_ajustado),
    confint(modelo_ajustado)
  )
)

OR_ajustado
##                    OR     2.5 %     97.5 %
## (Intercept) 2.5891029 0.2419222 30.6159150
## smokeYes    2.8374937 1.3237817  6.3195287
## age         0.9704002 0.9030621  1.0403580
## lwt         0.9825130 0.9682452  0.9954047
## raceBlack   3.6006075 1.2706104 10.4957734
## raceOther   2.2658344 0.9724586  5.4573304
## htYes       6.3523223 1.6939674 27.0431266
## uiYes       2.5701524 1.0449218  6.3699230

El Odds Ratio asociado con smoke en el primer modelo corresponde a un Odds Ratio crudo, mientras que el obtenido en el segundo modelo corresponde a un Odds Ratio ajustado por las demás variables incluidas.

5.4 Evaluación de interacción

En algunos estudios puede existir modificación de efecto o interacción, donde la asociación entre una variable independiente y el outcome cambia dependiendo del nivel de otra variable.

Una interacción puede incorporarse al modelo utilizando el operador *. Como ejemplo ilustrativo, se puede evaluar una posible interacción entre smoke y race.

# Modelo con término de interacción
modelo_interaccion <- glm(
  low ~ smoke * race + age + lwt + ht + ui,
  data = data,
  family = binomial(link = "logit")
)

# Mostrar resultados
summary(modelo_interaccion)
## 
## Call:
## glm(formula = low ~ smoke * race + age + lwt + ht + ui, family = binomial(link = "logit"), 
##     data = data)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)  
## (Intercept)         0.496429   1.366458   0.363   0.7164  
## smokeYes            1.531586   0.628371   2.437   0.0148 *
## raceBlack           1.475502   0.818642   1.802   0.0715 .
## raceOther           1.341786   0.628884   2.134   0.0329 *
## age                -0.028933   0.037216  -0.777   0.4369  
## lwt                -0.017113   0.007155  -2.392   0.0168 *
## htYes               1.779843   0.702125   2.535   0.0112 *
## uiYes               1.022735   0.465341   2.198   0.0280 *
## smokeYes:raceBlack -0.115570   1.109892  -0.104   0.9171  
## smokeYes:raceOther -1.400067   0.944329  -1.483   0.1382  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 230.87  on 183  degrees of freedom
## Residual deviance: 195.96  on 174  degrees of freedom
## AIC: 215.96
## 
## Number of Fisher Scoring iterations: 5

El operador * incluye automáticamente los efectos principales de ambas variables y sus términos de interacción.

La presencia e interpretación de una interacción debe evaluarse considerando los coeficientes correspondientes, la incertidumbre de las estimaciones, el contexto de la investigación y la plausibilidad de la interacción.

5.5 Comparación de modelos

Cuando los modelos son apropiadamente comparables, se pueden utilizar diferentes herramientas para evaluar si la incorporación de variables adicionales mejora el ajuste.

5.5.1 Comparación mediante AIC

# Comparar AIC de los modelos
AIC(
  modelo_crudo,
  modelo_ajustado,
  modelo_interaccion
)

5.5.2 Prueba de razón de verosimilitudes

También se pueden comparar modelos anidados mediante una prueba de razón de verosimilitudes.

# Comparar modelos anidados
anova(
  modelo_ajustado,
  modelo_interaccion,
  test = "Chisq"
)

Esta comparación permite evaluar si la incorporación de los términos adicionales produce una mejora estadísticamente detectable en el ajuste del modelo.

5.6 Probabilidades predichas (modelo ajustado)

Además de los Odds Ratios, la regresión logística permite obtener probabilidades predichas para cada observación.

# Obtener probabilidades predichas del modelo ajustado
data$probabilidad_predicha <- predict(
  modelo_ajustado,
  type = "response"
)

# Mostrar algunas observaciones
head(
  data %>%
    select(id, low, probabilidad_predicha)
)

Las probabilidades predichas se encuentran entre 0 y 1 y representan la probabilidad estimada del evento de interés según las características incluidas en el modelo.

5.7 Regresión lineal múltiple

La regresión lineal se utiliza cuando la variable dependiente es numérica continua. Permite evaluar la relación entre una variable dependiente y una o más variables independientes.

Cuando se incluye una sola variable independiente se denomina regresión lineal simple, mientras que cuando se incluyen múltiples variables independientes se denomina regresión lineal múltiple.

En este ejemplo se utilizará el peso al nacer (bwt) como variable dependiente.

5.7.1 Regresión lineal simple

Inicialmente se puede evaluar la relación entre el peso materno (lwt) y el peso al nacer (bwt).

# Crear modelo de regresión lineal simple
modelo_lineal_simple <- lm(
  bwt ~ lwt,
  data = data
)

# Mostrar resultados
summary(modelo_lineal_simple)
## 
## Call:
## lm(formula = bwt ~ lwt, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2185.01  -496.60    -1.66   515.97  2082.63 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2359.384    233.809   10.09   <2e-16 ***
## lwt            4.455      1.747    2.55   0.0116 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 726 on 182 degrees of freedom
## Multiple R-squared:  0.03448,    Adjusted R-squared:  0.02918 
## F-statistic:   6.5 on 1 and 182 DF,  p-value: 0.01161

El coeficiente asociado con lwt representa el cambio promedio estimado en bwt por cada aumento de una unidad en lwt.

5.7.2 Visualización de la relación

Antes de interpretar un modelo lineal, es útil visualizar la relación entre las variables.

# Scatterplot
plot(
  data$lwt,
  data$bwt,
  xlab = "Peso materno",
  ylab = "Peso al nacer",
  main = "Relación entre peso materno y peso al nacer"
)

# Añadir línea de regresión
abline(modelo_lineal_simple)

5.7.3 Regresión lineal múltiple

Posteriormente se pueden incorporar múltiples variables independientes para evaluar sus asociaciones con el peso al nacer mientras se controla por las demás variables incluidas.

# Crear modelo de regresión lineal múltiple
modelo_lineal_multiple <- lm(
  bwt ~ age + lwt + race + smoke + ht + ui,
  data = data
)

# Mostrar resultados
summary(modelo_lineal_multiple)
## 
## Call:
## lm(formula = bwt ~ age + lwt + race + smoke + ht + ui, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1801.42  -446.15    53.08   452.03  1702.30 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2890.659    319.683   9.042 2.76e-16 ***
## age           -3.271      9.537  -0.343 0.732026    
## lwt            4.425      1.737   2.548 0.011702 *  
## raceBlack   -479.626    153.045  -3.134 0.002021 ** 
## raceOther   -346.928    115.781  -2.996 0.003126 ** 
## smokeYes    -369.183    105.835  -3.488 0.000614 ***
## htYes       -586.380    202.181  -2.900 0.004204 ** 
## uiYes       -548.592    138.413  -3.963 0.000107 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 652.6 on 176 degrees of freedom
## Multiple R-squared:  0.2455, Adjusted R-squared:  0.2155 
## F-statistic: 8.182 on 7 and 176 DF,  p-value: 1.28e-08

5.7.4 Intervalos de confianza

Se pueden calcular intervalos de confianza del 95% para los coeficientes del modelo.

# Intervalos de confianza del 95%
confint(
  modelo_lineal_multiple,
  level = 0.95
)
##                    2.5 %      97.5 %
## (Intercept) 2259.7540563 3521.563345
## age          -22.0931048   15.551034
## lwt            0.9970722    7.852336
## raceBlack   -781.6667640 -177.585718
## raceOther   -575.4247291 -118.430416
## smokeYes    -578.0524202 -160.313261
## htYes       -985.3918150 -187.367689
## uiYes       -821.7542332 -275.429256

5.7.5 Multicolinealidad

# Evaluar multicolinealidad
vif(modelo_lineal_multiple)
##           GVIF Df GVIF^(1/(2*Df))
## age   1.094327  1        1.046101
## lwt   1.222509  1        1.105671
## race  1.324136  2        1.072712
## smoke 1.146849  1        1.070910
## ht    1.076715  1        1.037649
## ui    1.036393  1        1.018034

5.7.6 Supuestos del modelo lineal

Entre los principales aspectos que deben evaluarse en una regresión lineal se encuentran la linealidad, independencia de los errores, homocedasticidad, distribución aproximadamente normal de los residuos y presencia de observaciones influyentes.

Los gráficos diagnósticos de R permiten realizar una evaluación inicial de varios de estos aspectos.

# Gráficos diagnósticos del modelo
par(mfrow = c(2, 2))

plot(modelo_lineal_multiple)

par(mfrow = c(1, 1))

Los gráficos permiten evaluar:

  1. Residuals vs Fitted: ayuda a evaluar linealidad y patrones en los residuos.
  2. Normal Q-Q: permite evaluar si la distribución de los residuos se aproxima razonablemente a una distribución normal.
  3. Scale-Location: ayuda a evaluar la homocedasticidad o constancia de la variabilidad de los residuos.
  4. Residuals vs Leverage: ayuda a identificar observaciones potencialmente influyentes.

Normalidad de los residuos

La normalidad relevante en la regresión lineal se refiere a los residuos del modelo, no necesariamente a que cada variable independiente tenga una distribución normal.

# Extraer residuos
residuos <- residuals(modelo_lineal_multiple)

# Histograma de residuos
hist(
  residuos,
  main = "Distribución de los residuos",
  xlab = "Residuos"
)

# Q-Q plot
qqnorm(residuos)
qqline(residuos)

# Shapiro-Wilk como herramienta complementaria
shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.99516, p-value = 0.8197

La prueba de Shapiro-Wilk puede utilizarse como herramienta complementaria, pero no debe ser el único criterio para determinar si el supuesto es razonable.

Homocedasticidad

La homocedasticidad implica que la variabilidad de los residuos sea aproximadamente constante a través de los valores predichos.

Además de examinar los gráficos Residuals vs Fitted y Scale-Location, se puede utilizar una prueba formal como Breusch-Pagan.

# Prueba de Breusch-Pagan
bptest(modelo_lineal_multiple)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_lineal_multiple
## BP = 9.2269, df = 7, p-value = 0.2368

Un valor p > 0.05 indica que no existe evidencia estadísticamente significativa para rechazar la hipótesis de homocedasticidad. La prueba debe interpretarse junto con los gráficos diagnósticos.

5.7.7 Observaciones influyentes

# Calcular distancia de Cook
cook_lineal <- cooks.distance(modelo_lineal_multiple)

# Visualizar
plot(
  cook_lineal,
  type = "h",
  main = "Distancia de Cook",
  xlab = "Observación",
  ylab = "Distancia de Cook"
)

# Mostrar observaciones con valores más altos
head(
  sort(cook_lineal, decreasing = TRUE),
  10
)
##        226         11         10        202        197          4         13 
## 0.11423822 0.06768312 0.06263714 0.05588187 0.04588598 0.04315898 0.04210075 
##        210        187        188 
## 0.03304067 0.03158834 0.03098419

Las observaciones potencialmente influyentes deben investigarse antes de considerar cualquier modificación o exclusión de los datos.

5.7.8 Evaluación general del modelo

# Resumen del modelo
summary(modelo_lineal_multiple)
## 
## Call:
## lm(formula = bwt ~ age + lwt + race + smoke + ht + ui, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1801.42  -446.15    53.08   452.03  1702.30 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2890.659    319.683   9.042 2.76e-16 ***
## age           -3.271      9.537  -0.343 0.732026    
## lwt            4.425      1.737   2.548 0.011702 *  
## raceBlack   -479.626    153.045  -3.134 0.002021 ** 
## raceOther   -346.928    115.781  -2.996 0.003126 ** 
## smokeYes    -369.183    105.835  -3.488 0.000614 ***
## htYes       -586.380    202.181  -2.900 0.004204 ** 
## uiYes       -548.592    138.413  -3.963 0.000107 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 652.6 on 176 degrees of freedom
## Multiple R-squared:  0.2455, Adjusted R-squared:  0.2155 
## F-statistic: 8.182 on 7 and 176 DF,  p-value: 1.28e-08

El resultado permite evaluar los coeficientes estimados, errores estándar, valores p, R² y R² ajustado.

El R² representa la proporción de variabilidad observada en la variable dependiente que es explicada por las variables incluidas en el modelo. El R² ajustado incorpora una penalización por el número de predictores y resulta particularmente útil al evaluar modelos con múltiples variables independientes.

Conclusión

Este documento presenta un flujo de trabajo estructurado y reproducible para el manejo y análisis estadístico de datos en R, incluyendo procesos de limpieza y control de calidad, análisis descriptivos, análisis bivariados y modelos multivariables.

La selección de cada método estadístico debe realizarse considerando el tipo de variables, los objetivos del análisis y los supuestos correspondientes. Además, la interpretación de los resultados debe considerar no solo la significancia estadística, sino también la magnitud y dirección de las asociaciones, los intervalos de confianza y el contexto del estudio.

La integración del código, análisis y documentación en un mismo flujo de trabajo facilita la reproducibilidad, transparencia y reutilización del proceso analítico en futuros proyectos.