Regresión Lineal
Diseño Experimental
Edimer David Jaramillo
2023-02-28
1 Bibliotecas
library(tidyverse) # Manipulación y visualización de datos
library(agridat) # Base de datos de ejemplo
library(ggpubr) # Gráfico cuantil-cuantil
library(yardstick) # Métricas de error: R-cuadrado y para el RMSE
library(plotly) # Gráficos interactivos
library(rsample)2 Datos de ejemplo
datos <- hernandez.nitrogen
datos %>% head()3 Distribuciones
4 ¿Existe normalidad?
ggqqplot(datos$yield)
datos %>%
ggplot(aes(x = yield)) +
geom_density()
5 Correlación
\[\rho_{(X,Y)} = \frac{Cov_{(X,Y)}}{\sigma_X\times\sigma_Y} = \frac{\sum_{i=1}^{n}(X_i-\mu_X)(Y_i-\mu_Y)}{\sigma_X\times\sigma_Y}\]
- Cálculo manual:
covarianza <- cov(x = datos$nitro, y = datos$yield)
desv_nitro <- sd(datos$nitro)
desv_yield <- sd(datos$yield)
covarianza / (desv_nitro * desv_yield)## [1] 0.7159447- Con la función “cor”:
cor(x = datos$nitro, y = datos$yield, method = "pearson")## [1] 0.71594476 Dispersión
datos %>%
ggplot(aes(x = nitro, y = yield)) +
geom_point()
7 Regresión lineal simple
- El propósito principal del análisis de regresión es estimar la función de regresión poblacional con base en la función de regresión muestral.
\[Y_i = \beta_0\ +\ \beta_1X\] \[\hat{Y_i} = \hat{\beta_0}\ + \hat{\beta_1}X_i\ + \hat{\epsilon_i}\]
- Función matemática:
\[Y_i = E(Y|X_i)\ +\ \epsilon_i\\\]
- Asumiendo que \(E(Y|X_i)\) es lineal en \(X_i\):
\[Y_i = E(Y|X_i)\ +\ \epsilon_i\\ Y_i = \beta_0 +\ \beta_1X_i\ +\ \epsilon_i\]
modelo <- lm(yield ~ nitro, data = datos)
summary(modelo)##
## Call:
## lm(formula = yield ~ nitro, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1752 -0.8803 0.1086 1.1148 3.4049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.872032 0.251752 35.24 <2e-16 ***
## nitro 0.023434 0.001974 11.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.674 on 134 degrees of freedom
## Multiple R-squared: 0.5126, Adjusted R-squared: 0.5089
## F-statistic: 140.9 on 1 and 134 DF, p-value: < 2.2e-168 Coeficientes
\[Y_i = 8.87203185 + 0.02343431 \times Dosis\]
8.1 Estimación puntual
modelo %>% coef()## (Intercept) nitro
## 8.87203185 0.023434318.2 Intervalos
modelo %>% confint(level = 0.95)## 2.5 % 97.5 %
## (Intercept) 8.37411051 9.36995320
## nitro 0.01952984 0.027338779 Predicciones
9.1 Sin intervalo
- ¿Son correctas estas predicciones?
nuevos_datos <- data.frame(nitro = c(122.4, 12.5, 33.8))
predict(modelo, newdata = nuevos_datos)## 1 2 3
## 11.740391 9.164961 9.6641119.2 Intervalo de confianza
predict(modelo, newdata = nuevos_datos, interval = "confidence")## fit lwr upr
## 1 11.740391 11.448177 12.032605
## 2 9.164961 8.706284 9.623637
## 3 9.664111 9.267389 10.06083410 ¿Residuales del modelo?
10.1 Supuestos matemáticos
- Normalidad de los residuales con valor medio igual a cero: \(E(\epsilon_i|X_i) = 0\).
- Homocedasticidad o varianza constante de los errores \(\epsilon_i\).
- Linealidad en los parámetros.
- Independencia de los errores (autocorrelación): \(cov(\epsilon_i, \epsilon_j)=0\).
\[\epsilon_i\ \overset{\text{i.i.d.}}\sim\ N(\mu = 0,\ \sigma)\]
10.2 Residuales
- ¿Real vs Predicho?
yield_real <- datos$yield
yield_pred <- modelo$fitted.values
residuales <- yield_real - yield_pred
residuales %>% head()## 1 2 3 4 5 6
## -5.1751819 -4.9939019 -4.4833319 -4.2326419 -2.7606746 -0.8994646- Podemos obtenerlos del modelo a través de la función “residuals”:
residuals(modelo) %>% head()## 1 2 3 4 5 6
## -5.1751819 -4.9939019 -4.4833319 -4.2326419 -2.7606746 -0.899464610.3 Estimados vs Reales
df_real_estimado <- bind_cols(real = yield_real, estimado = yield_pred)
df_real_estimado %>% head()- Graficamos los reales vs los estimados:
df_real_estimado %>%
ggplot(aes(x = real, y = estimado)) +
geom_point()
10.4 Normalidad de residuales
ggqqplot(residuales)
- También podemos graficar la densidad de los residuales:
residuales %>%
enframe(value = "residual") %>%
ggplot(aes(x = residual)) +
geom_density()
- ¿Cuál es el promedio de los residuales?
residuales %>% mean()## [1] 1.232164e-1710.5 Homocedasticidad
df_real_estimado %>%
mutate(residual = residuales) %>%
ggplot(aes(x = estimado, y = residual)) +
geom_point() +
geom_hline(yintercept = 0, color = "red")
11 ¿Desempeño del modelo?
11.1 RMSE
\[ RMSE = \sqrt{\frac{\sum_{i=1}^{n}(y_i - \hat{y_i})^2}{n}} \]
Donde \(y_i\) es el valor real y \(\hat{y_i}\) es el valor predicho o estimado por el modelo.
sqrt(sum((yield_real - yield_pred)^2) / length(yield_real))## [1] 1.662047- Podemos extraer el RMSE con la función “rmse_vec” de la biblioteca yardstick:
rmse_vec(truth = yield_real, estimate = yield_pred)## [1] 1.66204711.2 R-Cuadrado
\[ R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y})^2}{\sum_{i = 1}^{n} (y_i - \bar{y})^2} \]
1 - ( sum((yield_real - yield_pred)^2) / sum((yield_real - mean(yield_real))^2) )## [1] 0.5125768- Podemos extraer el coeficiente de determinación (r-cuadrado) con la función “rsq_vec” de la biblioteca yardstick:
rsq_vec(truth = yield_real, estimate = yield_pred)## [1] 0.5125768- El coeficiente de determinación es igual al coeficiente de correlación de pearson al cuadrado:
cor(x = datos$nitro, y = datos$yield, method = "pearson")^2## [1] 0.512576812 LM con ggplot2
datos %>%
ggplot(aes(x = nitro, y = yield)) +
geom_point() +
geom_smooth(method = "lm")
13 ¿Cómo mejorar el modelo?
13.1 Transformación logarítmica
datos %>%
ggplot(aes(x = nitro, y = log(yield))) +
geom_point() +
geom_smooth(method = "lm")
13.2 Modelo con logaritmos
modelo_log <- lm(log(yield) ~ nitro, data = datos)
modelo_log %>% summary()##
## Call:
## lm(formula = log(yield) ~ nitro, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.84331 -0.07030 0.03508 0.11883 0.28701
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1507960 0.0283609 75.84 <2e-16 ***
## nitro 0.0023708 0.0002224 10.66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1886 on 134 degrees of freedom
## Multiple R-squared: 0.4589, Adjusted R-squared: 0.4549
## F-statistic: 113.6 on 1 and 134 DF, p-value: < 2.2e-16- La inversa del logaritmo natural es la exponencial:
exp(2.1507960) # B0## [1] 8.591695exp(0.0023708) # B1## [1] 1.00237413.3 RMSE Y R-Cuadrado
yield_pred_log <- modelo_log$fitted.values %>% exp()
yield_real_log <- yield_real
rmse_vec(truth = yield_real_log, estimate = yield_pred_log)## [1] 1.78757rsq_vec(truth = yield_real_log, estimate = yield_pred_log)## [1] 0.461094413.4 Estimados vs Reales
df_real_estimado_log <- bind_cols(real = yield_real_log, estimado = yield_pred_log)
df_real_estimado_log %>% head()- Graficamos los reales vs los estimados:
df_real_estimado_log %>%
ggplot(aes(x = real, y = estimado)) +
geom_point()
13.5 Normalidad de residuales
ggqqplot(residuals(modelo_log))
13.6 Homocedasticidad
df_real_estimado_log %>%
mutate(residual = residuals(modelo_log)) %>%
ggplot(aes(x = estimado, y = residual)) +
geom_point() +
geom_hline(yintercept = 0, color = "red")
14 ¿Cómo comparar dos modelos?
14.1 Métricas
- RMSE
- R-Cuadrado
- Criterios de información estadística:
- \(LogLik\): logaritmo de la verosimilitud
- \(AIC = -2 \times logLik + k + n\)
- \(BIC = AIC,\ con\ k = log(n)\)
14.2 Predicciones simuladas
- Vamos a crear un vector de 100 valores hipotéticos de posibles dosis de nitrogeno entre 0 y 300 kg/ha:
dosis_nitrogeno <- seq(from = 0, to = 300, length.out = 100)
dosis_nitrogeno %>% head()## [1] 0.000000 3.030303 6.060606 9.090909 12.121212 15.151515- Incorporamos estos nuevos datos a un dataframe para luego hacer predicciones:
nuevas_dosis <- data.frame(nitro = dosis_nitrogeno)
nuevas_dosis %>% head()- Ahora hagamos la predicción con ambos modelos, el modelo lineal original y el modelo lineal con transformación logarítmica:
pred_original <- predict(modelo, newdata = nuevas_dosis)
pred_logaritmos <- predict(modelo_log, newdata = nuevas_dosis) %>% exp()- Incorporamos las predicciones en la tabla de “nuevas_dosis” y graficamos las predicciones realizadas por ambos modelos:
ggplotly(
nuevas_dosis %>%
mutate(pred_orig = pred_original,
pred_log = pred_logaritmos) %>%
pivot_longer(cols = -nitro, names_to = "modelo", values_to = "pred_kg") %>%
ggplot(aes(x = nitro, y = pred_kg, color = modelo)) +
geom_line()
)15 Solución manual
\[ (X^TX)^{-1}X^TY \]
- Primero conformamos la matriz X:
matriz_x <- cbind(1, datos$nitro)
matriz_x %>% head()## [,1] [,2]
## [1,] 1 0.0
## [2,] 1 0.0
## [3,] 1 0.0
## [4,] 1 0.0
## [5,] 1 33.6
## [6,] 1 33.6- Cuando ejecutamos “lm” en R automáticamente se construye la matriz anterior, denominada matriz de diseño:
model.matrix(yield ~ nitro, data = datos) %>% head()## (Intercept) nitro
## 1 1 0.0
## 2 1 0.0
## 3 1 0.0
## 4 1 0.0
## 5 1 33.6
## 6 1 33.6- Definimos el vector Y:
vector_y <- datos$yield
vector_y %>% head()## [1] 3.69685 3.87813 4.38870 4.63939 6.89875 8.75996- Transponemos la matriz X:
x_transpuesta <- matriz_x %>% t()
x_transpuesta %>% head()## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1 1 1 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
## [2,] 0 0 0 0 33.6 33.6 33.6 33.6 67.2 67.2 67.2 67.2 100.8 100.8
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## [1,] 1.0 1.0 1.0 1.0 1.0 1.0 1 1 1 1 1.0 1.0
## [2,] 100.8 100.8 134.4 134.4 134.4 134.4 168 168 168 168 201.6 201.6
## [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## [1,] 1.0 1.0 1 1 1 1 1.0 1.0 1.0 1.0 1.0 1.0
## [2,] 201.6 201.6 0 0 0 0 33.6 33.6 33.6 33.6 67.2 67.2
## [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,] 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1
## [2,] 67.2 67.2 100.8 100.8 100.8 100.8 134.4 134.4 134.4 134.4 168 168
## [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62]
## [1,] 1 1 1.0 1.0 1.0 1.0 1 1 1 1 1.0 1.0
## [2,] 168 168 201.6 201.6 201.6 201.6 0 0 0 0 33.6 33.6
## [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74]
## [1,] 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
## [2,] 33.6 33.6 67.2 67.2 67.2 67.2 100.8 100.8 100.8 100.8 134.4 134.4
## [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86]
## [1,] 1.0 1.0 1 1 1 1 1.0 1.0 1.0 1.0 1 1
## [2,] 134.4 134.4 168 168 168 168 201.6 201.6 201.6 201.6 0 0
## [,87] [,88] [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98]
## [1,] 1 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
## [2,] 0 0 44.8 44.8 44.8 44.8 89.6 89.6 89.6 89.6 134.4 134.4
## [,99] [,100] [,101] [,102] [,103] [,104] [,105] [,106] [,107] [,108]
## [1,] 1.0 1.0 1.0 1.0 1.0 1.0 1 1 1 1
## [2,] 134.4 134.4 179.2 179.2 179.2 179.2 224 224 224 224
## [,109] [,110] [,111] [,112] [,113] [,114] [,115] [,116] [,117] [,118]
## [1,] 1.0 1.0 1.0 1.0 1 1 1 1 1.0 1.0
## [2,] 268.8 268.8 268.8 268.8 0 0 0 0 33.6 33.6
## [,119] [,120] [,121] [,122] [,123] [,124] [,125] [,126] [,127] [,128]
## [1,] 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
## [2,] 33.6 33.6 67.2 67.2 67.2 67.2 100.8 100.8 100.8 100.8
## [,129] [,130] [,131] [,132] [,133] [,134] [,135] [,136]
## [1,] 1.0 1.0 1.0 1.0 1 1 1 1
## [2,] 134.4 134.4 134.4 134.4 168 168 168 168- Calculamos la multiplicación de \(X^TX\):
x_transx <- x_transpuesta %*% matriz_x
x_transx %>% head()## [,1] [,2]
## [1,] 136.0 14246.4
## [2,] 14246.4 2211758.1- Obtenemos la inversa de \(X^TX\):
inversa_x_tranx <- x_transx %>% solve()
inversa_x_tranx## [,1] [,2]
## [1,] 0.0226060556 -1.456104e-04
## [2,] -0.0001456104 1.390036e-06- Finalmente multiplicamos la transpuesta de X por Y \((X^TY)\):
transx_y <- x_transpuesta %*% vector_y
transx_y## [,1]
## [1,] 1540.451
## [2,] 178225.535- Obtenemos \((X^TX)^{-1}X^TY\):
betas <- inversa_x_tranx %*% transx_y
betas## [,1]
## [1,] 8.87203185
## [2,] 0.0234343116 ¿Validación cruzada?
16.1 Train - Test
- Podemos usar la función “initial_split” para indicar cómo queremos hacer la partición. Luego usamos las funciones “training” y “testing” para obtener los dos dataframe que usaremos para estimar los parámetros del modelo.
set.seed(2023)
particion <- initial_split(data = datos, prop = 0.7, strata = "yield")
datos_train <- training(particion)
datos_test <- testing(particion)- ¿Cómo quedan divididos los datos?
particion## <Training/Testing/Total>
## <92/44/136>- ¿Son parecidos los promedios de la variable “yield” en ambos conjuntos de datos”?
datos_train %>%
summarise(mean(yield))datos_test %>%
summarise(mean(yield))16.2 K-Fold
k_fold <- vfold_cv(data = datos_train, v = 5, strata = "yield")
k_fold- Imprimimos el primer conjunto de datos de los 5 existentes:
k_fold$splits[[1]] %>% training()