Ajuste de modelos con R

Bibliotecas R

Code

library(tidyverse)  # Manipulación de datos y gráficosos
library(janitor)    # Edición automática de nombreses
library(yardstick)  # Métricas para evaluar modeloss
library(ggeffects)  # Gráficos de efectos estimados en modelo lineall
library(FactoMineR) # Métodos multivariados  
library(factoextra) # Gráficos para métodos multivariadosos

Preguntas pendientes

Distancia vs tiempo

• Problema 1: si una persona está interesada en estimar (o predecir) el tiempo que le tomaría llegar a la universidad desde un lugar ubicado a una distancia de 7.83 km, ¿Qué modelo le recomendaría a esta persona?
• Problema 2: la misma persona del problema 1 ahora nos menciona que se movilizará hasta la universidad en motocicleta, ¿cambia en algo su modelo? ¿Cuántos modelos son posibles?

Análisis multivariado

• ¿Cómo podemos analizar todas las variables (de frutas y vegetales) juntas? ¿Es posible?

Datos

Encuesta

Code

encuesta <- read_csv("EncuestaDepurada.csv")
encuesta %>% head()

Frutas y vegetales

Code

frutas <- read_csv("fruits.csv", na = "-") %>%
  clean_names() %>%
  mutate(grupo = "Frutas") %>%
  select(-c(energy_kcal_k_j))

vegetales <- read_csv("vegetables.csv", na = "-") %>%
  clean_names() %>%
  mutate(grupo = "Vegetales") %>%
  select(-c(energy_kcal_k_j))

frutas_vegetales <- bind_rows(frutas, vegetales) %>% 
  filter(!is.na(sugars_g)) %>% 
  filter(!is.na(vitamin_a_iu)) %>% 
  filter(!is.na(vitamin_b5_mg)) %>% 
  filter(!is.na(vitamin_e_mg))

frutas_vegetales %>% head()

Modelo lineal

Modelo 1

• Problema 1: si una persona está interesada en estimar (o predecir) el tiempo que le tomaría llegar a la universidad desde un lugar ubicado a una distancia de 7.83 km, ¿Qué modelo le recomendaría a esta persona?

\[y = mx + b\] \[y = \beta_1x + \beta_0\] \[\hat{y} = \hat{\beta_1}x + \hat{\beta_0} + \hat{\epsilon} \\ \]

Gráfico de dispersión

Code

encuesta %>% 
  filter(tiempo_casa_u < 750) %>% 
  ggplot(aes(x = distancia_casa_u, y = tiempo_casa_u)) +
  geom_point(size = 2) 

Ajuste

Code

# Datos para modelo
encuesta_mod <- encuesta %>% 
  filter(tiempo_casa_u < 750)

# Ajuste del modelo
modelo1 <- lm(tiempo_casa_u ~ distancia_casa_u,
              data = encuesta)

modelo2 <- lm(tiempo_casa_u ~ distancia_casa_u,
              data = encuesta_mod)

• ¿El resultado de los modelos es diferente?

Code

modelo1 %>% summary()

Call:
lm(formula = tiempo_casa_u ~ distancia_casa_u, data = encuesta)

Residuals:
   Min     1Q Median     3Q    Max 
-57.48 -50.70 -30.31 -13.70 838.70 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)
(Intercept)       62.7513    39.2581   1.598    0.122
distancia_casa_u  -0.9055     3.3508  -0.270    0.789

Residual standard error: 165.9 on 27 degrees of freedom
Multiple R-squared:  0.002698,  Adjusted R-squared:  -0.03424 
F-statistic: 0.07304 on 1 and 27 DF,  p-value: 0.789

Code

modelo2 %>% summary()

Call:
lm(formula = tiempo_casa_u ~ distancia_casa_u, data = encuesta_mod)

Residuals:
    Min      1Q  Median      3Q     Max 
-13.241 -11.400  -5.928   9.694  25.146 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       17.5930     3.0734   5.724 5.03e-06 ***
distancia_casa_u   1.1306     0.2578   4.385  0.00017 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.68 on 26 degrees of freedom
Multiple R-squared:  0.4251,    Adjusted R-squared:  0.403 
F-statistic: 19.23 on 1 and 26 DF,  p-value: 0.0001704

¿Mejor modelo?

• ¿Valores predichos vs valores reales? Raíz del cuadrado medio del error (RMSE):

\[RMSE = \sqrt{\frac{\sum_{i = 1}^{n}(y_i-\hat{y_i})^2}{n}}\]

Code

reales_mod1 <- encuesta$tiempo_casa_u
reales_mod2 <- encuesta_mod$tiempo_casa_u

predichos_mod1 <- predict(modelo1, newdata = encuesta)
predichos_mod2 <- predict(modelo2, newdata = encuesta_mod)

rmse_mod1 <- rmse_vec(truth = reales_mod1, estimate = predichos_mod1)
rmse_mod2 <- rmse_vec(truth = reales_mod2, estimate = predichos_mod2)

• RMSE del modelo 1:

Code

rmse_mod1

[1] 160.0359

• RMSE del modelo 2:

Code

rmse_mod2

[1] 12.21476

Predicciones

\[\hat{y} = \hat{\beta_1}x + \hat{\beta_0} + \hat{\epsilon} \\ \] \[\hat{y} = 1.1306x + 17.5930 \\ \]

Code

datos_nuevos <- data.frame(distancia_casa_u = c(7.83))

predict(modelo2, newdata = datos_nuevos)

       1 
26.44556 

• Manualmente:

Code

(1.1306 * 7.83) + 17.5930

[1] 26.4456

Modelo 2

• Problema 2: la misma persona del problema 1 ahora nos menciona que se movilizará hasta la universidad en motocicleta, ¿cambia en algo su modelo? ¿Cuántos modelos son posibles?

\[\hat{y} = \hat{\beta_0} + \hat{\beta_1}x_1 + \hat{\beta_2}x_2 + ... +\hat{\beta_k}x_k + \hat{\epsilon} \\ \]

Gráfico de dispersión

Code

encuesta %>% 
  filter(tiempo_casa_u < 750) %>% 
  ggplot(aes(x = distancia_casa_u, y = tiempo_casa_u, color = medio_transporte)) +
  geom_point(size = 2) 

Ajuste

Code

# Ajuste del modelo
modelo3 <- lm(tiempo_casa_u ~ distancia_casa_u + medio_transporte,
              data = encuesta_mod)

• ¿El resultado de los modelos es diferente?

Code

modelo3 %>% summary()

Call:
lm(formula = tiempo_casa_u ~ distancia_casa_u + medio_transporte, 
    data = encuesta_mod)

Residuals:
    Min      1Q  Median      3Q     Max 
-18.867  -5.421  -0.265   3.134  31.805 

Coefficients:
                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                         14.7719     4.8083   3.072 0.005225 ** 
distancia_casa_u                     1.0234     0.2306   4.438 0.000173 ***
medio_transporteMotocicleta         -3.6239     6.0051  -0.603 0.551861    
medio_transporteTransporte público  11.7717     6.0244   1.954 0.062451 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.73 on 24 degrees of freedom
Multiple R-squared:  0.6195,    Adjusted R-squared:  0.5719 
F-statistic: 13.02 on 3 and 24 DF,  p-value: 2.99e-05

¿Mejor modelo?

• ¿Valores predichos vs valores reales? Raíz del cuadrado medio del error (RMSE):

\[RMSE = \sqrt{\frac{\sum_{i = 1}^{n}(y_i-\hat{y_i})^2}{n}}\]

Code

reales_mod3 <- encuesta_mod$tiempo_casa_u

predichos_mod3 <- predict(modelo3, newdata = encuesta_mod)

rmse_mod3 <- rmse_vec(truth = reales_mod3, estimate = predichos_mod3)

• RMSE de los modelos 1, 2 y 3:

Code

rmse_mod1

[1] 160.0359

Code

rmse_mod2

[1] 12.21476

Code

rmse_mod3

[1] 9.937715

Predicciones

\[\hat{y} = \hat{\beta_0} + \hat{\beta_1}x_1 + \hat{\beta_2}x_2 + ... +\hat{\beta_k}x_k + \hat{\epsilon} \\ \]

Code

coefficients(modelo3)

                       (Intercept)                   distancia_casa_u 
                         14.771884                           1.023430 
       medio_transporteMotocicleta medio_transporteTransporte público 
                         -3.623867                          11.771737 

Code

datos_nuevos2 <- data.frame(distancia_casa_u = 7.83, medio_transporte = "Motocicleta")

predict(modelo3, newdata = datos_nuevos2)

       1 
19.16147 

• Manualmente:

Code

14.771884 + ((1.023430 * 7.83) + (-3.623867 * 1) + (11.771737 * 0))                         

[1] 19.16147

Estimación con modelo final

• Efecto marginal de la variable distancia:

Code

graficos_modelo <-modelo3 %>% 
  ggpredict() %>% 
  plot()

graficos_modelo$distancia_casa_u

• Efecto marginal de la variable medio de transporte:

Code

graficos_modelo$medio_transporte

Análisis de componentes principales (PCA)

Idea intuitiva

Ajuste y resumen

Code

acp <- PCA(X = frutas_vegetales %>% select(-name),
           quali.sup = 21,
           graph = FALSE)

summary(acp)

Call:
PCA(X = frutas_vegetales %>% select(-name), quali.sup = 21, graph = FALSE) 


Eigenvalues
                       Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6   Dim.7
Variance               7.537   3.235   1.821   1.381   1.254   0.893   0.781
% of var.             37.685  16.176   9.105   6.907   6.268   4.465   3.905
Cumulative % of var.  37.685  53.861  62.966  69.872  76.140  80.605  84.510
                       Dim.8   Dim.9  Dim.10  Dim.11  Dim.12  Dim.13  Dim.14
Variance               0.660   0.488   0.405   0.394   0.309   0.233   0.209
% of var.              3.300   2.441   2.025   1.969   1.546   1.165   1.045
Cumulative % of var.  87.810  90.251  92.277  94.246  95.792  96.957  98.002
                      Dim.15  Dim.16  Dim.17  Dim.18  Dim.19  Dim.20
Variance               0.139   0.111   0.063   0.046   0.030   0.010
% of var.              0.697   0.554   0.317   0.230   0.149   0.052
Cumulative % of var.  98.699  99.253  99.570  99.799  99.948 100.000

Individuals (the 10 first)
                    Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
1               |  3.231 | -2.603  0.737  0.649 | -1.322  0.443  0.167 |  0.066
2               |  1.908 | -1.200  0.157  0.396 | -0.805  0.164  0.178 | -0.281
3               |  7.032 |  3.607  1.415  0.263 | -1.787  0.809  0.065 | -2.267
4               |  3.236 | -0.176  0.003  0.003 | -1.556  0.614  0.231 |  1.334
5               |  2.326 | -0.744  0.060  0.102 | -0.847  0.182  0.133 | -0.417
6               |  2.732 | -1.874  0.382  0.471 | -1.408  0.502  0.266 | -0.026
7               |  2.764 | -2.038  0.452  0.544 | -1.182  0.354  0.183 | -0.895
8               |  2.452 | -0.060  0.000  0.001 | -1.289  0.421  0.276 |  0.873
9               |  2.517 | -1.579  0.271  0.394 | -1.460  0.540  0.337 |  0.637
10              |  2.885 | -1.620  0.285  0.315 | -1.226  0.381  0.180 | -0.494
                   ctr   cos2  
1                0.002  0.000 |
2                0.035  0.022 |
3                2.312  0.104 |
4                0.801  0.170 |
5                0.078  0.032 |
6                0.000  0.000 |
7                0.360  0.105 |
8                0.343  0.127 |
9                0.182  0.064 |
10               0.110  0.029 |

Variables (the 10 first)
                   Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
water_g         | -0.773  7.933  0.598 |  0.526  8.543  0.276 | -0.222  2.703
protein_g       |  0.901 10.759  0.811 | -0.158  0.772  0.025 | -0.024  0.032
total_fat_g     |  0.751  7.474  0.563 | -0.280  2.427  0.079 | -0.500 13.712
carbohydrates_g |  0.296  1.160  0.087 | -0.486  7.315  0.237 |  0.606 20.153
fiber_g         |  0.629  5.256  0.396 | -0.289  2.587  0.084 |  0.305  5.107
sugars_g        | -0.045  0.027  0.002 | -0.523  8.452  0.273 |  0.445 10.885
calcium_mg      |  0.472  2.958  0.223 |  0.622 11.971  0.387 |  0.240  3.166
iron_mg         |  0.702  6.536  0.493 |  0.415  5.332  0.172 |  0.302  5.000
magnessium_mg   |  0.881 10.289  0.775 |  0.223  1.543  0.050 | -0.019  0.019
phosphorus_mg   |  0.904 10.851  0.818 | -0.140  0.608  0.020 |  0.064  0.227
                  cos2  
water_g          0.049 |
protein_g        0.001 |
total_fat_g      0.250 |
carbohydrates_g  0.367 |
fiber_g          0.093 |
sugars_g         0.198 |
calcium_mg       0.058 |
iron_mg          0.091 |
magnessium_mg    0.000 |
phosphorus_mg    0.004 |

Supplementary categories
                    Dist    Dim.1   cos2 v.test    Dim.2   cos2 v.test    Dim.3
Frutas          |  1.812 | -1.138  0.395 -3.068 | -1.219  0.452 -5.013 |  0.116
Vegetales       |  0.820 |  0.515  0.395  3.068 |  0.551  0.452  5.013 | -0.053
                  cos2 v.test  
Frutas           0.004  0.638 |
Vegetales        0.004 -0.638 |

Gráfico de variables (1 y 2)

• Todas:

Code

fviz_pca_var(X = acp, axes = c(1, 2))

• 10 variables de mayor contribución:

Code

fviz_pca_var(X = acp, axes = c(1, 2), select.var = list(contrib = 10))

Gráfico de individuos

Code

fviz_pca_ind(X = acp, axes = c(1, 2))

Biplot

Code

fviz_pca_biplot(X = acp, axes = c(1, 2), select.var = list(contrib = 10),
                col.ind = frutas_vegetales$grupo)

Explor (App)

Code

library(explor)
explor(acp)

Análisis de clúster

K-Means idea intuitiva

K = 2

Code

set.seed(1)
kmeans_k2 <- kmeans(x = frutas_vegetales %>% select(where(is.numeric)),
                    centers = 2)

kmeans_k2

K-means clustering with 2 clusters of sizes 103, 19

Cluster means:
   water_g protein_g total_fat_g carbohydrates_g  fiber_g sugars_g calcium_mg
1 84.64553  2.249515   1.0338835       12.126990 2.855340 5.332136   30.40777
2 90.20421  1.755789   0.2726316        6.513684 2.280526 1.732632   71.42105
    iron_mg magnessium_mg phosphorus_mg potassium_mg sodium_g vitamin_a_iu
1 0.7029126      20.29126      46.88350     254.8447 15.49515     632.0971
2 1.4947368      39.00000      40.31579     393.4211 73.10526    9785.4211
  vitamin_c_mg vitamin_b1_mg vitamin_b2_mg viatmin_b3_mg vitamin_b5_mg
1     27.38155    0.06896117    0.05618447     0.6859806     0.2730485
2     19.00526    0.08068421    0.16868421     0.6124737     0.2167895
  vitamin_b6_mg vitamin_e_mg
1     0.1314951    0.5194175
2     0.1428421    1.3447368

Clustering vector:
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 2 2 1 1 1 1 1 1 1 1 1 1 1 2 2 2 1 2 2 1 2 1 2 2 1 2 1 2 2 1 1 1 1 1 1 1
 [75] 2 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 2 2 1 1 1 1 2 1 1 1

Within cluster sum of squares by cluster:
[1] 125979382 232730820
 (between_SS / total_SS =  78.9 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"      

¿Cuál es el mejor K?

Code

set.seed(1)
fviz_nbclust(x = frutas_vegetales %>% select(where(is.numeric)),
             method = "wss",
             FUNcluster = kmeans,
             k.max = 10)

K = 3

Code

set.seed(1)
kmeans_k3 <- kmeans(x = frutas_vegetales %>% select(where(is.numeric)),
                    centers = 3)

kmeans_k3

K-means clustering with 3 clusters of sizes 99, 15, 8

Cluster means:
   water_g protein_g total_fat_g carbohydrates_g  fiber_g sugars_g calcium_mg
1 84.29869  2.260101   1.0616162       12.475556 2.914141 5.487677   28.61616
2 92.39933  1.918667   0.2973333        4.098667 1.835333 1.224000   72.13333
3 87.60125  1.566250   0.2637500        9.535000 2.675000 2.561250   71.75000
    iron_mg magnessium_mg phosphorus_mg potassium_mg sodium_g vitamin_a_iu
1 0.6965657      20.29293      46.78788     253.2424 14.80808     469.3333
2 1.4406667      38.53333      43.73333     409.4000 79.20000    6602.5333
3 1.2787500      30.50000      38.37500     314.0000 41.37500   13190.8750
  vitamin_c_mg vitamin_b1_mg vitamin_b2_mg viatmin_b3_mg vitamin_b5_mg
1     27.37980    0.06871717    0.05392929     0.6760404     0.2749899
2     21.77333    0.07566667    0.19340000     0.6156000     0.1930667
3     18.02500    0.08725000    0.09400000     0.7663750     0.2653750
  vitamin_b6_mg vitamin_e_mg
1     0.1305960    0.5008081
2     0.1309333    1.1826667
3     0.1706250    1.4662500

Clustering vector:
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 2 2 1 2 1 1 2 1 1 1 1 1 1 2 2 2 2 3 2 1 3 1 2 2 1 2 1 2 3 1 1 1 1 1 1 1
 [75] 2 2 1 1 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 3 3 1 1 1 1 3 1 1 1

Within cluster sum of squares by cluster:
[1] 58240314 32104412 61118463
 (between_SS / total_SS =  91.1 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"      

Grupos

Code

df_cluster_k3 <-
  frutas_vegetales %>%
  mutate(cluster = as.factor(kmeans_k3$cluster))

df_cluster_k3 %>% head()

• ¿Cómo quedaron asignadas las frutas y vegetales con el análisis de clúster?

Code

df_cluster_k3  %>%
  count(grupo, cluster)

PCA + K-Means

Code

fviz_pca_biplot(X = acp, axes = c(1, 2), select.var = list(contrib = 10),
                col.ind = df_cluster_k3$cluster)

Interpolación espacial

Datos de ejemplo

Code

# Bibliotecas
library(sp)
library(sf)
library(raster)
library(leaflet)
library(gstat)

# Datos para ejemplo
data("meuse")
data_zinc <- meuse

# Formato sf
zinc_sf <- data_zinc %>% st_as_sf(coords = c(1, 2), crs = "+init=epsg:28992")

# Coordenadas
# x: longitud
# y: latitud
coordinates(data_zinc) <- ~ x + y

Gráfico con ggplot2

Code

zinc_sf %>% 
  ggplot() +
  geom_sf(aes(color = zinc)) +
  theme_bw() +
  scale_color_viridis_c()

Gráfico con leaflet

Code

leaflet(as_Spatial(st_transform(zinc_sf, crs = 4326))) %>%
  addTiles() %>%
  addCircles()

Variogramas

Modelo nulo

Code

variograma_nulo <- variogram(object = log(zinc) ~ 1, data = data_zinc) %>% 
  mutate(Modelo = "Nulo")

variograma_nulo %>% 
  ggplot(aes(x = dist, y = gamma)) +
  geom_point()

Modelo distancias

Code

variograma_dist <- variogram(log(zinc) ~ dist, data = data_zinc) %>%
  mutate(Modelo = "Distancia")

variograma_dist %>% 
  ggplot(aes(x = dist, y = gamma)) +
  geom_point()

Modelo espacial

Code

variograma_esp <- variogram(log(zinc) ~ x + y, data = data_zinc) %>%
  mutate(Modelo = "Espacial")

variograma_esp %>% 
  ggplot(aes(x = dist, y = gamma)) +
  geom_point()

Comparando variogramas

Code

total_variogramas <- bind_rows(variograma_nulo, variograma_dist, variograma_esp)

total_variogramas %>% 
  ggplot(aes(x = dist, y = gamma, color = Modelo)) +
  geom_point()

Regresión Kriging

• Ajuste de variogramas:
  • psill: meseta parcial. Estabilización de la semivarianza.
  • range: separación o distancia entre pares de puntos en la cual ya no hay dependencia espacial.
  • nugget: semivarianza esperada a la separación o distancia de 0m
  • Los valores de psill, range y nugget sirven como valores iniciales para el ajuste del modelo.

Code

# Biblioteca necesaria para ajustar variogramas
library(gstat)

# Ajuste de variorama
ajuste_var <- fit.variogram(object = variograma_esp,
                            model = vgm(
                              model = "Sph",
                              psill = 1,
                              range =  700,
                              nugget =  1
                            ))

# Grid de datos a predecir
data("meuse.grid")

# Conversión a sff
meusegrid_sf <- meuse.grid %>% 
  st_as_sf(coords = c(1, 2), crs = "+init=epsg:28992")

prediccion_spat <- krige(log(zinc) ~ 1, zinc_sf, meusegrid_sf, model = ajuste_var)

[using ordinary kriging]

Code

prediccion_spat %>% 
  ggplot() +
  geom_sf(aes(color = exp(var1.pred))) +
  scale_color_viridis_c(name = "Zinc (predicho)")

Anexo

• Mapa universidad:

Code

library(leaflet)
leaflet() %>% 
  addTiles() %>%
  setView(lng=-75.20082, lat=7.99156, zoom = 16) %>% 
  addProviderTiles(provider = providers$Esri.WorldImagery, group = "Satellite") %>% 
  addLayersControl(baseGroups = c("Street View", "Satellite")) %>% 
  addMarkers(lng=-75.20082, lat=7.99156, popup = "Belmont Park")