#DATA MANIPULATION ASSIGNMENT ii VICTORIA PRICE 2315903

############################################
# TASK 3: DISPLAY SINGLE POLLUTANT GRAPH


#Creating the function 
displayPollutantGraph<- function(pollutant, startMonth, startYear, 
                                 endMonth, endYear) {
  
 #Confirming the data format and filtering needed for the task 
  PortTalbotAQ$Date <- as.Date(as.character(PortTalbotAQ$Date), format="%Y-%m-%d")
  start_date <- as.Date(paste(startYear, startMonth, "01", sep="-"))
  end_date   <- as.Date(paste(endYear, endMonth, "01", sep="-"))
  
  
  df <- PortTalbotAQ[
    PortTalbotAQ$Date >= start_date &
      PortTalbotAQ$Date <= end_date &
      format(PortTalbotAQ$Date, "%d") == "01" &
      substr(PortTalbotAQ$Time, 1, 2) == "12",
  
  ]
  
 
  if (nrow(df) == 0) {
    stop("No data available for this pollutant in the given range.")
  }
  

  x <- df$Date
  y <- df[[pollutant]]
  
 
  plot(
    x, y,
    type = "l",
    col = "purple",
    xlab = "Date",
    ylab = paste(pollutant, "Level"),
    main = paste(pollutant, "Levels (Noon, 1st of Month)")
  )
  
}

############################################
# TASK 4: COMPARE MULTIPLE POLLUTANTS


#Creating the function
comparePollutants <- function(pollutants) {
  
  PortTalbotAQ$Date <- as.Date(PortTalbotAQ$Date, format = "%Y-%m-%d")
  
  df <- PortTalbotAQ[
    format(PortTalbotAQ$Date, "%d") == "01" &
      substr(PortTalbotAQ$Time, 1, 2) == "12" &
      PortTalbotAQ$Date >= as.Date("2015-01-01") &
      PortTalbotAQ$Date <= as.Date("2024-12-31"),
  ]
  
  x <- df$Date
  colours <- c("purple", "blue", "green", "orange")
  
#plotting the first pollutant
  y <- df[[pollutants[1]]]
  keep <- !is.na(x) & !is.na(y)
  ys <- lowess(as.numeric(x[keep]), y[keep], f = 0.2)
  
  if (sum(keep) < 2) {
    stop(paste("Not enough data to plot", pollutants[1]))
  }
  
  plot(
    x[keep], ys$y,
    type = "l",
    col = colours[1],
    ylim = c(0, 50), #Fixed scale for the graph
    xlab = "Date",
    ylab = "Pollutant Level",
    main = "Smoothed Variation in Pollutant Levels (Noon, 1st of Month)"
  )
  
#Loop for the other pollutants
  for (i in 2:length(pollutants)) {
    y <- df[[pollutants[i]]]
    keep <- !is.na(x) & !is.na(y)
    ys <- lowess(as.numeric(x[keep]), y[keep], f = 0.2)
    lines(x[keep], ys$y, col = colours[i])
    if (sum(keep) < 2) next
  }
  
  legend("topright", legend = pollutants, col = colours, lty = 1)
}


############################################
# TASK 5: OZONE–NITROGEN DIOXIDE CORRELATION



O3NOCorrelation <- function() {
 
  corr <- cor(PortTalbotAQ$O3, PortTalbotAQ$NO2, use = "complete.obs")
  print(paste("Correlation:", round(corr, 3)))
  
#Creating the plot WITH added correlation line 
  plot(PortTalbotAQ$NO2, PortTalbotAQ$O3,
       xlab = "NO2", ylab = "O3",
       main = "O3 vs NO2", pch = 19, col = "purple")
  abline(lm(O3 ~ NO2, data = PortTalbotAQ), col = "orange")
}

############################################
# TASK 6: T-TEST – OZONE AND WIND SPEED


#Hypothesis

# H0- mean ozone levels are equal when wind speed < 5 m/s
#     and when wind speed ≥ 5 m/s
# H1- mean ozone levels differ between the two wind speed groups

#Creating the function
O3WindTTest <- function() {
  

  O3_low  <- PortTalbotAQ$O3[PortTalbotAQ$WindSpeed < 5]
  O3_high <- PortTalbotAQ$O3[PortTalbotAQ$WindSpeed >= 5]
  
 
  test <- t.test(O3_low, O3_high)
  
  # If p-value < 0.05, reject the null hypothesis
  print(test)
}

# If the p-value is less than 0.05, the null hypothesis is rejected,
# meaning wind speed has a significant effect on ozone levels.
# If the p-value is greater than or equal to 0.05, the null hypothesis
# is not rejected, meaning no significant difference is detected.

O3WindBoxplot <- function() {
  
 
  O3_low  <- PortTalbotAQ$O3[PortTalbotAQ$WindSpeed < 5]
  O3_high <- PortTalbotAQ$O3[PortTalbotAQ$WindSpeed >= 5]
  
 
  boxplot(O3_low, O3_high,
          names = c("< 5 m s⁻¹", "≥ 5 m s⁻¹"),
          ylab = "Ozone (O3)",
          main = "Ozone Levels Under Different Wind Speeds",
          col = c("lightblue", "lightgreen"))
}




#################################MAIN SCRIPT###################################

############################################
# TASK 1: IMPORT DATA AND RENAME HEADERS


#Import Data File
PortTalbotAQ<-read.csv("Port Talbot Hourly AQ CSV.csv", header=TRUE, stringsAsFactors = FALSE)

#Print Data File
PortTalbotAQ

#Renaming headers to shorter yet understandable titles
colnames(PortTalbotAQ) <- c(
  "Date",
  "Time",
  "O3",
  "O3_status",
  "NO",
  "NO_status",
  "NO2",
  "NO2_status",
  "NOx",
  "NOx_status",
  "SO2",
  "SO2status",
  "CO",
  "COstatus",
  "PM10",
  "PM10status",
  "PM2_5",
  "PM2_5_status",
  "WindDir",
  "WindDirstatus",
  "WindSpeed",
  "WindSpeedstatus",
  "Temp",
  "Tempstatus",
  "Extra"
)

#Display new DF with new titles
PortTalbotAQ


############################################
# TASK 2: MISSING DATA PERCENTAGES

pollutantNames <- c("O3", "NO", "NO2", "NOx", "SO2", "CO", "PM10", "PM2_5")

#Converting 'No Data' entries to NA
PortTalbotAQ[PortTalbotAQ == "No data"] <- NA

#Converting all data to numeric
PortTalbotAQ$O3    <- as.numeric(PortTalbotAQ$O3)
PortTalbotAQ$NO    <- as.numeric(PortTalbotAQ$NO)
PortTalbotAQ$NO2   <- as.numeric(PortTalbotAQ$NO2)
PortTalbotAQ$NOx   <- as.numeric(PortTalbotAQ$NOx)
PortTalbotAQ$SO2   <- as.numeric(PortTalbotAQ$SO2)
PortTalbotAQ$CO    <- as.numeric(PortTalbotAQ$CO)
PortTalbotAQ$PM10  <- as.numeric(PortTalbotAQ$PM10)
PortTalbotAQ$PM2_5 <- as.numeric(PortTalbotAQ$PM2_5)

#Confirming that the conversion was correct
class(PortTalbotAQ$O3)

#Calculate the percentage of NA results for each pollutant
invalidData <- data.frame(
  Pollutant = pollutantNames,
  PercentMissing = c(
    mean(is.na(PortTalbotAQ$O3))    * 100,
    mean(is.na(PortTalbotAQ$NO))    * 100,
    mean(is.na(PortTalbotAQ$NO2))   * 100,
    mean(is.na(PortTalbotAQ$NOx))   * 100,
    mean(is.na(PortTalbotAQ$SO2))   * 100,
    mean(is.na(PortTalbotAQ$CO))    * 100,
    mean(is.na(PortTalbotAQ$PM10))  * 100,
    mean(is.na(PortTalbotAQ$PM2_5)) * 100
  )
)

#Print percentage of missing data for each pollutant 
invalidData

#Calling the function and outputting a graph for required specification
displayPollutantGraph("NO2", 7, 2015, 11, 2024)

pollutants <-c("NOx", "NO2", "SO2", "PM10")

comparePollutants(c("NOx", "NO2", "SO2", "PM10"))

O3NOCorrelation()

# If p-value < 0.05, reject the null hypothesis
O3WindTTest()

O3WindBoxplot()

############################################
# TASK 7: ANOVA – NO2 ACROSS THREE PERIODS


#Hypothesis

#H0- Mean NO2 levels are equal across all three periods
#H1- At least one period has a different mean NO2 level


anovaData <- PortTalbotAQ[, c("Date", "Time", "NO2")]



anovaData$Date <- as.Date(as.character(anovaData$Date), format = "%Y-%m-%d")

#removing NA entries
anovaData <- anovaData[!is.na(anovaData$Date), ]
sum(is.na(anovaData$Date))

anovaData$Year  <- format(anovaData$Date, "%Y")
anovaData$Month <- format(anovaData$Date, "%m")

# Filter for 17:00 November–December for certain years

anovaData <- anovaData[
  substr(as.character(anovaData$Time), 1, 2) == "17" &
    anovaData$Month %in% c("11", "12") &
    anovaData$Year %in% c("2016", "2018", "2020"),
]


anovaData$Period <- factor(anovaData$Year,
                           levels = c("2016", "2018", "2020"),
                           labels = c("Nov–Dec 2016",
                                      "Nov–Dec 2018",
                                      "Nov–Dec 2020"))

anovaData


table(anovaData$Period)

# Perform ANOVA
resaov <- aov(NO2 ~ Period, data = anovaData)
summary(resaov)

# If the p-value for Period is less than 0.05, the null hypothesis is rejected,
# indicating that NO2 levels differ significantly between at least two periods.
# If the p-value is greater than or equal to 0.05, the null hypothesis is not
# rejected and no significant difference is found.


#The Tukey test is an additional and more exact result

# Tukey test is used to identify which
# specific periods differ from ANOVA

TukeyHSD(resaov)

# Period pairs with adjusted p-values below 0.05 are considered
# significantly different.