Regresión Logística en R: Tutorial Completo con Código (glm)

La regresion logistica es una de las tecnicas estadisticas mas importantes en psicologia clinica e investigacion en salud. A diferencia de la regresion multiple paso a paso, que predice valores continuos, la regresion logistica nos permite predecir la probabilidad de que ocurra un evento binario: si un paciente abandonara el tratamiento, si cumple criterios diagnosticos, si respondera a una intervencion. En este tutorial completo aprenderemos a implementarla en R paso a paso, desde la preparacion de datos hasta el reporte en formato APA 7.

Que es la regresion logistica?

La regresion logistica binaria es un modelo de regresion generalizado (GLM) que se utiliza cuando la variable dependiente es dicotomica (0/1, si/no, presente/ausente). Mientras que la regresion lineal modela directamente el valor esperado de Y, la regresion logistica modela el logaritmo de los odds (log-odds o logit) de que Y = 1.

La diferencia fundamental con la regresion lineal es que en la logistica los coeficientes (B) representan el cambio en log-odds por cada unidad de incremento en el predictor, no un cambio directo en la variable dependiente. Para obtener una interpretacion mas intuitiva, exponenciamos los coeficientes para obtener los Odds Ratios (OR), una medida clave en la odds ratio e investigacion clinica.

En psicologia clinica, la regresion logistica es ampliamente utilizada para predecir diagnosticos (por ejemplo, presencia o ausencia de un trastorno), respuesta al tratamiento (mejoria clinicamente significativa si/no), abandono terapeutico y factores de riesgo asociados a condiciones clinicas.

Cuando usarla?

Utiliza regresion logistica binaria cuando:

• Tu variable dependiente es dicotomica (dos categorias mutuamente excluyentes).
• Quieres estimar la probabilidad de pertenecer a una categoria en funcion de uno o mas predictores.
• Tus predictores pueden ser continuos, categoricos o una mezcla de ambos.
• Buscas identificar que variables son factores de riesgo o proteccion significativos.

Escenarios tipicos en psicologia: predecir si un paciente abandonara la terapia, si un individuo desarrollara TEPT tras una experiencia traumatica, si un tratamiento farmacologico producira respuesta clinica, o si un estudiante sera identificado con dificultades de aprendizaje.

Ejemplo practico: el dataset

Trabajaremos con un ejemplo clinico realista. Supongamos que un equipo de investigacion recogio datos de N = 280 pacientes adultos con trastorno de ansiedad generalizada (TAG) que iniciaron un programa de terapia cognitivo-conductual grupal de 12 sesiones. El objetivo es identificar que factores predicen el abandono del tratamiento (definido como asistir a menos de 8 de las 12 sesiones).

Primero, crearemos los datos simulados en R para que puedas seguir el tutorial completo:

# Crear datos simulados reproducibles
set.seed(2026)
n <- 280

age        <- round(rnorm(n, mean = 38, sd = 11))
bai        <- round(rnorm(n, mean = 28, sd = 9))
depression <- rbinom(n, 1, prob = 0.40)
wai        <- round(rnorm(n, mean = 55, sd = 12))

# Generar variable dependiente con relaciones conocidas
logit_p <- -1.5 + 0.02 * age + 0.07 * bai + 0.80 * depression - 0.05 * wai
prob    <- plogis(logit_p)
dropout <- rbinom(n, 1, prob)

# Crear data frame
datos <- data.frame(dropout, age, bai, depression, wai)

# Asegurar rangos realistas
datos$age <- pmax(18, pmin(65, datos$age))
datos$bai <- pmax(0, pmin(63, datos$bai))
datos$wai <- pmax(12, pmin(84, datos$wai))
datos$depression <- factor(datos$depression, levels = c(0, 1),
                           labels = c("No", "Si"))

Paso 1: Exploracion y preparacion de datos

Antes de ajustar cualquier modelo, debemos conocer la estructura de nuestros datos, verificar la distribucion de la variable dependiente y buscar datos faltantes.

# Estructura del data frame
str(datos)

# Resumen descriptivo
summary(datos)

# Distribucion de la variable dependiente
table(datos$dropout)
prop.table(table(datos$dropout))

# Verificar datos faltantes
colSums(is.na(datos))

# Descriptivos por grupo (abandono vs. no abandono)
aggregate(cbind(age, bai, wai) ~ dropout, data = datos, FUN = mean)
aggregate(cbind(age, bai, wai) ~ dropout, data = datos, FUN = sd)

Es crucial verificar el balance de clases de la variable dependiente. Si una categoria tiene menos del 10-15% de los casos, podriamos necesitar tecnicas especiales (sobremuestreo, submuestreo o pesos). En nuestro ejemplo esperamos aproximadamente un 30-35% de abandono, lo cual es razonable.

Paso 2: Estimar el modelo

En R, la regresion logistica se ajusta con la funcion glm() especificando family = binomial. Por defecto, R usa la funcion de enlace logit.

# Ajustar el modelo de regresion logistica
modelo <- glm(dropout ~ age + bai + depression + wai,
              data = datos,
              family = binomial(link = "logit"))

# Ver el resumen completo
summary(modelo)

La formula dropout ~ age + bai + depression + wai indica que queremos predecir dropout a partir de los cuatro predictores de forma aditiva. El resumen del modelo muestra:

• Estimate: Coeficiente B en escala log-odds. Un valor positivo indica que al aumentar el predictor, aumenta la probabilidad de dropout = 1.
• Std. Error: Error estandar del coeficiente.
• z value: Estadistico de Wald (Estimate / Std. Error). Prueba si el coeficiente es significativamente diferente de cero.
• Pr(>|z|): Valor p asociado al test de Wald.

Al final del resumen, R muestra la devianza nula (modelo sin predictores, solo intercepto) y la devianza residual (modelo con predictores). La diferencia entre ambas indica cuanta variacion explican los predictores. Tambien se reporta el AIC (Akaike Information Criterion), que es util para comparar modelos: un AIC menor indica mejor ajuste relativo.

Paso 3: Interpretar los Odds Ratios

Los coeficientes en log-odds son dificiles de interpretar directamente. Exponenciandolos obtenemos los Odds Ratios (OR), que son mucho mas intuitivos clinicamente. Si necesitas verificar tus calculos rapidamente, puedes usar nuestra calculadora de odds ratio.

# Odds Ratios
OR <- exp(coef(modelo))
round(OR, 3)

# Intervalos de confianza del 95% para los OR
IC <- exp(confint(modelo))
round(IC, 3)

# Tabla combinada
tabla_or <- data.frame(
  B       = round(coef(modelo), 3),
  OR      = round(OR, 3),
  IC_inf  = round(IC[, 1], 3),
  IC_sup  = round(IC[, 2], 3),
  p       = round(summary(modelo)$coefficients[, 4], 4)
)
print(tabla_or)

La tabla de resultados esperada se veria aproximadamente asi:

Como interpretar cada OR:

• age (OR = 1.02): Cada anio adicional de edad multiplica los odds de abandono por 1.02 (un 2% mas), pero el efecto no es estadisticamente significativo (p = .102). La edad no es un predictor relevante en este modelo.
• bai (OR = 1.07): Cada punto adicional en el BAI multiplica los odds de abandono por 1.07. Es decir, mayor ansiedad basal se asocia significativamente con mayor probabilidad de abandonar el tratamiento (p < .001).
• depressionSi (OR = 2.23): Los pacientes con depresion comorbida tienen 2.23 veces mas odds de abandonar el tratamiento que quienes no presentan depresion (p = .003). Este es un factor de riesgo clinicamente importante.
• wai (OR = 0.95): Cada punto adicional en alianza terapeutica multiplica los odds de abandono por 0.95, es decir, los reduce en un 5%. Mayor alianza terapeutica es un factor protector significativo (p = .001).

Paso 4: Evaluar el ajuste del modelo

A diferencia de la regresion lineal, no existe un R-cuadrado directo en regresion logistica. Usamos pseudo-R-cuadrado y otras medidas de ajuste.

# Pseudo-R cuadrado de Nagelkerke
install.packages("DescTools")  # solo la primera vez
library(DescTools)
PseudoR2(modelo, which = "Nagelkerke")

# Test de Hosmer-Lemeshow
install.packages("ResourceSelection")
library(ResourceSelection)
hoslem.test(datos$dropout, fitted(modelo), g = 10)

# Matriz de confusion (punto de corte = 0.5)
predicciones <- ifelse(fitted(modelo) > 0.5, 1, 0)
confusion    <- table(Observado = datos$dropout, Predicho = predicciones)
print(confusion)

# Metricas de clasificacion
accuracy    <- sum(diag(confusion)) / sum(confusion)
sensitivity <- confusion[2, 2] / sum(confusion[2, ])
specificity <- confusion[1, 1] / sum(confusion[1, ])

cat("Accuracy:", round(accuracy, 3), "\n")
cat("Sensibilidad:", round(sensitivity, 3), "\n")
cat("Especificidad:", round(specificity, 3), "\n")

# AIC del modelo
AIC(modelo)

El pseudo-R-cuadrado de Nagelkerke oscila entre 0 y 1 y se interpreta de forma analoga al R-cuadrado de la regresion lineal, aunque no es exactamente equivalente. Valores entre .20 y .40 se consideran aceptables en ciencias sociales. El test de Hosmer-Lemeshow evalua si las probabilidades predichas coinciden con las frecuencias observadas: un resultado no significativo (p > .05) indica buen ajuste. La sensibilidad es la proporcion de verdaderos abandonos correctamente identificados, y la especificidad es la proporcion de no-abandonos correctamente identificados.

Paso 5: Curva ROC y AUC

La curva ROC (Receiver Operating Characteristic) es la herramienta estandar para evaluar la capacidad discriminativa de un modelo logistico. Representa la sensibilidad frente a 1 - especificidad para todos los posibles puntos de corte.

# Instalar y cargar pROC
install.packages("pROC")
library(pROC)

# Calcular la curva ROC
roc_obj <- roc(datos$dropout, fitted(modelo))

# Graficar la curva ROC
plot(roc_obj,
     col = "#2d9cdb",
     lwd = 2,
     main = "Curva ROC - Modelo de Abandono",
     print.auc = TRUE,
     auc.polygon = TRUE,
     auc.polygon.col = "#2d9cdb22")

# Valor del AUC con intervalo de confianza
auc(roc_obj)
ci.auc(roc_obj)

# Punto de corte optimo (maximiza sensibilidad + especificidad)
corte_optimo <- coords(roc_obj, "best", ret = c("threshold", "sensitivity", "specificity"))
print(corte_optimo)

El punto de corte optimo no siempre es 0.5. La funcion coords() con el criterio "best" utiliza el indice de Youden (maximiza sensibilidad + especificidad - 1) para encontrar el umbral que mejor separa ambos grupos. En contextos clinicos, puede ser preferible priorizar la sensibilidad (detectar todos los abandonos posibles) sobre la especificidad.

Paso 6: Verificar supuestos

Aunque la regresion logistica tiene menos supuestos que la lineal, hay varios que debemos verificar:

1. Linealidad del logit

Para las variables continuas, la relacion entre el predictor y el log-odds debe ser lineal. Esto se comprueba con el test de Box-Tidwell, que incluye la interaccion entre cada predictor continuo y su logaritmo natural:

# Test de Box-Tidwell para linealidad del logit
# Crear los terminos de interaccion con el log
datos$age_log <- datos$age * log(datos$age)
datos$bai_log <- datos$bai * log(datos$bai + 1)  # +1 para evitar log(0)
datos$wai_log <- datos$wai * log(datos$wai)

modelo_bt <- glm(dropout ~ age + bai + depression + wai +
                           age_log + bai_log + wai_log,
                 data = datos, family = binomial)
summary(modelo_bt)

# Si los terminos *_log NO son significativos,
# el supuesto de linealidad se cumple

2. Ausencia de multicolinealidad

# VIF (Variance Inflation Factor)
library(car)
vif(modelo)

# VIF > 5 indica multicolinealidad problematica
# VIF > 10 indica multicolinealidad severa

3. Independencia de observaciones

Cada observacion debe ser independiente de las demas. Esto es un requisito de disenio, no se verifica estadisticamente de forma simple. En nuestro ejemplo, cada paciente aparece una sola vez, por lo que se cumple. Si tuvieramos medidas repetidas o datos agrupados (pacientes dentro de terapeutas), necesitariamos un modelo multinivel.

4. Tamano muestral adecuado

Como ya mencionamos, la regla EPV (eventos por variable) requiere al menos 10 casos de la categoria menos frecuente por cada predictor. Verificamos:

# Verificar EPV
n_eventos   <- min(table(datos$dropout))
n_predictores <- 4
epv <- n_eventos / n_predictores
cat("Eventos por variable:", epv, "\n")
cat("Criterio cumplido:", epv >= 10, "\n")

Paso 7: Modelo paso a paso (opcional)

En un enfoque exploratorio, podemos usar seleccion por pasos (stepwise) basada en AIC para identificar el mejor subconjunto de predictores. Sin embargo, esta tecnica tiene limitaciones importantes y no debe sustituir a la seleccion guiada por teoria.

# Modelo completo
modelo_full <- glm(dropout ~ age + bai + depression + wai,
                   data = datos, family = binomial)

# Seleccion hacia atras con AIC
modelo_step <- step(modelo_full, direction = "backward", trace = 1)
summary(modelo_step)

# Comparar modelos con test de razon de verosimilitud
modelo_reducido <- glm(dropout ~ bai + depression + wai,
                       data = datos, family = binomial)
anova(modelo_reducido, modelo_full, test = "Chisq")

El test de razon de verosimilitud compara el modelo reducido (sin edad) contra el completo. Si el resultado no es significativo (p > .05), el modelo mas simple es preferible porque explica practicamente lo mismo con menos parametros.

Como reportar los resultados en APA 7

Errores frecuentes

Estos son los errores mas comunes al realizar e interpretar una regresion logistica:

Codigo completo reproducible

A continuacion se presenta el script completo que puedes copiar y pegar directamente en R o RStudio:

# =============================================================
# REGRESION LOGISTICA EN R: TUTORIAL COMPLETO
# Prediccion de abandono terapeutico en pacientes con TAG
# =============================================================

# --- 0. Instalar y cargar paquetes ---
# install.packages(c("DescTools", "ResourceSelection", "pROC", "car"))
library(DescTools)
library(ResourceSelection)
library(pROC)
library(car)

# --- 1. Crear datos simulados ---
set.seed(2026)
n <- 280

age        <- round(rnorm(n, mean = 38, sd = 11))
bai        <- round(rnorm(n, mean = 28, sd = 9))
depression <- rbinom(n, 1, prob = 0.40)
wai        <- round(rnorm(n, mean = 55, sd = 12))

logit_p <- -1.5 + 0.02 * age + 0.07 * bai + 0.80 * depression - 0.05 * wai
prob    <- plogis(logit_p)
dropout <- rbinom(n, 1, prob)

datos <- data.frame(dropout, age, bai, depression, wai)
datos$age <- pmax(18, pmin(65, datos$age))
datos$bai <- pmax(0, pmin(63, datos$bai))
datos$wai <- pmax(12, pmin(84, datos$wai))
datos$depression <- factor(datos$depression, levels = c(0, 1),
                           labels = c("No", "Si"))

# --- 2. Exploracion de datos ---
str(datos)
summary(datos)
table(datos$dropout)
prop.table(table(datos$dropout))
colSums(is.na(datos))
aggregate(cbind(age, bai, wai) ~ dropout, data = datos, FUN = mean)

# --- 3. Ajustar el modelo ---
modelo <- glm(dropout ~ age + bai + depression + wai,
              data = datos, family = binomial(link = "logit"))
summary(modelo)

# --- 4. Odds Ratios e intervalos de confianza ---
OR <- exp(coef(modelo))
IC <- exp(confint(modelo))
tabla_or <- data.frame(
  B      = round(coef(modelo), 3),
  OR     = round(OR, 3),
  IC_inf = round(IC[, 1], 3),
  IC_sup = round(IC[, 2], 3),
  p      = round(summary(modelo)$coefficients[, 4], 4)
)
print(tabla_or)

# --- 5. Evaluacion del ajuste ---
PseudoR2(modelo, which = "Nagelkerke")
hoslem.test(datos$dropout, fitted(modelo), g = 10)

predicciones <- ifelse(fitted(modelo) > 0.5, 1, 0)
confusion    <- table(Observado = datos$dropout, Predicho = predicciones)
print(confusion)

accuracy    <- sum(diag(confusion)) / sum(confusion)
sensitivity <- confusion[2, 2] / sum(confusion[2, ])
specificity <- confusion[1, 1] / sum(confusion[1, ])
cat("Accuracy:", round(accuracy, 3), "\n")
cat("Sensibilidad:", round(sensitivity, 3), "\n")
cat("Especificidad:", round(specificity, 3), "\n")

# --- 6. Curva ROC y AUC ---
roc_obj <- roc(datos$dropout, fitted(modelo))
plot(roc_obj, col = "#2d9cdb", lwd = 2,
     main = "Curva ROC - Modelo de Abandono",
     print.auc = TRUE, auc.polygon = TRUE,
     auc.polygon.col = "#2d9cdb22")
auc(roc_obj)
ci.auc(roc_obj)
corte_optimo <- coords(roc_obj, "best",
                       ret = c("threshold", "sensitivity", "specificity"))
print(corte_optimo)

# --- 7. Verificar supuestos ---
# Multicolinealidad
vif(modelo)

# Linealidad del logit (Box-Tidwell)
datos$age_log <- datos$age * log(datos$age)
datos$bai_log <- datos$bai * log(datos$bai + 1)
datos$wai_log <- datos$wai * log(datos$wai)
modelo_bt <- glm(dropout ~ age + bai + depression + wai +
                           age_log + bai_log + wai_log,
                 data = datos, family = binomial)
summary(modelo_bt)

# EPV
n_eventos <- min(table(datos$dropout))
cat("EPV:", n_eventos / 4, "\n")

# --- 8. Modelo paso a paso (opcional) ---
modelo_step <- step(modelo, direction = "backward", trace = 1)
summary(modelo_step)

# Comparar modelo completo vs. reducido
modelo_reducido <- glm(dropout ~ bai + depression + wai,
                       data = datos, family = binomial)
anova(modelo_reducido, modelo, test = "Chisq")