library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
# initializing pollution data
pollutionData <- read.csv("/home/nlc62/csds133/nlc62/24sp-csds133-nlc62/5-FinalProject/pollution_2000_2023.csv")
pollutionData$State <- factor(pollutionData$State)
pollutionData$County <- factor(pollutionData$County)
#pollutionData$City <- factor(pollutionData$City)
pollutionData$Date <- as.POSIXct(pollutionData$Date, format = "%Y-%m-%d")
pollutionData$Year <- as.numeric(format(pollutionData$Date, "%Y"))
pollutionData <- pollutionData %>% filter(Year >= 2016) %>% filter(Year < 2023)
# initializing weather data
weatherData <- read.csv("/home/nlc62/csds133/nlc62/24sp-csds133-nlc62/5-FinalProject/WeatherEvents_Jan2016-Dec2022.csv")
weatherData$Severity <- factor(weatherData$Severity)
weatherData$Type <- factor(weatherData$Type)
weatherData$AirportCode <- factor(weatherData$AirportCode)
weatherData$TimeZone <- factor(weatherData$TimeZone)
weatherData$City <- factor(weatherData$City)
weatherData$County <- factor(weatherData$County)
weatherData$State <- factor(weatherData$State)
weatherData$ZipCode <- factor(weatherData$ZipCode)
# exploring the data
#levels(weatherData$Severity) # severity levels are "Light", "Moderate", "Heavy", "Severe", "Other", and "UNK" (unknown)
#levels(weatherData$Type) # types of weather are "Cold", "Fog", "Hail", "Precipitation", "Rain", "Snow", and "Storm"
#weatherData %>% filter(Severity == "Other") %>% filter(Type == "Cold" |Type == "Fog" |Type == "Precipitation" | Type =="Rain" | Type =="Snow" | Type =="Storm")
# "Other" severity is only used for hail
#weatherData %>% filter(Severity=="UNK")
#weatherData %>% filter(Severity == "UNK") %>% filter(Type == "Cold" |Type == "Fog" |Type == "Hail" | Type =="Rain" | Type =="Snow" | Type =="Storm")
# "UNK" severity is only used for precipitation
# check what states are included in each data set
#levels(weatherData$State)
#levels(pollutionData$State)
# filtering by state but not necessary since filtering by cities
#weatherData <- weatherData %>% filter(State!="MT") %>% filter(State!="NE") %>% filter(State!="WV")
#pollutionData <- pollutionData %>% filter(State!="Alaska") %>% filter(State!="District of Columbia") %>% filter(State!="Hawaii")
# filtering so that only cities that exist between both data sets are used
pollutionData$StateAbb <- state.abb[match(pollutionData$State, state.name)]
pollutionData$CityAndState <- paste(pollutionData$City, pollutionData$StateAbb)
weatherData$CityAndState <- paste(weatherData$City, weatherData$State)
pollutionCities <- unique(pollutionData$CityAndState)
weatherCities <- unique(weatherData$CityAndState)
overlapCitiesWeatherData <- weatherData %>% filter(weatherData$CityAndState %in% pollutionCities)
overlapCitiesPollutionData <- pollutionData %>% filter(pollutionData$CityAndState %in% weatherCities)
# ordering by city
overlapCitiesPollutionData <- overlapCitiesPollutionData[order(overlapCitiesPollutionData$City),]
overlapCitiesWeatherData <- overlapCitiesWeatherData[order(overlapCitiesWeatherData$City),]
overlapCitiesPollutionData$City <- factor(overlapCitiesPollutionData$City)
# ensure no repeats (don't want to have certain datapoints overrepresented)
overlapCitiesPollutionData <- distinct(overlapCitiesPollutionData)
# checking if cities were filtered correctly
uniqueCities1 <- unique(overlapCitiesPollutionData$City)
uniqueCities2 <- unique(overlapCitiesWeatherData$City)
print(uniqueCities1 %in% uniqueCities2)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE
# adding columns to separate the start/end dates and times, making data type date
overlapCitiesWeatherData <- overlapCitiesWeatherData %>% separate(StartTime.UTC.,c("startDate", "startTime"), sep = " ", remove = FALSE)
overlapCitiesWeatherData <- overlapCitiesWeatherData %>% separate(EndTime.UTC.,c("endDate", "endTime"), sep = " ", remove = FALSE)
overlapCitiesWeatherData$startDate <- as.POSIXct(overlapCitiesWeatherData$startDate, format="%Y-%m-%d")
overlapCitiesWeatherData$endDate <- as.POSIXct(overlapCitiesWeatherData$endDate, format="%Y-%m-%d")
# calculating amount of time event lasted
overlapCitiesWeatherData$timeEventLasted <- difftime(overlapCitiesWeatherData$EndTime.UTC., overlapCitiesWeatherData$StartTime.UTC., units = "mins")
LosAngelesPollutionData <- overlapCitiesPollutionData %>% filter(City=='Los Angeles') %>% filter(State=="California")
LosAngelesWeatherData <- overlapCitiesWeatherData %>% filter(City == 'Los Angeles') %>% filter(State=="CA")
LosAngelesWeatherData$startDate <- as.POSIXct(LosAngelesWeatherData$startDate, format="%Y-%m-%d")
LosAngelesWeatherData$endDate <- as.POSIXct(LosAngelesWeatherData$endDate, format="%Y-%m-%d")
OriginalLAWeatherData <-overlapCitiesWeatherData %>% filter(City == 'Los Angeles')
# Reset row names
row.names(LosAngelesWeatherData) <- NULL
# loop to separate weather events that span two days, correct time of event
k=nrow(LosAngelesWeatherData)
LosAngelesWeatherData$StartTime.UTC. <- as.POSIXct(LosAngelesWeatherData$StartTime.UTC., format = "%Y-%m-%d %H:%M:%S")
# HERE I need help figuring out how to make this work. I'm trying to make it such that if a weather event spans between two days (ex: 23:40 first day to 03:27 the second day, I want to split these up into their own weather events so that we can look at each day specifically, and change the time of the event to reflect this.) I think it also would be good for events with precipitation to try to split the amount of precipitation by the same proportion as the time is split between the two days.
for(i in 1:k) {
if(as.Date(LosAngelesWeatherData$EndTime.UTC.[i]) > as.Date(LosAngelesWeatherData$StartTime.UTC.[i])) {
LosAngelesWeatherData <- rbind(LosAngelesWeatherData, LosAngelesWeatherData[i,])
row.names(LosAngelesWeatherData) <- NULL
date_only <- as.Date(LosAngelesWeatherData$StartTime.UTC.[i])
next_day <- date_only + 1
midnight_next_day <- as.POSIXct(paste(next_day, "00:00:00"), format = "%Y-%m-%d %H:%M:%S")
time_diff <- difftime(midnight_next_day, LosAngelesWeatherData$StartTime.UTC.[i], units = "mins")
LosAngelesWeatherData$timeEventLasted[nrow(LosAngelesWeatherData)] <- LosAngelesWeatherData$timeEventLasted[i]-time_diff
LosAngelesWeatherData$timeEventLasted[i] <- time_diff
LosAngelesWeatherData$StartTime.UTC.[nrow(LosAngelesWeatherData)] <- as.POSIXct(paste(LosAngelesWeatherData$endDate[i], "00:00:01"))#, format = "%Y-%m-%d %H:%M:%S")
LosAngelesWeatherData$startDate[nrow(LosAngelesWeatherData)] <- as.POSIXct(LosAngelesWeatherData$endDate[i], format = "%Y-%m-%d")
LosAngelesWeatherData$startTime[nrow(LosAngelesWeatherData)] <- "00:00:00"
LosAngelesWeatherData$EndTime.UTC.[i] <- (paste((LosAngelesWeatherData$startDate[i]),"23:59:59"))#, format = "%Y-%m-%d %H:%M:%S")
LosAngelesWeatherData$endDate[i] <- as.POSIXct(LosAngelesWeatherData$startDate[i], format = "%Y-%m-%d")
LosAngelesWeatherData$endTime[i] <- "23:59:59"
}
}
# planning how to calculate times
# Extract the date portion from the datetime column
#date_only <- as.Date(LosAngelesWeatherData$StartTime.UTC.[27])
# Add one day to the date portion to get the date of the next day
#next_day <- date_only + 1
# Combine the next day's date with midnight time to get midnight of the next day
#midnight_next_day <- as.POSIXct(paste(next_day, "00:00:00"), format = "%Y-%m-%d %H:%M:%S")
# Calculate the time difference from midnight of the next day for each row
#time_diff <- difftime(midnight_next_day, LosAngelesWeatherData$StartTime.UTC.[27], units = "mins")
#print(time_diff)
#print(midnight_next_day)
#print(LosAngelesWeatherData$StartTime.UTC.[27])
#print(difftime(next_day,LosAngelesWeatherData$StartTime.UTC.[27], units = "mins"))
#print(difftime(as.POSIXct(paste(as.Date(LosAngelesWeatherData$StartTime.UTC.[27]+1), "00:00:00")),LosAngelesWeatherData$StartTime.UTC.[27], units = "mins"))
#print(difftime(as.POSIXct(LosAngelesWeatherData$endTime[5855], format = "%H:%M:%D"), as.POSIXct(LosAngelesWeatherData$startTime[5855],format = "%H:%M:%D"), units = "hours"))
#print(difftime(LosAngelesWeatherData$endTime[5855], LosAngelesWeatherData$startTime[5855], units = "minutes"))
#LosAngelesWeatherData$timeEventLasted <- difftime(LosAngelesWeatherData$EndTime.UTC., LosAngelesWeatherData$StartTime.UTC., units = "mins")
# graph for Los Angeles
# loop to put true/false if weather event occurred on a given day in LA
LosAngelesWeatherData <- LosAngelesWeatherData[order(LosAngelesWeatherData$startDate),]
row.names(LosAngelesWeatherData) <- NULL
LosAngelesPollutionData$weatherEvents <- NA
LosAngelesPollutionData <- LosAngelesPollutionData[!duplicated(LosAngelesPollutionData$Date), ]
c
## function (...) .Primitive("c")
j <- 1
i <- 1
for(i in 1:nrow(LosAngelesPollutionData)) {
while((as.numeric(difftime(LosAngelesPollutionData$Date[i],LosAngelesWeatherData$startDate[j], units="days"))) >0 && j <= nrow(LosAngelesWeatherData)) {
j <- j+1
}
if(as.Date(LosAngelesPollutionData$Date[i])==as.Date(LosAngelesWeatherData$startDate[j]) && j <= nrow(LosAngelesWeatherData)) {
LosAngelesPollutionData$weatherEvents[i] <- TRUE
}
else {
LosAngelesPollutionData$weatherEvents[i] <- FALSE
}
}
# making separate data sets for days with vs without weather events
LosAngelesDataNoWeatherEvents <- LosAngelesPollutionData %>% filter(weatherEvents==FALSE)
LosAngelesDataNoWeatherEvents <- LosAngelesDataNoWeatherEvents[!duplicated(LosAngelesDataNoWeatherEvents$Date), ]
LosAngelesDataWeatherEvents <- LosAngelesPollutionData %>% filter(weatherEvents==TRUE)
LosAngelesDataWeatherEvents <- LosAngelesDataWeatherEvents[!duplicated(LosAngelesDataWeatherEvents$Date), ]
#LosAngelesDataWeatherEvents$Date <- as.Date(LosAngelesDataWeatherEvents$Date, format="%Y-%m-%d")
#LosAngelesWeatherData$startDate <- as.Date(LosAngelesWeatherData$startDate, format="%Y-%m-%d")
# violin plot and histogram: no weather events, SO2 AQI
p <- ggplot(LosAngelesDataNoWeatherEvents, aes(x=City, y=SO2.AQI)) +
geom_violin()
p

hist(LosAngelesDataNoWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in LA on Days with No Recorded Weather")

# violin plot and histogram: no weather events, O3 AQI
p <- ggplot(LosAngelesDataNoWeatherEvents, aes(x=City, y=O3.AQI)) +
geom_violin()
p

hist(LosAngelesDataNoWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in LA on Days with No Recorded Weather")

# violin plot and histogram: recorded weather events, SO2 AQI
p <- ggplot(LosAngelesDataWeatherEvents, aes(x=City, y=SO2.AQI)) +
geom_violin()
p

hist(LosAngelesDataWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in LA on Days with Recorded Weather")

# violin plot and histogram: recorded weather events, O3 AQI
p <- ggplot(LosAngelesDataWeatherEvents, aes(x=City, y=O3.AQI)) +
geom_violin()
p

hist(LosAngelesDataWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in LA on Days with Recorded Weather")

# cleaning data for the loop (order, duplicates for pollution)
row.names(LosAngelesDataWeatherEvents) <- NULL
LosAngelesDataWeatherEvents <- LosAngelesDataWeatherEvents[!duplicated(LosAngelesDataWeatherEvents$Date), ]
LosAngelesDataWeatherEvents <- LosAngelesDataWeatherEvents[order(LosAngelesDataWeatherEvents$Date),]
LosAngelesWeatherData <- LosAngelesWeatherData[order(LosAngelesWeatherData$startDate),]
# calculate total time of weather events for the day
LosAngelesDataWeatherEvents$timeLasted <- NA
j <- 1
i <- 1
for(i in 1:nrow(LosAngelesDataWeatherEvents)) {
p <- 0
while((as.numeric(difftime(as.Date(LosAngelesDataWeatherEvents$Date[i]),as.Date(LosAngelesWeatherData$startDate[j]), units="days")))>0 && j <= nrow(LosAngelesWeatherData)) {
j <- j+1
}
while((as.numeric(difftime(as.Date(LosAngelesDataWeatherEvents$Date[i]),as.Date(LosAngelesWeatherData$startDate[j]), units="days")))==0 && j <= nrow(LosAngelesWeatherData)) {
p <- p + LosAngelesWeatherData$timeEventLasted[j]
j <- j+1
}
LosAngelesDataWeatherEvents$timeLasted[i] <- p
p <- 0
}
# plotting boxplots for O3 AQI for all cities
boxplot(O3.AQI~City,
data=overlapCitiesPollutionData,
main="O3 AQI for Different Cities",
xlab="City",
ylab="O3 AQI",
col="pink",
border="purple", las = 2
)

# boxplots for SO2 AQI for all cities shared between weather and pollution datasets
boxplot(SO2.AQI~City,
data=overlapCitiesPollutionData,
main="SO2 AQI for Different Cities",
xlab="City",
ylab="SO2 AQI",
col="pink",
border="purple", las = 2
)

# boxplots for NO2 AQI for all cities shared between weather and pollution datasets
boxplot(NO2.AQI~City,
data=overlapCitiesPollutionData,
main="NO2 AQI for Different Cities",
xlab="City",
ylab="NO2 AQI",
col="pink",
border="purple", las = 2
)

# boxplots for CO AQI for all cities shared between weather and pollution datasets
boxplot(CO.AQI~City,
data=overlapCitiesPollutionData,
main="CO AQI for Different Cities",
xlab="City",
ylab="CO AQI",
col="pink",
border="purple", las = 2
)

ggplot(overlapCitiesPollutionData, aes(x = O3.AQI, y = SO2.AQI, color=City)) +
geom_point() +
labs(x = "O3 AQI", y = "SO2 AQI",
title = "Relationship between SO2 AQI and O3 AQI",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0))

ggplot(overlapCitiesPollutionData, aes(x = O3.AQI, y = CO.AQI, color=City)) +
geom_point() +
labs(x = "O3 AQI", y = "CO AQI",
title = "Relationship between CO AQI and O3 AQI",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0))

ggplot(overlapCitiesPollutionData, aes(x = NO2.AQI, y = CO.AQI, color=City)) +
geom_point() +
labs(x = "NO2 AQI", y = "CO AQI",
title = "Relationship between NO2 AQI and O3 AQI",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0))

# Plot heatmap
library(reshape2)
##
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
##
## smiths
pollutionMeasures <- overlapCitiesPollutionData[, c("O3.AQI", "SO2.AQI", "NO2.AQI", "CO.AQI")]
correlation_matrix <- cor(pollutionMeasures)
correlation_long <- melt(correlation_matrix)
ggplot(data = correlation_long, aes(x = Var1, y = Var2, fill = value)) +
geom_tile() +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, vjust = 1,
size = 12, hjust = 1)) +
labs(x = "Variables", y = "Variables")

# graph time of weather events vs SO2 AQI
summary(LosAngelesDataWeatherEvents)
## X Date Address
## Min. :438394 Min. :2016-01-01 00:00:00.00 Length:1283
## 1st Qu.:470028 1st Qu.:2017-11-15 12:00:00.00 Class :character
## Median :530951 Median :2019-08-02 00:00:00.00 Mode :character
## Mean :533139 Mean :2019-07-20 18:44:31.24
## 3rd Qu.:593778 3rd Qu.:2021-04-24 12:00:00.00
## Max. :626077 Max. :2022-12-31 00:00:00.00
##
## State County City O3.Mean
## California:1283 Los Angeles :1283 Los Angeles:1283 Min. :0.00075
## Alabama : 0 Ada : 0 Albuquerque: 0 1st Qu.:0.01961
## Alaska : 0 Adair : 0 Baton Rouge: 0 Median :0.02829
## Arizona : 0 Adams : 0 Birmingham : 0 Mean :0.02722
## Arkansas : 0 Alameda : 0 Bismarck : 0 3rd Qu.:0.03471
## Colorado : 0 Alexandria City: 0 Cheyenne : 0 Max. :0.06577
## (Other) : 0 (Other) : 0 (Other) : 0
## O3.1st.Max.Value O3.1st.Max.Hour O3.AQI CO.Mean
## Min. :0.00200 Min. : 7.00 Min. : 2.00 Min. :0.1458
## 1st Qu.:0.03300 1st Qu.:10.00 1st Qu.: 31.00 1st Qu.:0.2896
## Median :0.04200 Median :10.00 Median : 39.00 Median :0.3917
## Mean :0.04198 Mean :10.48 Mean : 41.74 Mean :0.4596
## 3rd Qu.:0.05100 3rd Qu.:11.00 3rd Qu.: 47.00 3rd Qu.:0.5675
## Max. :0.11800 Max. :23.00 Max. :214.00 Max. :1.3958
##
## CO.1st.Max.Value CO.1st.Max.Hour CO.AQI SO2.Mean
## Min. :0.2000 Min. : 0.000 Min. : 2.000 Min. :-0.4667
## 1st Qu.:0.4000 1st Qu.: 1.000 1st Qu.: 5.000 1st Qu.: 0.1167
## Median :0.5000 Median : 6.000 Median : 6.000 Median : 0.2625
## Mean :0.6266 Mean : 6.712 Mean : 7.185 Mean : 0.3028
## 3rd Qu.:0.8000 3rd Qu.: 9.000 3rd Qu.: 9.000 3rd Qu.: 0.4521
## Max. :1.8000 Max. :23.000 Max. :20.000 Max. : 1.4917
##
## SO2.1st.Max.Value SO2.1st.Max.Hour SO2.AQI NO2.Mean
## Min. :-0.3000 Min. : 0.00 Min. : 0.0000 Min. : 2.729
## 1st Qu.: 0.3000 1st Qu.: 7.00 1st Qu.: 0.0000 1st Qu.:11.825
## Median : 0.6000 Median :11.00 Median : 0.0000 Median :16.279
## Mean : 0.7158 Mean :10.37 Mean : 0.3297 Mean :18.486
## 3rd Qu.: 0.9000 3rd Qu.:13.00 3rd Qu.: 0.0000 3rd Qu.:24.036
## Max. :17.9000 Max. :23.00 Max. :24.0000 Max. :55.604
##
## NO2.1st.Max.Value NO2.1st.Max.Hour NO2.AQI Year
## Min. : 5.10 Min. : 0.00 Min. : 5.00 Min. :2016
## 1st Qu.:22.95 1st Qu.: 6.00 1st Qu.:21.50 1st Qu.:2017
## Median :30.10 Median : 9.00 Median :28.00 Median :2019
## Mean :31.86 Mean :11.22 Mean :29.67 Mean :2019
## 3rd Qu.:39.80 3rd Qu.:20.00 3rd Qu.:37.00 3rd Qu.:2021
## Max. :80.00 Max. :23.00 Max. :79.00 Max. :2022
##
## StateAbb CityAndState weatherEvents timeLasted
## Length:1283 Length:1283 Mode:logical Min. : 3.0
## Class :character Class :character TRUE:1283 1st Qu.: 113.5
## Mode :character Mode :character Median : 264.0
## Mean : 432.1
## 3rd Qu.: 625.5
## Max. :2981.0
##
plot(LosAngelesDataWeatherEvents$timeLasted, LosAngelesDataWeatherEvents$SO2.AQI)

LosAngelesDataNoWeatherEvents$timeLasted <- 0
# scatter plot of SO2 AQI values vs time of weather events
mergedData <- rbind(LosAngelesDataNoWeatherEvents, LosAngelesDataWeatherEvents)
ggplot(mergedData, aes(x = timeLasted, y = SO2.AQI, color=weatherEvents)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Time of Weather Events in a Day", y = "SO2 AQI",
title = "Daily SO2 AQI in Los Angeles on Days with vs without Weather Events",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(mergedData, aes(x = Date, y = SO2.AQI, color=weatherEvents)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI from 2016-2022 in Los Angeles",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# boxplot of SO2 AQI for LA
boxplot(SO2.AQI~weatherEvents,
data=LosAngelesPollutionData,
main="SO2 AQI for LA on Days With vs Without Weather Events",
xlab="Weather Events Occured",
ylab="SO2 AQI",
col="purple",
border="blue", las = 2
)

# boxplot of O3 AQI for LA
boxplot(O3.AQI~weatherEvents,
data=LosAngelesPollutionData,
main="O3 AQI for LA on Days With vs Without Weather Events",
xlab="Weather Events Occured",
ylab="O3 AQI",
col="purple",
border="blue", las = 2
)

# adding fog as variable to los angeles pollution dataset
LosAngelesFogData <- LosAngelesWeatherData %>% filter(Type=="Fog")
fogDays <- unique(LosAngelesFogData$startDate)
LAFogDaysPollution <- LosAngelesPollutionData %>% filter(Date %in% fogDays)
LosAngelesPollutionData$Fog <- ifelse(LosAngelesPollutionData$Date %in% fogDays, "Foggy", "No fog")
# graphing pollution levels on foggy vs no fog days
ggplot(LosAngelesPollutionData, aes(x = Fog, y = SO2.AQI, color=Fog)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI on Foggy and Non-Foggy Days in Los Angeles",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(breaks = c("Foggy", "Not foggy"),
values=c("cyan2", "darkslategray"))

ggplot(LosAngelesPollutionData, aes(x = Fog, y = O3.AQI, color=Fog)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "O3 AQI",
title = "O3 AQI on Foggy and Non-Foggy Days in Los Angeles",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(breaks = c("Foggy", "Not foggy"),
values=c("cyan2", "darkslategray"))

library(ggstatsplot)
## You can cite this package as:
## Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
## Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ lubridate 1.9.3 ✔ stringr 1.5.1
## ✔ purrr 1.0.2 ✔ tibble 3.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# violin plot with statistical analysis of fog vs no fog for LA
plt <- ggbetweenstats(
data = LosAngelesPollutionData,
x = Fog,
y = SO2.AQI
) + labs(
x = "Presence of Fog",
y = "SO2 AQI",
title = "SO2 AQI Compared Between Foggy and Non-Foggy in Los Angeles"
)
plt

plt2 <- ggbetweenstats(
data = LosAngelesPollutionData,
x = Fog,
y = O3.AQI
) + labs(
x = "Presence of Fog",
y = "O3 AQI",
title = "O3 AQI Compared Between Foggy and Non-Foggy in Los Angeles"
)
plt2

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(LosAngelesPollutionData, aes(x = Date, y = SO2.AQI, color=Fog)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI from 2016-2022 in Los Angeles",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("Foggy", "No fog"),
values=c("cyan2", "darkslategray"), labels = c("Foggy", "No fog"))

# scatter plot of O3 AQI values over time with distinction between days with vs without weather events
ggplot(LosAngelesPollutionData, aes(x = Date, y = O3.AQI, color=Fog)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "O3 AQI from 2016-2022 in Los Angeles",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("Foggy", "No fog"),
values=c("cyan2", "darkslategray"), labels = c("Foggy", "No fog"))

DenverPollutionData <- overlapCitiesPollutionData %>% filter(City=='Denver')%>% filter(State=='Colorado')
DenverWeatherData <- overlapCitiesWeatherData %>% filter(City == 'Denver') %>% filter(State=='CO')
DenverWeatherData$startDate <- as.POSIXct(DenverWeatherData$startDate, format="%Y-%m-%d")
DenverWeatherData$endDate <- as.POSIXct(DenverWeatherData$endDate, format="%Y-%m-%d")
# Reset row names
row.names(DenverWeatherData) <- NULL
# loop to separate weather events that span two days, correct time of event
k=nrow(DenverWeatherData)
DenverWeatherData$StartTime.UTC. <- as.POSIXct(DenverWeatherData$StartTime.UTC., format = "%Y-%m-%d %H:%M:%S")
for(i in 1:k) {
if(as.Date(DenverWeatherData$EndTime.UTC.[i]) > as.Date(DenverWeatherData$StartTime.UTC.[i])) {
DenverWeatherData <- rbind(DenverWeatherData, DenverWeatherData[i,])
row.names(DenverWeatherData) <- NULL
date_only <- as.Date(DenverWeatherData$StartTime.UTC.[i])
next_day <- date_only + 1
midnight_next_day <- as.POSIXct(paste(next_day, "00:00:00"), format = "%Y-%m-%d %H:%M:%S")
time_diff <- difftime(midnight_next_day, DenverWeatherData$StartTime.UTC.[i], units = "mins")
DenverWeatherData$timeEventLasted[nrow(DenverWeatherData)] <- DenverWeatherData$timeEventLasted[i]-time_diff
DenverWeatherData$timeEventLasted[i] <- time_diff
DenverWeatherData$StartTime.UTC.[nrow(DenverWeatherData)] <- as.POSIXct(paste(DenverWeatherData$endDate[i], "00:00:01"))#, format = "%Y-%m-%d %H:%M:%S")
DenverWeatherData$startDate[nrow(DenverWeatherData)] <- as.POSIXct(DenverWeatherData$endDate[i], format = "%Y-%m-%d")
DenverWeatherData$startTime[nrow(DenverWeatherData)] <- "00:00:00"
DenverWeatherData$EndTime.UTC.[i] <- (paste((DenverWeatherData$startDate[i]),"23:59:59"))#, format = "%Y-%m-%d %H:%M:%S")
DenverWeatherData$endDate[i] <- as.POSIXct(DenverWeatherData$startDate[i], format = "%Y-%m-%d")
DenverWeatherData$endTime[i] <- "23:59:59"
}
}
# loop to put true/false if weather event occurred on a given day in LA
DenverWeatherData <- DenverWeatherData[order(DenverWeatherData$startDate),]
row.names(DenverWeatherData) <- NULL
DenverPollutionData$weatherEvents <- NA
DenverPollutionData <- DenverPollutionData[!duplicated(DenverPollutionData$Date), ]
j <- 1
i <- 1
for(i in 1:nrow(DenverPollutionData)) {
while((as.numeric(difftime(DenverPollutionData$Date[i],DenverWeatherData$startDate[j], units="days"))) >0 && j <= nrow(DenverWeatherData)) {
j <- j+1
}
if(as.Date(DenverPollutionData$Date[i])==as.Date(DenverWeatherData$startDate[j]) && j <= nrow(DenverWeatherData)) {
DenverPollutionData$weatherEvents[i] <- TRUE
}
else {
DenverPollutionData$weatherEvents[i] <- FALSE
}
}
# making separate data sets for days with vs without weather events
DenverDataNoWeatherEvents <- DenverPollutionData %>% filter(weatherEvents==FALSE)
DenverDataNoWeatherEvents <- DenverDataNoWeatherEvents[!duplicated(DenverDataNoWeatherEvents$Date), ]
DenverDataWeatherEvents <- DenverPollutionData %>% filter(weatherEvents==TRUE)
DenverDataWeatherEvents <- DenverDataWeatherEvents[!duplicated(DenverDataWeatherEvents$Date), ]
# cleaning data for the loop (order, duplicates for pollution)
row.names(DenverDataWeatherEvents) <- NULL
DenverDataWeatherEvents <- DenverDataWeatherEvents[!duplicated(DenverDataWeatherEvents$Date), ]
DenverDataWeatherEvents <- DenverDataWeatherEvents[order(DenverDataWeatherEvents$Date),]
DenverWeatherData <- DenverWeatherData[order(DenverWeatherData$startDate),]
# calculate total time of weather events for the day
DenverDataWeatherEvents$timeLasted <- NA
j <- 1
i <- 1
for(i in 1:nrow(DenverDataWeatherEvents)) {
p <- 0
while((as.numeric(difftime(as.Date(DenverDataWeatherEvents$Date[i]),as.Date(DenverWeatherData$startDate[j]), units="days")))>0 && j <= nrow(DenverWeatherData)) {
j <- j+1
}
while((as.numeric(difftime(as.Date(DenverDataWeatherEvents$Date[i]),as.Date(DenverWeatherData$startDate[j]), units="days")))==0 && j <= nrow(DenverWeatherData)) {
p <- p + DenverWeatherData$timeEventLasted[j]
j <- j+1
}
DenverDataWeatherEvents$timeLasted[i] <- p
p <- 0
}
# violin plot and histogram: no weather events, SO2 AQI
p <- ggplot(DenverDataNoWeatherEvents, aes(x=City, y=SO2.AQI)) +
geom_violin()
p

hist(DenverDataNoWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in Denver on Days with No Recorded Weather")

# violin plot and histogram: no weather events, O3 AQI
p <- ggplot(DenverDataNoWeatherEvents, aes(x=City, y=O3.AQI)) +
geom_violin()
p

hist(DenverDataNoWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in Denver on Days with No Recorded Weather")

# violin plot and histogram: recorded weather events, SO2 AQI
p <- ggplot(DenverDataWeatherEvents, aes(x=City, y=SO2.AQI)) +
geom_violin()
p

hist(DenverDataWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in Denver on Days with Recorded Weather")

# violin plot and histogram: recorded weather events, O3 AQI
p <- ggplot(DenverDataWeatherEvents, aes(x=City, y=O3.AQI)) +
geom_violin()
p

hist(DenverDataWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in Denver on Days with Recorded Weather")

ClevelandPollutionData <- overlapCitiesPollutionData %>% filter(City=='Cleveland') %>% filter(State=='Ohio')
ClevelandWeatherData <- overlapCitiesWeatherData %>% filter(City == 'Cleveland') %>% filter(State=="OH")
ClevelandWeatherData$startDate <- as.POSIXct(ClevelandWeatherData$startDate, format="%Y-%m-%d")
ClevelandWeatherData$endDate <- as.POSIXct(ClevelandWeatherData$endDate, format="%Y-%m-%d")
# Reset row names
row.names(ClevelandWeatherData) <- NULL
# loop to separate weather events that span two days, correct time of event
k=nrow(ClevelandWeatherData)
ClevelandWeatherData$StartTime.UTC. <- as.POSIXct(ClevelandWeatherData$StartTime.UTC., format = "%Y-%m-%d %H:%M:%S")
for(i in 1:k) {
if(as.Date(ClevelandWeatherData$EndTime.UTC.[i]) > as.Date(ClevelandWeatherData$StartTime.UTC.[i])) {
ClevelandWeatherData <- rbind(ClevelandWeatherData, ClevelandWeatherData[i,])
row.names(ClevelandWeatherData) <- NULL
date_only <- as.Date(ClevelandWeatherData$StartTime.UTC.[i])
next_day <- date_only + 1
midnight_next_day <- as.POSIXct(paste(next_day, "00:00:00"), format = "%Y-%m-%d %H:%M:%S")
time_diff <- difftime(midnight_next_day, ClevelandWeatherData$StartTime.UTC.[i], units = "mins")
ClevelandWeatherData$timeEventLasted[nrow(ClevelandWeatherData)] <- ClevelandWeatherData$timeEventLasted[i]-time_diff
ClevelandWeatherData$timeEventLasted[i] <- time_diff
ClevelandWeatherData$StartTime.UTC.[nrow(ClevelandWeatherData)] <- as.POSIXct(paste(ClevelandWeatherData$endDate[i], "00:00:00"))#, format = "%Y-%m-%d %H:%M:%S") # need to be 00 or 01??
ClevelandWeatherData$startDate[nrow(ClevelandWeatherData)] <- as.POSIXct(ClevelandWeatherData$endDate[i], format = "%Y-%m-%d")
ClevelandWeatherData$startTime[nrow(ClevelandWeatherData)] <- "00:00:00"
ClevelandWeatherData$EndTime.UTC.[i] <- (paste((ClevelandWeatherData$startDate[i]),"23:59:59"))#, format = "%Y-%m-%d %H:%M:%S")
ClevelandWeatherData$endDate[i] <- as.POSIXct(ClevelandWeatherData$startDate[i], format = "%Y-%m-%d")
ClevelandWeatherData$endTime[i] <- "23:59:59"
}
}
# loop to put true/false if weather event occurred on a given day in Cleveland
ClevelandWeatherData <- ClevelandWeatherData[order(ClevelandWeatherData$startDate),]
row.names(ClevelandWeatherData) <- NULL
ClevelandPollutionData$weatherEvents <- NA
ClevelandPollutionData <- ClevelandPollutionData[!duplicated(ClevelandPollutionData$Date), ]
c
## function (...) .Primitive("c")
j <- 1
i <- 1
for(i in 1:nrow(ClevelandPollutionData)) {
while((as.numeric(difftime(ClevelandPollutionData$Date[i],ClevelandWeatherData$startDate[j], units="days"))) >0 && j <= nrow(ClevelandWeatherData)) {
j <- j+1
}
if(as.Date(ClevelandPollutionData$Date[i])==as.Date(ClevelandWeatherData$startDate[j]) && j <= nrow(ClevelandWeatherData)) {
ClevelandPollutionData$weatherEvents[i] <- TRUE
}
else {
ClevelandPollutionData$weatherEvents[i] <- FALSE
}
}
# making separate data sets for days with vs without weather events
ClevelandDataNoWeatherEvents <- ClevelandPollutionData %>% filter(weatherEvents==FALSE)
ClevelandDataNoWeatherEvents <- ClevelandDataNoWeatherEvents[!duplicated(ClevelandDataNoWeatherEvents$Date), ]
ClevelandDataWeatherEvents <- ClevelandPollutionData %>% filter(weatherEvents==TRUE)
ClevelandDataWeatherEvents <- ClevelandDataWeatherEvents[!duplicated(ClevelandDataWeatherEvents$Date), ]
# cleaning data for the loop (order, duplicates for pollution)
row.names(ClevelandDataWeatherEvents) <- NULL
ClevelandDataWeatherEvents <- ClevelandDataWeatherEvents[!duplicated(ClevelandDataWeatherEvents$Date), ]
ClevelandDataWeatherEvents <- ClevelandDataWeatherEvents[order(ClevelandDataWeatherEvents$Date),]
ClevelandWeatherData <- ClevelandWeatherData[order(ClevelandWeatherData$startDate),]
# calculate total time of weather events for the day
ClevelandDataWeatherEvents$timeLasted <- NA
j <- 1
i <- 1
for(i in 1:nrow(ClevelandDataWeatherEvents)) {
p <- 0
while((as.numeric(difftime(as.Date(ClevelandDataWeatherEvents$Date[i]),as.Date(ClevelandWeatherData$startDate[j]), units="days")))>0 && j <= nrow(ClevelandWeatherData)) {
j <- j+1
}
while((as.numeric(difftime(as.Date(ClevelandDataWeatherEvents$Date[i]),as.Date(ClevelandWeatherData$startDate[j]), units="days")))==0 && j <= nrow(ClevelandWeatherData)) {
p <- p + ClevelandWeatherData$timeEventLasted[j]
j <- j+1
}
ClevelandDataWeatherEvents$timeLasted[i] <- p
p <- 0
}
# violin plot and histogram: no weather events, SO2 AQI
p <- ggplot(ClevelandDataNoWeatherEvents, aes(x=City, y=SO2.AQI)) +
geom_violin()
p

hist(ClevelandDataNoWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in Cleveland on Days with No Recorded Weather")

# violin plot and histogram: no weather events, O3 AQI
p <- ggplot(ClevelandDataNoWeatherEvents, aes(x=City, y=O3.AQI)) +
geom_violin()
p

hist(ClevelandDataNoWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in Cleveland on Days with No Recorded Weather")

# violin plot and histogram: recorded weather events, SO2 AQI
p <- ggplot(ClevelandDataWeatherEvents, aes(x=City, y=SO2.AQI)) +
geom_violin()
p

hist(ClevelandDataWeatherEvents$SO2.AQI, breaks=20, main="SO2 AQI in Cleveland on Days with Recorded Weather")

# violin plot and histogram: recorded weather events, O3 AQI
p <- ggplot(ClevelandDataWeatherEvents, aes(x=City, y=O3.AQI)) +
geom_violin()
p

hist(ClevelandDataWeatherEvents$O3.AQI, breaks=20, main="O3 AQI in Cleveland on Days with Recorded Weather")

# graph time of weather events vs SO2 AQI
plot(ClevelandDataWeatherEvents$timeLasted, ClevelandDataWeatherEvents$SO2.AQI)

ClevelandDataNoWeatherEvents$timeLasted <- 0
# scatter plot of SO2 AQI values vs time of weather events
mergedClevelandData <- rbind(ClevelandDataNoWeatherEvents, ClevelandDataWeatherEvents)
ggplot(mergedClevelandData, aes(x = timeLasted, y = SO2.AQI, color=weatherEvents)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Time of Weather Events in a Day", y = "SO2 AQI",
title = "Daily SO2 AQI in Cleveland on Days with vs without Weather Events",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(mergedClevelandData, aes(x = Date, y = SO2.AQI, color=weatherEvents)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI from 2016-2022 in Cleveland",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# boxplot of SO2 AQI for Cleveland
boxplot(SO2.AQI~weatherEvents,
data=ClevelandPollutionData,
main="SO2 AQI for Cleveland on Days With vs Without Weather Events",
xlab="Weather Events Occured",
ylab="SO2 AQI",
col="purple",
border="blue", las = 2
)

# boxplot of O3 AQI for Cleveland
boxplot(O3.AQI~weatherEvents,
data=ClevelandPollutionData,
main="O3 AQI for Cleveland on Days With vs Without Weather Events",
xlab="Weather Events Occured",
ylab="O3 AQI",
col="purple",
border="blue", las = 2
)

# adding fog as variable to Cleveland pollution dataset
ClevelandFogData <- ClevelandWeatherData %>% filter(Type=="Fog")
fogDays <- unique(ClevelandFogData$startDate)
ClevelandFogDaysPollution <- ClevelandPollutionData %>% filter(Date %in% fogDays)
ClevelandPollutionData$Fog <- ifelse(ClevelandPollutionData$Date %in% fogDays, "Foggy", "No fog")
# graphing pollution levels on foggy vs no fog days
ggplot(ClevelandPollutionData, aes(x = Fog, y = SO2.AQI, color=Fog)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI on Foggy and Non-Foggy Days in Cleveland",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(breaks = c("Foggy", "No fog"),
values=c("cyan2", "darkslategray"))

library(ggstatsplot)
library(tidyverse)
# violin plot with statistical analysis of fog vs no fog for Cleveland
plt <- ggbetweenstats(
data = ClevelandPollutionData,
x = Fog,
y = SO2.AQI
) + labs(
x = "Presence of Fog",
y = "SO2 AQI",
title = "SO2 AQI Compared Between Foggy and Non-Foggy in Cleveland"
)
plt

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(ClevelandPollutionData, aes(x = Date, y = SO2.AQI, color=Fog)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI from 2016-2022 in Cleveland",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("Foggy", "No fog"),
values=c("cyan2", "darkslategray"), labels = c("Foggy", "No fog"))

summary(ClevelandDataWeatherEvents$timeLasted)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.0 180.0 580.5 864.5 1287.0 4971.0
# graph time of weather events vs SO2 AQI
plot(DenverDataWeatherEvents$timeLasted, DenverDataWeatherEvents$SO2.AQI)

DenverDataNoWeatherEvents$timeLasted <- 0
# scatter plot of SO2 AQI values vs time of weather events
mergedDenverData <- rbind(DenverDataNoWeatherEvents, DenverDataWeatherEvents)
ggplot(mergedDenverData, aes(x = timeLasted, y = SO2.AQI, color=weatherEvents)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Time of Weather Events in a Day", y = "SO2 AQI",
title = "Daily SO2 AQI in Denver on Days with vs without Weather Events",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(mergedDenverData, aes(x = Date, y = SO2.AQI, color=weatherEvents)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI from 2016-2022 in Denver",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("FALSE", "TRUE"),
values=c("deeppink", "cornflowerblue"), labels = c("No weather events", "Weather events"))

# boxplot of SO2 AQI for Denver
boxplot(SO2.AQI~weatherEvents,
data=DenverPollutionData,
main="SO2 AQI for Denver on Days With vs Without Weather Events",
xlab="Weather Events Occured",
ylab="SO2 AQI",
col="purple",
border="blue", las = 2
)

# boxplot of O3 AQI for Denver
boxplot(O3.AQI~weatherEvents,
data=DenverPollutionData,
main="O3 AQI for Denver on Days With vs Without Weather Events",
xlab="Weather Events Occured",
ylab="O3 AQI",
col="purple",
border="blue", las = 2
)

# adding fog as variable to Denver pollution dataset
DenverFogData <- DenverWeatherData %>% filter(Type=="Fog")
fogDays <- unique(DenverFogData$startDate)
DenverFogDaysPollution <- DenverPollutionData %>% filter(Date %in% fogDays)
DenverPollutionData$Fog <- ifelse(DenverPollutionData$Date %in% fogDays, "Foggy", "No fog")
# graphing pollution levels on foggy vs no fog days
ggplot(DenverPollutionData, aes(x = Fog, y = SO2.AQI, color=Fog)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI on Foggy and Non-Foggy Days in Denver",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(breaks = c("Foggy", "No fog"),
values=c("cyan2", "darkslategray"))

library(ggstatsplot)
library(tidyverse)
# violin plot with statistical analysis of fog vs no fog for LA
plt <- ggbetweenstats(
data = DenverPollutionData,
x = Fog,
y = SO2.AQI
) + labs(
x = "Presence of Fog",
y = "SO2 AQI",
title = "SO2 AQI Compared Between Foggy and Non-Foggy in Denver"
)
plt

# scatter plot of SO2 AQI values over time with distinction between days with vs without weather events
ggplot(DenverPollutionData, aes(x = Date, y = SO2.AQI, color=Fog)) +
geom_point() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(x = "Date", y = "SO2 AQI",
title = "SO2 AQI from 2016-2022 in Denver",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(name = "Legend", breaks = c("Foggy", "No fog"),
values=c("cyan2", "darkslategray"), labels = c("Foggy", "No fog"))

pltDen <- ggbetweenstats(
data = DenverPollutionData,
x = weatherEvents,
y = SO2.AQI
) + labs(
x = "Presence of Weather Events",
y = "SO2 AQI",
title = "SO2 AQI Between Days With and Without Weather in Denver"
)
plt

pltLA <- ggbetweenstats(
data = LosAngelesPollutionData,
x = weatherEvents,
y = SO2.AQI
) + labs(
x = "Presence of Weather Events",
y = "SO2 AQI",
title = "SO2 AQI Between Days With and Without Weather in LA"
)
pltLA

pltCle <- ggbetweenstats(
data = ClevelandPollutionData,
x = weatherEvents,
y = SO2.AQI
) + labs(
x = "Presence of Weather Events",
y = "SO2 AQI",
title = "SO2 AQI Between Days With and Without Weather in Cleveland"
)
pltCle

CleLADenPollution <- rbind(ClevelandPollutionData, LosAngelesPollutionData, DenverPollutionData)
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin,
draw_group = function(self, data, ..., draw_quantiles = NULL) {
data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
grp <- data[1, "group"]
newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
1))
quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
aesthetics$alpha <- rep(1, nrow(quantiles))
both <- cbind(quantiles, aesthetics)
quantile_grob <- GeomPath$draw_panel(both, ...)
ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
}
else {
ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
}
})
geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ...,
draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE,
show.legend = NA, inherit.aes = TRUE) {
layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin,
position = position, show.legend = show.legend, inherit.aes = inherit.aes,
params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}
ggplot(CleLADenPollution, aes(City, SO2.AQI, fill = weatherEvents)) + geom_split_violin()

ggplot(data = CleLADenPollution, aes (x = weatherEvents, y = SO2.AQI, fill=weatherEvents)) + geom_violin(position = "dodge") + facet_wrap(~City)

t.test(DenverDataNoWeatherEvents$SO2.AQI, DenverDataWeatherEvents$SO2.AQI, alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: DenverDataNoWeatherEvents$SO2.AQI and DenverDataWeatherEvents$SO2.AQI
## t = 6.0243, df = 1864, p-value = 2.042e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.5138864 1.0099953
## sample estimates:
## mean of x mean of y
## 2.119004 1.357063
t.test(LosAngelesDataNoWeatherEvents$SO2.AQI, LosAngelesDataWeatherEvents$SO2.AQI, alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: LosAngelesDataNoWeatherEvents$SO2.AQI and LosAngelesDataWeatherEvents$SO2.AQI
## t = 1.3113, df = 2390.9, p-value = 0.1899
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.02453863 0.12360448
## sample estimates:
## mean of x mean of y
## 0.379229 0.329696
t.test(ClevelandDataNoWeatherEvents$SO2.AQI, ClevelandDataWeatherEvents$SO2.AQI, alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: ClevelandDataNoWeatherEvents$SO2.AQI and ClevelandDataWeatherEvents$SO2.AQI
## t = 8.1233, df = 1074.5, p-value = 1.237e-15
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2.731518 4.471355
## sample estimates:
## mean of x mean of y
## 6.885204 3.283768
DenverFogPollution <- DenverPollutionData %>% filter(Fog == "Foggy")
DenverNoFogPollution <- DenverPollutionData %>% filter(Fog == "No fog")
t.test(DenverFogPollution$SO2.AQI, DenverNoFogPollution$SO2.AQI, alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: DenverFogPollution$SO2.AQI and DenverNoFogPollution$SO2.AQI
## t = -2.7745, df = 560.04, p-value = 0.005713
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.6451957 -0.1103281
## sample estimates:
## mean of x mean of y
## 1.520124 1.897886
LANoFogPollution <- LosAngelesPollutionData %>% filter(Fog == "No fog")
t.test(LAFogDaysPollution$SO2.AQI, LANoFogPollution$SO2.AQI, alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: LAFogDaysPollution$SO2.AQI and LANoFogPollution$SO2.AQI
## t = 1.4167, df = 1568.5, p-value = 0.1568
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.02289013 0.14194291
## sample estimates:
## mean of x mean of y
## 0.3900293 0.3305029
ClevelandNoFog <- ClevelandPollutionData %>% filter(Fog=="No fog")
t.test(ClevelandFogDaysPollution$SO2.AQI, ClevelandNoFog$SO2.AQI, alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: ClevelandFogDaysPollution$SO2.AQI and ClevelandNoFog$SO2.AQI
## t = -3.4819, df = 1399.5, p-value = 0.0005131
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.0662502 -0.5770584
## sample estimates:
## mean of x mean of y
## 3.478571 4.800226
plot1 <- ggplot(data = LosAngelesPollutionData, aes(x=Fog, y=SO2.AQI)) + geom_boxplot(aes(fill=Fog)) + labs(x = "", y = "SO2 AQI",
title = "Los Angeles") + theme(legend.position="none") + ylim(0,70)
plot2 <- ggplot(data = DenverPollutionData, aes(x=Fog, y=SO2.AQI)) + geom_boxplot(aes(fill=Fog)) + labs(x = "", y = "SO2 AQI",
title = "Denver") + theme(legend.position="none")
plot3 <- ggplot(data = ClevelandPollutionData, aes(x=Fog, y=SO2.AQI)) + geom_boxplot(aes(fill=Fog)) + labs(x = "", y = "SO2 AQI",
title = "Cleveland") + theme(legend.position="none")
require(gridExtra)
## Loading required package: gridExtra
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
grid.arrange(plot1, plot2, plot3, ncol=3)

LAFogDaysPollution$Fog <- "Foggy"
DenverFogDaysPollution$Fog <- "Foggy"
ClevelandFogDaysPollution$Fog <- "Foggy"
mergedFogData <- rbind(LAFogDaysPollution, DenverFogDaysPollution, ClevelandFogDaysPollution)
mergedNoFogData <- rbind(LANoFogPollution, DenverNoFogPollution, ClevelandNoFog)
mergedFogAndNoFogData <- rbind(mergedFogData, mergedNoFogData)
t.test(mergedFogData$SO2.AQI, mergedNoFogData$SO2.AQI, alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: mergedFogData$SO2.AQI and mergedNoFogData$SO2.AQI
## t = -5.5553, df = 4081.8, p-value = 2.948e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.040657 -0.497734
## sample estimates:
## mean of x mean of y
## 1.625122 2.394318
ggplot(mergedFogAndNoFogData, aes(x = Fog, y = SO2.AQI, color=Fog)) +
geom_violin() +
scale_y_continuous(label = function(x) {return(paste(x, "AQI"))}) +
labs(y = "SO2 AQI",
title = "SO2 AQI on Foggy and Non-Foggy Days in Los Angeles and Denver",
caption = "Data: Kaggle") +
theme(plot.title = element_text(size = 14, face = "bold.italic"),
plot.caption = element_text(hjust = 0)) +
scale_color_manual(breaks = c("Foggy", "No fog"),
values=c("cyan2", "darkslategray"))

plt <- ggbetweenstats(
data = mergedFogAndNoFogData,
x = Fog,
y = SO2.AQI
) + labs(
x = "Presence of Fog",
y = "SO2 AQI",
title = "SO2 AQI Compared Between Foggy and Non-Foggy in Los Angeles, Denver, and Cleveland"
)
plt

ggplot(pollutionData, aes(x = Date, y = SO2.AQI, color=City)) +
geom_point() + theme(legend.position="none")

ggplot(overlapCitiesPollutionData, aes(x = Date, y = SO2.AQI)) +
geom_line(position = "dodge", stat = "identity", aes(col = City)) +
facet_wrap(~City) +
theme(legend.position = "none") +
labs(title = "SO2 AQI vs Date", x = "Date", y = "SO2 AQI", caption = "source: Kaggle")
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`

ggplot(pollutionData, aes(x = Date, y = O3.AQI, color=City)) +
geom_point() + theme(legend.position="none")

ggplot(overlapCitiesPollutionData, aes(x = Date, y = O3.AQI)) +
geom_line(position = "dodge", stat = "identity", aes(col = City)) +
facet_wrap(~City) +
theme(legend.position = "none") +
labs(title = "O3 AQI vs Date", x = "Date", y = "O3 AQI", caption = "source: Kaggle")
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`

ggplot(pollutionData, aes(x = Date, y = CO.AQI, color=City)) +
geom_point() + theme(legend.position="none")

ggplot(overlapCitiesPollutionData, aes(x = Date, y = CO.AQI)) +
geom_line(position = "dodge", stat = "identity", aes(col = City)) +
facet_wrap(~City) +
theme(legend.position = "none") +
labs(title = "CO AQI vs Date", x = "Date", y = "CO AQI", caption = "source: Kaggle")
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`

ggplot(pollutionData, aes(x = Date, y = NO2.AQI, color=City)) +
geom_point() + theme(legend.position="none")

ggplot(overlapCitiesPollutionData, aes(x = Date, y = NO2.AQI)) +
geom_line(position = "dodge", stat = "identity", aes(col = City)) +
facet_wrap(~City) +
theme(legend.position = "none") +
labs(title = "NO2 AQI vs Date", x = "Date", y = "NO2 AQI", caption = "source: Kaggle")
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`
