-
Notifications
You must be signed in to change notification settings - Fork 0
/
.Rhistory
229 lines (229 loc) · 7.67 KB
/
.Rhistory
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
summary(air_mad)
#| message: false
#| warning: false
#| code-summary: Lectura del conjuto de de datos `air_mad.RDS`
library(tidyverse)
air_mad <- readr::read_rds("data/air_mad.RDS")
class(air_mad)
#| code-summary: Visualización de las primeras filas del dataset
head(air_mad)
#| code-summary: Dimensión del conjunto de datos
dim(air_mad)
#| code-summary: Dimensión del conjunto de datos
summary(air_mad)
#| code-summary: Estructura del conjunto de datos
str(air_mad)
#| code-summary: Representación de un gráfico de lineas con la evolución de temporal de los distintos contaminantes en la ciudad de Madrid (frecuencia semanal)
#| message: false
#| warning: false
#| label: fig-lineas-madrid
#| fig-cap: Evolución de contaminantes en Madrid (2013-2023)
library(lubridate)
air_mad %>%
filter(V == "V") %>%
group_by(semana = floor_date(fecha, unit = "week"), nom_mag) %>%
summarise(media_estaciones = mean(valor, na.rm = TRUE)) %>%
ggplot(aes(x = semana, y = media_estaciones)) +
geom_line(aes(color = nom_mag)) +
geom_smooth(size = 0.5, color = "red") +
labs(
x = NULL, y = "(µg/m3)", title = "Evolución de contaminantes en Madrid",
subtitle = "Concentración media semanal en las estaciones de medición",
caption = "Fuente: Portal de datos abiertos del Ayuntamiento de Madrid"
) +
theme_minimal() +
theme(legend.position = "none") +
facet_wrap(~nom_mag, scales = "free_y", ncol = 2)
#| message: false
#| warning: false
#| code-summary: Representación de un gráfico de lineas con la evolución del NO2 desde 2018 hasta la actualidad según las zonas definidas de calidad del aire en Madrid
plot_air_mad <- air_mad %>%
filter(V == "V" & nom_abv == "NO2" & fecha >= "2018-01-01") %>%
group_by(fecha, nom_mag, zona) %>%
summarise(media_estaciones = mean(valor, na.rm = TRUE)) %>%
ggplot(aes(x = fecha, y = media_estaciones)) +
geom_line(aes(color = zona)) +
geom_smooth(size = 0.5, color = "red") +
labs(
x = NULL, y = "(µg/m3)", title = "Evolución de NO2 en Madrid",
subtitle = "Concentración por zonas en las estaciones de medición",
caption = "Fuente: Portal de datos abiertos del Ayuntamiento de Madrid"
) +
theme_minimal() +
theme(legend.position = "none") +
facet_wrap(~zona, ncol = 1)
plot_air_mad
#| code-summary: Selección del contaminante NO2 y creación de un nuevo objeto llamado `no2`.
no2 <- air_mad %>%
filter(nom_abv == "NO2")
#| code-summary: Conteo de los datos validados disponibles para el NO2.
no2 %>%
count(V)
#| code-summary: Filtrado de los datos validados disponibles. Se "machaca" el objeto `no2` con el nuevo dataset.
no2 <- no2 %>%
filter(V == "V")
#| code-summary: Exploración de datos faltanes con la función is.na()
no2 %>%
count(is.na(valor))
#| message: false
#| warning: false
#| code-summary: Manipulación del data frame para su análisis. Se agrupa por zona y fecha y se calcula la media.
no2_zona <- no2 %>%
group_by(zona, fecha) %>%
summarise(valor_zona = mean(valor, na.rm = T))
#| message: false
#| warning: false
#| code-summary: Manipulación del data frame para la detección de valores atípicos
anomalias <- lapply(unique(no2_zona$zona), function(x) {
aux_data <- no2_zona %>%
filter(zona == x) %>%
select(zona, fecha, valor_zona)
aux_data <- aux_data %>%
group_by(month(fecha), year(fecha)) %>%
mutate(anomaly = !between(
valor_zona,
quantile(aux_data$valor_zona, probs = c(.25)) - 1.5 * IQR(valor_zona),
quantile(aux_data$valor_zona, probs = c(.75)) + 1.5 * IQR(valor_zona)
)) %>%
ungroup() %>%
select(zona, fecha, anomaly)
})
anomalias <- do.call(rbind, anomalias)
no2_zona <- left_join(no2_zona, anomalias, by = c("zona", "fecha"))
#| code-summary: Representación de ouliers en la zona de caliad del aire interior de la M30 en el periodo considerado.
no2_zona %>%
filter(zona == "Interior M30") %>%
ggplot(aes(fecha, valor_zona, color = anomaly)) +
geom_point() +
scale_x_date() +
scale_color_manual(values = c("black", "red")) +
labs(title = "Anomalías NO2 por zona de interes",
subtitle = "Interior M30",
x = NULL,
y = "(µg/m3)") +
theme_minimal() +
theme(legend.position = "bottom")
#| message: false
#| warning: false
#| code-summary: Manipulación de datos, creación de una nueva variable dicotómica que determina si el valor es inferior o igual a 40 µg/m3 (límite legal).
no2_zona <- no2_zona %>%
mutate(objetivo_anual = valor_zona <= 40)
no2_zona %>%
ggplot(aes(zona, fill = objetivo_anual)) +
geom_bar(position = "fill") +
scale_y_continuous(labels = scales::percent) +
labs(title = "Cumplimiento del tope de 40 µg/m3 de NO2",
subtitle = "porcentaje de cumplimiento por zona",
x = NULL,
y = "%",
fill = "Cumplimiento del objetivo") +
theme(legend.position = "bottom")
#| message: false
#| warning: false
library(zoo)
no2_zona <- no2_zona %>%
group_by(zona) %>%
mutate(valor_zona_anual = rollmean(valor_zona, k = 365, fill = NA, align = 'right'))
no2_zona %>%
filter(fecha >= "2014-01-01") %>%
ggplot(aes(fecha, valor_zona_anual, color = zona)) +
geom_line() +
scale_x_date() +
geom_hline(yintercept = 40, color = "red") +
labs(title = "Cumplimiento del tope de 40 µg/m3 de NO2",
subtitle = "media móvil anual por zona",
x = NULL,
y = "µg/m3 de NO2",
color = "Zona") +
theme(legend.position = "bottom")
#| echo: false
write_rds(no2_zona, "data/no2_zona.RDS")
#| code-summary: Se filtra la zona de calidad del aire con mayores niveles de contaminación, la zona "Interior M30" y se crea el subconjunto `no2_zona_m30`
no2_zona_m30 <- no2_zona %>%
filter(zona == "Interior M30" & anomaly == FALSE) %>%
ungroup() %>%
select(fecha, valor_zona)
#| code-summary: se realiza un gráfico de lineas para ver la evolución del NO2 en el periodo consideraro (2013-2023)
library(timetk)
no2_zona_m30 %>%
plot_time_series(fecha, valor_zona)
#| message: false
#| warning: false
#| code-summary: generación de valiables para la predicción y el horizonte temporal para dicha predicción
# el horizonte
h <- 21
# función auxilar para hacer los lags
lag_transformer <- function(data){
data %>%
tk_augment_lags(valor_zona, .lags = 1:h) %>%
ungroup()
}
# Extensión de los datos en h
no2_zona_m30_ext <- no2_zona_m30 %>%
future_frame(
.date_var = "fecha",
.length_out = h,
.bind_data = TRUE
) %>%
ungroup()
# Generación de variables
no2_zona_m30_ext <- no2_zona_m30_ext %>%
mutate(
fecha_num = normalize_vec(as.numeric(fecha)),
dow = wday(fecha, week_start = 1),
month = month(fecha),
quarter = quarter(fecha),
year = year(fecha)
) %>%
ungroup()
# Añadir lags
no2_zona_m30_ext <- no2_zona_m30_ext %>%
lag_transformer()
#| code-summary: se crea el conjunto de datos de entrenamiento `train_data` y de predicción `future_data`
# Datos de entrenamiento
train_data <- no2_zona_m30_ext %>%
drop_na()
# Datos para la prediccción
future_data <- no2_zona_m30_ext %>%
filter(is.na(valor_zona))
#| code-summary: Estructura de los datos de entrenamiento y validación
str(train_data)
#| code-summary: Estructura de los datos de predicción
str(future_data)
#| message: false
#| warning: false
#| code-summary: Creación de 4 subconjnotos (folders) para llevar a cabo el entrenamiento y, posteriomente, la evaluación de modelos de ML propuestos.
resamples_tscv <- time_series_cv(
data = train_data,
assess = h,
initial = 12 * 4 * h,
skip = 12 * 3 * h,
slice_limit = 4
)
resamples_tscv
resamples_tscv %>%
tk_time_series_cv_plan() %>%
plot_time_series_cv_plan(fecha, valor_zona, .facet_ncol = 2, .interactive = FALSE)
knitr::include_graphics(here::here("img/modelo.png"))
#| message: false
#| warning: false
# Tidymodeling
library(tidymodels)
library(modeltime)
library(modeltime.resample)
# Base Models
library(earth)
library(glmnet)
library(xgboost)
library(lightgbm)
# Core Packages
library(tidyverse)
library(lubridate)
#library(timetk)
library(bonsai)
set.seed(160191)
names(train_data)
names(train_data)
shiny::runApp('air_app')
runApp('air_app')
runApp('air_app')