Tutorial R: Regresión lineal en R

Regresión simple y múltiple con lm()

La regresión lineal es una técnica estadística fundamental que nos permite modelar la relación entre variables. En R, esta técnica se implementa principalmente a través de la función lm() (linear model), una herramienta versátil y potente para análisis predictivo.

La función lm() nos permite crear modelos que describen cómo una variable dependiente (también llamada variable respuesta) se relaciona con una o más variables independientes (o predictoras). Dependiendo del número de variables predictoras, trabajaremos con regresión simple o múltiple.

Regresión lineal simple

La regresión lineal simple modela la relación entre dos variables: una dependiente y una independiente. La ecuación básica es:

y = β₀ + β₁x + ε

Donde:

• y es la variable dependiente
• x es la variable independiente
• β₀ es el intercepto (valor de y cuando x=0)
• β₁ es la pendiente (cambio en y por cada unidad de cambio en x)
• ε representa el error aleatorio

Veamos cómo implementar esto en R con un ejemplo sencillo. Primero, creemos algunos datos simulados:

# Crear datos de ejemplo
set.seed(123)  # Para reproducibilidad
x <- 1:30
y <- 2 + 3*x + rnorm(30, mean=0, sd=5)
datos <- data.frame(x=x, y=y)

# Visualizar los datos
plot(datos$x, datos$y, 
     main="Relación entre X e Y",
     xlab="Variable X", 
     ylab="Variable Y",
     pch=19, col="blue")

Ahora, ajustemos un modelo de regresión lineal simple:

# Ajustar modelo de regresión simple
modelo_simple <- lm(y ~ x, data=datos)

# Ver resumen del modelo
summary(modelo_simple)

El resultado del summary() nos mostrará información detallada sobre el modelo, incluyendo los coeficientes estimados, errores estándar, valores p y estadísticas de ajuste.

Podemos visualizar la línea de regresión sobre nuestros datos:

# Graficar los datos con la línea de regresión
plot(datos$x, datos$y, 
     main="Regresión Lineal Simple",
     xlab="Variable X", 
     ylab="Variable Y",
     pch=19, col="blue")
abline(modelo_simple, col="red", lwd=2)

Sintaxis básica de lm()

La función lm() tiene una sintaxis muy intuitiva:

lm(formula, data, subset, weights, ...)

Donde:

• formula: especifica la relación entre variables usando la notación y ~ x
• data: el data frame que contiene las variables
• subset: opcionalmente, permite seleccionar un subconjunto de observaciones
• weights: permite especificar pesos para las observaciones

La fórmula es el componente más importante y sigue una sintaxis especial:

variable_dependiente ~ variable_independiente

El símbolo ~ se lee como "es modelado por" o "en función de".

Regresión lineal múltiple

Cuando tenemos más de una variable predictora, utilizamos la regresión lineal múltiple. La ecuación se expande a:

y = β₀ + β₁x₁ + β₂x₂ + ... + βₙxₙ + ε

Donde cada βᵢ representa el efecto de la variable xᵢ sobre y, manteniendo constantes las demás variables.

Veamos un ejemplo con datos simulados:

# Crear datos para regresión múltiple
set.seed(456)
n <- 100
x1 <- rnorm(n, mean=10, sd=2)
x2 <- rnorm(n, mean=20, sd=5)
x3 <- rnorm(n, mean=15, sd=3)
y <- 5 + 2*x1 - 1.5*x2 + 0.5*x3 + rnorm(n, mean=0, sd=3)

datos_multi <- data.frame(y=y, x1=x1, x2=x2, x3=x3)

# Ajustar modelo de regresión múltiple
modelo_multiple <- lm(y ~ x1 + x2 + x3, data=datos_multi)

# Ver resumen del modelo
summary(modelo_multiple)

En la fórmula de regresión múltiple, simplemente añadimos más variables con el operador +:

y ~ x1 + x2 + x3

Fórmulas avanzadas en lm()

R ofrece una sintaxis flexible para especificar modelos más complejos:

• Incluir todas las variables de un data frame:

# Usar todas las variables excepto y como predictores
modelo <- lm(y ~ ., data=datos_multi)

• Excluir variables específicas:

# Excluir x3 del modelo
modelo <- lm(y ~ . - x3, data=datos_multi)

• Interacciones entre variables (cuando el efecto de una variable depende del valor de otra):

# Incluir interacción entre x1 y x2
modelo_interaccion <- lm(y ~ x1 + x2 + x3 + x1:x2, data=datos_multi)

# Forma abreviada para incluir variables e interacción
modelo_interaccion2 <- lm(y ~ x1*x2 + x3, data=datos_multi)

• Transformaciones de variables:

# Incluir término cuadrático
modelo_cuadratico <- lm(y ~ x1 + I(x1^2) + x2 + x3, data=datos_multi)

# Transformación logarítmica
modelo_log <- lm(log(y) ~ x1 + x2 + x3, data=datos_multi)

Ejemplo práctico con datos reales

Usemos el conjunto de datos mtcars incluido en R para un ejemplo más realista:

# Cargar datos
data(mtcars)
head(mtcars)

# Regresión simple: Consumo en función de potencia
modelo_simple <- lm(mpg ~ hp, data=mtcars)
summary(modelo_simple)

# Visualizar
plot(mtcars$hp, mtcars$mpg, 
     main="Consumo vs. Potencia",
     xlab="Potencia (hp)", 
     ylab="Consumo (mpg)",
     pch=19)
abline(modelo_simple, col="red", lwd=2)

Ahora, un modelo múltiple para predecir el consumo de combustible:

# Regresión múltiple: Consumo en función de potencia, peso y cilindros
modelo_multiple <- lm(mpg ~ hp + wt + cyl, data=mtcars)
summary(modelo_multiple)

Extracción de componentes del modelo

Después de ajustar un modelo, podemos extraer diversos componentes:

# Coeficientes del modelo
coef(modelo_multiple)

# Valores ajustados (predicciones para los datos de entrenamiento)
fitted_values <- fitted(modelo_multiple)

# Residuos
residuos <- residuals(modelo_multiple)

# Graficar residuos vs valores ajustados
plot(fitted_values, residuos,
     main="Residuos vs Valores Ajustados",
     xlab="Valores Ajustados", 
     ylab="Residuos",
     pch=19)
abline(h=0, col="red", lty=2)

Comparación de modelos

Podemos comparar modelos para determinar cuál es más adecuado:

# Modelo simple
modelo1 <- lm(mpg ~ hp, data=mtcars)

# Modelo con más variables
modelo2 <- lm(mpg ~ hp + wt, data=mtcars)

# Modelo completo
modelo3 <- lm(mpg ~ hp + wt + cyl, data=mtcars)

# Comparar con ANOVA
anova(modelo1, modelo2, modelo3)

# Comparar con AIC (menor es mejor)
AIC(modelo1, modelo2, modelo3)

La función lm() es el punto de partida para el análisis de regresión en R, proporcionando una base sólida para modelar relaciones lineales entre variables. Dominar su uso nos permite construir modelos predictivos efectivos y comprender mejor las relaciones entre nuestras variables de interés.

Interpretación de coeficientes y diagnósticos

Una vez que hemos ajustado un modelo de regresión lineal con la función lm(), el siguiente paso crucial es interpretar correctamente los resultados y evaluar la calidad del modelo. Esta interpretación nos permite entender qué nos dice realmente el modelo sobre las relaciones entre nuestras variables.

Interpretación de los coeficientes

Los coeficientes de regresión representan la esencia de nuestro modelo. Para examinarlos, utilizamos la función summary():

# Continuando con el ejemplo de mtcars
modelo <- lm(mpg ~ hp + wt, data=mtcars)
resumen <- summary(modelo)
resumen

La salida incluirá una tabla de coeficientes con esta estructura:

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
hp          -0.03177    0.00903  -3.519  0.00145 ** 
wt          -3.87783    0.63273  -6.129 1.12e-06 ***

Cada coeficiente tiene varios componentes importantes:

• Estimate (Estimación): El valor del coeficiente que indica el cambio esperado en la variable dependiente por cada unidad de cambio en la variable independiente, manteniendo constantes las demás variables.

• Std. Error (Error estándar): Mide la precisión de nuestra estimación. Valores más pequeños indican estimaciones más precisas.

• t value (Valor t): Es el resultado de dividir la estimación entre su error estándar. Valores absolutos grandes sugieren que el coeficiente es significativamente diferente de cero.

• Pr(>|t|) (Valor p): Indica la probabilidad de observar un valor t tan extremo si el coeficiente fuera realmente cero. Valores pequeños (típicamente < 0.05) sugieren que la variable es estadísticamente significativa.

• Símbolos de significancia: Los asteriscos proporcionan una referencia visual rápida del nivel de significancia (*** para p < 0.001, ** para p < 0.01, * para p < 0.05).

Interpretemos los coeficientes de nuestro ejemplo:

# Extraer coeficientes
coef(modelo)

• Intercepto (37.22727): Representa el valor esperado de mpg cuando todas las variables predictoras son cero. En este caso, un vehículo teórico con 0 caballos de potencia y 0 peso tendría un consumo de aproximadamente 37.23 millas por galón.

• hp (-0.03177): Por cada caballo de potencia adicional, el consumo disminuye en aproximadamente 0.032 millas por galón, manteniendo el peso constante.

• wt (-3.87783): Por cada 1000 libras adicionales de peso (la unidad de wt en mtcars), el consumo disminuye en aproximadamente 3.88 millas por galón, manteniendo la potencia constante.

Estadísticas de ajuste global

El resumen también proporciona estadísticas sobre la calidad general del modelo:

# Estadísticas de ajuste
resumen$r.squared        # R-cuadrado
resumen$adj.r.squared    # R-cuadrado ajustado
resumen$sigma            # Error estándar residual

• R-cuadrado: Indica la proporción de varianza en la variable dependiente explicada por el modelo. Varía de 0 a 1, donde valores más altos indican mejor ajuste.

• R-cuadrado ajustado: Similar al R-cuadrado, pero penaliza la adición de variables predictoras que no mejoran sustancialmente el modelo. Es más útil para comparar modelos con diferentes números de predictores.

• Error estándar residual: Mide la dispersión típica de los residuos alrededor de la línea de regresión. Valores más pequeños indican mejor ajuste.

• Estadístico F y su valor p: Evalúa si al menos una de las variables predictoras tiene una relación significativa con la variable dependiente. Un valor p pequeño sugiere que el modelo es estadísticamente significativo.

Diagnósticos del modelo

Para evaluar si nuestro modelo cumple con los supuestos de la regresión lineal, R ofrece varias herramientas de diagnóstico:

# Gráficos de diagnóstico básicos
par(mfrow=c(2,2))  # Dividir la ventana gráfica en 2x2
plot(modelo)

Este comando genera cuatro gráficos diagnósticos clave:

1. Residuos vs Valores ajustados: Permite verificar la linealidad y homocedasticidad (varianza constante). Idealmente, los puntos deberían distribuirse aleatoriamente alrededor de la línea horizontal en y=0, sin patrones evidentes.

2. Q-Q normal: Evalúa la normalidad de los residuos. Los puntos deberían seguir aproximadamente la línea diagonal.

3. Scale-Location: Otra forma de verificar la homocedasticidad. La línea roja debería ser aproximadamente horizontal con puntos distribuidos uniformemente.

4. Residuos vs Leverage: Ayuda a identificar observaciones influyentes que podrían estar distorsionando el modelo.

Podemos examinar estos aspectos individualmente:

# Verificar normalidad de residuos
shapiro.test(residuals(modelo))

# Identificar valores atípicos potenciales
outliers <- which(abs(rstandard(modelo)) > 2)
mtcars[outliers,]

Multicolinealidad

En regresión múltiple, es importante verificar si existe multicolinealidad (correlación alta entre predictores), que puede afectar la interpretación de los coeficientes:

# Instalar si es necesario
# install.packages("car")
library(car)

# Factor de inflación de varianza (VIF)
vif(modelo)

Valores VIF superiores a 5 o 10 suelen indicar problemas de multicolinealidad.

Análisis de observaciones influyentes

Las observaciones influyentes pueden afectar desproporcionadamente nuestro modelo:

# Calcular distancias de Cook
cooks_d <- cooks.distance(modelo)

# Identificar observaciones potencialmente influyentes
influential <- which(cooks_d > 4/(nrow(mtcars)-length(coef(modelo))))
mtcars[influential,]

# Visualizar
plot(cooks_d, type="h", main="Distancias de Cook", ylab="Distancia de Cook", xlab="Observación")
abline(h=4/(nrow(mtcars)-length(coef(modelo))), col="red", lty=2)

Ejemplo práctico de interpretación

Veamos un ejemplo completo con el conjunto de datos mtcars:

# Ajustar modelo
modelo <- lm(mpg ~ hp + wt + am, data=mtcars)
summary(modelo)

Interpretación de los resultados:

• Si el coeficiente de hp es -0.025 y es significativo (p < 0.05), podemos decir: "Por cada caballo de potencia adicional, el consumo de combustible disminuye en promedio 0.025 millas por galón, manteniendo constantes el peso y el tipo de transmisión."

• Si el coeficiente de am (transmisión manual vs. automática) es 2.5 y es significativo, podemos decir: "Los vehículos con transmisión manual tienen, en promedio, 2.5 millas por galón más de rendimiento que los automáticos, manteniendo constantes la potencia y el peso."

• Si el R-cuadrado ajustado es 0.84, podemos decir: "Este modelo explica aproximadamente el 84% de la variabilidad en el consumo de combustible."

Diagnóstico de problemas comunes

Basándonos en los diagnósticos, podemos identificar y abordar problemas:

• Patrón curvo en Residuos vs Valores ajustados: Sugiere que la relación no es lineal. Podríamos considerar transformaciones o términos cuadráticos:

# Añadir término cuadrático
modelo_mejorado <- lm(mpg ~ hp + I(hp^2) + wt, data=mtcars)

• Embudo en Residuos vs Valores ajustados: Indica heterocedasticidad. Podríamos considerar transformaciones logarítmicas:

# Transformación logarítmica
modelo_log <- lm(log(mpg) ~ hp + wt, data=mtcars)

• Desviaciones en el gráfico Q-Q: Sugieren no normalidad. Podríamos considerar transformaciones o métodos robustos:

# Regresión robusta
# install.packages("MASS")
library(MASS)
modelo_robusto <- rlm(mpg ~ hp + wt, data=mtcars)

La interpretación adecuada de coeficientes y diagnósticos es fundamental para construir modelos de regresión confiables y útiles. Estos análisis nos permiten no solo entender las relaciones entre variables, sino también evaluar críticamente la calidad y limitaciones de nuestros modelos.

Predicciones y validación del modelo

Una vez que hemos construido nuestro modelo de regresión lineal y analizado sus coeficientes, el siguiente paso crucial es evaluar su capacidad predictiva. Un modelo estadísticamente significativo no necesariamente es bueno para hacer predicciones. En esta sección, aprenderemos a utilizar nuestros modelos para realizar predicciones y a validar su rendimiento.

Realizando predicciones con el modelo

La función predict() es la herramienta principal para generar predicciones a partir de un modelo de regresión en R. Esta función toma un modelo ajustado y un conjunto de datos con las variables predictoras, y devuelve los valores estimados de la variable dependiente.

# Continuando con nuestro modelo de mtcars
modelo <- lm(mpg ~ hp + wt, data=mtcars)

# Predicción para los datos originales
predicciones <- predict(modelo)
head(predicciones)

Para hacer predicciones con nuevos datos, creamos un data frame con las mismas variables predictoras que usamos en el modelo:

# Crear nuevos datos para predicción
nuevos_datos <- data.frame(
  hp = c(100, 150, 200),
  wt = c(2.5, 3.0, 3.5)
)

# Realizar predicciones
predicciones_nuevas <- predict(modelo, newdata=nuevos_datos)
predicciones_nuevas

También podemos obtener intervalos de confianza o predicción para nuestras estimaciones:

# Intervalos de confianza (para la media)
predicciones_ic <- predict(modelo, newdata=nuevos_datos, interval="confidence", level=0.95)
predicciones_ic

# Intervalos de predicción (para observaciones individuales)
predicciones_ip <- predict(modelo, newdata=nuevos_datos, interval="prediction", level=0.95)
predicciones_ip

La diferencia entre estos intervalos es importante:

• El intervalo de confianza estima el rango donde se encuentra la media verdadera de la variable dependiente.
• El intervalo de predicción estima el rango donde se encontrará una nueva observación individual.

Visualizando las predicciones

Podemos visualizar nuestras predicciones junto con los datos originales:

# Para regresión simple
modelo_simple <- lm(mpg ~ hp, data=mtcars)

# Crear secuencia de valores para hp
nuevos_hp <- data.frame(hp = seq(min(mtcars$hp), max(mtcars$hp), length.out=100))

# Predecir con intervalos
pred_con_intervalos <- predict(modelo_simple, newdata=nuevos_hp, 
                              interval="prediction", level=0.95)

# Combinar con los valores de hp
datos_grafico <- cbind(nuevos_hp, pred_con_intervalos)

# Graficar
plot(mtcars$hp, mtcars$mpg, 
     main="Predicciones con intervalos",
     xlab="Potencia (hp)", 
     ylab="Consumo (mpg)",
     pch=19, col="blue")

# Añadir línea de predicción
lines(nuevos_hp$hp, datos_grafico[,"fit"], col="red", lwd=2)

# Añadir bandas de predicción
lines(nuevos_hp$hp, datos_grafico[,"lwr"], col="darkgray", lty=2)
lines(nuevos_hp$hp, datos_grafico[,"upr"], col="darkgray", lty=2)

Validación del modelo

Para evaluar adecuadamente el rendimiento predictivo de un modelo, necesitamos ir más allá del análisis de los datos de entrenamiento. Existen varias técnicas para validar modelos:

1. División en conjuntos de entrenamiento y prueba

La técnica más básica consiste en dividir nuestros datos en dos conjuntos:

# Establecer semilla para reproducibilidad
set.seed(123)

# Crear índices para la división (70% entrenamiento, 30% prueba)
n <- nrow(mtcars)
indices_entrenamiento <- sample(1:n, size=round(0.7*n))
datos_entrenamiento <- mtcars[indices_entrenamiento, ]
datos_prueba <- mtcars[-indices_entrenamiento, ]

# Ajustar modelo con datos de entrenamiento
modelo_entrenamiento <- lm(mpg ~ hp + wt, data=datos_entrenamiento)

# Predecir en datos de prueba
predicciones_prueba <- predict(modelo_entrenamiento, newdata=datos_prueba)

# Calcular error de predicción
errores <- datos_prueba$mpg - predicciones_prueba

2. Métricas de evaluación

Para cuantificar el rendimiento del modelo, utilizamos diversas métricas de error:

# Error Cuadrático Medio (MSE)
mse <- mean(errores^2)

# Raíz del Error Cuadrático Medio (RMSE)
rmse <- sqrt(mse)

# Error Absoluto Medio (MAE)
mae <- mean(abs(errores))

# Error Porcentual Absoluto Medio (MAPE)
mape <- mean(abs(errores/datos_prueba$mpg)) * 100

# Mostrar métricas
cat("MSE:", mse, "\n")
cat("RMSE:", rmse, "\n")
cat("MAE:", mae, "\n")
cat("MAPE:", mape, "%\n")

Estas métricas tienen interpretaciones diferentes:

• RMSE: Error típico en las unidades de la variable dependiente.
• MAE: Error promedio en valor absoluto.
• MAPE: Error porcentual promedio (útil para comparar entre diferentes escalas).

3. Validación cruzada

La validación cruzada es una técnica más robusta que utiliza múltiples divisiones de los datos:

# Validación cruzada k-fold (requiere el paquete caret)
# install.packages("caret")
library(caret)

# Configurar control de entrenamiento (5 folds)
control <- trainControl(method="cv", number=5)

# Entrenar modelo con validación cruzada
modelo_cv <- train(mpg ~ hp + wt, data=mtcars, 
                  method="lm", 
                  trControl=control)

# Ver resultados
print(modelo_cv)
modelo_cv$results

La validación cruzada divide los datos en k subconjuntos (folds), entrena el modelo k veces usando k-1 subconjuntos y evalúa en el subconjunto restante. Esto proporciona una estimación más estable del rendimiento del modelo.

4. Comparación de modelos mediante validación

Podemos usar la validación para comparar diferentes modelos y seleccionar el mejor:

# Definir modelos a comparar
modelo1 <- train(mpg ~ hp, data=mtcars, 
                method="lm", trControl=control)

modelo2 <- train(mpg ~ hp + wt, data=mtcars, 
                method="lm", trControl=control)

modelo3 <- train(mpg ~ hp + wt + cyl, data=mtcars, 
                method="lm", trControl=control)

# Comparar resultados
resultados <- resamples(list(Modelo1=modelo1, Modelo2=modelo2, Modelo3=modelo3))
summary(resultados)

# Visualizar comparación
dotplot(resultados)

Diagnóstico de sobreajuste y subajuste

Un aspecto crucial de la validación es detectar problemas de sobreajuste (overfitting) o subajuste (underfitting):

• Sobreajuste: El modelo se ajusta demasiado a los datos de entrenamiento y no generaliza bien a nuevos datos.
• Subajuste: El modelo es demasiado simple y no captura la estructura de los datos.

Podemos detectar estos problemas comparando el rendimiento en los datos de entrenamiento y prueba:

# Calcular RMSE en entrenamiento
predicciones_entrenamiento <- predict(modelo_entrenamiento, newdata=datos_entrenamiento)
errores_entrenamiento <- datos_entrenamiento$mpg - predicciones_entrenamiento
rmse_entrenamiento <- sqrt(mean(errores_entrenamiento^2))

# Comparar con RMSE en prueba
cat("RMSE en entrenamiento:", rmse_entrenamiento, "\n")
cat("RMSE en prueba:", rmse, "\n")

Si el RMSE en entrenamiento es mucho menor que en prueba, podría indicar sobreajuste.

Validación de supuestos para predicción

Para que las predicciones sean confiables, debemos verificar que se cumplan los supuestos del modelo:

# Comprobar normalidad de residuos para intervalos de predicción válidos
shapiro.test(residuals(modelo))

# Comprobar homocedasticidad
plot(fitted(modelo), residuals(modelo),
     xlab="Valores ajustados", ylab="Residuos",
     main="Residuos vs Valores ajustados")
abline(h=0, col="red", lty=2)

Ejemplo práctico: Predicción del consumo de combustible

Veamos un ejemplo completo de predicción y validación:

# Cargar datos
data(mtcars)

# Dividir datos
set.seed(123)
indices <- sample(1:nrow(mtcars), size=0.7*nrow(mtcars))
train_data <- mtcars[indices, ]
test_data <- mtcars[-indices, ]

# Ajustar modelo
modelo_final <- lm(mpg ~ hp + wt + am, data=train_data)
summary(modelo_final)

# Predecir en datos de prueba
predicciones <- predict(modelo_final, newdata=test_data)

# Evaluar rendimiento
errores <- test_data$mpg - predicciones
rmse <- sqrt(mean(errores^2))
cat("RMSE:", rmse, "\n")

# Crear un data frame con valores reales y predicciones
resultados <- data.frame(
  Actual = test_data$mpg,
  Predicción = predicciones,
  Error = errores
)
print(resultados)

# Visualizar predicciones vs valores reales
plot(test_data$mpg, predicciones,
     xlab="Valores reales", ylab="Predicciones",
     main="Predicciones vs Valores reales")
abline(0, 1, col="red", lty=2)  # Línea de referencia perfecta

Consideraciones finales

Al realizar predicciones y validar modelos, es importante tener en cuenta:

• La calidad de los datos de entrada afecta directamente la calidad de las predicciones.
• Las predicciones son más confiables dentro del rango de los datos utilizados para entrenar el modelo.
• Los intervalos de predicción son generalmente más útiles que las predicciones puntuales, ya que capturan la incertidumbre.
• La validación cruzada proporciona una evaluación más robusta que una simple división en entrenamiento/prueba.
• El modelo más complejo no siempre es el mejor para predicción; a veces un modelo más simple generaliza mejor.

La capacidad de generar predicciones precisas y validar adecuadamente los modelos es lo que hace que la regresión lineal sea una herramienta tan valiosa en el análisis de datos y la toma de decisiones basada en evidencia.