--- title: "Airpollution" author: "Nisarg Jani" date: "2025-11-06" output: word_document: default html_document: default --- I write the code for each Question in separate section ```{r} # 1 – Create PortTalbotAQ data frame with renamed headers # Read the CSV file (assume the file "Port Talbot Hourly AQ.csv" # is in the same folder as this script) PortTalbotAQ <- read.csv("Port_Talbot_Hourly_AQ.csv", header = TRUE, stringsAsFactors = FALSE) # Rename headers to shorter but understandable names colnames(PortTalbotAQ) <- c( "Date", # Date of measurement "Time", # Time of measurement "O3", # Ozone "NO", # Nitric oxide "NO2", # Nitrogen dioxide "NOx", # Nitrogen oxides "SO2", # Sulphur dioxide "CO", # Carbon monoxide "PM10", # Particulate matter <10μm "PM2_5", # Particulate matter <2.5μm "WindSpeed", # Wind speed "WindDirection", # Wind direction "Temperature", # Temperature "Units" # Units of measurement ) # Display the resulting data frame (first few rows) print("PortTalbotAQ data frame:") head(PortTalbotAQ) #In this task, I imported the raw air-quality dataset using read.csv() and renamed the column headers to shorter, clearer labels. This step was important because the original column names were long and difficult to reference. Cleaning and renaming the headers made the data easier to handle in later analysis and ensured consistency when applying functions across pollutants. ``` ```{r} # Task 2 – Create invalidData data frame showing % of "No data" readings per pollutant # Convert all "No data" strings to proper NA values PortTalbotAQ[PortTalbotAQ == "No data"] <- NA # List of pollutant columns to analyse pollutants <- c("O3", "NO", "NO2", "NOx", "SO2", "CO", "PM10", "PM2_5") # Calculate the percentage of missing data for each pollutant invalid_percent <- sapply(PortTalbotAQ[pollutants], function(x) { mean(is.na(x)) * 100 }) # Create a clean data frame of results invalidData <- data.frame( Pollutant = pollutants, Missing_Percentage = round(invalid_percent, 2) ) # Display the resulting data frame print("Percentage of missing readings per pollutant:") print(invalidData) #In this task, I calculated how much data was missing for each pollutant. I first replaced all “No data” text entries with proper NA values so R could recognise them as missing. Then, I used sapply() with mean(is.na(x)) to find the percentage of missing readings for each pollutant column. This helped identify which pollutants had the most gaps in measurement, ensuring better data quality awareness before analysis. ``` ```{r} # Task 3 – Function to plot variation in pollutant levels # measured at noon on the 1st day of each month between two dates library(ggplot2) library(lubridate) displayPollutantGraph <- function(pollutant, start_month, start_year, end_month, end_year) { df <- PortTalbotAQ # Fix invalid "24:00:00" times -> "00:00:00" (next day) df$Time[df$Time == "24:00:00"] <- "00:00:00" # Combine date and time df$DateTime <- suppressWarnings(dmy_hms(paste(df$Date, df$Time))) if (all(is.na(df$DateTime))) { df$DateTime <- suppressWarnings(ymd_hms(paste(df$Date, df$Time))) } # Extract components df$Year <- year(df$DateTime) df$Month <- month(df$DateTime) df$Day <- day(df$DateTime) df$Hour <- hour(df$DateTime) # Convert pollutant to numeric df[[pollutant]] <- as.numeric(df[[pollutant]]) # Filter: noon (12:00:00), first day of each month, within range subsetData <- subset( df, Day == 1 & Hour == 12 & ((Year > start_year | (Year == start_year & Month >= start_month)) & (Year < end_year | (Year == end_year & Month <= end_month))) ) if (nrow(subsetData) == 0) { warning(paste("No matching data found for", pollutant, "at noon on 1st of month")) return(NULL) } # Plot ggplot(subsetData, aes(x = DateTime, y = get(pollutant))) + geom_line(color = "steelblue", linewidth = 1) + geom_point(color = "darkred", size = 2) + labs( title = paste("Variation in", pollutant, "levels at noon on 1st of each month"), x = "Date", y = paste(pollutant, "(µg/m³)") ) + theme_minimal(base_size = 12) } # Required example: displayPollutantGraph("NO2", 7, 2015, 11, 2024) #I created a function to visualise how pollutant levels change over time. The function extracts data recorded at noon on the 1st of each month within a chosen date range. I used lubridate to handle date-time conversions and ggplot2 for plotting because it allows for clean and readable visualisations. The graph helps identify monthly trends in pollutant concentration, making it easier to observe long-term changes or seasonal patterns. ``` ```{r} # Task 4 – Compare multiple pollutants on one graph # Variation in pollutant levels measured at noon on the 1st of each month # over the full 10-year period library(ggplot2) library(lubridate) library(dplyr) library(tidyr) comparePollutants <- function(pollutantVector) { df <- PortTalbotAQ # --- CLEAN COLUMN NAMES --- # Replace blanks or NA names with generic "X1", "X2", etc. colnames(df) <- make.names(colnames(df), unique = TRUE) # Keep only pollutant columns that actually exist valid_pollutants <- pollutantVector[pollutantVector %in% colnames(df)] if (length(valid_pollutants) == 0) { stop("None of the pollutants in pollutantVector were found in the dataset.") } # Fix invalid "24:00:00" times (convert to 00:00:00) df$Time[df$Time == "24:00:00"] <- "00:00:00" # Combine and parse Date + Time safely df$DateTime <- suppressWarnings(dmy_hms(paste(df$Date, df$Time))) if (all(is.na(df$DateTime))) { df$DateTime <- suppressWarnings(ymd_hms(paste(df$Date, df$Time))) } # Extract date parts df <- df %>% mutate( Year = year(DateTime), Month = month(DateTime), Day = day(DateTime), Hour = hour(DateTime) ) # Convert valid pollutants to numeric df[valid_pollutants] <- lapply(df[valid_pollutants], as.numeric) # Filter: 1st day of month, 12:00 (noon) subsetData <- df %>% filter(Day == 1, Hour == 12) # If no rows found, warn and exit if (nrow(subsetData) == 0) { warning("No matching data found for noon on 1st of month.") return(NULL) } # Create clean DateLabel column subsetData <- subsetData %>% mutate(DateLabel = as.Date(DateTime)) %>% select(DateLabel, all_of(valid_pollutants)) # Reshape to long format longData <- subsetData %>% pivot_longer( cols = all_of(valid_pollutants), names_to = "Pollutant", values_to = "Value" ) # --- PLOT --- ggplot(longData, aes(x = DateLabel, y = Value, color = Pollutant)) + geom_line(linewidth = 1) + geom_point(size = 1.5) + labs( title = "Comparison of Pollutant Levels at Noon on the 1st of Each Month (2015–2024)", x = "Date", y = "Concentration (µg/m³)", color = "Pollutant" ) + theme_minimal(base_size = 12) } comparePollutants(c("NO2", "O3", "SO2", "PM10")) #I built a function to compare multiple pollutants on the same graph. I used dplyr for filtering and tidyr to reshape the data into a long format suitable for plotting. The function selects values recorded at noon on the 1st of each month and plots each pollutant as a separate line using ggplot2. This makes it easy to visually compare how different pollutants vary over time and observe any correlations or seasonal trends between them. ``` ```{r} # Task 5 – Correlation between O3 and NO2 + scatterplot (final version) library(ggplot2) library(dplyr) library(rlang) O3NOCorrelation <- function() { df <- PortTalbotAQ colnames(df) <- make.names(colnames(df), unique = TRUE) # Identify columns nm <- tolower(colnames(df)) o3_col <- colnames(df)[which(nm %in% c("o3", "ozone"))[1]] no2_col <- colnames(df)[which(nm %in% c("no2", "nitrogendioxide", "nitrogen.dioxide"))[1]] if (is.na(o3_col) || is.na(no2_col)) stop("Could not find O3 and/or NO2 columns.") # Clean and convert to numeric for (v in c(o3_col, no2_col)) { df[[v]][df[[v]] %in% c("No data", "")] <- NA df[[v]] <- as.numeric(df[[v]]) } clean <- df %>% select(all_of(c(o3_col, no2_col))) %>% filter(!is.na(.data[[o3_col]]), !is.na(.data[[no2_col]])) if (nrow(clean) == 0) stop("No overlapping O3/NO2 data after cleaning.") # Pearson correlation ct <- cor.test(clean[[o3_col]], clean[[no2_col]], method = "pearson") message(sprintf( "Pearson correlation (O3 vs NO2): r = %.3f, 95%% CI [%.3f, %.3f], p = %.3g, n = %d", unname(ct$estimate), ct$conf.int[1], ct$conf.int[2], ct$p.value, nrow(clean) )) # Build the plot p <- ggplot(clean, aes(x = !!sym(no2_col), y = !!sym(o3_col))) + geom_point(alpha = 0.6, na.rm = TRUE, color = "green") + geom_smooth(method = "lm", se = TRUE, color = "red") + labs( title = "Ozone (O3) vs Nitrogen Dioxide (NO2) — Entire Period", x = "NO2 (µg/m³)", y = "O3 (µg/m³)" ) + theme_minimal(base_size = 12) # Print the plot so it appears in the output print(p) invisible(ct) } O3NOCorrelation() #I explored the relationship between ozone (O₃) and nitrogen dioxide (NO₂) using a Pearson correlation test. I cleaned the data by converting “No data” values to NA and kept only valid numeric records. Then, I used cor.test() to calculate the correlation coefficient and plotted the results using ggplot2 with a regression line. This scatterplot visually shows whether O₃ and NO₂ have a positive or negative relationship, helping to understand how changes in one pollutant relate to the other. ``` ```{r} # Task 6 – Independent t-test comparing ozone levels for low vs high wind speeds # H0: Mean O3 level is the same when wind speed <5 m/s and ≥5 m/s # H1: Mean O3 levels differ between these two wind-speed categories library(dplyr) O3WindSpeed_ttest <- function() { df <- PortTalbotAQ colnames(df) <- make.names(colnames(df), unique = TRUE) # Identify relevant columns (tolerant to naming) nm <- tolower(colnames(df)) o3_col <- colnames(df)[which(nm %in% c("o3", "ozone"))[1]] ws_col <- colnames(df)[which(nm %in% c("windspeed", "wind.speed", "ws"))[1]] if (is.na(o3_col) || is.na(ws_col)) stop("Could not find O3 and/or WindSpeed columns.") # Clean "No data" and convert to numeric for (v in c(o3_col, ws_col)) { df[[v]][df[[v]] %in% c("No data", "")] <- NA df[[v]] <- as.numeric(df[[v]]) } # Remove missing values df_clean <- df %>% filter(!is.na(.data[[o3_col]]), !is.na(.data[[ws_col]])) # Create wind speed groups df_clean <- df_clean %>% mutate(WindGroup = ifelse(.data[[ws_col]] < 5, "Low (<5 m/s)", "High (≥5 m/s)")) # Run independent samples t-test t_res <- t.test( formula = as.formula(paste(o3_col, "~ WindGroup")), data = df_clean, var.equal = FALSE, # Welch's t-test (unequal variances) alternative = "two.sided" ) # Print results neatly message(sprintf( "t = %.3f, df = %.1f, p = %.4g; Mean(O3)<5 = %.2f, Mean(O3)≥5 = %.2f", t_res$statistic, t_res$parameter, t_res$p.value, mean(df_clean[[o3_col]][df_clean$WindGroup == 'Low (<5 m/s)'], na.rm = TRUE), mean(df_clean[[o3_col]][df_clean$WindGroup == 'High (≥5 m/s)'], na.rm = TRUE) )) # Return full test output for inspection return(t_res) } # Run Task 6 O3WindSpeed_ttest() #I used a t-test to compare ozone (O₃) levels under low and high wind speeds. I grouped the data into two categories: wind speeds below 5 m/s and those equal to or above 5 m/s. After cleaning and converting the data, I applied Welch’s t-test (t.test() with var.equal = FALSE) because it does not assume equal variances. This test checks whether there is a statistically significant difference in mean ozone levels between the two wind-speed groups. ``` ```{r} # Task 7 – ANOVA on NO2 at 17:00 for Nov–Dec in 2016, 2018, 2020 # If significant, perform Tukey HSD post-hoc tests. library(dplyr) library(lubridate) library(ggplot2) NO2_ANOVA_17h <- function() { df <- PortTalbotAQ # Make names valid in case of blanks/dupes colnames(df) <- make.names(colnames(df), unique = TRUE) # Find the NO2 column (tolerant to naming) nm <- tolower(colnames(df)) no2_col <- colnames(df)[which(nm %in% c("no2", "nitrogendioxide", "nitrogen.dioxide"))[1]] if (is.na(no2_col)) stop("Could not find a NO2 column.") # Clean NO2: convert "No data"/"" -> NA, then numeric df[[no2_col]][df[[no2_col]] %in% c("No data", "")] <- NA df[[no2_col]] <- as.numeric(df[[no2_col]]) # Fix invalid "24:00:00" times before parsing if ("Time" %in% colnames(df)) { df$Time[df$Time == "24:00:00"] <- "00:00:00" } # Parse DateTime robustly (assumes Date + Time columns exist) # Your file uses dd/mm/YYYY HH:MM:SS; we try a couple formats just in case df$DateTime <- suppressWarnings(dmy_hms(paste(df$Date, df$Time))) if (all(is.na(df$DateTime))) { df$DateTime <- suppressWarnings(ymd_hms(paste(df$Date, df$Time))) } # Derive components df <- df %>% mutate( Year = year(DateTime), Month = month(DateTime), Day = day(DateTime), Hour = hour(DateTime) ) # Filter: 17:00 on each day, Nov–Dec, for 2016, 2018, 2020 target_years <- c(2016, 2018, 2020) sub <- df %>% filter( Year %in% target_years, Month %in% c(11, 12), Hour == 17, !is.na(.data[[no2_col]]) ) %>% mutate( Period = factor(paste0(Year, "_NovDec"), levels = paste0(target_years, "_NovDec")) ) if (nrow(sub) == 0) stop("No NO2 data found at 17:00 for Nov–Dec in 2016/2018/2020.") # Quick descriptive stats per period (useful in your write-up) desc <- sub %>% group_by(Period) %>% summarise( n = n(), mean_NO2 = mean(.data[[no2_col]], na.rm = TRUE), sd_NO2 = sd(.data[[no2_col]], na.rm = TRUE), .groups = "drop" ) print(desc) # One-way ANOVA fit <- aov(formula = as.formula(paste(no2_col, "~ Period")), data = sub) cat("\nANOVA summary:\n") print(summary(fit)) # If significant, do Tukey HSD post-hoc to locate pairwise differences p_val <- summary(fit)[[1]][["Pr(>F)"]][1] if (!is.na(p_val) && p_val < 0.05) { cat("\nPost-hoc (Tukey HSD) because ANOVA is significant:\n") tuk <- TukeyHSD(fit, "Period") print(tuk) } else { cat("\nANOVA not significant at alpha = 0.05; no post-hoc performed.\n") } # Optional: visualise distributions (helps marker) p <- ggplot(sub, aes(x = Period, y = .data[[no2_col]], fill = Period)) + geom_boxplot(outlier.alpha = 0.5) + labs( title = "NO2 at 17:00 — Nov–Dec (2016 vs 2018 vs 2020)", x = "Period", y = "NO2 (µg/m³)" ) + theme_minimal(base_size = 12) + theme(legend.position = "none") print(p) invisible(list(desc = desc, anova = fit)) } # Run Task 7 NO2_ANOVA_17h() #I performed an ANOVA test to check if nitrogen dioxide (NO₂) levels at 17:00 during November and December differed across the years 2016, 2018, and 2020. I used aov() for the one-way ANOVA because it compares mean values between multiple groups. If the ANOVA was significant, I applied the Tukey HSD post-hoc test to see which years showed the actual differences. I also included a boxplot to visualise the variation between years. This approach helps identify whether air quality improved or worsened over time. ```