#install and load package 
install.packages('readxl')
install.packages('dplyr')
install.packages("ggplot2") 
install.packages("lubridate") 
install.packages('janitor')
install.packages("rlang")

library('readxl')
library('dplyr')
library('ggplot2') 
library('lubridate')
library('janitor')
library('rlang')

#importing the dataset
aqpt <- read_excel('Port Talbot Hourly AQ.xlsx', skip = 10)
aqpt <- aqpt[-c(87649:87652), ]
colnames(aqpt) #displays the headings
newcol <- c("Date", "Time", "Ozone", "Status", "NO", "Status", "NO2", "Status", "NO_to_NO2", "Status", "SO2", "Status", "CO", "Status", "PM10_per_hour", "Status", "PM2point5_per_hour", "Status", "Wind_Direction", "Status", "Wind_Speed", "Status", "Temp", "Status") #sets variables for the row
names(aqpt) <- newcol #renames the row
aqpt <- clean_names(aqpt)
aqpt #shows dataset
#change date/time to date and time
aqpt <- aqpt |>
  mutate(
    dt = ymd_hms(paste(date, time))
  )
#clean dataframe
aqpt <- aqpt |>
  mutate(across(where(is.character), ~na_if(.x, "No data")))
aqpt <- aqpt |> 
  mutate(across(c("ozone", "no", "no2", "no_to_no2", "so2", "co", "pm10_per_hour", "pm2point5_per_hour", "wind_direction", "wind_speed", "temp"), as.numeric))



pollutecol <- c("ozone", "no", "no2", "so2", "co", "pm10_per_hour", "pm2point5_per_hour") #collects the specific columns
#function to group NA
missing <- function(x){
  x == "No data" | is.na(x)
}

#dataframe to display the percentage
invdata <- data.frame(
  PercentageMissing = sapply(aqpt[pollutecol], function(x){
    mean(missing(x) * 100)
  })
)

invdata #shows dataframe

#question 3 function for displaying pollutant
displaypolluteg <- function(pollute, mstart, mend, ystart, yend){
  #filtering to specified timeframe
  df <- aqpt |>
    filter(
      hour(dt) == 12,
      minute(dt) == 0,
      second(dt) == 0,
      day(dt) == 1
      )
  
  #declaring start and end dates
  dstart <- ymd(paste(ystart, mstart, '01', sep = '-'))
  dend <-  ymd(paste(yend, mend, '01', sep = '-'))
  
  df <- df |>
    filter(dt >= dstart,
           dt <= dend
           )
  
  #plot line and scatter graph
  
ggplot(df, aes(x = dt, y= .data[[pollute]])) +
  geom_line(color = 'green') +
  geom_point(colour = 'black')+
labs(
  title = 'Variation in selected pollutant starting at 1st of noon each month',
  x = 'Date',
  y = pollute
) +
  theme_minimal()
}

#display selected data using function
displaypolluteg("no2", 7, 11, 2015, 2024)

#question 4 function to compare pollutants in the selected timeframe
comppol <- function(pollutevect){
  df <- aqpt |>
    filter(
      hour(dt) == 12,
      minute(dt) == 0,
      second(dt) == 0,
      day(dt) == 1
    )
#displays the graph for easy reading
  dfl <- df |>
    select(dt, all_of(pollutevect)) |>
    tidyr::pivot_longer(
      cols = all_of(pollutevect),
      names_to = 'pollutant',
      values_to = 'value'
    )
  #create a multi-series graph
  ggplot(dfl, aes(x = dt, y = value, colour = pollutant)) +
    geom_line() +
    labs(
      title = "comparison of chosen pollutants",
      x = 'date',
      y = 'concentration'
    ) +
    theme_minimal()
}
#comparing the 4 pollutants in question 4
comppol(c('no2','no','so2', 'co'))

#question 5 function to compare Nitrogen Dioxide and Ozone
ozcor <- function() {
  df <- aqpt |>
    filter(!is.na(ozone), !is.na(no2))
  
  #plot
  plot(df$no2, 
       df$ozone,
       xlab = 'NO2',
       ylab = 'Ozone',
       main = 'Ozone and Nitrogen Dioxide comparison',
       pch = 16,
       col = 'lightblue'
  )
  #creates the correlation dots using selected data
  corval <- cor(df$no2, df$ozone)
  paste('correlation between Ozone and Nitrogen Dioxide', round(corval, 3))
  
  invisible(corval)
}

#displays the graph of the correlation
ozcor()

# question 6 t-test for ozone vs wind speed

#hypotesis
#H0 - ozone values when wind < 5m/s is equal to ozone values when wind greater or equal to 5 m/s
#H1 - ozone values when wind >= 5m/s is equal to ozone values when wind greater or equal to 5 m/s


lwind <- aqpt$ozone[aqpt$wind_speed < 5]
hwind <- aqpt$ozone[aqpt$wind_speed >= 5]

#t test for two sample 
tres <- t.test(lwind, hwind)

#displays the t-test results
tres

#P-values < 0.05 which rejects H0. There is no significant difference between the two values

# question 7 ANOVA test for Nitrogen dioxide between periods of times and dates

#hypothesis
# H0 - Nitrogen Dixoide is equal between the three periods of dates at 17:00
# H1 - at least one period at 17:00 has a different Nitrogen Dioxide mean

#filters between 1st nov to 31st dec at 17:00 for the selected years
getpdata <- function(y){
  aqpt |>
  filter(
    year(dt) == y,
    month(dt) %in% c(11,12),
    day (dt) >= 1,
    day (dt) <= 31,
    hour(dt) == 17,
    minute(dt) == 0
  ) |>
    mutate(period = paste0("novdec", y))
}

#defines sections of the years in the data
data16 <- getpdata(2016)
data18 <- getpdata(2018)
data20 <- getpdata(2020)

#collates rows and removes NA
anovat <- bind_rows(data16, data18, data20) |>
  filter(!is.na(no2))

#converts to factor
anovat$period <- factor(anovat$period,
                        levels = c('novdec2016', 'novdec2018', 'novdec2020'))

#ANOVA test
manova <- aov(no2 ~ period, data = anovat)
summary(manova)

#ANOVA test shows Nitrogen Dioxide has significant difference between Nov-Dec 2016, 2018, 2020 at 17:00
#Reject H0 

#Perform a POST-HOC test: Tukey method
phoc <- TukeyHSD(manova)
phoc