TFG CIENCIAS AMBIENTALES

Código utilizado para los 3 contaminantes del estudio del Trabajo de Fin de Grado de Álvaro Rodríguez Castrillo basado en el análisis de la contaminación atmosférica de la ciudad de Granada.

DIÓXIDO DE AZUFRE

#librerias y paquetes utilizados ########
search()
library(ggplot2)
library(highcharter)
library(modelr)
library(extrafont)
library(tidyverse)
library(lubridate)
library(scales)
library(evaluate)
library(knitr)
library(datasets)
library(data.table)
library(dplyr)
library(ggpubr)
library(rstatix) 
library(car)
library(carData)
library(nptest)
library(parallel)
library(PASWR2)
library(ggplot2)
library(lattice)
library(nortest)
library(datawizard)
library(pgirmess)
library(rmarkdown)


##preparación de los datos#####

GR_NORTE_19 <- GRANADA_NORTE_19 [complete.cases(GRANADA_NORTE_19),]
GRANADA_NORTE_20$MIN<- NULL
GR_NORTE_20 <- GRANADA_NORTE_20 [complete.cases(GRANADA_NORTE_20),]
#filtrado de datos para generar data frame solo de la fecha de la pandemia.
PANDEMIA <- GR_NORTE_20 %>% filter(MES %in% c(3,4,5,6))

NOPANDEMIA <- GR_NORTE_19 %>% filter(MES %in% c(3,4,5,6))


##estudio estadistico descriptivo####

tapply(PANDEMIA$SO2,PANDEMIA$MES,mean)
tapply(PANDEMIA$SO2,PANDEMIA$MES,sd)
tapply(PANDEMIA$SO2,PANDEMIA$MES,var)
tapply(PANDEMIA$SO2,list(PANDEMIA$MES,PANDEMIA$DIA),mean)

mean(PANDEMIA$SO2)
sd(PANDEMIA$SO2)
var(PANDEMIA$SO2)
coef_var(x= PANDEMIA$SO2, na.rm=T)


tapply(NOPANDEMIA$SO2,NOPANDEMIA$MES,mean)
tapply(NOPANDEMIA$SO2,NOPANDEMIA$MES,sd)
tapply(NOPANDEMIA$SO2,NOPANDEMIA$MES,var)
tapply(NOPANDEMIA$SO2,list(NOPANDEMIA$MES,NOPANDEMIA$DIA),mean)

mean(NOPANDEMIA$SO2)
sd(NOPANDEMIA$SO2)
var(NOPANDEMIA$SO2)
coef_var(x= NOPANDEMIA$SO2,na.rm=T)
##coef de variacion####
#formula
coef_var <- function(x, na.rm = FALSE) {
  sd(x, na.rm=na.rm) / mean(x, na.rm=na.rm)
}

coef_var(x= PANDEMIA$SO2, na.rm=T)
coef_var(x= NOPANDEMIA$SO2,na.rm=T)
tapply(PANDEMIA$SO2,PANDEMIA$MES,coef_var)
tapply(NOPANDEMIA$SO2,NOPANDEMIA$MES,coef_var)

#data wizard#
describe_distribution(round(PANDEMIA$SO2, 2))
describe_distribution(NOPANDEMIA$SO2)

#quantile
quantile(PANDEMIA$SO2, prob=c(0,0.25,0.5,0.75,1))
quantile(NOPANDEMIA$SO2, prob=c(0,0.25,0.5,0.75,1))
##tratamiento de datos por día para 2019 ####
daily19 <- NOPANDEMIA %>%
  mutate(date= make_date(ANNO,MES,DIA)) %>%
  group_by(date) %>%
  summarise(n=SO2)

daily19 %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> daily19
daily19$días <- as.factor(daily19$días)
daily19$días = factor(daily19$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(daily19$días)


meanSO219 <- daily19 %>% group_by(date) %>% summarise(n = mean(n))
meanSO219 %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> meanSO219
meanSO219$días <- as.factor(meanSO219$días)
meanSO219$días = factor(meanSO219$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(meanSO219$días)

#media para cada dia de la semana 7 valores en total
mediasemanal19 <- meanSO219 %>%  group_by(días) %>% summarise(n= mean(n))



####tratamiento de datos por dia para 2020####
daily <- PANDEMIA %>%
  mutate(date= make_date(ANNO,MES,DIA)) %>%
  group_by(date) %>%
  summarise(n=SO2)

daily %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> daily
daily$días <- as.factor(daily$días)
daily$días = factor(daily$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(daily$días)


meanSO2 <- daily %>% group_by(date) %>% summarise(n = mean(n))
meanSO2 %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> meanSO2
meanSO2$días <- as.factor(meanSO2$días)
meanSO2$días = factor(meanSO2$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(meanSO2$días)
mediasemanal20 <- meanSO2 %>%  group_by(días) %>% summarise(n= mean(n))

#data frame con los dos años juntos, valores de SO2#####
SO2test <- rbind(meanSO2, meanSO219)
date <- as.Date(SO2test$date, "%Y-%m-%d")
date
year <- as.numeric(format(date, "%Y"))
year

SO2test1 <- data.frame(
  year= rep(c("2020", "2019"), each=120),
  concentracion= c(SO2test$n[year==2020], SO2test$n[year==2019])
)

print(SO2test1)

###Plots 2020#####
boxplot20 <- meanSO2 %>%
  ggplot(aes(días, n, color= días)) +
  geom_boxplot() + 
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  labs(title = "Concentración media diaria de SO2 en µg/m3",
       subtitle= "Durante el confinamiento") +
  xlab("días") +
  ylab("Concentración media diaria de SO2 en µg/m3") +
  scale_color_discrete("Días de la semana") +
  ylim(0,10)
boxplot20

boxplot19 <- meanSO219 %>%
  ggplot(aes(días, n, color= días)) +
  geom_boxplot() + 
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  labs(title = "Concentración media diaria de SO2 en µg/m3",
       subtitle= "El mismo período del confinamiento pero en 2019") +
  xlab("días") +
  ylab("Concentración media diaria de SO2 en µg/m3") +
  scale_color_discrete("Días de la semana") +
  ylim(0,10)
boxplot19 


##representacion de los dos boxplots juntos##

ggarrange(boxplot19,boxplot20, ncol=2)
ggarrange(boxplot19, boxplot20, nrow=2)


## grafico de barras para la concentracion media cada dia de la semana
barplot20 <- mediasemanal20 %>% 
  ggplot(aes(x= días, y= n, fill= días)) +
  geom_bar(stat= "identity") +
  ylim(0,10) + 
  labs(title= "Concentración media diaria de SO2 ",
       subtitle= "Durante el confinamiento",
       x = "días",
       y= "Concentración media diaria de SO2 en µg/m3",
       col= "días") +
  theme_bw()

barplot19 <- mediasemanal19 %>% 
  ggplot(aes(x= días, y= n, fill= días)) +
  geom_bar(stat= "identity") +
  ylim(0,10) +
  labs(title= "Concentración media diaria de SO2 ",
       subtitle= "El mismo período del confinamiento pero en 2019",
       x = "días",
       y= "Concentración media diaria de SO2 en µg/m3",
       color= "Días") +
  theme_bw()
barplot19
##representacion de los dos barplots juntos##

ggarrange(barplot19,barplot20, ncol=2)
ggarrange(barplot19, barplot20, nrow=2)

##concentracion media diaria 1 linea
total20 <- meanSO2 %>%
  ggplot(aes(x=date, y=n)) +
  geom_point(mapping = aes(col=n)) + 
  geom_line() + 
  ggtitle("Concentración media diaria de SO2 en µg/m3") +
  labs(x= "Meses 2020", y= "SO2 en  µg/m3", subtitle = "Durante el confinamiento") +
  scale_y_continuous(limits = c(0,10)) +
  scale_color_continuous(name = "SO2 µg/m3") + 
  theme(plot.title = element_text(vjust = 0.5, hjust = 0.5)) +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  scale_x_date(date_breaks = "5 days",
               labels = label_date_short(format = c(NA, "%B", "%d", "%H:%M")),
               expand = c(0.005,0.005))
total20

total19 <- meanSO219 %>%
  ggplot(aes(x=date, y=n)) +
  geom_point(mapping = aes(col=n)) + 
  geom_line() + 
  ggtitle("Concentración media diaria de SO2 en µg/m3") +
  labs(x= "Meses 2019", y= "SO2 en  µg/m3", subtitle= "El mismo período del confinamiento pero en 2019") +
  scale_y_continuous(limits = c(0,10)) +
  scale_color_continuous(name = "SO2 µg/m3") + 
  theme(plot.title = element_text(vjust = 0.5, hjust = 0.5)) +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  scale_x_date(date_breaks = "5 days",
               labels = label_date_short(format = c(NA, "%B", "%d", "%H:%M")),
               expand = c(0.005,0.005))

total19 


##representacion de los dos lineas medias juntos##

ggarrange(total19,total20, ncol=2)
ggarrange(total19, total20, nrow=2)



#boxplot los dos años juntos
ggboxplot(SO2test1, x="year", y= "concentracion",
          color= "year",
          ylab= "Concentración media diaria de SO2 en µg/m3", xlab= "Año") +
  stat_summary(fun = mean, geom = "point", col="red") +
  stat_summary(fun= mean, geom= "text", col= "red",
               vjust= -1 , aes(label= paste("Mean:", round(..y.., digits = 2))))

#histograma los dos años juntos
##USAR ESTE##
ggplot(data= SO2test1, aes(concentracion)) +
  geom_histogram(bindwith= 10, aes(fill= ..count..), col="black") +
  labs(x= "Concentración media diaria de SO2 en µg/m3",
       y= "Frencuencia") +
  scale_color_continuous(name = "SO2 µg/m3") +
  facet_grid(year  ~.)
##USAR ESTE ##


hist(SO2test1$concentracion[year==2019],
     main= "Histograma de frecuencias de la concentración media de SO2 en µg/m3  en 2019",
     col= "lightblue",
     xlab= "Concentración SO2")
qqnorm(SO2test1$concentracion[year==2019])
qqline(SO2test1$concentracion[year==2019], col = 2)

#para 2020
hist(SO2test1$concentracion[year==2020],
     main= "Histograma de frecuencias de la concentración media de SO2 en µg/m3 en 2020 ",
     col= "lightblue",
     xlab= "Concentración SO2")

qqnorm(SO2test1$concentracion[year==2020])
qqline(SO2test1$concentracion[year==2020], col = 2)
## 2 años juntos

hist(SO2test1$concentracion,
     main= "Histograma de frecuencias de la concentración media de SO2 en µg/m3 en 2019 y 2020",
     col= "lightblue",
     xlab= "Concentración SO2")

##test de hipotesis####
#para comprobar normalidad# 
lillie.test(SO2test1$concentracion[year==2019])
lillie.test(SO2test1$concentracion[year==2020])
lillie.test(SO2test1$concentracion)$p.value

# para comprobar diferencias significativas entre n
wilcox.test(concentracion~year, data=SO2test1)
# para comprobar diferencias significativas entre días

kruskal.test(n~días, data= meanSO2) 
kruskal.test(n~días, data= meanSO219) 
# para 2020 
pairwise.wilcox.test(x= meanSO2$n, g= meanSO2$días, p.adjust.method = "bonferroni", exact= FALSE)
kruskalmc(n~días, data=meanSO2)
#para 2019
pairwise.wilcox.test(x= meanSO219$n, g= meanSO219$días, p.adjust.method = "bonferroni", exact= FALSE)
kruskalmc(n~días, data=meanSO219)

DIÓXIDO DE NITRÓGENO

##estudio estadistico descriptivo####

tapply(PANDEMIA$NO2,PANDEMIA$MES,mean)
tapply(PANDEMIA$NO2,PANDEMIA$MES,sd)
tapply(PANDEMIA$NO2,PANDEMIA$MES,var)
tapply(PANDEMIA$NO2,list(PANDEMIA$MES,PANDEMIA$DIA),mean)

mean(PANDEMIA$NO2)
sd(PANDEMIA$NO2)
var(PANDEMIA$NO2)
coef_var(x= PANDEMIA$NO2,na.rm=T)

tapply(NOPANDEMIA$NO2,NOPANDEMIA$MES,mean)
tapply(NOPANDEMIA$NO2,NOPANDEMIA$MES,sd)
tapply(NOPANDEMIA$NO2,NOPANDEMIA$MES,var)
tapply(NOPANDEMIA$NO2,list(NOPANDEMIA$MES,NOPANDEMIA$DIA),mean)

mean(NOPANDEMIA$NO2)
sd(NOPANDEMIA$NO2)
var(NOPANDEMIA$NO2)
coef_var(x= NOPANDEMIA$NO2,na.rm=T)
##coef de variacion####
#formula
coef_var <- function(x, na.rm = FALSE) {
  sd(x, na.rm=na.rm) / mean(x, na.rm=na.rm)
}

coef_var(x= PANDEMIA$NO2, na.rm=T)
coef_var(x= NOPANDEMIA$NO2,na.rm=T)
tapply(PANDEMIA$NO2,PANDEMIA$MES,coef_var)
tapply(NOPANDEMIA$NO2,NOPANDEMIA$MES,coef_var)

#data wizard#
describe_distribution(PANDEMIA$NO2)
describe_distribution(NOPANDEMIA$NO2)
##tratamiento de datos por día para 2019 ####
daily19 <- NOPANDEMIA %>%
  mutate(date= make_date(ANNO,MES,DIA)) %>%
  group_by(date) %>%
  summarise(n=NO2)

daily19 %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> daily19
daily19$días <- as.factor(daily$días)
daily19$días = factor(daily19$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(daily19$días)


meanno219 <- daily19 %>% group_by(date) %>% summarise(n = mean(n))
meanno219 %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> meanno219
meanno219$días <- as.factor(meanno219$días)
meanno219$días = factor(meanno219$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(meanno219$días)

#media para cada dia de la semana 7 valores en total
mediasemanal19 <- meanno219 %>%  group_by(días) %>% summarise(n= mean(n))

#data frame con los dos años juntos, valores de NO2
no2test <- rbind(meanno2, meanno219)
date <- as.Date(no2test$date, "%Y-%m-%d")
date
year <- as.numeric(format(date, "%Y"))
year

no2test1 <- data.frame(
  year= rep(c("2020", "2019"), each=120),
  concentracion= c(no2test$n[year==2020], no2test$n[year==2019])
)

print(NO2test1)

####tratamiento de datos por dia para 2020####
daily <- PANDEMIA %>%
  mutate(date= make_date(ANNO,MES,DIA)) %>%
  group_by(date) %>%
  summarise(n=NO2)

daily %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> daily
daily$días <- as.factor(daily$días)
daily$días = factor(daily$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(daily$días)


meanno2 <- daily %>% group_by(date) %>% summarise(n = mean(n))
meanno2 %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> meanno2
meanno2$días <- as.factor(meanno2$días)
meanno2$días = factor(meanno2$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(meanno2$días)
mediasemanal20 <- meanno2 %>%  group_by(días) %>% summarise(n= mean(n))
###Plots 2020#####
boxplot20 <- meanno2 %>%
  ggplot(aes(días, n, color= días)) +
  geom_boxplot()  +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  labs(title = "Diagrama de cajas y bigotes de la concentración media diaria NO2",
       subtitle= "Durante el confinamiento") +
  xlab("Días de la semana") +
  ylab("Concentración media diaria NO2 µg/m3") +
  scale_color_discrete("Días de la semana") +
  ylim(0,60)
boxplot20

boxplot19 <- meanno219 %>%
  ggplot(aes(días, n, color= días)) +
  geom_boxplot() + 
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  labs(title = "Diagrama de cajas y bigotes de la concentración media diaria NO2",
  subtitle = "El mismo período del confinamiento pero en 2019") +
  xlab("Días de la semana") +
  ylab("Concentración media diaria NO2 µg/m3") +
  scale_color_discrete("Días de la semana") +
  ylim(0,60)
boxplot19 


##representacion de los dos boxplots juntos##

ggarrange(boxplot19,boxplot20, ncol=2)
ggarrange(boxplot19, boxplot20, nrow=2)


## grafico de barras para la concentracion media cada dia de la semana
barplot20 <- mediasemanal20 %>% 
  ggplot(aes(x= días, y= n, fill= días)) +
  geom_bar(stat= "identity") +
  ylim(0,50) + 
  labs(title= "Concentración media diaria de NO2 ",
       subtitle= "Durante el confinamiento",
       x = "días",
       y= "Concentración media diaria de NO2 en µg/m3",
       col= "días") +
  theme_bw()
barplot19 <- mediasemanal19 %>% 
  ggplot(aes(x= días, y= n, fill= días)) +
  geom_bar(stat= "identity") +
  ylim(0,50) +
  labs(title= "Concentración media diaria de NO2 ",
       subtitle= "El mismo período del confinamiento pero en 2019",
       x = "días",
       y= "Concentración media diaria de NO2 en µg/m3",
       col= "días") +
  theme_bw()

##representacion de los dos barplots juntos##

ggarrange(barplot19,barplot20, ncol=2)
ggarrange(barplot19, barplot20, nrow=2)

##concentracion media diaria 1 linea
total20 <- meanno2 %>%
  ggplot(aes(x=date, y=n)) +
  geom_point(mapping = aes(col=n)) + 
  geom_line() + 
  ggtitle("Concentración media diaria de NO2 en µg/m3", subtitle= "Durante el confinamiento") +
  labs(x= "Meses 2020", y= "Concentración media diaria NO2 en µg/m3") +
  scale_y_continuous(limits = c(0,90)) +
  scale_color_continuous(name = "NO2 µg/m3") + 
  theme(plot.title = element_text(vjust = 0.5,hjust = 0.5)) +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  scale_x_date(date_breaks = "5 days",
               labels = label_date_short(format = c(NA, "%B", "%d", "%H:%M")),
               expand = c(0.005,0.005))
total20

total19 <- meanno219 %>%
  ggplot(aes(x=date, y=n)) +
  geom_point(mapping = aes(col=n)) + 
  geom_line() + 
  ggtitle("Concentración media diaria de NO2 en µg/m3", subtitle= "El mismo período del confinamiento pero en 2019") +
  labs(x= "Meses 2020", y= "Concentración media diaria NO2 en µg/m3") +
  scale_y_continuous(limits = c(0,90)) +
  scale_color_continuous(name = "NO2 µg/m3") + 
  theme(plot.title = element_text(vjust = 0.5,hjust = 0.5)) +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  scale_x_date(date_breaks = "5 days",
               labels = label_date_short(format = c(NA, "%B", "%d", "%H:%M")),
               expand = c(0.005,0.005))

total19 


##representacion de los dos lineas medias juntos##

ggarrange(total19,total20, ncol=2)
ggarrange(total19, total20, nrow=2)



#boxplot los dos años juntos
ggboxplot(no2test1, x="year", y= "concentracion",
          color= "year",
          ylab= "Concentración NO2", xlab= "Año") +
  stat_summary(fun = mean, geom = "point", col="red") +
  stat_summary(fun= mean, geom= "text", col= "red",
               vjust= -2, aes(label= paste("Mean:", round(..y.., digits = 2))))

#histograma los dos años juntos
ggplot(data= no2test1, aes(concentracion)) +
  geom_histogram(bindwith= 10, aes(fill= ..count..), col="black") +
  labs(x= "Concentración media diaria de NO2 en µg/m3",
       y= "Frencuencia") +
  scale_color_continuous(name = "NO2 µg/m3") +
  facet_grid(year  ~.)

hist(no2test1$concentracion[year==2019],
     main= "Histograma de frecuencias de la concentración media de NO2 en µg/m3  en 2019",
     col= "lightblue",
     xlab= "Concentración NO2")
qqnorm(no2test1$concentracion[year==2019])
qqline(no2test1$concentracion[year==2019], col = 2)

#para 2020
hist(no2test1$concentracion[year==2020],
     main= "Histograma de frecuencias de la concentración media de NO2 en µg/m3 en 2020 ",
     col= "lightblue",
     xlab= "Concentración NO2")

qqnorm(no2test1$concentracion[year==2020])
qqline(no2test1$concentracion[year==2020], col = 2)
## 2 años juntos

hist(no2test1$concentracion,
     main= "Histograma de frecuencias de la concentración media de NO2 en µg/m3 en 2019 y 2020",
     col= "lightblue",
     xlab= "Concentración NO2")

##test de hipotesis####
#para comprobar normalidad# 
lillie.test(no2test1$concentracion[year==2019])
lillie.test(no2test1$concentracion[year==2020])
lillie.test(no2test1$concentracion)

# para comprobar diferencias significativas entre n
wilcox.test(concentracion~year, data=no2test1)$p.value

# para comprobar diferencias significativas entre días

kruskal.test(n~días, data= meanno2)  

kruskal.test(n~días, data= meanno219) 
# para 2020 
pairwise.wilcox.test(x= meanno2$n, g= meanno2$días, p.adjust.method = "bonferroni", exact= FALSE)
kruskalmc(n~días, data=meanno2)
#para 2019
pairwise.wilcox.test(x= meanno219$n, g= meanno219$días, p.adjust.method = "bonferroni", exact= FALSE)
kruskalmc(n~días, data=meanno219)

PARTÍCULAS EN SUSPENSIÓN

##estudio estadistico descriptivo####

tapply(PANDEMIA$PART,PANDEMIA$MES,mean)
tapply(PANDEMIA$PART,PANDEMIA$MES,sd)
tapply(PANDEMIA$PART,PANDEMIA$MES,var)
tapply(PANDEMIA$PART,list(PANDEMIA$MES,PANDEMIA$DIA),mean)

mean(PANDEMIA$PART)
sd(PANDEMIA$PART)
var(PANDEMIA$PART)
coef_var(x= PANDEMIA$PART, na.rm=T)

tapply(NOPANDEMIA$PART,NOPANDEMIA$MES,mean)
tapply(NOPANDEMIA$PART,NOPANDEMIA$MES,sd)
tapply(NOPANDEMIA$PART,NOPANDEMIA$MES,var)
tapply(NOPANDEMIA$PART,list(NOPANDEMIA$MES,NOPANDEMIA$DIA),mean)

mean(NOPANDEMIA$PART)
sd(NOPANDEMIA$PART)
var(NOPANDEMIA$PART)
coef_var(x= NOPANDEMIA$PART, na.rm=T)

##coef de variacion####
#formula
coef_var <- function(x, na.rm = FALSE) {
  sd(x, na.rm=na.rm) / mean(x, na.rm=na.rm)
}

coef_var(x= PANDEMIA$PART, na.rm=T)
coef_var(x= NOPANDEMIA$PART,na.rm=T)
tapply(PANDEMIA$PART,PANDEMIA$MES,coef_var)
tapply(NOPANDEMIA$PART,NOPANDEMIA$MES,coef_var)

#data wizard#
describe_distribution(PANDEMIA$PART)
describe_distribution(NOPANDEMIA$PART)
##tratamiento de datos por día para 2019 ####
daily19 <- NOPANDEMIA %>%
  mutate(date= make_date(ANNO,MES,DIA)) %>%
  group_by(date) %>%
  summarise(n=PART)

daily19 %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> daily19
daily19$días <- as.factor(daily19$días)
daily19$días = factor(daily19$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(daily19$días)


meanpart19 <- daily19 %>% group_by(date) %>% summarise(n = mean(n))
meanpart19 %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> meanpart19
meanpart19$días <- as.factor(meanpart19$días)
meanpart19$días = factor(meanpart19$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(meanpart19$días)

#media para cada dia de la semana 7 valores en total
mediasemanal19 <- meanpart19 %>%  group_by(días) %>% summarise(n= mean(n))



####tratamiento de datos por dia para 2020####
daily <- PANDEMIA %>%
  mutate(date= make_date(ANNO,MES,DIA)) %>%
  group_by(date) %>%
  summarise(n=PART)

daily %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> daily
daily$días <- as.factor(daily$días)
daily$días = factor(daily$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(daily$días)


meanpart <- daily %>% group_by(date) %>% summarise(n = mean(n))
meanpart %>%
  mutate(días= weekdays(date, abbreviate = FALSE)) -> meanpart
meanpart$días <- as.factor(meanpart$días)
meanpart$días = factor(meanpart$días, levels=c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))
str(meanpart$días)
mediasemanal20 <- meanpart %>%  group_by(días) %>% summarise(n= mean(n))

#data frame con los dos años juntos, valores de PART#####
parttest <- rbind(meanpart, meanpart19)
date <- as.Date(parttest$date, "%Y-%m-%d")
date
year <- as.numeric(format(date, "%Y"))
year

parttest1 <- data.frame(
  year= rep(c("2020", "2019"), each=120),
  concentracion= c(parttest$n[year==2020], parttest$n[year==2019])
)

print(parttest1)

###Plots 2020#####
boxplot20 <- meanpart %>%
  ggplot(aes(días, n, color= días)) +
  geom_boxplot() + 
  theme(plot.title = element_text(vjust = 0.2, hjust = 0.2)) +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  labs(title = "Concentración media diaria de partículas en suspensión en 2020") +
  xlab("días") +
  ylab("Concentración media diaria de partículas en suspensión en µg/m3") +
  scale_color_discrete("Días de la semana") +
  ylim(0,40)
boxplot20

boxplot19 <- meanpart19 %>%
  ggplot(aes(días, n, color= días)) +
  geom_boxplot() + 
  theme(plot.title = element_text(vjust = 0.2, hjust = 0.2)) +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  labs(title = "Concentración media diaria de partículas en suspensión en 2019") +
  xlab("días") +
  ylab("Concentración media diaria de partículas en suspensión en µg/m3") +
  scale_color_discrete("Días de la semana") +
  ylim(0,40)
boxplot19 


##representacion de los dos boxplots juntos##

ggarrange(boxplot19,boxplot20, ncol=2)
ggarrange(boxplot19, boxplot20, nrow=2)


## grafico de barras para la concentracion media cada dia de la semana
barplot20 <- mediasemanal20 %>% 
  ggplot(aes(x= días, y= n, fill= días)) +
  geom_bar(stat= "identity") +
  ylim(0,30) + 
  labs(title= "Concentración media diaria de partículas en suspensión ",
       subtitle= "Durante el confinamiento",
       x = "días",
       y= "Concentración media diaria de partículas en suspensión en µg/m3",
       col= "días") +
  theme_bw()

barplot19 <- mediasemanal19 %>% 
  ggplot(aes(x= días, y= n, fill= días)) +
  geom_bar(stat= "identity") +
  ylim(0,30) +
  labs(title= "Concentración media diaria de partículas en suspensión ",
       subtitle= "El mismo período del confinamiento pero en 2019",
       x = "días",
       y= "Concentración media diaria de partículas en suspensión en µg/m3",
       color= "Días") +
  theme_bw()
barplot19
##representacion de los dos barplots juntos##

ggarrange(barplot19,barplot20, ncol=2)
ggarrange(barplot19, barplot20, nrow=2)

##concentracion media diaria 1 linea
total20 <- meanpart %>%
  ggplot(aes(x=date, y=n)) +
  geom_point(mapping = aes(col=n)) + 
  geom_line() + 
  ggtitle("Concentración media diaria de partículas en suspensión en µg/m3") +
  labs(x= "Meses 2020", y= "Partículas en suspensión en  µg/m3", subtitle = "Durante el confinamiento") +
  scale_y_continuous(limits = c(0,90)) +
  scale_color_continuous(name = "Partículas en suspensión µg/m3") + 
  theme(plot.title = element_text(vjust = 0.5, hjust = 0.5)) +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  scale_x_date(date_breaks = "5 days",
               labels = label_date_short(format = c(NA, "%B", "%d", "%H:%M")),
               expand = c(0.005,0.005))
total20

total19 <- meanpart19 %>%
  ggplot(aes(x=date, y=n)) +
  geom_point(mapping = aes(col=n)) + 
  geom_line() + 
  ggtitle("Concentración media diaria de partículas en suspensión en µg/m3") +
  labs(x= "Meses 2019", y= "Partículas en suspensión en  µg/m3", subtitle= "El mismo período del confinamiento pero en 2019") +
  scale_y_continuous(limits = c(0,90)) +
  scale_color_continuous(name = "Partículas en suspensión µg/m3") + 
  theme(plot.title = element_text(vjust = 0.5, hjust = 0.5)) +
  theme(axis.title.x = element_text(face="bold", vjust=0.5)) +
  theme(axis.title.y = element_text(face="bold", vjust=0.5)) +
  scale_x_date(date_breaks = "5 days",
               labels = label_date_short(format = c(NA, "%B", "%d", "%H:%M")),
               expand = c(0.005,0.005))

total19 


##representacion de los dos lineas medias juntos##

ggarrange(total19,total20, ncol=2)
ggarrange(total19, total20, nrow=2)



#boxplot los dos años juntos
ggboxplot(parttest1, x="year", y= "concentracion",
          color= "year",
          ylab= "Concentración media diaria de partículas en suspensión en µg/m3", xlab= "Año") +
  stat_summary(fun = mean, geom = "point", col="red") +
  stat_summary(fun= mean, geom= "text", col= "red",
               vjust= -1 , aes(label= paste("Mean:", round(..y.., digits = 2))))

#histograma los dos años juntos
##USAR ESTE##
ggplot(data= parttest1, aes(concentracion)) +
  geom_histogram(bindwith= 10, aes(fill= ..count..), col="black") +
  labs(x= "Concentración media diaria de partículas en suspensión en µg/m3",
       y= "Frencuencia") +
  scale_color_continuous(name = "Partículas en suspensión µg/m3") +
  facet_grid(year  ~.)
##USAR ESTE ##


hist(parttest1$concentracion[year==2019],
     main= "Histograma de frecuencias de la concentración media de SO2 en µg/m3  en 2019",
     col= "lightblue",
     xlab= "Concentración SO2")
qqnorm(parttest1$concentracion[year==2019])
qqline(parttest1$concentracion[year==2019], col = 2)

#para 2020
hist(parttest1$concentracion[year==2020],
     main= "Histograma de frecuencias de la concentración media de SO2 en µg/m3 en 2020 ",
     col= "lightblue",
     xlab= "Concentración SO2")

qqnorm(parttest1$concentracion[year==2020])
qqline(parttest1$concentracion[year==2020], col = 2)
## 2 años juntos

hist(parttest1$concentracion,
     main= "Histograma de frecuencias de la concentración media de SO2 en µg/m3 en 2019 y 2020",
     col= "lightblue",
     xlab= "Concentración SO2")

##test de hipotesis####
#para comprobar normalidad# 
lillie.test(parttest1$concentracion[year==2019])
lillie.test(parttest1$concentracion[year==2020])
lillie.test(parttest1$concentracion)

# para comprobar diferencias significativas entre n
wilcox.test(concentracion~year, data=parttest1)
# para comprobar diferencias significativas entre días

kruskal.test(n~días, data= meanpart) 
kruskal.test(n~días, data= meanpart19) 
# para 2020 
pairwise.wilcox.test(x= meanpart$n, g= meanpart$días, p.adjust.method = "bonferroni", exact= FALSE)
kruskalmc(n~días, data=meanpart)

#para 2019
pairwise.wilcox.test(x= meanpart19$n, g= meanpart19$días, p.adjust.method = "bonferroni", exact= FALSE)
kruskalmc(n~días, data=meanpart19)