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)     # Remuestreo - Validación cruzada
library(ggeffects)   # Gráficos de estimaciones de modelos

2 Datos de ejemplo

datos <- hernandez.nitrogen
datos %>% head()

3 Exploratorio

3.1 Promedios

datos %>% 
  group_by(site) %>% 
  summarise(promedio = mean(yield))

3.2 Sitios

datos %>% 
  ggplot(aes(x = site, y = yield)) +
  geom_boxplot()

3.3 Dosis y Sitios (1)

datos %>% 
  ggplot(aes(x = nitro, y = yield, color = site)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

3.4 Dosis y Sitios (2)

datos %>% 
  ggplot(aes(x = nitro, y = yield)) +
  facet_wrap(~site, scales = "free") +
  geom_point() +
  geom_smooth(method = "lm")

4 Regresión lineal múltiple

\[Y_i = \beta_0\ +\ \beta_1X_1 + \beta_{2}X_2\]

  • ¿Por qué no aparece el site1?
modelo1 <- lm(yield ~ nitro, data = datos)
modelo2 <- lm(yield ~ site, data = datos)
modelo3 <- lm(yield ~ nitro + site, data = datos)
summary(modelo3)
## 
## Call:
## lm(formula = yield ~ nitro + site, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4665 -0.9888  0.0308  0.9561  3.2666 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.163309   0.306519  23.370  < 2e-16 ***
## nitro       0.024889   0.001649  15.095  < 2e-16 ***
## siteS2      1.812934   0.364226   4.977 2.00e-06 ***
## siteS3      2.072701   0.364226   5.691 8.00e-08 ***
## siteS4      1.075991   0.368415   2.921  0.00412 ** 
## siteS5      3.030779   0.380109   7.973 6.92e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.363 on 130 degrees of freedom
## Multiple R-squared:  0.6867, Adjusted R-squared:  0.6747 
## F-statistic:    57 on 5 and 130 DF,  p-value: < 2.2e-16

5 ¿Variables dummy o indicadoras?

mtx_modelo <- model.matrix(yield ~ nitro + site, data = datos)
mtx_modelo %>% head()
##   (Intercept) nitro siteS2 siteS3 siteS4 siteS5
## 1           1   0.0      0      0      0      0
## 2           1   0.0      0      0      0      0
## 3           1   0.0      0      0      0      0
## 4           1   0.0      0      0      0      0
## 5           1  33.6      0      0      0      0
## 6           1  33.6      0      0      0      0

5.1 Ecuación

La ecucación del modelo de regresión lineal múltiple con variables categóricas puede ser representado de la siguiente manera:

\[ Y_i = 7.163309 +\ (0.024889 \times Nitro_i) + (1.812934 \times SiteS2_i) + \\ (2.072701 \times SiteS3_i) + (1.075991 \times SiteS4_i) + (3.030779 \times SiteS5_i) \]

5.2 Predicción manual

  • Ejemplo de predicción para una dosis de nitrógeno de 30 kg/ha para el sitio 5.
7.163309 + (0.024889 * 30) + (1.812934 * 0 ) + (2.072701 * 0) + (1.075991 * 0) + (3.030779 * 1)
## [1] 10.94076

5.3 Predicción con R

predict(modelo3, newdata = data.frame(nitro = 30, site = "S5"))
##        1 
## 10.94075

5.4 Residuales

  • ¿Real vs Predicho?
yield_real <- datos$yield
yield_pred <- modelo3$fitted.values
residuales <- residuals(modelo3)
residuales %>% head()
##          1          2          3          4          5          6 
## -3.4664585 -3.2851785 -2.7746085 -2.5239185 -1.1008243  0.7603857

5.5 Estimados vs Reales

df_real_estimado <- bind_cols(real = yield_real, estimado = yield_pred)

df_real_estimado %>% 
  ggplot(aes(x = real, y = estimado)) +
  geom_point()

5.6 Normalidad de residuales

ggqqplot(residuales)

  • También podemos graficar la densidad de los residuales:
residuales %>% 
  enframe(value = "residual") %>% 
  ggplot(aes(x = residual)) +
  geom_density()

5.7 Homocedasticidad

df_real_estimado %>% 
  mutate(residual = residuales) %>% 
  ggplot(aes(x = estimado, y = residual)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red")

5.8 ¿Cuál es mejor modelo?

predichos_mod1 <- modelo1$fitted.values
predichos_mod2 <- modelo2$fitted.values
predichos_mod3 <- modelo3$fitted.values

rmse_vec(truth = yield_real, estimate = predichos_mod1)
## [1] 1.662047
rmse_vec(truth = yield_real, estimate = predichos_mod2)
## [1] 2.210637
rmse_vec(truth = yield_real, estimate = predichos_mod3)
## [1] 1.332407

6 Estimaciones del modelo

  • Podemos usar la biblioteca ggeffects para graficar las estimaciones de los modelos a través de la función ggpredict

6.1 Modelo 1

modelo1 %>%
  ggpredict(terms = c("nitro")) %>%
  plot()

6.2 Modelo 2

modelo2 %>%
  ggpredict(terms = c("site")) %>%
  plot()

6.3 Modelo 3

modelo3 %>%
  ggpredict(terms = c("nitro",  "site")) %>%
  plot(facet = FALSE)

7 Regresión polinomial

\[Y_i = \beta_0\ +\ \beta_1X_i + \beta_{2}X_i^2 + \beta_{3}X_i^3 + ... + \beta_{n}X_i^p\]

  • Podemos agregar un polinomio de \(p\) grado con nuestras variables numéricas predictoras de tres formas:
    • Generando la variable manualmente
    • Usando la función I(x^p)
    • Usando la función poly(x, raw = TRUE)

7.1 p=2

# Nota: los resultados son equivalentes
modelo4 <- lm(yield ~ nitro + I(nitro^2) + site, data = datos)

modelo4.2 <- lm(yield ~ nitro + nitro2 + site,
              data = datos %>% mutate(nitro2 = nitro^2))

modelo4.3 <- lm(yield ~ poly(nitro, degree = 2, raw = TRUE) + site,
              data = datos)

# Nota: los resultados son equivalentes
modelo4 %>% 
  summary()
## 
## Call:
## lm(formula = yield ~ nitro + I(nitro^2) + site, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.24723 -0.69886  0.00713  0.76629  2.53071 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.9440772  0.2698833  22.025  < 2e-16 ***
## nitro        0.0607110  0.0040100  15.140  < 2e-16 ***
## I(nitro^2)  -0.0001630  0.0000173  -9.421 2.34e-16 ***
## siteS2       1.8129336  0.2814192   6.442 2.14e-09 ***
## siteS3       2.0727014  0.2814192   7.365 1.86e-11 ***
## siteS4       1.7325330  0.2930622   5.912 2.85e-08 ***
## siteS5       2.9273633  0.2938964   9.961  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.053 on 129 degrees of freedom
## Multiple R-squared:  0.8144, Adjusted R-squared:  0.8058 
## F-statistic: 94.36 on 6 and 129 DF,  p-value: < 2.2e-16
modelo4.2 %>% 
  summary()
## 
## Call:
## lm(formula = yield ~ nitro + nitro2 + site, data = datos %>% 
##     mutate(nitro2 = nitro^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.24723 -0.69886  0.00713  0.76629  2.53071 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.9440772  0.2698833  22.025  < 2e-16 ***
## nitro        0.0607110  0.0040100  15.140  < 2e-16 ***
## nitro2      -0.0001630  0.0000173  -9.421 2.34e-16 ***
## siteS2       1.8129336  0.2814192   6.442 2.14e-09 ***
## siteS3       2.0727014  0.2814192   7.365 1.86e-11 ***
## siteS4       1.7325330  0.2930622   5.912 2.85e-08 ***
## siteS5       2.9273633  0.2938964   9.961  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.053 on 129 degrees of freedom
## Multiple R-squared:  0.8144, Adjusted R-squared:  0.8058 
## F-statistic: 94.36 on 6 and 129 DF,  p-value: < 2.2e-16
modelo4.3 %>% 
  summary()
## 
## Call:
## lm(formula = yield ~ poly(nitro, degree = 2, raw = TRUE) + site, 
##     data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.24723 -0.69886  0.00713  0.76629  2.53071 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                           5.9440772  0.2698833  22.025  < 2e-16 ***
## poly(nitro, degree = 2, raw = TRUE)1  0.0607110  0.0040100  15.140  < 2e-16 ***
## poly(nitro, degree = 2, raw = TRUE)2 -0.0001630  0.0000173  -9.421 2.34e-16 ***
## siteS2                                1.8129336  0.2814192   6.442 2.14e-09 ***
## siteS3                                2.0727014  0.2814192   7.365 1.86e-11 ***
## siteS4                                1.7325330  0.2930622   5.912 2.85e-08 ***
## siteS5                                2.9273633  0.2938964   9.961  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.053 on 129 degrees of freedom
## Multiple R-squared:  0.8144, Adjusted R-squared:  0.8058 
## F-statistic: 94.36 on 6 and 129 DF,  p-value: < 2.2e-16
  • Grafiquemos la estimación del modelo:
modelo4 %>%
  ggpredict(terms = c("nitro",  "site")) %>%
  plot(facet = FALSE)

7.2 p=3

modelo5 <- lm(yield ~ nitro + I(nitro^2) + + I(nitro^3) + site, data = datos)

8 ¿Cuál es el mejor modelo?

predichos_mod4 <- modelo4$fitted.values
predichos_mod5 <- modelo5$fitted.values

rmse_vec(truth = yield_real, estimate = predichos_mod1)
## [1] 1.662047
rmse_vec(truth = yield_real, estimate = predichos_mod2)
## [1] 2.210637
rmse_vec(truth = yield_real, estimate = predichos_mod3)
## [1] 1.332407
rmse_vec(truth = yield_real, estimate = predichos_mod4)
## [1] 1.025518
rmse_vec(truth = yield_real, estimate = predichos_mod5)
## [1] 0.9871166

8.1 Residuales

  • ¿Real vs Predicho?
yield_pred_poli <- modelo5$fitted.values
residuales_poli <- residuals(modelo5)
residuales_poli %>% head()
##          1          2          3          4          5          6 
## -1.9817800 -1.8005000 -1.2899300 -1.0392400 -1.1295138  0.7316962

8.2 Estimados vs Reales

df_real_estimado_poli <- bind_cols(real = yield_real, estimado = yield_pred_poli)

df_real_estimado_poli %>% 
  ggplot(aes(x = real, y = estimado)) +
  geom_point()

8.3 Normalidad de residuales

ggqqplot(residuales_poli)

  • También podemos graficar la densidad de los residuales:
residuales_poli %>% 
  enframe(value = "residual") %>% 
  ggplot(aes(x = residual)) +
  geom_density()

8.4 Homocedasticidad

df_real_estimado_poli %>% 
  mutate(residual = residuales) %>% 
  ggplot(aes(x = estimado, y = residual)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red")

8.5 Estimaciones

modelo5 %>%
  ggpredict(terms = c("nitro",  "site")) %>%
  plot(facet = FALSE)