El análisis de regresión lineal simple es una técnica estadística que se utiliza para modelar la relación entre una variable dependiente y una variable independiente. En este documento se presenta un ejemplo de cómo realizar un análisis de regresión lineal simple en R.
La ecuación de regresión lineal simple es:
\[Y = \alpha + \beta X + \epsilon\]
Donde:
\(Y\) es la variable dependiente.
\(X\) es la variable independiente.
\(\alpha\) es la intersección de la recta de regresión con el eje Y.
\(\beta\) es la pendiente de la recta de regresión.
\(\epsilon\) es el error aleatorio o residual.
Estimar la ecuación de regresión mediante estimadores de \(\alpha\) (a) y \(\beta\) (b) a partir de una muestra y generar intervalos de confianza para los mismos.
Estimar que tanto está la variable dependiente controlada por la variable independiente, en otras palabras que tanto la variación en y se explica por la variación en x.
Usar la ecuación de regresión para predecir valores de y a partir de valores de x.
Linealidad: La relación entre las variables es lineal.
Homocedasticidad: La varianza de los residuales es constante para los valores de x. Esto es equivalente a la homogeneidad de varianza en los grupos del ANOVA.
Normalidad: Los residuales siguen una distribución normal.
Independencia: Los residuales son independientes entre sí.
Para estimar los parámetros de la regresión, se utiliza el método de mínimos cuadrados. El objetivo es encontrar los valores de \(\alpha\) y \(\beta\) que minimizan la suma de los cuadrados de los residuos:
\[\sum_{i=1}^{n} (Y_i - \hat{Y}_i)^2\]
Donde:
Para calcular los coeficientes de regresión manualmente, se utilizan las siguientes fórmulas:
\[\beta = \frac{\sum{(X_i - \bar{X})(Y_i - \bar{Y})}}{\sum{(X_i - \bar{X})^2}}\]
\[\alpha = \bar{Y} - \beta \bar{X}\]
Donde:
Vamos a usar los datos de la longitud (cm) del ala derecha de gorriones de diferentes edades (días). En este caso, podemos asumir que la variable dependiente es la longitud del ala y la independiente, la edad de las aves.
edad <- c(3.0,4.0,5.0,6.0,8.0,9.0,10.0,11.0,12.0,14.0,15.0,16.0,17.0)
l.ala <- c(1.4,1.5,2.2,2.4,3.1,3.2,3.2,3.9,4.1,4.7,4.5,5.2,5.0)
sparrow <- data.frame(edad,l.ala)
#tabla de datos
knitr::kable(sparrow, col.names = c("Edad (días)", "Longitud del ala (cm)"), caption = 'Tabla 1. Longitud del ala derecha de gorriones de diferente edad')
Edad (días) | Longitud del ala (cm) |
---|---|
3 | 1.4 |
4 | 1.5 |
5 | 2.2 |
6 | 2.4 |
8 | 3.1 |
9 | 3.2 |
10 | 3.2 |
11 | 3.9 |
12 | 4.1 |
14 | 4.7 |
15 | 4.5 |
16 | 5.2 |
17 | 5.0 |
Gráfica de dispersión
# Gráfica de dispersión
library(ggplot2)
ggplot(sparrow, aes(x = edad, y = l.ala)) +
geom_point() +
labs(x = "Edad (días)", y = "Longitud del ala (cm)")
Figura 1. Gráfica de dispersión de la longitud del ala derecha de gorriones en función de la edad.
Cálculo de los coeficientes de regresión con R
Para calcular los coeficientes de regresión con R, se utiliza la
función lm()
.
# Ajuste del modelo de regresión lineal
modelo <- lm(l.ala ~ edad, data = sparrow)
# Coeficientes de regresión
coeficientes <- coef(modelo)
summary(modelo)
##
## Call:
## lm(formula = l.ala ~ edad, data = sparrow)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30699 -0.21538 0.06553 0.16324 0.22507
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.71309 0.14790 4.821 0.000535 ***
## edad 0.27023 0.01349 20.027 5.27e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2184 on 11 degrees of freedom
## Multiple R-squared: 0.9733, Adjusted R-squared: 0.9709
## F-statistic: 401.1 on 1 and 11 DF, p-value: 5.267e-10
Los coeficientes de regresión son \(\alpha = 0.7130945\) y \(\beta = 0.270229\).
Gráfica de la recta de regresión con la ecuación y el intervalo de confianza
# Gráfica de la recta de regresión con ecuación
eq <- substitute(
italic(y) == a + b %.% italic(x),
list(
a = round(as.numeric(coeficientes[1]), 3),
b = round(as.numeric(coeficientes[2]), 3)
)
)
# plot
ggplot(sparrow, aes(x = edad, y = l.ala)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE, color = "red") +
annotate("text", x = 5.5, y = 5, label = as.character(as.expression(eq)), parse = TRUE) +
labs(x = "Edad (días)", y = "Longitud del ala (cm)")
## `geom_smooth()` using formula = 'y ~ x'
Figura 2. Gráfica de dispersión de la longitud del ala derecha de gorriones en función de la edad con la recta de regresión, intervalo de confianza (95%) y la ecuación de la recta.
Intervalo de confianza para la pendiente (\(\beta\))
# Intervalo de confianza para la pendiente
confint(modelo, "edad", level = 0.95)
## 2.5 % 97.5 %
## edad 0.2405308 0.2999272
Prueba de hipótesis para la pendiente
Para realizar una prueba de hipótesis para la pendiente, se plantean las siguientes hipótesis:
# Prueba de hipótesis para la pendiente
summary(modelo)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7130945 0.14790445 4.821319 5.347485e-04
## edad 0.2702290 0.01349312 20.027169 5.267053e-10
El valor p para la pendiente es 5.2670535^{-10}.
Prueba del supuesto de heterocedasticidad
Para probar el supuesto de homocedasticidad, se puede realizar un análisis de los residuales.
# Análisis de los residuales
residuos <- residuals(modelo)
plot(modelo, which = 1)
# Prueba de Breusch-Pagan
library(lmtest)
bptest(modelo)
##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 1.6349, df = 1, p-value = 0.201
Aunque el gráfico de residuos vs. valores ajustados muestra cierta dispersión, la prueba de Breusch-Pagan no rechaza la hipótesis nula de homocedasticidad.
Predicción de la longitud del ala para una edad dada con intervalo de confianza
Para predecir la longitud del ala para una edad dada, se utiliza la
función predict()
.
# Predicción de la longitud del ala para una edad dada
nueva_edad <- 7
prediccion <- predict(modelo, newdata = data.frame(edad = nueva_edad), interval = "confidence", level = 0.95)
prediccion
## fit lwr upr
## 1 2.604698 2.444344 2.765051
Predicción de la edad para una longitud del ala dada con intervalo de confianza
Para predecir la edad para una longitud del ala dada, se utiliza la
función predict()
con el modelo invertido.
# Predicción de la edad para una longitud del ala dada
nueva_l_ala <- c(4.5, 3)
prediccion_inversa <- predict(lm(edad ~ l.ala, data = sparrow), newdata = data.frame(l.ala = nueva_l_ala), interval = "confidence", level = 0.95)
prediccion_inversa
## fit lwr upr
## 1 13.906551 13.257517 14.55558
## 2 8.503874 7.990108 9.01764
Si los supuestos de normalidad y homocedasticidad no se cumplen, se puede utilizar otras alternativas que incluyen:
Transformaciones de las variables (logarítmica, cuadrática, etc.).
Regresión polinómica.
Modelos no-paramétricos como la regresión de Loess.
# Regresión polinómica
modelo_polinomico <- lm(l.ala ~ poly(edad, 2), data = sparrow)
# Gráfica de la regresión polinómica
ggplot(sparrow, aes(x = edad, y = l.ala)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = TRUE, color = "green") +
labs(x = "Edad (días)", y = "Longitud del ala (cm)")
Figura 3. Gráfica de dispersión de la longitud del ala derecha de gorriones en función de la edad con la regresión polinómica y el intervalo de confianza (95%).
Coeficientes de la regresión polinómica:
# Coeficientes de la regresión polinómica
coeficientes_polinomico <- coef(modelo_polinomico)
coeficientes_polinomico
## (Intercept) poly(edad, 2)1 poly(edad, 2)2
## 3.4153846 4.3740386 -0.3683301
# p-value
p_value_polinomico <- summary(modelo_polinomico)$coefficients[2,4]
p_value_polinomico
## [1] 7.801604e-10
# R^2
r2_polinomico <- summary(modelo_polinomico)$r.squared
r2_polinomico
## [1] 0.9802084
# Regresión de Loess
modelo_loess <- loess(l.ala ~ edad, data = sparrow)
# Gráfica de la regresión de Loess
ggplot(sparrow, aes(x = edad, y = l.ala)) +
geom_point() +
geom_smooth(method = "loess", se = TRUE, color = "blue") +
labs(x = "Edad (días)", y = "Longitud del ala (cm)")
## `geom_smooth()` using formula = 'y ~ x'
Figura 4. Gráfica de dispersión de la longitud del ala derecha de gorriones en función de la edad con la regresión de Loess y el intervalo de confianza (95%).
La regresión múltiple es una extensión de la regresión lineal simple que permite modelar la relación entre una variable dependiente y múltiples variables independientes. La ecuación de regresión múltiple es: \[Y = \alpha + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_n X_n + \epsilon\] Donde: - \(Y\) es la variable dependiente. - \(X_1, X_2, ..., X_n\) son las variables independientes. - \(\alpha\) es la intersección de la recta de regresión con el eje Y. - \(\beta_1, \beta_2, ..., \beta_n\) son las pendientes de la recta de regresión. - \(\epsilon\) es el error aleatorio o residual.
# Cargar los datos
library(MASS)
library(readxl)
small_cities <- read_excel("death_small_cities.xlsx", sheet = "data")
# modelo
modelo_multiplo <- lm(death1K ~ income1K + doctor100K + hosp100K + density, data = small_cities)
# Resumen del modelo
summary(modelo_multiplo)
##
## Call:
## lm(formula = death1K ~ income1K + doctor100K + hosp100K + density,
## data = small_cities)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6404 -0.7904 0.3053 0.9164 2.7906
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.2662552 2.0201467 6.072 1.95e-07 ***
## income1K -0.3302302 0.2345518 -1.408 0.1656
## doctor100K 0.0073916 0.0069336 1.066 0.2917
## hosp100K 0.0005837 0.0007219 0.809 0.4228
## density -0.0094629 0.0048868 -1.936 0.0587 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.601 on 48 degrees of freedom
## Multiple R-squared: 0.1437, Adjusted R-squared: 0.07235
## F-statistic: 2.014 on 4 and 48 DF, p-value: 0.1075
# Selección de mejor predictor
library(MASS)
modelo_seleccion <- stepAIC(modelo_multiplo, direction = "both")
## Start: AIC=54.65
## death1K ~ income1K + doctor100K + hosp100K + density
##
## Df Sum of Sq RSS AIC
## - hosp100K 1 1.6763 124.75 53.369
## - doctor100K 1 2.9139 125.99 53.892
## <none> 123.07 54.652
## - income1K 1 5.0825 128.16 54.797
## - density 1 9.6144 132.69 56.639
##
## Step: AIC=53.37
## death1K ~ income1K + doctor100K + density
##
## Df Sum of Sq RSS AIC
## <none> 124.75 53.369
## - doctor100K 1 5.1882 129.94 53.529
## - income1K 1 6.1544 130.91 53.921
## + hosp100K 1 1.6763 123.07 54.652
## - density 1 8.3192 133.07 54.791
# Resumen del modelo
summary(modelo_seleccion)
##
## Call:
## lm(formula = death1K ~ income1K + doctor100K + density, data = small_cities)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7852 -0.6675 0.4010 0.9348 2.7495
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.565900 1.978836 6.350 6.74e-08 ***
## income1K -0.359140 0.230990 -1.555 0.1264
## doctor100K 0.009284 0.006504 1.428 0.1598
## density -0.008580 0.004746 -1.808 0.0768 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.596 on 49 degrees of freedom
## Multiple R-squared: 0.132, Adjusted R-squared: 0.0789
## F-statistic: 2.485 on 3 and 49 DF, p-value: 0.07163