Tasta ahora, hemos trabajado con el Modelo Lineal Clásico (LM), una herramienta fundamental en estadística. En un modelo lineal, asumimos que nuestra variable de respuesta (o dependiente) es continua y que los errores o residuales siguen una distribución normal.
Pensemos en un ejemplo clásico en comportamiento alimentario: queremos predecir el Índice de Masa Corporal (IMC) de un paciente basándonos en su ingesta calórica diaria y sus niveles de ansiedad. La estructura del modelo es:
\[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \epsilon_i\]
Donde \(Y_i\) (el IMC) es continuo,
y \(\epsilon_i\) representa el error
aleatorio normalmente distribuido. Cuando nuestros datos cumplen estos
supuestos, la función lm() en R es perfecta.
En la investigación del comportamiento alimentario, la realidad clínica a menudo no se ajusta a una “distribución normal” ni a variables continuas. ¿Qué pasa cuando queremos investigar preguntas como las siguientes?
Variables binarias (Sí/No): Queremos predecir la presencia o ausencia de un patrón de alimentación disfuncional (ej. 1 = Sí presenta el patrón, 0 = No lo presenta) a partir de ciertos rasgos apetitivos.
Variables de conteo: Queremos modelar el número de episodios de atracón que un paciente reporta en una semana (ej. 0, 1, 2, 5 episodios).
Si intentamos trazar una línea recta clásica a través de datos que solo son ceros y unos (presencia/ausencia), la línea predecirá valores absurdos (como probabilidades negativas o mayores al 100%).
Como su nombre lo indica, un GLM es una ampliación o generalización del modelo lineal clásico. Los GLM nos permiten utilizar el mismo marco lógico de la regresión lineal, pero flexibilizando las reglas para poder trabajar con variables de respuesta que no son normales, como conteos, proporciones o categorías binarias.
Para entender cómo logran esto, debemos ver los tres componentes fundamentales de cualquier GLM:
En lugar de estar atados exclusivamente a la distribución Normal (Campana de Gauss), los GLM nos permiten elegir una distribución de la “familia exponencial” que se ajuste a la naturaleza de nuestra variable.
¿Es presencia/ausencia de un trastorno? Usamos la familia Binomial.
¿Es el conteo de episodios de purga? Usamos la familia Poisson.
¿Es el IMC continuo? Usamos la familia Gaussiana (Normal).
Esta es la parte matemática que ya conoces. Es la combinación de nuestras variables independientes (predictores) multiplicadas por sus coeficientes. Se suele representar con la letra griega eta (\(\eta\)): \[\eta_i = \beta_0 + \beta_1 X_{1i} + \dots + \beta_k X_{ki}\]
Esta es la “magia” del GLM. Es una transformación matemática, denotada como \(g(\mu)\), que sirve como puente entre el componente aleatorio y el predictor lineal. La función de enlace transforma la escala de los datos para que el modelo pueda calcular las predicciones correctamente sin salirse de los límites lógicos (por ejemplo, evitando que una probabilidad predicha sea menor a 0).
En regresión logística (datos binarios), la función de enlace es el Logit.
En regresión de Poisson (conteos), la función de enlace es el Log. \[g(\mu_i) = \eta_i = \beta_0 + \beta_1 X_{1i} + \dots + \beta_k X_{ki}\]
La belleza conceptual de los GLM es que los modelos LM clásicos son simplemente un caso especial dentro del universo de los GLM.
Si a un GLM le indicas que su distribución es Gaussiana (Normal) y que su función de enlace es la Identidad (es decir, que no transforme nada), el resultado matemático es exactamente idéntico a un modelo lineal estándar.
Por lo tanto, aprender a usar la función glm() en R (en
lugar de solo lm()) te abre la puerta a modelar casi
cualquier tipo de fenómeno clínico, conductual o nutricional, ajustando
únicamente los parámetros de family y
link.
Nombre del dataset:
behavior_food_glm
| Variable | Tipo | Descripción |
|---|---|---|
| edad | continua | edad del participante |
| sexo | binaria | 0 = mujer, 1 = hombre |
| estres | escala | nivel de estrés (1–10) |
| sueno | continua | horas de sueño |
| actividad | continua | minutos de actividad física por día |
| impulsividad | escala | impulsividad alimentaria (1–10) |
| IMC | continua | índice de masa corporal |
| obesidad | binaria | 1 = IMC ≥ 30 |
| snacks_semana | conteo | número de snacks consumidos por semana |
| episodios_emocionales | conteo | episodios de comer emocional |
| tiempo_redes | continua positiva | minutos diarios en redes sociales |
set.seed(123)
n <- 200
edad <- round(rnorm(n, 35, 10))
sexo <- rbinom(n, 1, 0.5)
estres <- round(runif(n, 1, 10),1)
sueno <- round(rnorm(n, 7, 1.2),1)
actividad <- round(rnorm(n, 45, 20))
impulsividad <- round(runif(n, 1, 10),1)
IMC <- round(22 +
0.4*estres +
0.3*impulsividad -
0.2*actividad/10 +
rnorm(n,0,2),1)
obesidad <- ifelse(IMC >= 30,1,0)
snacks_semana <- rpois(n,
lambda = exp(1 +
0.08*estres +
0.06*impulsividad))
episodios_emocionales <- rpois(n,
lambda = exp(0.5 +
0.12*estres +
0.08*impulsividad))
tiempo_redes <- rgamma(n,
shape = 5,
scale = 20 + 2*estres)
datos <- data.frame(
edad,
sexo,
estres,
sueno,
actividad,
impulsividad,
IMC,
obesidad,
snacks_semana,
episodios_emocionales,
tiempo_redes
)
head(datos)
## edad sexo estres sueno actividad impulsividad IMC obesidad snacks_semana
## 1 29 1 3.1 6.9 66 6.8 21.9 0 7
## 2 33 0 7.2 5.6 44 2.1 23.0 0 7
## 3 51 1 3.0 6.2 44 3.3 23.9 0 6
## 4 36 1 3.9 7.0 15 8.4 29.1 0 2
## 5 36 0 2.6 7.8 61 8.2 26.4 0 3
## 6 52 0 8.2 5.0 41 1.4 23.6 0 6
## episodios_emocionales tiempo_redes
## 1 6 217.88300
## 2 5 89.09039
## 3 1 159.86682
## 4 1 84.79594
## 5 7 40.85726
## 6 2 115.55849
A continuación, exploraremos cómo la naturaleza de nuestra variable de respuesta dicta el tipo de modelo que debemos utilizar.
Uso: Variable de respuesta continua y teóricamente normal.
Ejemplo en clase: Queremos evaluar cómo el estrés, la actividad física y la edad explican las variaciones en el Índice de Masa Corporal (IMC). Asumimos que el IMC es una métrica continua que puede tomar cualquier valor dentro de un rango razonable.
modelo_lm <- lm(IMC ~ estres + actividad + edad, data = datos)
summary(modelo_lm)
##
## Call:
## lm(formula = IMC ~ estres + actividad + edad, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9207 -1.4433 -0.1031 1.2254 5.3009
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.978851 0.810792 29.575 < 2e-16 ***
## estres 0.268310 0.060392 4.443 1.48e-05 ***
## actividad -0.011303 0.007718 -1.464 0.145
## edad 0.005123 0.016763 0.306 0.760
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.225 on 196 degrees of freedom
## Multiple R-squared: 0.1048, Adjusted R-squared: 0.09114
## F-statistic: 7.652 on 3 and 196 DF, p-value: 7.309e-05
Guía de Interpretación:
Los coeficientes nos hablan de cambios absolutos:
Un coeficiente positivo en estres indica que, a
mayor estrés, mayor es el IMC esperado (manteniendo lo demás
constante).
Un coeficiente negativo en actividad indica el
efecto protector: a mayor actividad, menor IMC.
Uso: Variable de respuesta binaria o dicotómica (0 o 1).
Ejemplo en clase: Ahora no nos interesa el IMC
exacto, sino clasificar si un paciente tiene obesidad
clínica (1 = Sí, 0 = No). Un modelo lineal clásico podría
predecir probabilidades sin sentido (ej. \(-15\%\) o \(120\%\)). Al usar
family = binomial, la función de enlace logit
comprime las predicciones para que vivan estrictamente entre \(0\) y \(1\) (probabilidad).
modelo_logit <- glm(obesidad ~ estres + impulsividad + actividad,
family = binomial,
data = datos)
summary(modelo_logit)
##
## Call:
## glm(formula = obesidad ~ estres + impulsividad + actividad, family = binomial,
## data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.0514402 2.0065835 -1.521 0.128
## estres 0.0003926 0.1984367 0.002 0.998
## impulsividad -0.0144320 0.1968880 -0.073 0.942
## actividad -0.0186707 0.0246006 -0.759 0.448
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 39.216 on 199 degrees of freedom
## Residual deviance: 38.624 on 196 degrees of freedom
## AIC: 46.624
##
## Number of Fisher Scoring iterations: 7
Interpretación mediante Odds Ratios (OR):
En regresión logística, los coeficientes crudos (los log-odds) son difíciles de interpretar clínicamente. Por ello, los exponenciamos (calculando \(e^{\beta}\)) para obtener los Odds Ratios (Razones de Momios).
Un OR \(> 1\) indica que la variable aumenta el riesgo de presentar obesidad.
Un OR \(< 1\) indica un factor protector.
# Cálculo de los Odds Ratios
exp(coef(modelo_logit))
## (Intercept) estres impulsividad actividad
## 0.04729077 1.00039271 0.98567163 0.98150251
Uso: Variables de conteo discreto (números enteros \(\ge 0\)). No puedes tener eventos negativos ni “fracciones de evento” (ej. no puedes consumir -3 snacks o 2.5 snacks enteros en ocasiones separadas).
Ejemplo en clase: Queremos predecir el número de snacks ultraprocesados consumidos a la semana basándonos en el estrés y la impulsividad.
modelo_poisson <- glm(snacks_semana ~ estres + impulsividad,
family = poisson,
data = datos)
summary(modelo_poisson)
##
## Call:
## glm(formula = snacks_semana ~ estres + impulsividad, family = poisson,
## data = datos)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.98781 0.10369 9.527 < 2e-16 ***
## estres 0.09290 0.01113 8.347 < 2e-16 ***
## impulsividad 0.05145 0.01071 4.806 1.54e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 305.73 on 199 degrees of freedom
## Residual deviance: 223.33 on 197 degrees of freedom
## AIC: 940.01
##
## Number of Fisher Scoring iterations: 4
Guía de Interpretación:
El modelo Poisson asume que la media y la varianza de los datos son iguales. Los coeficientes exponenciados representan razones de tasas (Incident Rate Ratios). Si el coeficiente exponencial del estrés es \(1.08\), significa que por cada punto extra de estrés, el número de snacks consumidos aumenta en un \(8\%\).
Uso: Variables continuas y estrictamente positivas, que a menudo presentan asimetría positiva (una “cola larga” hacia la derecha).
Ejemplo en clase: El tiempo invertido en redes sociales tiene un límite inferior estricto (\(0\) minutos) pero no tiene un límite superior fijo. La mayoría de las personas lo usa moderadamente, pero unos pocos pacientes reportan un uso extremo. La distribución Normal no ajusta bien aquí, pero la distribución Gamma sí.
modelo_gamma <- glm(tiempo_redes ~ estres + sueno,
family = Gamma(link="log"),
data = datos)
summary(modelo_gamma)
##
## Call:
## glm(formula = tiempo_redes ~ estres + sueno, family = Gamma(link = "log"),
## data = datos)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.67569 0.21271 21.982 < 2e-16 ***
## estres 0.05514 0.01240 4.448 1.45e-05 ***
## sueno 0.01080 0.02806 0.385 0.701
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.2106307)
##
## Null deviance: 46.125 on 199 degrees of freedom
## Residual deviance: 42.062 on 197 degrees of freedom
## AIC: 2249.5
##
## Number of Fisher Scoring iterations: 5
Uso: Variables de conteo donde ocurre sobredispersión (la varianza es mucho mayor que la media).
Ejemplo en clase: Los episodios de comer emocional son conteos. Podríamos intentar usar Poisson. Sin embargo, en psicología clínica es común que gran parte de la muestra reporte \(0\) o \(1\) episodio, mientras que un subgrupo clínico reporta números muy altos (ej. \(15\) o \(20\) episodios). Esto “rompe” el supuesto de la regresión de Poisson (donde media \(=\) varianza). El modelo Binomial Negativo introduce un parámetro extra que relaja este supuesto y corrige el error.
Nota: Para esto usamos la función glm.nb de la
paquetería MASS.
library(MASS)
modelo_nb <- glm.nb(episodios_emocionales ~ estres + impulsividad,
data = datos)
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
summary(modelo_nb)
##
## Call:
## glm.nb(formula = episodios_emocionales ~ estres + impulsividad,
## data = datos, init.theta = 74789.86771, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.45362 0.11675 3.885 0.000102 ***
## estres 0.11677 0.01226 9.523 < 2e-16 ***
## impulsividad 0.08399 0.01186 7.080 1.44e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(74789.87) family taken to be 1)
##
## Null deviance: 299.55 on 199 degrees of freedom
## Residual deviance: 175.31 on 197 degrees of freedom
## AIC: 852.52
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 74790
## Std. Err.: 869484
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -844.521
El paquete ggplot2 es extremadamente útil para graficar
modelos GLM, ya que nos permite visualizar la curva de predicción no
lineal de vuelta en la escala original de los datos. Aquí graficamos la
predicción de conteo (Poisson) del número de snacks en función del
estrés.
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
ggplot(datos, aes(x = estres, y = snacks_semana)) +
geom_point(alpha = 0.5, color = "darkorange") +
geom_smooth(method = "glm",
method.args = list(family = "poisson"),
color = "darkblue") +
theme_minimal() +
labs(title = "Predicción de Consumo de Snacks mediante Modelo Poisson",
x = "Nivel de Estrés Reportado (1-10)",
y = "Snacks Ultraprocesados por Semana")
## `geom_smooth()` using formula = 'y ~ x'