predicciones.models

---
title: "Predicciones"
author: "Juan Daniel Bula"
date: "9/12/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
library(dummies)
library(stringr)
library(readxl)
library(sf)
library(dplyr)
library(lubridate)
library(ggplot2)
library(GGally)
library(car)
library(MLmetrics)
library(wordcloud)
library(gplots)
library(R.utils)
library(tm)
library(DescTools)
library(raster)
library(mclust)
library(rgdal)
library(raster)
library(geosphere)
library(NbClust)
library(factoextra)
library(vegan)
library(qpcR)
library(leaflet)
```

```{r}
#Lectura de la Base de datos con la cual se trabajará en este proyecto
base_final <- read.csv("incidentes_final.csv", encoding="UTF-8")
```



# 1 - Entrenamiento de un Modelo Predictivo

Tenemos en cuenta el entrenamiento para los años 2014 hasta el AÑO 2017, para validacion los años 2018 y 2019.

### Modelo Lineal


```{r }
#Modelo lineal
base_final$CLASE <- as.factor(as.character(base_final$CLASE))
datos_vl <- subset(base_final, (AÑO == '2018'))
base_final01 <- subset(base_final, (AÑO != '2018'))
base_final02 <- subset(base_final01, (AÑO != '2019'))
base_final03 <- subset(base_final02, (AÑO != '2020'))
datos_ml1 <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   DISENO) %>% dplyr::count(name ="NRO_ACCID") 


ml1 <- lm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA+DISENO, data = datos_ml1)
promedio <- mean(datos_lm1$NRO_ACCID)
TSS <- sum((datos_lm1$NRO_ACCID - promedio)^2)
RSS <- RSS(lm1)
r2 <- 1-RSS/TSS
RSS2 <- anova(lm1)[4, 2]
r2 <- 1-RSS/TSS
```

En este modelo se va a observar el **MSE** (Error Cuadr?tico Medio) y el **$R^2$** (Coeficiente de Determinaci?n) para determinar la potencia del modelo para predecir.

### Predicci?n y Evaluaci?n para los datos de Entrenamiento

```{r }
ml1_data <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   DISENO) %>% dplyr::count(name = "NRO_ACCID") 
ml1_tr <- ml1_data[,-c(5)]
predicted <- round(predict(ml1, newdata=ml1_tr))
actual <- ml1_data$NRO_ACCID
ml1_mse <- MSE(predicted, actual) # MSE
ml1_mae <- MAE(predicted, actual) # MAE
ml1_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", ml1_mse, ml1_mae, ml1_r2)
```
### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2018

```{r }
ml1_2018 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                  DISENO) %>% dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm1, newdata=lm1_2018))
actual <- lm1_2018$NRO_ACCID
lm1_mse <- MSE(predicted, actual) # MSE
lm1_mae <- MAE(predicted, actual) # MAE
lm1_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", lm1_mse, lm1_mae, lm1_r2)
```
La diferencia del MSE entre los datos de entrenamiento y validaci?n del 2018 es del 34.48%, que al ser mayor que el 15% indica un posible sobreajuste. Adem?s se puede apreciar que el $R^2$ de los datos de validaci?n para el AÑO 2018 predice un 72.86%, sin embargo disminuy? un 16.82% en cuanto al $R^2$ para los datos de entrenamiento. As? que luego de ?sto se decide igualmente ver qu? sucede con el mismo modelo validando con el AÑO 2019.

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2019

```{r }
datos_vl02 <- subset(base_final, (AÑO == '2019'))
ml1_2019 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                  DISENO) %>% dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(ml1, newdata=ml1_2019))
actual <- ml1_2019$NRO_ACCID
ml1_mse <- MSE(predicted, actual) # MSE
ml1_mae <- MAE(predicted, actual) # MAE
ml1_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", ml1_mse, ml1_mae, ml1_r2)
```

La diferencia del MSE entre los datos de entrenamiento y validaci?n del 2019 es del 41.29%, que al ser mayor que el 15% indica un posible sobreajuste. Adem?s se puede evidenciar que el $R^2$ de los datos de validaci?n para el AÑO 2019 predice un 76.53%, sin embargo aunque mejor? con respecto a la validaci?n del AÑO 2018, disminuy? un 13.15% en cuanto al $R^2$ para los datos de entrenamiento. Por tanto, al obtener estos resultados validando con los AÑOs 2018 y 2019 se decide buscar un nuevo modelo cambiando las variables.

## Modelo Lineal con disminuci?n de Variables

Para este nuevo modelo se decide utilizar un modelo lineal ?nicamente con las variables "Festividad" y "D?a Semana". Es decir, se omite en este caso "Dise?o" para observar qu? cambios pueden ocurrir en el modelo.

```{r}
datos_ml2 <- base_final03 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
lm2 <- lm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA, data = datos_ml2)
```

### Predicci?n y Evaluaci?n para los datos de Entrenamiento

```{r}
lm2_tr <- base_final03 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm2, newdata=lm2_tr))
actual <- lm2_tr$NRO_ACCID
lm2_mse <- MSE(predicted, actual) # MSE
lm2_mae <- MAE(predicted, actual) # MAE
lm2_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", lm2_mse, lm2_mae, lm2_r2)
```
### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2018

```{r}
lm2_2018 <- datos_vl %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm2, newdata=lm2_2018))
actual <- lm2_2018$NRO_ACCID
lm2_mse <- MSE(predicted, actual) # MSE
lm2_mae <- MAE(predicted, actual) # MAE
lm2_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", lm2_mse, lm2_mae, lm2_r2)
```
Con este nuevo modelo disminuyendo variables no se obtuvieron resultados positivos, ya que el R2 disminuy? notablemente tanto en el entrenamiento, como en la validaci?n con el AÑO 2018 y la diferencia entre el MSE de datos de entrenamiento y validaci?n para el AÑO 2018 fue de 57.89%, lo cual indica que aument?, indicando as? una alta variabilidad en cuanto a las predicciones y de esa forma, una baja variabilidad explicada por este nuevo modelo.

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2019

```{r }
lm2_2019 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm2, newdata=lm2_2019))
actual <- lm2_2019$NRO_ACCID
lm2_mse <- MSE(predicted, actual) # MSE
lm2_mae <- MAE(predicted, actual) # MAE
lm2_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", lm2_mse, lm2_mae, lm2_r2)
```

Para este modelo la diferencia del MSE entre los datos de entrenamiento y validaci?n del 2019 es del 91.02%, que claramente indica un sobreajuste. Y se puede observar que el $R^2$ de los datos de validaci?n para el AÑO 2019 predice un 10.07%. As? que se concluye que este modelo no sirve para predecir seg?n los resultados obtenidos tanto en el entrenamiento, como para la validaci?n en los AÑOs de 2018 y 2019. As? que se decide utilizar un modelo lineal generalizado.

## Modelo Lineal Generalizado

```{r}
datos_lm3 <- base_final03 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
lm3 <- glm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA, family = "poisson", data = datos_lm3) # Modelo lineal generalizado, con familia poisson
```

### Predicci?n y Evaluaci?n para los datos de Entrenamiento

```{r}
lm3_tr <- base_final03 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
lm3_tr_1 <- lm3_tr[,-4]
predicted <- round(predict(lm3, newdata=lm3_tr_1, type="response"))
actual <- lm3_tr$NRO_ACCID
lm3_mse <- MSE(predicted, actual) # MSE
lm3_mae <- MAE(predicted, actual) # MAE
lm3_r2 <- R2_Score(predicted, actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm3_mse, lm3_mae, lm3_r2)
```

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2018

```{r}
lm3_2018 <- datos_vl %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm3, newdata=lm3_2018, type="response")) 
actual <- lm3_2018$NRO_ACCID
lm3_mse <- MSE(predicted, actual) # MSE
lm3_mae <- MAE(predicted, actual) # MAE
lm3_r2 <- R2_Score(predicted, actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm3_mse, lm3_mae, lm3_r2)
```

En este modelo lineal generalizado con la familia de distribuci?n Poisson se obtuvo un $R^2$ en la etapa de entrenamiento de 7.21%, y para la etapa de validaci?n para el AÑO 2018 el $R^2$ fue de 2.88%,lo cual indica que dicho modelo no sirve para predecir seg?n los resultados obtenidos. Adem?s la diferencia entre el MSE de entrenamiento y validaci?n para el AÑO 2018 fue de 57.89%, que indica un sobreajuste. Sin embargo se procede a validar igualmente con el AÑO 2019.

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2019

```{r}
lm3_2019 <- datos_vl02 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm3, newdata=lm3_2019, type="response")) 
actual <- lm3_2019$NRO_ACCID
lm3_mse <- MSE(predicted, actual) # MSE
lm3_mae <- MAE(predicted, actual) # MAE
lm3_r2 <- R2_Score(predicted, actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm3_mse, lm3_mae, lm3_r2)
```

Para este modelo la diferencia del MSE entre los datos de entrenamiento y validaci?n del 2019 es del 90.91%, que claramente indica un sobreajuste. Y se puede observar que el $R^2$ de los datos de validaci?n para el AÑO 2019 predice un 10.1%. As? que se evidencia que este modelo no sirve para predecir seg?n los resultados obtenidos tanto en el entrenamiento, como para la validaci?n en los AÑOs de 2018 y 2019. As? que se decide utilizar el modelo lineal generalizado adicionando otra variable.

## Modelo lineal generalizado con Adici?n de la variable Clase

```{r}
datos_lm4 <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   CLASE) %>% dplyr::count(name = "NRO_ACCID")
lm4 <- glm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA+CLASE, family = "poisson", 
           data = datos_lm4)
```

### Predicci?n y Evaluaci?n para los datos de Entrenamiento

```{r}
datos_lm4_p <- datos_lm4[,-5]
y_train <- round(predict(lm4, newdata= datos_lm4_p, type="response"))
y_actual <- datos_lm4$NRO_ACCID
lm4_tmse <- MSE(y_train, y_actual)
lm4_tmae <-  MAE(y_train, y_actual)
lm4_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm4_tmse, lm4_tmae, lm4_r2)
```

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2018

```{r}
datos_lm4_v1 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm4_v2 <- datos_lm4_v1[,-5]
y_train <- round(predict(lm4, newdata= datos_lm4_v2, type="response"))
y_actual <- datos_lm4_v1$NRO_ACCID
lm4_tmse <- MSE(y_train, y_actual)
lm4_tmae <-  MAE(y_train, y_actual)
lm4_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm4_tmse, lm4_tmae, lm4_r2)
```

En este modelo lineal generalizado con la familia de distribuci?n Poisson con adici?n de la variable Clase se obtuvo un $R^2$ en la etapa de entrenamiento de 88.6%,y para la etapa de validaci?n para el AÑO 2018 el $R^2$ fue de 89.51%,lo cual indica que dicho modelo sirve para predecir seg?n los resultados obtenidos. Adem?s la diferencia entre el MSE de entrenamiento y validaci?n para el AÑO 2018 fue de 10.03%, lo que indica que al ser menor del 15% no hay problemas de sobreentrenamiento. Luego se procede a validar con el AÑO 2019. 

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2019

```{r}
datos_lm4_v1 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm4_v2 <- datos_lm4_v1[,-5]
y_train <- round(predict(lm4, newdata= datos_lm4_v2, type="response"))
y_actual <- datos_lm4_v1$NRO_ACCID
lm4_tmse <- MSE(y_train, y_actual)
lm4_tmae <-  MAE(y_train, y_actual)
lm4_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm4_tmse, lm4_tmae, lm4_r2)
```

En este modelo la diferencia del MSE entre los datos de entrenamiento y validaci?n del 2019 es del 0.94%, que fue menor al valor obtenido con la validaci?n del 2018 (10.03%), que indica claramente que no hay problemas de sobreentrenamiento. Adem?s se puede observar que el $R^2$ de los datos de validaci?n para el año 2019 predice un 88.63%, evidenciando as? que este modelo lineal generalizado con la adici?n de la variable Clase es un buen candidato para predecir la accidentalidad en Medell?n. Sin embargo, se adicionar? otra variable para ver si se obtiene un mejor modelo.

## Modelo lineal generalizado con Adici?n de la variable Dise?o

```{r}
datos_lm5 <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   CLASE, DISENO) %>% dplyr::count(name = "NRO_ACCID")
lm5 <- glm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA+CLASE+DISENO, family = "poisson", 
           data = datos_lm5)
```

### Predicci?n y Evaluaci?n para los datos de Entrenamiento

```{r}
datos_lm5_p <- datos_lm5[,-6]
y_train <- round(predict(lm5, newdata= datos_lm5_p, type="response"))
y_actual <- datos_lm5$NRO_ACCID
lm5_tmse01 <- MSE(y_train, y_actual)
lm5_tmae <-  MAE(y_train, y_actual)
lm5_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm5_tmse01, lm5_tmae, lm5_r2)
```
### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2018

```{r}
datos_lm5_v1 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                      CLASE, DISENO) %>% dplyr::count(name = "NRO_ACCID")
datos_lm5_v2 <- datos_lm5_v1[,-6]
y_train <- round(predict(lm5, newdata= datos_lm5_v2, type="response"))
y_actual <- datos_lm5_v1$NRO_ACCID
lm5_tmse02 <- MSE(y_train, y_actual)
lm5_tmae <-  MAE(y_train, y_actual)
lm5_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm5_tmse02, lm5_tmae, lm5_r2)
```
En este modelo lineal generalizado con la familia de distribuci?n Poisson con adici?n de la variable Dise?o se obtuvo un $R^2$ en la etapa de entrenamiento de 86.56%, y para la etapa de validaci?n para el AÑO 2018 el $R^2$ fue de 81.93%,lo cual indica que dicho modelo sirve para predecir seg?n los resultados obtenidos. La diferencia entre el MSE de entrenamiento y validaci?n para el AÑO 2018 fue de 4.92%, lo que indica que al ser menor del 15% no hay problemas de sobreentrenamiento. As?, se procede a validar con el AÑO 2019. 

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2019

```{r}
datos_lm5_v1 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                      CLASE, DISENO) %>% dplyr::count(name = "NRO_ACCID")
datos_lm5_v2 <- datos_lm5_v1[,-6]
y_train <- round(predict(lm5, newdata= datos_lm5_v2, type="response"))
y_actual <- datos_lm5_v1$NRO_ACCID
lm5_tmse02 <- MSE(y_train, y_actual)
lm5_tmae <-  MAE(y_train, y_actual)
lm5_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm5_tmse02, lm5_tmae, lm5_r2)
```

En este modelo la diferencia del MSE entre los datos de entrenamiento y validaci?n del 2019 es del 3.4%, que indica que no hay problemas de sobreentrenamiento. Tambi?n se puede observar que el $R^2$ de los datos de validaci?n para el año 2019 predice un 82.91%, que evidencia que este modelo lineal generalizado con la adici?n de la variable Dise?o es buen candidato para predecir la accidentalidad en Medell?n.

Este modelo presenta el mejor MSE, pero no es un modelo viable ya que no es posible obtener la variable DISENO para realizar predicciones, según las pautas dadas.

Finalmente, se decide trabajar con el modelo lineal generalizado que posee la variable CLASE, ya que como se analiz? anteriormente, se observa que tanto para los datos de entrenamiento, como para los de validaci?n (2018, 2019) se pueden obtener buenas predicciones.

Luego, con el modelo elegido se procedi? a realizar la etapa de entrenamiento y validaci?n de forma semanal y mensual ya que dicho modelo para la predicci?n diaria tuvo buenos resultados, tanto para el MSE, como para el $R^2$, para as? observar si tambi?n con dicho modelo se pueden obtener buenas predicciones no solamente de manera diaria.



## Mensual

## Modelo lineal generalizado con Adici?n de la variable Clase

```{r}
datos_lm7 <- base_final03 %>% group_by(FECHA, FESTIVIDAD, MES, 
                                   CLASE) %>% dplyr::count(name = "NRO_ACCID")
lm7 <- glm(NRO_ACCID ~ FESTIVIDAD+MES+CLASE, family = "poisson", 
           data = datos_lm7)
```

### Predicci?n y Evaluaci?n para los datos de Entrenamiento

```{r}
datos_lm7_p <- datos_lm7[,-5]
y_train <- round(predict(lm7, newdata= datos_lm7_p, type="response"))
y_actual <- datos_lm7$NRO_ACCID
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm7_tmse, lm7_tmae, lm7_r2)
```

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2018

```{r}
datos_lm7_v1 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, MES, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm7_v2 <- datos_lm7_v1[,-5]
y_train <- round(predict(lm7, newdata= datos_lm7_v2, type="response"))
y_actual <- datos_lm7_v1$NRO_ACCID
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm7_tmse, lm7_tmae, lm7_r2)
```

En este modelo lineal generalizado con la familia de distribuci?n Poisson con adici?n de la variable Clase se obtuvo un $R^2$ en la etapa de entrenamiento de 88.84%,y para la etapa de validaci?n para el AÑO 2018 el $R^2$ fue de 88.81%,lo cual indica que dicho modelo sirve para predecir seg?n los resultados obtenidos. Adem?s la diferencia entre el MSE de entrenamiento y validaci?n para el AÑO 2018 fue de 2.19%, lo que indica que al ser menor del 15% no hay problemas de sobreentrenamiento. Luego se procede a validar con el AÑO 2019. 

### Predicci?n y Evaluaci?n para los datos de Validaci?n en el AÑO 2019

```{r}
datos_lm7_v1 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, MES, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm7_v2 <- datos_lm7_v1[,-5]
y_train <- round(predict(lm7, newdata= datos_lm7_v2, type="response"))
y_actual <- datos_lm7_v1$NRO_ACCID
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm7_tmse, lm7_tmae, lm7_r2)
```

En este modelo la diferencia del MSE entre los datos de entrenamiento y validaci?n del 2019 es del  1.3%, que fue menor al valor obtenido con la validaci?n del 2018 (2.19%), que indica claramente que no hay problemas de sobreentrenamiento. Adem?s se puede observar que el $R^2$ de los datos de validaci?n para el AÑO 2019 predice un 88.83%, evidenciando as? que este modelo lineal generalizado con la adici?n de la variable Clase es un buen candidato para predecir la accidentalidad en Medell?n.

## Predicciones

Ahora, se muestra gr?ficamente la potencia de la predicci?n para el AÑO 2020, seg?n los datos reales de 2020 que se poseen desde la fecha 2020-01-01, hasta la fecha 2020-08-31.

Para ello se cre? una base que contuviese los datos del AÑO 2020 real y el del AÑO 2020 con las predicciones; como en este AÑO se tom? desde 2020-01-01 hasta 2020-12-31 porque es el predictivo, entonces para esta comparaci?n se decide que va hasta 2020-08-31, para que ?sta posea la misma longitud de los datos que se tienen reales.

Igualmente, se dise?? una funci?n para contener la relaci?n entre la desviaci?n est?ndar y la ra?z cuadrada de la longitud. Para luego con ayuda de esto y junto a las especificaciones del gr?fico, observar la cantidad de accidentes reales y predichos mediante boxplots, desde el mes 01 hasta el mes 08.

```{r }
#Se cre? una funci?n que contuviera la relaci?n entre desviaci?n est?ndar y la ra?z cuadrada de la longitud
sem <- function(x, na.rm = FALSE) {
  out <- sd(x, na.rm = na.rm)/sqrt(length(x))
  return(out)}
#Se a?ade la base de datos especificada
datos_real_pred <- read_excel("Base_2020_real_predict.xlsx")
#Se hacen los agrupamientos necesarios y se omite la varible Incendio para que en ambos AÑOs est?n las mismas 'Clases'
datos_real_pred_001 <- datos_real_pred[datos_real_pred$CLASE!="Incendio", ]
datos_real_pred_001$AÑO<-as.integer(datos_real_pred_001$AÑO)
datos_real_pred_01 <- datos_real_pred_001 %>% group_by(AÑO, MES, CLASE, NRO_ACCID) %>% 
  count() %>% mutate(NRO_ACCID=as.integer(NRO_ACCID))
#Gr?fico Comparativo de AÑOs reales y predictivos del 2020
datos_real_pred_01 %>% 
  group_by(MES, AÑO) %>% 
  summarize(Número_de_Accidentes = mean(NRO_ACCID),
            se = sem(NRO_ACCID)) %>% 
  ggplot(aes(x = MES, 
             y = Número_de_Accidentes, 
             group = AÑO,  
             color = AÑO)) +
  geom_point() +
  geom_errorbar(aes(ymin = Número_de_Accidentes - se, 
                    ymax = Número_de_Accidentes + se)) +
  geom_line()
```

Nota: Se aclara que en esta gr?fica los AÑOs que se comparan son dos (2020 con datos reales y 2020 con los datos predictivos). En ese orden, se hace referencia a 'AÑO 1' como el AÑO 2020 con los datos reales que se poseen, y de la misma forma, se hace referencia a 'AÑO 2', como el AÑO 2020 con las predicciones. Por tanto, se aconseja hacer caso omiso a los valores de AÑOs 1.25, 1.5 y 1.75.

En la gr?fica se observa que el AÑO 1 (representado con la l?nea y boxplots color negro) tiene menos variabilidad ya que son los valores reales que se poseen de la cantidad de accidentes durante el tiempo que se especific? anteriormente, e igualmente se puede ver su promedio de Números de accidentes a trav?s de los 8 meses, representado por el punto de los diagramas de Cajas y Bigotes.

Luego, para el AÑO 2 (representado con la l?nea y boxplots color azul) se nota m?s variabilidad, lo cual tiene sentido ya que hace referencia a los valores predictivos seg?n el modelo empleado, e igualmente se puede ver su promedio de Números de accidentes a trav?s de los 8 meses, representado por el punto que se observa para cada Boxplot.

Se nota claramente que para los dos primeros meses del AÑO 2020 el promedio de Número de accidentes del 'AÑO 2', que contiene los valores predictivos, son muy similares a los que se observan del AÑO 1, con los valores reales. Luego, para los dem?s meses del AÑO 2020 se puede apreciar que el AÑO con los datos reales (AÑO 1) su promedio de accidentes baj?, observ?ndose as? que el AÑO con las predicciones qued? por encima de ?ste. Sin embargo, ?sto no significa que el AÑO 2 tenga malas predicciones, ni tampoco quiere decir que el modelo no ten?a potencia predictiva, ya que anteriormente se hab?a notado su capacidad predictiva en el estudio de la Accidentalidad de Medell?n seg?n su MSE y $R^2$, lo que realmente sucede es que se debe poner en contexto lo que suced?a en el pa?s en dichos meses.

Para el AÑO 2020 en Colombia, El Ministerio de Salud y Protecci?n Social confirma el primer caso de COVID-19 (enfermedad infecciosa, causada por el coronavirus) en el mes 3, a causa de ello el 8 de marzo de 2020, el presidente Iv?n Duque M?rquez dio a conocer las acciones de contenci?n en el pa?s. Por consiguiente, los ciudadanos adoptando dichas medidas y con el prop?sito de que no aumentara la tasa de contagios por coronavirus, empezaron a tener conciencia social y a cuidarse y as? mismo a cuidar a los dem?s. De hecho, desde febrero que ya se estaba especulando que el COVID-19 iba a llegar al pa?s, muchas personas ya se estaban cuidando.

As?, entre las medidas que se dieron, se encuentra la de cierre de v?as en marzo, que es por eso que gr?ficamente se not? un declive para los datos del promedio de accidentes a partir del mes 3 para los datos reales, 'AÑO 1'. Adem?s de que luego se implement? cuarentena obligatoria por el causal de pandemia, que fue conocido como 'El Aislamiento Preventivo Obligatorio en Colombia o confinamiento de Colombia de 2020' que se inici? el 25 de marzo de 2020, y que en un principio se hab?a decretado por 19 d?as, pero luego fue variando y extendi?ndose dicha cuarentena, seg?n la tasa de contagios.

Para el AÑO 2021 no se poseen valores reales, puesto que no tenemos datos con respecto a dicha informaci?n.

A continuaci?n, se presentan las predicciones diarias, semanales y mensuales que se obtuvieron para los AÑOs 2020 y 2021.

## Predicci?n Diaria del 2020

Se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones diarias obtenidas para el AÑO 2020:

```{r}
Base_prediccion <- read.csv("prediccion.csv", sep = ",", encoding = "UTF-8")
Base_prediccion <- Base_prediccion[,-1]
Base_prediccion_2020 <- subset(Base_prediccion, (AÑO != '2021'))
Base_prediccion_2020$FECHA <- as.Date(Base_prediccion_2020$FECHA)
Base_prediccion_2020$CLASE <- as.factor(Base_prediccion_2020$CLASE)
Base_prediccion_2020$DIA_SEMANA <- as.factor(Base_prediccion_2020$DIA_SEMANA)
Base_prediccion_2020$AÑO <- as.integer(Base_prediccion_2020$AÑO)
Base_prediccion_2020$FESTIVIDAD <- as.factor(Base_prediccion_2020$FESTIVIDAD)
prediccion_2020 <- predict(object = lm4, newdata = Base_prediccion_2020,
                          type = "response")
prediccion_diaria2020 <- Base_prediccion_2020 %>% 
  mutate(NRO_ACCID = round(prediccion_2020,0))
diario_20_02 <- prediccion_diaria2020 %>%
  group_by(FECHA, DIA_SEMANA, CLASE, FESTIVIDAD) %>%
  summarise(NRO_TOTAL_ACCID=NRO_ACCID)
head(diario_20_02, 10)
```

## Predicci?n Diaria del 2021

A continuaci?n, se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones diarias obtenidas para el AÑO 2021:

```{r}
Base_prediccion_2021 <- subset(Base_prediccion, (AÑO != '2020'))
Base_prediccion_2021$FECHA <- as.Date(Base_prediccion_2021$FECHA)
Base_prediccion_2021$CLASE <- as.factor(Base_prediccion_2021$CLASE)
Base_prediccion_2021$DIA_SEMANA <- as.factor(Base_prediccion_2021$DIA_SEMANA)
Base_prediccion_2021$AÑO <- as.integer(Base_prediccion_2021$AÑO)
Base_prediccion_2021$FESTIVIDAD <- as.factor(Base_prediccion_2021$FESTIVIDAD)
prediccion_2021 <- predict(object = lm4, newdata = Base_prediccion_2021,
                          type = "response")
prediccion_diaria2021 <- Base_prediccion_2021 %>% 
  mutate(NRO_ACCID = round(prediccion_2021,0))
diario_21_02 <- prediccion_diaria2021 %>%
  group_by(FECHA, DIA_SEMANA, CLASE, FESTIVIDAD) %>%
  summarise(NRO_TOTAL_ACCID=NRO_ACCID)
head(diario_21_02, 10)
```




## Predicci?n Mensual del 2020

De manera resumida, se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones mensuales que se obtuvieron para el AÑO 2020:

```{r}

Base_prediccion <- read.csv("prediccion.csv", sep = ",", encoding = "UTF-8")
Base_prediccion_2020 <- subset(Base_prediccion, (AÑO != '2021'))
Base_prediccion_2020$FECHA <- as.Date(Base_prediccion_2020$FECHA)
Base_prediccion_2020$CLASE <- as.factor(Base_prediccion_2020$CLASE)
Base_prediccion_2020$DIA_SEMANA <- as.factor(Base_prediccion_2020$DIA_SEMANA)
Base_prediccion_2020$AÑO <- as.integer(Base_prediccion_2020$AÑO)
Base_prediccion_2020$FESTIVIDAD <- as.factor(Base_prediccion_2020$FESTIVIDAD)
prediccion_2020 <- predict(object = lm4, newdata = Base_prediccion_2020,
                          type = "response")
prediccion_mensual2020 <- Base_prediccion_2020 %>% 
  mutate(NRO_ACCID = round(prediccion_2020,0))
#Se borraron columnas no necesarias
prediccion_mensual2020 <-  prediccion_mensual2020[,c(-1,-2,-3,-5,-7)]
#Agrupamiento por mes 2020
#Se agrup? por mes y si fue en d?a festivo o no
mensual_20 <- prediccion_mensual2020 %>% group_by(CLASE, MES, NRO_ACCID, FESTIVIDAD) %>% dplyr::summarize(total = n())
mensual_20 <- mutate(mensual_20, NRO_ACCID_TOTAL=NRO_ACCID*total)
mensual_20_02 <- mensual_20 %>%
  group_by(MES, CLASE, FESTIVIDAD) %>%
  summarise(NRO_TOTAL_ACCID=sum(NRO_ACCID_TOTAL))
head(mensual_20_02, 10)
```

## Predicci?n Mensual del 2021

Igualmente, se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones mensuales que se obtuvieron para el AÑO 2021:

```{r}
Base_prediccion_2021 <- subset(Base_prediccion, (AÑO != '2020'))
Base_prediccion_2021$FECHA <- as.Date(Base_prediccion_2021$FECHA)
Base_prediccion_2021$CLASE <- as.factor(Base_prediccion_2021$CLASE)
Base_prediccion_2021$DIA_SEMANA <- as.factor(Base_prediccion_2021$DIA_SEMANA)
Base_prediccion_2021$AÑO <- as.integer(Base_prediccion_2021$AÑO)
Base_prediccion_2021$FESTIVIDAD <- as.factor(Base_prediccion_2021$FESTIVIDAD)
prediccion_2021 <- predict(object = lm4, newdata = Base_prediccion_2021,
                          type = "response")
prediccion_mensual2021 <- Base_prediccion_2021 %>% 
  mutate(NRO_ACCID = round(prediccion_2021,0))
#Se borranron columnas no necesarias
prediccion_mensual2021 <-  prediccion_mensual2021[,c(-1,-2,-3,-5,-7)]
#Agrupamiento por mes 2021
#Se agrupó por mes y si fue en día festivo o no
mensual_21 <- prediccion_mensual2021 %>% group_by(CLASE, MES, NRO_ACCID, FESTIVIDAD) %>% dplyr::summarize(total = n())
mensual_21 <- mutate(mensual_21, NRO_ACCID_TOTAL=NRO_ACCID*total)
mensual_21_02 <- mensual_21 %>%
  group_by(MES, CLASE, FESTIVIDAD) %>%
  summarise(NRO_TOTAL_ACCID=sum(NRO_ACCID_TOTAL))
head(mensual_21_02, 10)
```

### 2 - Agrupamiento de los barrios de Medell?n de acuerdo a su accidentalidad

Antes de realizar la agrupaci?n de barrios seg?n la accidentalidad se decide mostrar primero un mapa de calor de accidentalidad en la ciudad de Medell?n entre los AÑOs de 2014 y 2019, ya que este mapa a su vez representar?a la historia de la accidentalidad. Lo anterior, teniendo presente que para este proyecto se pretende predecir para los AÑOs de 2020 y 2021.

### Mapa de Accidentes en la ciudad de Medell?n

Para la elaboraci?n del mapa calor seg?n la accidentalidad, se descarg? un archivo .shp del L?mite Barrio Vereda Catastral

```{r}
#lectura de .csv y .shp
catastral <- read.csv("Limite_Barrio_Vereda_Catastral.csv", encoding="UTF-8")
catastro <- read_sf("Limite_Barrio_Vereda_Catastral.shp")
barrio_vereda <- read.csv("Barrio_Vereda_2014.csv", encoding="UTF-8")
```

```{r}
#Mapa para todos los barrios, usando 'innerjoin' con el .shp de Limite_Barrio_Vereda_Catastral
Unido <- inner_join(catastral, base_final, by = c("COMUNA" = "NUMCOMUNA"))
nueva_base <- Unido %>% filter(AÑO >= 2014 & AÑO <= 2019) %>% 
  group_by(CODIGO) %>%
  dplyr::summarise(accidentes = n()) %>%
  dplyr::ungroup()
#Se realiz? la conversi?n de CODIGO a formato num?rico
catastro$CODIGO <- as.numeric(as.character(catastro$CODIGO))
#Se utiliz? 'inner join' para unir dos bases y para luego generar mapa
mapa <- inner_join(catastro, nueva_base, by = c("CODIGO" = "CODIGO"))
mypal <- colorNumeric(palette = c("#000000","#280100","#3D0201","#630201","#890100","#B00100","#DD0100","#F50201",
                                   "#FF5F5E","#FF7A79","#FF9796","#FEB1B0","#FDC9C8", "#FFE5E4"), domain = mapa$accidentes, reverse = T)
# Creaci?n del mapa
leaflet() %>% addPolygons(data = mapa, color = "#0A0A0A", opacity = 0.6, weight = 1, fillColor = ~mypal(mapa$accidentes),
                          fillOpacity = 0.6, label = ~NOMBRE_BAR,
                          highlightOptions = highlightOptions(color = "black", weight = 3, bringToFront = T, opacity = 1),
                          popup = paste("Barrio: ", mapa$NOMBRE_BAR, "<br>", "Accidentes: ", mapa$accidentes, "<br>")) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend(position = "bottomright", pal = mypal, values = mapa$accidentes, title = "Accidentes", opacity = 0.6)
```

```{r}
# Cantidad de Accidentes por Comuna
medellin_comuna <- base_final %>% filter(AÑO >= 2014 & AÑO <= 2019) %>% 
  group_by(COMUNA) %>% 
  dplyr::summarize(accidentes = n())
ggplot(data = medellin_comuna, aes(x = reorder(COMUNA,+accidentes), y = accidentes)) +
  geom_bar(stat = "identity", position = "dodge", fill = "blue", color = "black", alpha = 0.6) +
  geom_text(aes(y = accidentes, label = accidentes),
                position = position_dodge(width = 0.7), size = 3.5, vjust = 0.5, hjust = -0.1, col = "black") +
  xlab("Comuna") + 
  ylab("Total de Accidentes") +
  ggtitle("Total de Accidentes por Comuna entre los AÑOs 2014 y 2019") +
  ylim(c(0,50000)) +
  theme_minimal() +
  coord_flip()
```

Al analizar la accidentalidad por Barrio y Comuna, tanto en el mapa de accidentalidad como en el gr?fico, se observ? que La Candelaria es el lugar donde m?s se presentan accidentes entre los AÑOs 2014 y 2019, seguida por Laureles y Castilla. Tambi?n, se evidencia que Palmitas es el lugar donde ocurre menos accidentalidad en Medell?n durante los AÑOs 2014-2019. Adem?s en el mapa de calor se puede notar que la zona centro y las v?as principales son las m?s afectadas por los accidentes.

Luego se procede a crear una funci?n para calular distancias para datos geoespaciales y adem?s se crea un dendograma para as? emprender la b?squeda del $K$ ?ptimo para el agrupamiento (tomando los datos de 2014 hasta 2017, ya que ?stos fueron los AÑOs que se utilizaron para el modelo de predicci?n).

```{r}
#Haciendo uso de la librer?a 'geosphere', se cre? una funci?n para calcular las distancias para datos geoespaciales
geo.dist = function(df) {
  require(geosphere)
  d <- function(i,z){         # z[1:2] contain long, lat
    dist <- rep(0,nrow(z))
    dist[i:nrow(z)] <- distHaversine(z[i:nrow(z),1:2],z[i,1:2])
    return(dist)
  }
  dm <- do.call(cbind,lapply(1:nrow(df),d,df))
  return(as.dist(dm))
}
```

```{r }
#Se realiz? la conversi?n de la latitud y longitud al formato num?rico
base_final03$LATITUD <- as.numeric(as.character(base_final03$LATITUD))
base_final03$LONGITUD<- as.numeric(as.character(base_final03$LONGITUD))
```

```{r }
#Se cre? un nuevo dataset para el agrupamiento, seg?n longitud, latitud y barrio almacenado en 'df'
df <- data.frame(long = base_final03$LONGITUD, lat = base_final03$LATITUD, barrios = base_final03$BARRIO)
```

```{r}
#Se cre? con la funci?n 'geo.dist', una matriz de distancias
df1 <- df[1:1000, ]
d <- geo.dist(df1)
hc <- hclust(d)
plot(hc, main = "Dendograma", col = "#00AFBB")
df1$clust <- cutree(hc, k = 6)
head(df1,10)
```

### Mapa de agrupamiento seg?n latitud y longitud.

Para la realizaci?n del mapa de agrupamiento seg?n latitud y longitud, se descarg? un archivo .shp del L?mite Catastral de Comunas y Corregimientos.

```{r}
s <- shapefile("Limite_Catastral_de__Comunas_y_Corregimientos.shp")
map.df1 <- (s)
ggplot(map.df1)+
  geom_path(aes(x=long, y=lat, group=group))+
  geom_point(data=df1, aes(x=long, y=lat, color=factor(clust)), size=4)+
  scale_color_discrete("Cluster")+
  coord_fixed()
```

El anterior mapa muestra una posible agrupaci?n, seg?n las medidas geoespaciales de la latitud y la longitud de los accidentes. Sin embargo, este agrupamiento se utiliza como referencia porque para su creaci?n no se utiliz? ning?n m?todo para la elecci?n del $K$ ?ptimo.

As? que se procede a realizar una clusterizaci?n del Número de accidentes por Gravedad y Barrio haciendo uso del algoritmo "k means" para la b?squeda del $K$ ?ptimo.

### Clusterizaci?n con Número de accidentes por Gravedad y Barrio.

Con los datos preprocesados y el subconjunto de datos que se seleccionaron, se les realiz? un escalamiento y centrado de la base de datos.

```{r}
#Numero de accidentes por Barrio
datos_cluster <- base_final03 %>% group_by(BARRIO) %>% dplyr::count(name = "TOTAL_ACCIDENTES")
#Número de accidentes por barrio, seg?n gravedad almacenado en 'df'
df <- as.matrix(table(base_final03$BARRIO, base_final03$GRAVEDAD))
df <- data.frame(Con_heridos = df[,1], Con_muertos = df[,2], Solo_danos = df[,3])
#Escalamiento y centrado de la base de datos.
scaled_data = as.matrix(scale(df))
head(scaled_data, 10)
kmm = kmeans(scaled_data, 5, nstart = 50, iter.max = 15 )
```

Para este caso, con un k=5, se evidencia que el valor de Suma de Cuadrados entre grupos (SS between) sobre la Suma de Cuadrados Totales fue de aproximadamente 85.1% (0.851), que indica un buen ajuste porque es cercano a 1. Sin embargo, es mejor graficar el WSS contra el Número de cl?ster, ya que este Número se debe especificar de antemano.

Luego se procede a hallar el k ?ptimo.

### El Método del Codo

```{r}
#Se fij? una semilla y se realiz? el calculo y se grafico el WSS(total within - cluster sum of square) para k = 2 hasta k = 10
set.seed(2021022)
k.max <- 10
datos <- scaled_data
wss <- sapply(2:k.max, 
              function(k){kmeans(datos, k, nstart = 50, iter.max = 15 )$tot.withinss})
plot(2:k.max, wss, 
     type = "b", pch = 19, frame = FALSE,
     xlab = "Número de Clusters (k)",
     ylab = "WSS Total", 
     main = "M?todo del Codo", col="forestgreen")
```

```{r}
#Con k=5, se obtiene between_SS / total_SS =  85.1 %) almacenado en 'km'
km <- kmeans(datos, 5)
```

Seg?n la gr?fica del M?todo del Codo posiblemente el k=4 o k=5 ser?an buenos candidatos para el k ?ptimo, ya que presentan un cambio más suave en las pendientes en comparaci?n con k=2 o k=3. Igualmente al observar el between SS / total SS para k=5, mencionado anteriomente, se evidencia un 85.1 %, lo cual indica un buen ajuste. Adem?s como se grafic? el WSS contra el Número de cl?steres, se refleja que es un buen candidato.

Despu?s, se busca el $K$ ?ptimo haciendo uso del paquete NbClust, de la siguiente forma:

```{r }
nb <- NbClust(scaled_data, diss=NULL, distance = "euclidean", 
              min.nc=4, max.nc=8, method = "kmeans", 
              index = "all", alphaBeale = 0.1)
```
El cual sugiere que el $K$ ?ptimo es 4.


```{r }
hist(nb$Best.nc[1,], breaks = max(na.omit(nb$Best.nc[1,])), main = "Histograma del K ?ptimo ", xlab = "K", ylab = "Frecuencia", col="darkorchid4")
```

Seg?n el histograma, también indica que el $K$ ?ptimo es el 4.

### Resumen de M?todos

En este resumen se encuentra el método de la Silueta, el M?todo del Codo y Gap Statistic. Donde se obtuvieron los siguientes resultados:

### M?todo de la silueta.
```{r echo=FALSE, message=FALSE, warning=FALSE}
fviz_nbclust(scaled_data, kmeans, method = c("silhouette"))
```

### M?todo del Codo.
```{r echo=FALSE, message=FALSE, warning=FALSE}
fviz_nbclust(scaled_data, kmeans, method = "wss") +
    geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "M?todo del Codo")
```

###Gap Statistic
```{r }
set.seed(123)
fviz_nbclust(scaled_data, kmeans, nstart = 25,  method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap statistic method")
```

### Generaci?n de clusterizaci?n, seg?n el k ?ptimo.

Seg?n los diferentes m?todos $k = 4$, parec?a ser muy ?ptimo para la generaci?n de la clusterizaci?n.

Luego, se muestran las primeras 10 observaciones de los barrios ordenados alfab?ticamente, donde se observa el tipo de gravedad y el grupo al cual pertenecen.

```{r }
kmm = kmeans(scaled_data, 4, nstart = 50, iter.max = 15 )
df_clust <- data.frame(Con_heridos = df[,1], Con_muertos = df[,2], Solo_danos = df[,3], kmm$cluster)
head(df_clust, 10)
```

Con lo realizado anteriormente poseemos el agrupamiento seg?n el k ?ptimo igual a 4, es decir, ahora los barrios pertenecen a un tipo de cluster (enumerado del 1 al 4). Para el nombramiento de cada grupo se tom? como referencia el mapa de calor, quedando de la siguiente forma:

- Grupo 1: Accidentalidad Moderada
- Grupo 2: Accidentalidad Baja
- Grupo 3: Accidentalidad Alta
- Grupo 4: Accidentalidad Media-Alta

La informaci?n que se obtuvo de cada grupo es la siguiente:

## Grupo 1 - Accidentalidad Moderada

```{r}
#Accidentalidad Moderada
dfclust_clust1 <- df_clust[df_clust$kmm.cluster == 1, ]
dfclust_clust1$total <- rowSums(dfclust_clust1[,1:3])
sum(dfclust_clust1$Con_heridos)
sum(dfclust_clust1$Con_muertos)
sum(dfclust_clust1$Solo_danos)
sum(dfclust_clust1$total)
```

Para este grupo se observa que la cantidad total de accidentes son de 27895, la cantidad de accidentes con heridos es de 13145, la cantidad de accidentes con muertos 111 y la cantidad de accidentes con solo dAÑOs es de 14639.

## Grupo 2 - Accidentalidad Baja

```{r}
#Accidentalidad Baja
dfclust_clust2 <- df_clust[df_clust$kmm.cluster == 2, ]
dfclust_clust2$total <- rowSums(dfclust_clust2[,1:3])
sum(dfclust_clust2$Con_heridos)
sum(dfclust_clust2$Con_muertos)
sum(dfclust_clust2$Solo_danos)
sum(dfclust_clust2$total)
```

Para este grupo se observa que la cantidad total de accidentes son de 28622, la cantidad de accidentes con heridos es de 17817, la cantidad de accidentes con muertos 106 y la cantidad de accidentes con solo dAÑOs es de 10699.

## Grupo 3 - Accidentalidad Alta

```{r}
#Accidentalidad Alta
dfclust_clust3 <- df_clust[df_clust$kmm.cluster == 3, ]
dfclust_clust3$total <- rowSums(dfclust_clust3[,1:3])
sum(dfclust_clust3$Con_heridos)
sum(dfclust_clust3$Con_muertos)
sum(dfclust_clust3$Solo_danos)
sum(dfclust_clust3$total)
```

Para este grupo se observa que la cantidad total de accidentes son de 36247, la cantidad de accidentes con heridos es de 17147, la cantidad de accidentes con muertos 237 y la cantidad de accidentes con solo dAÑOs es de 18863.

## Grupo 4 - Accidentalidad Media-Alta

```{r}
#Accidentalidad Media-Alta
dfclust_clust4 <- df_clust[df_clust$kmm.cluster == 4, ]
dfclust_clust4$total <- rowSums(dfclust_clust4[,1:3])
sum(dfclust_clust4$Con_heridos)
sum(dfclust_clust4$Con_muertos)
sum(dfclust_clust4$Solo_danos)
sum(dfclust_clust4$total)
```

Para este grupo se observa que la cantidad total de accidentes son de 54841, la cantidad de accidentes con heridos es de 32887, la cantidad de accidentes con muertos 222 y la cantidad de accidentes con solo dAÑOs es de 21732.

Finalmente, se realiz? el mapa de agrupamiento, seg?n lo mostrado anteriormente.

## Mapa de Accidentalidad en la ciudad de Medell?n por Agrupamiento

```{r}
#Se vuelve a utlizar catastro para este mapa
#Se import? el archivo .xlsx basemapa
basemapa <- read_excel("basemapa.xlsx")
base_mapa <- data_frame(basemapa)
catastro$CODIGO <- as.numeric(as.character(catastro$CODIGO))
base_mapa$Codigo <- as.numeric(as.character(base_mapa$Codigo))
#Se utiliz? 'inner join' de nuevo para unir dos bases y para as? luego generar mapa
mapa02 <- inner_join(catastro, base_mapa, by = c("CODIGO" = "Codigo"))
colorgrupos <- c("#00FF66", "#CCFF00", "#FF0000", "#0066FF")
mapa02$colores <- ifelse(mapa02$kmm.cluster == "1", "#00FF66",
                            ifelse(mapa02$kmm.cluster == "2", "#CCFF00",
                                   ifelse(mapa02$kmm.cluster == "3", "#FF0000",
                                          ifelse(mapa02$kmm.cluster == "4", "#0066FF",0))))
#Mapa final
leaflet() %>% addPolygons(data = mapa02, opacity = 0.4, color = "#545454",weight = 1, fillColor = mapa02$colores,
                          fillOpacity = 0.4, label = ~NOMBRE_BAR,
                          highlightOptions = highlightOptions(color = "#262626", weight = 3, bringToFront = T, opacity = 1),
                          popup = paste("Barrio: ", mapa02$NOMBRE_BAR, "<br>", "Grupo: ", mapa02$kmm.cluster, "<br>", "Número de Accidentes con heridos: ", mapa02$Con_heridos, "<br>", "Número de Accidentes con muertos: ", mapa02$Con_muertos, "<br>", "Número de Accidentes con solo dAÑOs: ", mapa02$Solo_danos)) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend(position = "bottomright", colors = colorgrupos, labels = c("Grupo 1: Accidentalidad Moderada", "Grupo 2: Accidentalidad Baja", "Grupo 3: Accidentalidad Alta", "Grupo 4: Accidentalidad Media-Alta"))
```