~thrill_seeker/dataviz-challenge-2019

1a9978becc422221365ab3d32350842959d5ec8c — ThrillSeeker 5 months ago 8939ff0
Updating .Rmd files with present tense
2 files changed, 89 insertions(+), 91 deletions(-)

M dataviz_cat.Rmd
M dataviz_eng.Rmd
M dataviz_cat.Rmd => dataviz_cat.Rmd +45 -43
@@ 6,11 6,11 @@ output:

# Introducció

Aquest anàlisi va ser realitzat per a poder participar al [World Dataviz Challenge 2019](https://opendata-ajuntament.barcelona.cat/ca/data-viz-kobe-home).
Aquest anàlisi ha estat realitzat per a poder participar al [World Dataviz Challenge 2019](https://opendata-ajuntament.barcelona.cat/ca/data-viz-kobe-home).
Amb la recent preocupació envers temes medioambientals, i tenint en compte la recomanació explícita de la competició per a considerar aquest tema, l'elecció era prou clara.
[R](https://www.r-project.org/), un dels llenguatges de programació més estesos en l'anàlisi de dades (i el meu preferit), va ser l'escollit per a desenvolupar aquest anàlisi.
[R](https://www.r-project.org/), un dels llenguatges de programació més estesos en l'anàlisi de dades (i el meu preferit), ha estat l'escollit per a desenvolupar aquest anàlisi.

L'anàlisi seguí els següents passos:
L'anàlisi segueix els següents passos:
1. Ingesta: Ingestió de dades sobre contaminació.
2. Preprocessament: Processar les dades fins que estiguin estructurades per a l'anàlisi.
3. Exploració: Extreure patrons i evolucions de la contaminació utilitzant gràfiques.


@@ 22,14 22,14 @@ L'aplicació es pot trobar al [portal de shiny](x).

# Ingesta

En primer lloc, era necessari carregar les dades. Aquestes van ser extretes dels següents enllaços:
En primer lloc, és necessari carregar les dades. Aquestes han estat extretes dels següents enllaços:
- [Qualitat de l'aire](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/qualitat-aire-detall-bcn)
- [Estacions de mesura de contaminació](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/qualitat-aire-estacions-bcn)
- [Contaminants](https://opendata-ajuntament.barcelona.cat/data/ca/dataset/contaminants-estacions-mesura-qualitat-aire )

Tot i que existeix la possibilitat d'utilitzar una API directament per a aconseguir les dades, en aquest cas vaig preferir descarregar directament els fitxers en format .csv.
Tot i que existeix la possibilitat d'utilitzar una API directament per a aconseguir les dades, en aquest cas he preferit descarregar directament els fitxers en format .csv.

Els següents paquets van ser carregats pel desenvolupament.
Els següents paquets es carreguen pel desenvolupament.

```{r package-loading}
library(data.table)


@@ 44,10 44,12 @@ library(shinydashboard)
library(shinyWidgets)
```

De la mateixa manera que s'especifica a la URL que conté les dades sobre qualitat de l'aire, la estructura dels datasets havia canviat a partir de l'Abril del 2019.
Vaig decidir definir dues parts diferents d'ingesta: 
De la mateixa manera que s'especifica a la URL que conté les dades sobre qualitat de l'aire, l'estructura dels datasets va canviar a partir de l'Abril del 2019.
He decidit definir dues parts diferents d'ingesta: 
1. La *vella* que llegeix les dades a partir de juny de 2018 (primer arxiu disponible) fins a març de 2019 (l'últim arxiu amb la mateixa estructura).
2. La *nova* que llegeix les dades a partir d'abril de 2019 fins a setembre del 2019 (l'últim arxiu pel moment). Variables auxiliars van ser definides també.
2. La *nova* que llegeix les dades a partir d'abril de 2019 fins a setembre del 2019 (l'últim arxiu pel moment). 

Es defineixen també algunes variables auxiliars.

```{r variable-loading}
# Define variables


@@ 63,8 65,8 @@ old_files <- files_to_read[!(filter >= 201904)]
new_files <- files_to_read[(filter >= 201904)]
```

Then it was possible to load the data and inspect both structures.
For most of data processing, I tend to use [data.table package](https://rdatatable.gitlab.io/data.table/), which will be seen from now on.
Ara és possible carregar les dades i inspeccionar ambdues estructures.
Per la gran majoria de processament de dades, acostumo a utilitzar [el paquet data.table](https://rdatatable.gitlab.io/data.table/), el qual apareixerà a partir d'ara.

```{r data-loading}
# Read files


@@ 86,12 88,12 @@ str(new_quality_dt)

# Preprocessament

Sent el format completament diferent, i considerant que volia totes les combinacions d'estacions, contaminants i temps, l'enfoc va ser:
Sent el format completament diferent, i considerant que vull disposar de totes les combinacions d'estacions, contaminants i temps, l'enfoc és el següent:
1. Crear el dataset amb totes les combinacions desitjades, i separar-lo en dues parts (*vella* i *nova*).
2. Processar ambdós datasets separadament fins que es considerin estructurats, [tal com explica un gran gurú de R](https://vita.had.co.nz/papers/tidy-data.html).
3. Concatenar ambdós datasets per a crear el definitiu.

Així, el primer pas es va donar.
Així, es pot donar el primer pas.

```{r quality_dt-structure}
stations_dt <- file.path(input_path, 'Qualitat_Aire_Estacions.csv') %>% fread()


@@ 111,7 113,7 @@ quality_dt_old <- quality_dt[date < '2019-01-04']
quality_dt_new <- quality_dt[date >= '2019-01-04']
```

Després, les dades *velles* van ser processades.
Després, les dades *velles* són processades.

```{r old-quality-preprocessing}
# In order, get day, convert hours to integer, and filter by pollutant


@@ 141,7 143,7 @@ old_qual_pm10_dt <- na.omit(old_qual_pm10_dt)
old_quality_dt <- rbind(old_qual_o3_dt, old_qual_no2_dt, old_qual_pm10_dt)
```

Les dades *noves* van ser processades també.
Les dades *noves* són processades també.

```{r new-quality-preprocessing}
# In order get date, from wide (hour) to values


@@ 152,7 154,7 @@ new_quality_dt <- melt(new_quality_dt, id.vars=c('estacio', 'codi_contaminant', 
new_quality_dt[, time:=as.integer(gsub(pattern='h', replacement='', x=time)) -1]
```

Finalment, les dades van poder ser concatenades.
Finalment, les dades poden ser concatenades.

```{r dt-preprocessed}
# Join datasets


@@ 169,7 171,7 @@ quality_dt[, yearmonth:=format(date, '%Y%m')]
# Exploració

És mes fàcil entendre dades si aquestes estan agregades.
En primer lloc vaig fer una gràfica que mostrava valors mitjans per mes i contaminant, independentment de la estació on s'havia fet la mesura.
En primer lloc es pot crear una gràfica que mostra valors mitjans per mes i contaminant, independentment de la estació on s'hagi fet la mesura.

```{r monthly-eda-plot}
monthly_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, yearmonth)]


@@ 177,9 179,9 @@ m <- ggplot(data=monthly_dt) + aes(x=yearmonth, y=value, color=desc_contaminant,
plot(m)
```

Es podia veure que no hi havia dades ni per febrer ni març de 2019, i que també els nivells de O3 són superiors a la resta de contaminants, i amb una variabilitat superior també.
Es pot veure que no hi ha dades ni per febrer ni març de 2019, i que també els nivells de O3 són superiors a la resta de contaminants, i amb una variabilitat també superior.

Estava interessat amb l'evolució dels nivells de contaminació durant el dia, i per això també es va fer la gràfica de valors mitjans per hora i contaminant.
És interessant visualitzar l'evolució dels nivells de contaminació durant el dia, i per això també es crea una gràfica de valors mitjans per hora i contaminant.

```{r hourly-eda-plot}
hourly_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, time)]


@@ 187,15 189,15 @@ h <- ggplot(data=hourly_dt) + aes(x=time, y=value, color=desc_contaminant, group
plot(h)
```

Es podia veure que hi havia una clara diferència entre contaminants: sent els nivells de PM10 i NO2 són relativament baixos després del migdia, els nivells de O3 són molt alts en comparació.
També hi havia la pregunta de si hi havia diferències entre estacions, independentment de l'hora.
Es pot observar que hi ha una clara diferència entre contaminants: sent els nivells de PM10 i NO2 relativament baixos després del migdia, els nivells de O3 són molt alts en comparació.
També sorgeix la pregunta de si hi ha diferències entre estacions, independentment de l'hora.

```{r station-eda-plot}
stations_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, estacio)]
s <- ggplot(data=quality_dt) + aes(x=as.factor(estacio), y=value) + geom_boxplot() + facet_wrap(desc_contaminant~.)
plot(s)
```
En aquest cas hi havia valors extrems que no permetien veure les diferències. Es va comprovar que no hi havia cap patró, i es va procedir a excloure aquests valors per a la gràfica.
En aquest cas hi ha valors extrems que no permetien veure les diferències. Es comprova que no hi ha cap patró, i es procedeix a excloure aquests valors per a la gràfica.

```{r station-eda-plot-no-outliers}
quality_dt[value > 1000]


@@ 203,8 205,8 @@ s <- ggplot(data=quality_dt[value < 1000]) + aes(x=as.factor(estacio), y=value) 
plot(s)
```

No hi havia una clara diferència entre estacions, i per tant, entre districtes de Barcelona. El que va ser curiós és que els valors de PM10 estaven molt dispersats per sobre del percentil 75.
Havent vist estacions sense pràcticament dades, vaig decidir mantenir només les combinacions de estacions i contaminant que contenen valors.
No s'observa una clara diferència entre estacions, i per tant, entre districtes de Barcelona. El que és curiós és que els valors de PM10 estan molt dispersats per sobre del percentil 75.
Havent vist estacions sense pràcticament dades, decideixo mantenir només les combinacions de estacions i contaminant que contenen valors.

```{r exclude-combinations-without-values}
perc_na_dt <- quality_dt[, .(perc_nas=sum(is.na(value))/.N), by=.(desc_contaminant, estacio)]


@@ 212,7 214,7 @@ include_comb_dt <- perc_na_dt[perc_nas < 0.8, .(desc_contaminant, estacio)]
quality_dt <- merge(quality_dt, include_comb_dt, by=c('desc_contaminant', 'estacio'))
```

Considerant valors extrems, i tenint en compte que els mètodes de sèries temporals en general no gestionen valors no definits (NA), es va decidir intentar corregir ambdós problemes.
Considerant valors extrems, i tenint en compte que els mètodes de sèries temporals en general no gestionen valors no definits (NA), aprofito per a intentar corregir ambdós problemes amb la funció *tsclean*.

```{r remove-outliers-nas}
# Replace NA values and outliers - modelling purposes


@@ 222,7 224,7 @@ quality_dt[, imputed_values:=0]
quality_dt[is.na(original_value) | original_value != value, imputed_values:=1]
```

Després d'una búsqueda ràpida, vaig trobar [estàndards definits per la Organització Mundial de la Salut](https://apps.who.int/iris/bitstream/handle/10665/69477/WHO_SDE_PHE_OEH_06.02_eng.pdf) 
Després d'una cerca ràpida, he trobat [estàndards definits per la Organització Mundial de la Salut](https://apps.who.int/iris/bitstream/handle/10665/69477/WHO_SDE_PHE_OEH_06.02_eng.pdf) 
relatiu als 3 contaminants considerats. Així, podria ser molt útil visualitzar els nivells actuals en comparació amb els estàndars.

```{r standards-comparison}


@@ 256,18 258,18 @@ plot(no2_p)

# Modelització

Predir els nivells de contaminació pot proveïr d'informació als ajuntaments, amb la qual cosa podria mitigar nivells alts de predicció, per exemple.
Predir els nivells de contaminació pot proveïr d'informació als ajuntaments, amb la qual cosa es podria mitigar futurs nivells alts de contaminació, per exemple.

En general, predicció de sèries temporals pot no ser molt acurada per prediccions a llarg termini; per aquest motiu només els valors de octubre de 2019 són els que es prediran. 
En general, la predicció de sèries temporals pot no ser molt acurada a llarg termini; per aquest motiu només els valors d'octubre de 2019 són els que es prediran. 

Per simplificació i motius computacionals, dades diàries (mitjana de valors per dia) enlloc de dades horàries va ser utilitzat per a l'anàlisi.
Per simplificació i motius computacionals, s'utilitzaran dades diàries (mitjana de valors per dia) enlloc de dades horàries per a l'anàlisi.

Per tal de saber el potencial comportament del model sobre noves dades, va ser necessari de fer una prova a passat prèviament.
Això implica que podíem entrenar un model de juny del 2018 a novembre del 2018, predint valors de desembre de 2018, i així poder comparar-les amb la realitat.
El model que millor funcioni a passat seria l'utilitzat a futur.
Per tal de saber el potencial comportament del model sobre noves dades, és necessari fer una prova a passat prèviament.
En aquest cas, es podria entrenar un model de juny del 2018 a novembre del 2018, predint valors de desembre de 2018, i així poder comparar-los amb els valors reals.
El model que millor funcioni a passat serà l'utilitzat a futur.

En aquest cas, els dos models a comparar van ser: [ARIMA](https://en.wikipedia.org/wiki/Autoregressive_integrated_moving_average) i [Exponential Smoothing](https://en.wikipedia.org/wiki/Exponential_smoothing).
La mètrica per a comparar varis models era [SMAPE](https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error), la qual considera diferències relatives.
En aquest cas, els dos models a comparar són: [ARIMA](https://en.wikipedia.org/wiki/Autoregressive_integrated_moving_average) i [Exponential Smoothing](https://en.wikipedia.org/wiki/Exponential_smoothing).
La mètrica per a comparar els models serà [SMAPE](https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error), la qual considera diferències relatives.


```{r model-comparison}


@@ 292,7 294,7 @@ cat('Smape Autoarima', smape_autoarima)
cat('Smape ETS', smape_ets)
```

Tenint en compte que el model d'Autoarima va ser el de menor mètrica de les dues opcions, va ser el model escollit per a ser utilitzat amb les dades més recents.
Tenint en compte que el model d'Autoarima és el de menor mètrica, s'escull aquest per a ser utilitzat amb les dades més recents.

```{r predictions}
apply_model_dates <- list(train_ini=as.Date('2019/04/01'), train_end=as.Date('2019/09/30'), test_end=as.Date('2019/10/31'))


@@ 312,7 314,7 @@ daily_dt <- rbind(daily_dt, predictions_dt)

# Aplicació

L'últim pas va ser preparar el dataset que va ser utilitzat per a la aplicació interactiva. Variables importants eren les següents:
L'últim pas és preparar el conjunt de dades per a ser utilitzat en la aplicació interactiva. Considero que les variables més importants són les següents:
- Contaminant
- Districte
- Districte i coordenades de l'estació de contaminació


@@ 320,9 322,9 @@ L'últim pas va ser preparar el dataset que va ser utilitzat per a la aplicació
- Estàndars de la OMS
- Mètriques

Per motius d'experiència d'usuari, vaig pensar que seria millor visualitzar dades a nivell mensual. La pregunta llavors era, quines mètriques s'haurien de considerar per a aquest nou dataset.
Per motius d'experiència d'usuari, crec que seria millor visualitzar dades a nivell mensual. La pregunta llavors és, quines mètriques s'haurien de considerar per a aquest nou conjunt de dades.

En aquest case, les següents tres van ser les mètriques definides:
En aquest cas, es generen les següents mètriques:
- Mitjana: Valor mitjà, tenint en compte valors mitjans per dia (mitjana de mitjanes).
- Màxima: Valor màxim, tenint en compte valors mitjans per dia (màxim de mitjanes).
- Percentatge per sobre dels estàndars: Percentatge de dies el màxim dels quals està per sobre dels estàndars de la OMS.


@@ 351,7 353,7 @@ monthly_dt <- daily_dt[, .(mean_value=round(mean(value), 2), max_value=round(max
fwrite(monthly_dt, file.path(output_path, 'app_data.csv'))
```

[Shiny](https://shiny.rstudio.com/) és un paquet de R per a desenvolupar informes interactius a partir de codi R. Aquest va ser el paquet utilitzat per a la aplicació final.
[Shiny](https://shiny.rstudio.com/) és un paquet de R per a desenvolupar informes interactius a partir de codi R. Aquest serà el paquet utilitzat per a la aplicació final.
L'estructura genèrica d'una aplicació shiny és la següent:
- global.R: Carrega les dades i les processa si és necessari.
- server.R: Genera objectes d'interacció amb les dades.


@@ 359,14 361,14 @@ L'estructura genèrica d'una aplicació shiny és la següent:

El codi global.R és el que val la pena detallar.

En la aplicació, hi ha cercles que representen estacions, i ambdós radis i colors són parametritzables; per aquest motiu dues columnes auxiliars van ser generades.
A la aplicació, hi ha cercles que representen estacions, i ambdós radis i colors són parametritzables; per aquest motiu es generen dues columnes auxiliar.

La lògica de colors es va definir de la manera següent:
La lògica de colors es defineix de la manera següent:
- Si la mètrica *Max* és superior als estàndards, l'estació apareixerà de color vermell.
- Si la mètrica *Max* és superior o igual al 80% dels estàndards, l'estació apareixerà de color groc.
- Si la mètrica *Max* és inferior al 80% dels estàndars, l'estació apareixerà de color verd.

L'intuició dels radis és que els més grans indicaven valors més extrems (fossin aquests bons o dolents). La lògica concreta era la següent:
L'intuició dels radis és que els més grans indicaven valors més extrems (fossin aquests bons o dolents). La lògica concreta és la següent:
- Color verd:
	- El radi màxim es pot veure quan els valors van del 0 al 50% dels estàndars.
	- El radi decreixerà linealment quan els valors s'apropin al 80% dels estàndars, convergint al radi mínim.


@@ 409,4 411,4 @@ El principal objectiu d'aquest anàlisi ha estat mostrar el potencial de:
- Modelització de dades
- Eines de visualització

Usant-les apropiadament, és possible tenir un impacte positiu en la qualitat de vida de les persones.
Usant aquestes eines apropiadament, és possible disposar de la informació necessària per a generar un impacte positiu en la qualitat de vida de les persones.

M dataviz_eng.Rmd => dataviz_eng.Rmd +44 -48
@@ 6,11 6,11 @@ output:

# Introduction

This analysis was made in order to participate in the [World Dataviz Challenge 2019](https://opendata-ajuntament.barcelona.cat/en/data-viz-kobe-home).
This analysis has been made in order to participate in the [World Dataviz Challenge 2019](https://opendata-ajuntament.barcelona.cat/en/data-viz-kobe-home).
Given the growing environmental concerns nowadays, and the explicit recommendation of the challenge to address pollution as a topic, the choice was pretty clear.
[R](https://www.r-project.org/), one of the main programming languages used in data science (and my personal favourite), was the chosen one for this development.
[R](https://www.r-project.org/), one of the main programming languages used in data science (and my personal favourite), has been the chosen one for this development.

The analysis consisted in the following steps:
The analysis consists in the following steps:
1. Loading: Load pollution data.
2. Preprocessing: Preprocess the data until it is structured enough for analyze.
3. Exploration: Get insights from pollution trends and characteristics using plots.


@@ 22,14 22,14 @@ The application can be seen at [shiny portal](x).

# Loading

First, it was required to load the data. They were retrieved from the following links:
First, it is required to load the data. They can be retrieved from the following links:
- [Air data quality](https://opendata-ajuntament.barcelona.cat/data/en/dataset/qualitat-aire-detall-bcn)
- [Pollution measure stations](https://opendata-ajuntament.barcelona.cat/data/en/dataset/qualitat-aire-estacions-bcn)
- [Pollutants](https://opendata-ajuntament.barcelona.cat/data/en/dataset/contaminants-estacions-mesura-qualitat-aire )

Althought there is the possibility to use an API directly to retrieve the data, in this case I decided to just download the .csv files.
Althought there is the possibility to use an API directly to retrieve the data, in this case I prefer to just download the .csv files.

The required packages for development were loaded.
The required packages for development are loaded.

```{r package-loading}
library(data.table)


@@ 44,10 44,11 @@ library(shinydashboard)
library(shinyWidgets)
```

As specified in the Air Data Quality URL, dataset structure had changed completely since April 2019.
I decided then to define two different loading parts:
As specified in the Air Data Quality URL, dataset structure changed completely since April 2019. I decide then to define two different loading parts:
1. The *old* one that gets the data from the June 2018 (first file available) until March 2019 (the last with the same structure).
2. The *new* one that gets data from April 2019 to September 2019 (the last file at that moment). Auxiliary variables were also defined.
2. The *new* one that gets data from April 2019 to September 2019 (the last file at that moment).

Auxiliary variables will also be defined.

```{r variable-loading}
# Define variables


@@ 63,7 64,7 @@ old_files <- files_to_read[!(filter >= 201904)]
new_files <- files_to_read[(filter >= 201904)]
```

Then it was possible to load the data and inspect both structures.
Then it is possible to load the data and inspect both structures.
For most of data processing, I tend to use [data.table package](https://rdatatable.gitlab.io/data.table/), which will be seen from now on.

```{r data-loading}


@@ 86,12 87,12 @@ str(new_quality_dt)

# Preprocessing

Being the format completely different, and given that I wanted to have all combinations of station, pollutant, and times, the approach followed was:
Being the format completely different, and given that I want to have all combinations of station, pollutant, and times, the approach followed will be:
1. Create a dataset will all the desired combinations, and separate it into two parts (*old* and *new*).
2. Process both datasets separately for them to be considered tidy, [as explained by a great R guru](https://vita.had.co.nz/papers/tidy-data.html).
3. Concatenate both datasets to create the definite one.

Then, the first step was taken.
Then, the first step is taken.

```{r quality_dt-structure}
stations_dt <- file.path(input_path, 'Qualitat_Aire_Estacions.csv') %>% fread()


@@ 111,7 112,7 @@ quality_dt_old <- quality_dt[date < '2019-01-04']
quality_dt_new <- quality_dt[date >= '2019-01-04']
```

Afterwards, old structured data was processed.
Afterwards, old structured data is processed.

```{r old-quality-preprocessing}
# In order, get day, convert hours to integer, and filter by pollutant


@@ 141,7 142,7 @@ old_qual_pm10_dt <- na.omit(old_qual_pm10_dt)
old_quality_dt <- rbind(old_qual_o3_dt, old_qual_no2_dt, old_qual_pm10_dt)
```

New structured data was processed as well.
New structured data is processed as well.

```{r new-quality-preprocessing}
# In order get date, from wide (hour) to values


@@ 152,7 153,7 @@ new_quality_dt <- melt(new_quality_dt, id.vars=c('estacio', 'codi_contaminant', 
new_quality_dt[, time:=as.integer(gsub(pattern='h', replacement='', x=time)) -1]
```

Finally, data could be appended.
Finally, data can be appended.

```{r dt-preprocessed}
# Join datasets


@@ 169,7 170,7 @@ quality_dt[, yearmonth:=format(date, '%Y%m')]
# Exploration

It is easier to understand data if they are aggregated.
I plotted mean values by month and pollutant, regardless of the station in which they are measured.
Mean values by month and pollutant are plotted, regardless of the station in which they are measured.

```{r monthly-eda-plot}
monthly_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, yearmonth)]


@@ 177,8 178,8 @@ m <- ggplot(data=monthly_dt) + aes(x=yearmonth, y=value, color=desc_contaminant,
plot(m)
```

It could be seen that there were no data for both February 2019 and March 2019, and that O3 levels are higher than the rest, with higher variability also.
I was also interested in pollution levels evolution during the day, so I also plotted mean values by time and pollutant.
It can be seen that there are no data for both February 2019 and March 2019, and that O3 levels are higher than the rest, with higher variability also.
I am also interested in pollution levels evolution during the day, so mean values by time and pollutant are also plotted.

```{r hourly-eda-plot}
hourly_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, time)]


@@ 186,8 187,8 @@ h <- ggplot(data=hourly_dt) + aes(x=time, y=value, color=desc_contaminant, group
plot(h)
```

It could be seen that there is a clear difference between pollutants: whereas PM10 and NO2 levels are relatively low after noon, O3 levels are really high in comparison.
There was also the question regarding differences between stations, regardless of time.
It can be seen that there is a clear difference between pollutants: whereas PM10 and NO2 levels are relatively low after noon, O3 levels are really high in comparison.
The question regarding differences between stations arises as well, regardless of time.

```{r station-eda-plot}
stations_dt <- quality_dt[, .(value=mean(value, na.rm=TRUE)), by=.(desc_contaminant, estacio)]


@@ 195,8 196,7 @@ s <- ggplot(data=quality_dt) + aes(x=as.factor(estacio), y=value) + geom_boxplot
plot(s)
```

In this case there were outliers that did not permit to see differences. They were first check, to ensure that there was not
any strange pattern, and they were then removed for exploratory purposes.
In this case there are outliers that do not permit to see differences. They are first checked, to ensure that there was not any strange pattern, and they are removed for exploratory purposes.

```{r station-eda-plot-no-outliers}
quality_dt[value > 1000]


@@ 204,8 204,8 @@ s <- ggplot(data=quality_dt[value < 1000]) + aes(x=as.factor(estacio), y=value) 
plot(s)
```

There were no clear differences between stations, and therefore between Barcelona districts. What was curious was that PM10 values were really dispersed above the 75th percentile.
Having seen some stations without data, I decided to just keep combinations of stations and pollutants that contain values.
There seems not to be clear differences between stations, and therefore between Barcelona districts. What is curious is that PM10 values are really dispersed above the 75th percentile.
Having seen some stations without data, I will just keep combinations of stations and pollutants that contain values.

```{r exclude-combinations-without-values}
perc_na_dt <- quality_dt[, .(perc_nas=sum(is.na(value))/.N), by=.(desc_contaminant, estacio)]


@@ 213,8 213,7 @@ include_comb_dt <- perc_na_dt[perc_nas < 0.8, .(desc_contaminant, estacio)]
quality_dt <- merge(quality_dt, include_comb_dt, by=c('desc_contaminant', 'estacio'))
```

Considering outliers, and taking into account that time series forecasting methods in general do not handle NA values, it was decided to correct both problems.
In this case, it was done using a function called tsclean, from the forecast package.
Considering outliers, and taking into account that time series forecasting methods in general do not handle NA values, it will try to tackle both with *tsclean* function.

```{r remove-outliers-nas}
# Replace NA values and outliers - modelling purposes


@@ 224,7 223,7 @@ quality_dt[, imputed_values:=0]
quality_dt[is.na(original_value) | original_value != value, imputed_values:=1]
```

After a quick search, I found [standards defined by the World Health Organization](https://apps.who.int/iris/bitstream/handle/10665/69477/WHO_SDE_PHE_OEH_06.02_eng.pdf) 
After a quick search, I have found [standards defined by the World Health Organization](https://apps.who.int/iris/bitstream/handle/10665/69477/WHO_SDE_PHE_OEH_06.02_eng.pdf) 
related to all 3 pollutants considered: O3, NO2 and PM10. Then, it could be very useful to visualize their current levels compared to the standards defined.

```{r standards-comparison}


@@ 258,19 257,19 @@ plot(no2_p)

# Modelling

Predict future pollution levels can provide information to city councils, which can be used to mitigate high-predicted pollution levels, for example.
Predicting future pollution levels can provide information to city councils, which can be used to mitigate future high pollution levels, for example.

In general, time series forecasting may not be really accurate for long-term predictions; this is why just October 2019 values will be predicted.
Taking into account that there were no data for February and March 2019, the idea was to use from April 2019 to September 2019 to train the model.
Taking into account that there are no data for February and March 2019, the idea is to use data from April 2019 to September 2019 to train the model.

For simplification and computation reasons, daily data (mean values by day) instead of hourly data was used for the analysis.
For simplification and computation reasons, daily data (mean values by day) will be used for the analysis, instead of hourly data.

In order to know how good the model can behave on unseen data, it was necessary to backtest it previously. 
This means that we could train a model from June 2018 to November 2018, predict December 2018 values, and compare them with the reality.
The one that performs best for the past would be the one used for the future.
In order to know how good the model can behave on unseen data, it is necessary to backtest it. 
A model will be trained from June 2018 to November 2018; it will be used to predict December 2018 values, and its results will be compared with the actual values.
The one that performs best for the past will be the one used for the future.

In this case the two models to compare were: [ARIMA](https://en.wikipedia.org/wiki/Autoregressive_integrated_moving_average) and [Exponential Smoothing](https://en.wikipedia.org/wiki/Exponential_smoothing).
The metric to compare both models was [SMAPE](https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error), which considers relative difference.
In this case the two models to compare are: [ARIMA](https://en.wikipedia.org/wiki/Autoregressive_integrated_moving_average) and [Exponential Smoothing](https://en.wikipedia.org/wiki/Exponential_smoothing).
The metric to compare both models is [SMAPE](https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error), which considers relative differences.

```{r model-comparison}
gen_model_dates <- list(train_ini=as.Date('2018/06/01'), train_end=as.Date('2018/11/30'), test_end=as.Date('2018/12/31'))


@@ 294,7 293,7 @@ cat('Smape Autoarima', smape_autoarima)
cat('Smape ETS', smape_ets)
```

As Autoarima had the lowest smape of the two options, it was the chosen model to be applied to recent data.
As Autoarima has the lowest smape of the two options, it will be the model applied to recent data.

```{r predictions}
apply_model_dates <- list(train_ini=as.Date('2019/04/01'), train_end=as.Date('2019/09/30'), test_end=as.Date('2019/10/31'))


@@ 314,7 313,7 @@ daily_dt <- rbind(daily_dt, predictions_dt)

# Application

The last step was to prepare the dataset that would be used for the interactive application. Important features where the following:
The last step is to prepare the dataset that will be used for the interactive application. From my point of view, the most important features are the following:
- Pollutant
- District
- Pollution station district and coordinates


@@ 322,9 321,9 @@ The last step was to prepare the dataset that would be used for the interactive 
- WHO Standards
- Metrics

For user experience purposes, I thought it would be better to visualize monthly data. The question then was, which metrics should be taken into account for this new dataset.
For user experience purposes, I think it would be better to visualize monthly data. The question then is, which metrics should be taken into account for this new dataset.

In this case, three were the metrics chosen:
In this case, the metrics are:
- Mean: Mean value, taking into account mean values per day (mean of means).
- Max: Maximum value, taking into account mean values per day (max of means).
- Percentage above standards: Percentage of days whose maximum value is above the defined WHO standards.


@@ 353,25 352,22 @@ monthly_dt <- daily_dt[, .(mean_value=round(mean(value), 2), max_value=round(max
fwrite(monthly_dt, file.path(output_path, 'app_data.csv'))
```

**It is important to note that comparisons to WHO standards may not be completely fair. 
This is because WHO standards consider different time ranges for mean values, whereas this metrics are computed as monthly means, maximums and percentages from daily means.**

[Shiny](https://shiny.rstudio.com/) is an R package to develop interactive dashboards from R code. This was the package used for the final application.
[Shiny](https://shiny.rstudio.com/) is an R package to develop interactive dashboards from R code. This will be the package used for the final application.
The structure of a generic shiny application is the following:
- global.R: Loads the data and process them if necessary.
- server.R: Generates build objects that interact with the data.
- ui.R: Generates the user interface.

The code global.R is the one worth detailing.
The code global.R is the one worth being detailed.

In the application, there are circles representing stations, and both radius and colour are customizable, that is why two auxiliar columns will be generated.
In the application, there are circles representing stations, and both radius and colour are customizable; that is why two auxiliar columns will be generated.

Colour logic will be defined as follows:
- If *Max* metric is greater than standards, station colour will be red.
- If *Max* metric is greater or equal than 80% of standards, station colour will be yellow.
- If *Max* metric is lower than 80% of standards, station colour will be green.

Radius intuition is that bigger radius means more extreme values (either good or bad). The specific logic is the following:
Radius intuition is that bigger radius means more extreme values (either good or bad). The specific logic will be the following:
- Green colour:
	- Maximum radius can be seen when values range from 0 to 50% of standards.
	- Radius will linearly decrease as values approaches 80% of standards, converging then to minimum radius.


@@ 406,10 402,10 @@ The rest of code for shiny deployment is not detailed in this analysis. As has b

# Wrap up

The main objective of the analysis is to show the potential of using:
The main objective of the analysis has been to show the potential of using:
- Open source technologies
- Open data
- Data modelling
- Visualization tools

Using them appropriately, it is possible to have an impact on people's quality of life.
Using these tools appropriately, it is possible to provide information that can generate a positive impact on people's quality of life.