# ##############################################################################
# Task 1 Question 4 - Creating a function to compare multiple pollutants on the
# same graph (noon on the 1st of each month across 2015-2024)
# ##############################################################################
comparePollutants <- function(pollutants){
PortTalbotAQ$date <- as.Date(PortTalbotAQ$date, format ="%Y-%m-%d") # Ensuring the column is date format.
start_date <- as.Date("2015-01-01") # Specifying start date
end_date <- as.Date("2024-12-31") # specifying end date
keep_rows <- PortTalbotAQ$date >= start_date & # This section is setting the parameters of the requested conditions
PortTalbotAQ$date <= end_date &
format(PortTalbotAQ$date, "%d") == "01" &
PortTalbotAQ$time == "12:00:00"
df <- PortTalbotAQ[keep_rows, c("date", pollutants)] # Ensure the function only keeps date & requested pollutant columns
for (p in pollutants) {
df[[p]][df[[p]] == "No data"] <- NA
df[[p]] <- suppressWarnings(as.numeric(df[[p]]))
} # This section is converting "No data" to NA and each pollutant column to numeric
df <- df[rowSums(is.na(df[pollutants])) < length(pollutants), ] # Removing rows where all selected pollutants are missing
if (nrow(df) == 0) {
stop("No valid data found for the selected pollutants at noon on the 1st of the month (2015-2024).")
}
y_min <- min(sapply(pollutants, function(p) df[[p]]), na.rm =TRUE)
y_max <- max(sapply(pollutants, function(p) df[[p]]), na.rm =TRUE)
# The above determines the parameters for the selected pollutants
cols <- rep(c("blue", "red", "darkgreen", "purple"), length.out = length(pollutants)) # Adding colours to the plot lines to distinguish pollutants apart
plot( # Arranging plot similar to Task 1 Question 3.
df$date,
df[[pollutants[1]]],
type = "o",
pch = 16,
col = cols[1], # Adding colours in the plot
xlab ="Date",
ylab = "Pollutant level",
main = "Comparison of pollutant levels (Noon on 1st of month, 2015-2024)",
ylim = c(y_min, y_max) # This is being used as plotting many unlike Task 1 question 3 where we only plotted one.
)
if (length(pollutants) > 1) { # Adding the remaining pollutants as lines/points on the same graph
for (i in 2:length(pollutants)) {
lines(df$date,
df[[pollutants[i]]],
type = "o",
pch =16,
col = cols[i] # Adding colours into the legend
)
}
}
legend("topright",
legend = pollutants,
col = cols,
lty = 1,
pch = 16,
bty = "n"
) # lty = line type (1 = solid line), bty = box type (n = no box)
}
# Add 4 pollutants together
pollutants4 <- c("ozone", "nitrogen_dioxide", "PM10", "PM2_5")
comparePollutants(pollutants4) # Checking to see if worked
# ##############################################################################
# Task 1 Question 5 - Creating a function to calculate correlation between ozone
# and nitrogen dioxide and plot the data as a scatterplot.
# ##############################################################################
O3NOCorrelation <- function(p1,p2) {    # Creating the function
df <- PortTalbotAQ[, c(p1, p2)] # Creating a smaller data frame with only the two columns needed
df$ozone[df$ozone == "No data"] <- NA
df$nitrogen_dioxide[df$nitrogen_dioxide == "No data"] <- NA # Replace "No data" with NA for both variables
# Time to convert both columns to numeric (suppress warning caused by NA coercion)
df$ozone <- suppressWarnings(as.numeric(df$ozone))
df$nitrogen_dioxide <- suppressWarnings(as.numeric(df$nitrogen_dioxide))
df <- df[!is.na(df$ozone) & !is.na(df$nitrogen_dioxide), ] # Keeps only rows where both variables required are valid numbers
if (nrow(df) == 0) {
stop("No valid paired ozone and nitorgen_dioxide measurements available to calculate correlation.") # Safety check in cvase nothing remains
}
corr_value <- cor(df$ozone, df$nitrogen_dioxide) # Calculate Pearson correlation
plot(                    # Plot the paired data as a scatterplot
df$ozone,
df$nitrogen_dioxide,
pch = 16,
col = rgb(0, 0, 0, 0.2), # Helps show density patterns
xlab = "Ozone",
ylab = "Nitrogen dioxide",
main = paste("Ozone vs Nitrogen dioxide (Correlation =", round(corr_value, 3), ")") # rounded to 3 decimal places and rounded
)
return(round(corr_value, 3))
}
O3NOCorrelation("ozone", "nitrogen_dioxide") # Returns output as required
# rm(O3NOCorrelation)
# ##############################################################################
# Task 1 Question 6 - T-test: Ozone levels when wind speed < 5 m/s vs >= 5 m/s
# ##############################################################################
# Experimental hypotheses:
# H0 (null): The mean ozone level is the same when wind speed < 5 m/s and when wind speed >= 5 m/s.
# H1 (alternative): The mean ozone level is different between the two wind speed groups.
df <- PortTalbotAQ[, c("ozone", "wind_speed")] # Creating a smaller data frame with only required columns needed
df$ozone[df$ozone == "No data"] <- NA
df$wind_speed[df$wind_speed == "No data"] <- NA # Similar to previous code replacing No Data with NA for both data types
df$ozone <- suppressWarnings(as.numeric(df$ozone))
df$wind_speed <- suppressWarnings(as.numeric(df$wind_speed)) #Converting both columns to numeric (suppress warning from coercion)
df <- df[!is.na(df$ozone) & !is.na(df$wind_speed), ] # Keeps only rows where both ozone and wind_speed are valid
# Splitting ozone into two groups to assist further calculations (based on wind speed)
ozone_low_wind <- df$ozone[df$wind_speed < 5] # Wind speed below 5 m/s
ozone_high_wind <- df$ozone[df$wind_speed >= 5] # wind speed 5 m/s or higher
if (length(ozone_low_wind) == 0 | length(ozone_high_wind) == 0) {
stop("Not enough valid data in one or both wind speed groups to perform the t-test.") # This ensures both groups hold data and if no data is present this comment will be outputted.
}
ttest_result <- t.test(ozone_low_wind, ozone_high_wind, alternative = "two.sided") # Two-sample t-test (Welch)
ttest_result # Provides output of above Welch T-test
# The above worked and provided an output. Time to add the decision rule and conclusion
# If p-value < 0.05, reject H0 (significant difference in mean ozone between groups)
# If p-value >= 0.05, fail to reject H0 (no significant difference detected)
alpha <- 0.05
if (ttest_result$p.value < alpha) {
message ("Conclusion: p-value = ", round(ttest_result$p.value, 4),
" < ", alpha, ". Reject H0. There is a statistically significant difference in mean ozone levels between wind speed groups.") # Reject H0
} else {
message ("Conclusion: p-value = ", round(ttest_result$p.value, 4),
" >= ", alpha, ". Fail to reject H0. There is no statistically significant difference in mean ozone levels between wind speed groups.") # Fail to reject H0
}
# When the above has been run it has outputted the Reject H0 option which shows the code is working.
# ##############################################################################
# Task 1 Question 7 - ANOVA: Nitrogen dioxide levels across selected periods
# (Nov-Dec 2016, 2018, 2020 at 17:00)
# ##############################################################################
# Experimental Hypotheses again similar to Task 1 Question 6:
# H0 (null): Mean nitrogen dioxide levels are equal across all three periods required
# H1 (alternative): At least one period has a different mean nitrogen dioxide level
# As the question is requiring specific dates will intially ensure the date column is date format.
PortTalbotAQ$date <- as.Date(PortTalbotAQ$date, format = "%Y-%m-%d")
df <- PortTalbotAQ[, c("date", "time", "nitrogen_dioxide")] # Picking only relevant columns
df$nitrogen_dioxide[df$nitrogen_dioxide == "No data"] <- NA # Replacing "No data with NA like before
df$nitrogen_dioxide <- suppressWarnings(as.numeric(df$nitrogen_dioxide)) # Converting to numeric
df <- df[!is.na(df$nitrogen_dioxide), ] # Removing rows with missing values
df <- df[df$time == "17:00:00", ] # Setting the exact time in which the data need to be pulled
df$year <- as.numeric(format(df$date, "%Y")) # Creating year & month variables
df$month <- as.numeric(format(df$date, "%m"))
df <- df[df$month %in% c(11, 12), ] # Keeping only Nov & December data
df <- df[df$year %in% c(2016, 2018, 2020), ] # Keeping only the three years of interest
df$period <- factor(df$year,
levels = c(2016, 2018, 2020),
labels = c("Nov-Dec 2016", "Nov-Dec 2018", "Nov-Dec 2020")) # Creating a factor representing the three comparison periods
if (length(unique(df$period)) < 3) {
stop("Not enough data available for all three periods to perform ANOVA.") # Automated response if not enough available data to carry out request.
}
# Main code has now been added along with a safety check time to add the ANOVA code
anova_result <- aov(nitrogen_dioxide ~ period, data = df) # Run one-way ANOVA
summary(anova_result) # Displaying result
# Post-Hoc tests if significant
Alpha <- 0.05
anova_p <- summary(anova_result)[[1]]["period", "Pr(>F)"]
if (anova_p < alpha) {
message("ANOVA result: p-value = ", round(anova_p, 4),
" < ", alpha, ". Reject H0. At least one period differs significantly.")
# Post-Hoc explanation: ANOVA tells us a difference exists overall, but not which groups differ
# while Tukey HSD compares all pairs while controlling for multiple comparisons.
TukeyHSD(anova_result)
temp<-TukeyHSD(anova_result)
} else {
message("ANOVA result: p-value = ", round(anova_p, 4),
" >= ", alpha, ". Fail to reject H0. No signficant differences detected.")
}
# ANOVA results indicated a statistically significant difference in mean nitrogen dioxide
# levels across the three periods (p < 0.05), therefore the null hypothesis was rejected.
#
# A Tukey HSD post-hoc test was conducted to identify where the differences lay.
# The results showed that nitrogen dioxide levels in Nov-Dec 2016 were significantly
# higher than in both Nov-Dec 2018 and Nov-Dec 2020.
#
# No statistically significant difference was found between the Nov-Dec 2018 and
# Nov-Dec 2020 periods.
# ##############################################################################
install.packages("readxl")
View(temp)
#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)
#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", 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.
View(PortTalbotAQ)
View(invalidData)
View(displayPollutantGraph)
View(displayPollutantGraph)
View(comparePollutants)
View(O3NOCorrelation)
View(O3WindTTest)