Environmental Data Analysis Project
Shihab Shahriar; PSID:2159285
2023-12-11
This project was done as a part of the course “Environmental Data Analysis” instructed by Dr. Yuxuan Wang (Department of Earth and Atmospherics Sciences, University of Houston)
Study Site
Location: Dhaka, Bangladesh
Air Monitoring Station Coordinates: 23.7561, 90.3872
Rationale: One of most polluted cities in the world.
Study Site (Dhaka, Bangladesh)
Dataset
The dataset described here offers a comprehensive view of environmental and air quality metrics, as monitored by the Department of Environment (DoE) of the Government of Bangladesh. It presents daily average readings (temporal resolution of 24 hours) covering a substantial period from January 1, 2013, to May 31, 2019. The variables included in this dataset encompass a range of criteria air pollutants, such as Fine particulate matter (PM2.5), Coarse particulate matter (PM10), sulfur dioxide (SO2), nitrogen oxides (NOx), carbon monoxide (CO), and ozone (O3). In addition, the dataset provides detailed meteorological data, including wind speed (WS), wind direction (WD), mean temperature (T), relative humidity (RH), solar radiation (SR), barometric pressure (BP), and rainfall (R).
Now, let’s look at the datset
#Load the required libraries first
install_and_load <- function(package) {
if (!require(package, character.only = TRUE)) {
install.packages(package)
library(package, character.only = TRUE)
}
}
install_and_load("readr")
install_and_load("leaflet")
install_and_load("openair")
install_and_load("tidyverse")
install_and_load("pastecs")
install_and_load("psych")
install_and_load("ggplot2")
install_and_load("reshape2")
install_and_load("VIM")
install_and_load("lubridate")
install_and_load("corrr")
install_and_load("dplyr")
install_and_load("gridExtra")
install_and_load("Kendall")
install_and_load("car")
install_and_load("randomForest")
install_and_load("Metrics")
install_and_load("caret")
# Read the CSV file
data <- read_csv("dhaka.csv")
data<-as.data.frame(data)
names(data)## [1] "Date" "PM2.5" "WS" "WD" "T" "RH" "SR" "BP" "R"
## [10] "NOX" "SO2" "O3" "CO"# Print the number of rows and columns
cat("Number of Rows:", nrow(data), "\n")## Number of Rows: 2191cat("Number of Columns:", ncol(data), "\n")## Number of Columns: 13cat("Dimensions of the dataset:", dim(data), "\n")## Dimensions of the dataset: 2191 13print(head(data))## Date PM2.5 WS WD T RH SR BP R
## 1 1/1/13 274.5879 1.3795 206.5433 17.8800 69.8887 140.4457 1013.258 0.0583
## 2 1/2/13 292.1430 1.4104 190.9030 18.6961 67.1973 164.8315 1013.401 0.0467
## 3 1/3/13 297.3933 1.6614 196.1533 20.0829 68.4158 93.8048 1014.109 0.0482
## 4 1/4/13 266.0096 1.2790 285.2713 20.1200 70.4625 114.0026 1013.962 0.0226
## 5 1/5/13 189.6463 1.3413 263.0129 18.7043 62.1733 131.8348 1013.007 0.0200
## 6 1/6/13 205.8704 1.6908 263.6258 18.3696 58.7264 118.7886 1013.552 0.0200
## NOX SO2 O3 CO
## 1 127.3578 4.2444 3.88 1.9007
## 2 303.8237 29.5099 12.38 3.5839
## 3 385.6099 29.3615 10.54 3.3892
## 4 377.5306 29.0645 12.40 3.5381
## 5 264.2676 40.6711 24.28 3.1488
## 6 125.3342 27.1825 16.66 2.3244# Count and print the number of missing values in each column
cat("Number of Missing Values in Each Column:\n")## Number of Missing Values in Each Column:print(colSums(is.na(data)))## Date PM2.5 WS WD T RH SR BP R NOX SO2 O3 CO
## 0 20 0 0 0 0 0 0 0 10 12 7 12#Data statistics
data<-import(
file = "dhaka.csv",
file.type = "csv",
sep = ",",
header.at = 1,
data.at = 2,
date = "Date",
date.format = "%m/%d/%Y",
time = NULL,
time.format = NULL)## date1 date2 PM2.5 WS WD T RH SR
## "POSIXct" "POSIXt" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## BP R NOX SO2 O3 CO
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"#descriptive statistics
describe(data)#Time series
timePlot(data,
pollutant = c("PM2.5", "NOX", "SO2", "CO", "O3"),
y.relation = "free")
timePlot(data,
pollutant = c("WS", "WD", "T", "RH", "SR", "BP", "R"),
y.relation = "free")
#From the time series, it is clear that, there are multiple outliers across the dataset.
#Boxplots
# Remove the 'Date' column for the boxplot analysis
data_no_date <- data[, !(names(data) %in% c("date"))]
data_long <- melt(data_no_date, variable.name = "Variable", value.name = "Value")
ggplot(data_long, aes(x = Variable, y = Value)) +
geom_boxplot() +
facet_wrap(~ Variable, scales = "free") + theme_bw()+
theme(axis.text.x = element_text(hjust = 1)) +
labs(title = "Boxplots of Variables")
Data Preprocessing
From the descriptive statistics and figures it is clear that the data has missing values and outliers. Therefore, for the further analysis, these should be corrected.
#DATA PRPROCESSING
#Missing data correction
#For missing data correction I used KNN method, as similarly used by (Shahriar et al. 2021)
data_imputed <- kNN(data, k = 5)
head(data_imputed)#Outlier removal
remove_outliers <- function(df, column) {
if(column == "R") {
return(df)
}
Q1 <- quantile(df[[column]], 0.25, na.rm = TRUE)
Q3 <- quantile(df[[column]], 0.75, na.rm = TRUE)
IQR <- Q3 - Q1
df <- df %>%
dplyr::filter(df[[column]] >= (Q1 - 1.5 * IQR) & df[[column]] <= (Q3 + 1.5 * IQR))
return(df)
}
data_cleaned <- data_imputed
# Apply the function to each column
for(column in names(data_cleaned)) {
data_cleaned <- remove_outliers(data_cleaned, column)
}
#Now, lets check the boxplots after outlier removal
data_cleaned_no_date <- select(data_cleaned, -date)
data_long <- melt(data_cleaned_no_date, variable.name = "Variable", value.name = "Value")
ggplot(data_long, aes(x = Variable, y = Value)) +
geom_boxplot() +
facet_wrap(~ Variable, scales = "free") +
theme_bw() +
theme(axis.text.x = element_text(hjust = 1)) +
labs(title = "Boxplots of Variables After Outlier Removal")
#Now, lets check the time series after outlier removal
timePlot(data_cleaned,
pollutant = c("PM2.5", "NOX", "SO2", "CO", "O3"),
y.relation = "free")
timePlot(data_cleaned,
pollutant = c("WS", "WD", "T", "RH", "SR", "BP", "R"),
y.relation = "free")
Results
#Correlation among the variables
# Making sure the Date column is in the correct format
data_cleaned$Date <- as.Date(data_cleaned$date, format = "%Y-%m-%d") # Adjust the format if needed
# Adding the Season column
data_cleaned$Season <- ifelse(month(data_cleaned$Date) %in% c(3,4,5), "Spring",
ifelse(month(data_cleaned$Date) %in% c(6,7,8), "Summer",
ifelse(month(data_cleaned$Date) %in% c(9,10,11), "Autumn",
"Winter")))
# Splitting the data by Season
spring_data <- data_cleaned[data_cleaned$Season == "Spring",]
summer_data <- data_cleaned[data_cleaned$Season == "Summer",]
autumn_data <- data_cleaned[data_cleaned$Season == "Autumn",]
winter_data <- data_cleaned[data_cleaned$Season == "Winter",]
# Function to calculate and return correlation matrix
calculate_correlation <- function(data) {
numeric_data <- select(data, where(is.numeric))
return(cor(numeric_data, use = "complete.obs"))
}
# Calculate correlations for each season
correlation_spring <- calculate_correlation(spring_data)
correlation_summer <- calculate_correlation(summer_data)
correlation_autumn <- calculate_correlation(autumn_data)
correlation_winter <- calculate_correlation(winter_data)
# View the correlations
correlation_spring## PM2.5 WS WD T RH SR
## PM2.5 1.00000000 -0.468380830 0.47709485 -0.28125522 -0.63208349 0.093914706
## WS -0.46838083 1.000000000 -0.20221932 0.24283467 0.16085856 0.223655587
## WD 0.47709485 -0.202219322 1.00000000 -0.08659792 -0.44715592 0.161945783
## T -0.28125522 0.242834673 -0.08659792 1.00000000 -0.09597448 0.336317452
## RH -0.63208349 0.160858562 -0.44715592 -0.09597448 1.00000000 -0.488603795
## SR 0.09391471 0.223655587 0.16194578 0.33631745 -0.48860380 1.000000000
## BP 0.61908314 -0.349178641 0.34925971 -0.38852314 -0.47686327 0.046367233
## R -0.12431275 -0.003722512 -0.06419161 -0.13395553 0.26378492 -0.218521255
## NOX 0.52058831 -0.228488712 0.33169716 -0.36851112 -0.35750498 0.012357130
## SO2 0.43884067 -0.288378607 0.18633912 -0.20183064 -0.35277312 -0.027376470
## O3 0.19022502 -0.142750586 -0.06015740 -0.12260203 -0.08388494 0.008471684
## CO 0.45026364 -0.160786117 0.17822990 -0.13339467 -0.43186060 0.183374940
## BP R NOX SO2 O3 CO
## PM2.5 0.61908314 -0.124312747 0.52058831 0.43884067 0.190225021 0.4502636
## WS -0.34917864 -0.003722512 -0.22848871 -0.28837861 -0.142750586 -0.1607861
## WD 0.34925971 -0.064191613 0.33169716 0.18633912 -0.060157401 0.1782299
## T -0.38852314 -0.133955532 -0.36851112 -0.20183064 -0.122602034 -0.1333947
## RH -0.47686327 0.263784916 -0.35750498 -0.35277312 -0.083884938 -0.4318606
## SR 0.04636723 -0.218521255 0.01235713 -0.02737647 0.008471684 0.1833749
## BP 1.00000000 -0.116789774 0.43406895 0.37609806 0.078756937 0.3658932
## R -0.11678977 1.000000000 -0.11615384 -0.04554432 -0.032756569 -0.1267255
## NOX 0.43406895 -0.116153843 1.00000000 0.37266005 0.104588478 0.2820777
## SO2 0.37609806 -0.045544315 0.37266005 1.00000000 0.081771302 0.2755671
## O3 0.07875694 -0.032756569 0.10458848 0.08177130 1.000000000 0.1935076
## CO 0.36589323 -0.126725497 0.28207766 0.27556715 0.193507646 1.0000000correlation_summer## PM2.5 WS WD T RH SR
## PM2.5 1.00000000 -0.26452831 0.34433277 0.025304364 -0.147804379 0.01028049
## WS -0.26452831 1.00000000 -0.31183312 -0.044417664 -0.008718180 0.08818037
## WD 0.34433277 -0.31183312 1.00000000 0.157937595 -0.090479457 0.08504568
## T 0.02530436 -0.04441766 0.15793759 1.000000000 -0.630864314 0.55187692
## RH -0.14780438 -0.00871818 -0.09047946 -0.630864314 1.000000000 -0.78078341
## SR 0.01028049 0.08818037 0.08504568 0.551876918 -0.780783406 1.00000000
## BP 0.15026544 -0.19130361 0.32281270 0.048450593 -0.086062643 0.13392280
## R -0.01358759 0.09740938 -0.09840304 -0.246798416 0.209091243 -0.17439832
## NOX 0.12225675 -0.09032777 -0.02977606 -0.054497683 -0.049136096 0.02123323
## SO2 0.10004179 -0.06796833 0.15666881 -0.006226037 -0.001150890 -0.05121380
## O3 -0.02934832 0.08885273 -0.15904600 -0.037026689 -0.062807869 0.01245162
## CO 0.03153042 -0.11565130 0.19668125 0.098269148 0.006015649 -0.02654156
## BP R NOX SO2 O3 CO
## PM2.5 0.15026544 -0.01358759 0.12225675 0.100041794 -0.02934832 0.031530418
## WS -0.19130361 0.09740938 -0.09032777 -0.067968330 0.08885273 -0.115651304
## WD 0.32281270 -0.09840304 -0.02977606 0.156668808 -0.15904600 0.196681251
## T 0.04845059 -0.24679842 -0.05449768 -0.006226037 -0.03702669 0.098269148
## RH -0.08606264 0.20909124 -0.04913610 -0.001150890 -0.06280787 0.006015649
## SR 0.13392280 -0.17439832 0.02123323 -0.051213804 0.01245162 -0.026541563
## BP 1.00000000 -0.03613480 -0.02110713 -0.065849277 -0.21137717 0.038212052
## R -0.03613480 1.00000000 0.03448533 0.113715826 0.03589604 -0.039076089
## NOX -0.02110713 0.03448533 1.00000000 0.127921191 0.08199121 0.100756644
## SO2 -0.06584928 0.11371583 0.12792119 1.000000000 0.18428131 0.303847259
## O3 -0.21137717 0.03589604 0.08199121 0.184281307 1.00000000 0.054661516
## CO 0.03821205 -0.03907609 0.10075664 0.303847259 0.05466152 1.000000000correlation_autumn## PM2.5 WS WD T RH SR
## PM2.5 1.0000000 -0.33844665 0.389172836 -0.50994768 -0.5704325 -0.10563615
## WS -0.3384467 1.00000000 -0.351800413 0.12221633 0.1728196 0.12362025
## WD 0.3891728 -0.35180041 1.000000000 -0.21065012 -0.4681600 0.18505484
## T -0.5099477 0.12221633 -0.210650116 1.00000000 0.2455939 0.32686821
## RH -0.5704325 0.17281963 -0.468159951 0.24559389 1.0000000 -0.48681705
## SR -0.1056362 0.12362025 0.185054837 0.32686821 -0.4868170 1.00000000
## BP 0.6230176 -0.34343269 0.330598087 -0.61696746 -0.5229358 -0.05672371
## R -0.1796929 0.02835139 -0.214321734 -0.01013060 0.3447703 -0.25981446
## NOX 0.4336507 -0.14800202 0.176071512 -0.28804922 -0.4208904 0.02900015
## SO2 0.1912294 -0.05524058 0.007243236 -0.23615558 -0.1194771 -0.01568423
## O3 0.3147946 0.03888117 0.196168845 -0.46312445 -0.3068702 0.05887616
## CO 0.1981375 0.01551647 0.027896881 -0.00239113 -0.1521861 0.06868677
## BP R NOX SO2 O3 CO
## PM2.5 0.62301761 -0.17969290 0.43365073 0.191229434 0.31479455 0.19813751
## WS -0.34343269 0.02835139 -0.14800202 -0.055240576 0.03888117 0.01551647
## WD 0.33059809 -0.21432173 0.17607151 0.007243236 0.19616885 0.02789688
## T -0.61696746 -0.01013060 -0.28804922 -0.236155577 -0.46312445 -0.00239113
## RH -0.52293578 0.34477028 -0.42089043 -0.119477088 -0.30687022 -0.15218606
## SR -0.05672371 -0.25981446 0.02900015 -0.015684235 0.05887616 0.06868677
## BP 1.00000000 -0.10356063 0.45274727 0.238760752 0.37426870 0.06599921
## R -0.10356063 1.00000000 -0.15004174 -0.131951639 -0.02575476 -0.08424288
## NOX 0.45274727 -0.15004174 1.00000000 0.237952443 0.23472822 0.39575231
## SO2 0.23876075 -0.13195164 0.23795244 1.000000000 0.22330694 0.32564518
## O3 0.37426870 -0.02575476 0.23472822 0.223306940 1.00000000 0.21138418
## CO 0.06599921 -0.08424288 0.39575231 0.325645180 0.21138418 1.00000000correlation_winter## PM2.5 WS WD T RH
## PM2.5 1.00000000 -0.334736958 0.05656918 -0.384035429 0.114062500
## WS -0.33473696 1.000000000 -0.07795454 0.006317882 0.022393633
## WD 0.05656918 -0.077954542 1.00000000 -0.396428102 0.046425562
## T -0.38403543 0.006317882 -0.39642810 1.000000000 -0.393762324
## RH 0.11406250 0.022393633 0.04642556 -0.393762324 1.000000000
## SR -0.36387128 0.367563557 0.05889684 0.240235994 -0.318274878
## BP 0.24616893 -0.101037423 0.15903913 -0.427438966 0.166884821
## R -0.02330109 0.105068500 -0.01003532 -0.057764384 -0.009758111
## NOX 0.14725839 0.073240947 0.08075686 -0.142601287 0.026237654
## SO2 0.12233901 -0.263774918 0.05790841 0.070819659 -0.205859335
## O3 -0.11392196 0.134492268 0.01203551 0.082259461 -0.346915763
## CO 0.04802255 -0.096587859 0.04161272 0.022679386 0.070051237
## SR BP R NOX SO2 O3
## PM2.5 -0.36387128 0.24616893 -0.023301088 0.14725839 0.12233901 -0.11392196
## WS 0.36756356 -0.10103742 0.105068500 0.07324095 -0.26377492 0.13449227
## WD 0.05889684 0.15903913 -0.010035317 0.08075686 0.05790841 0.01203551
## T 0.24023599 -0.42743897 -0.057764384 -0.14260129 0.07081966 0.08225946
## RH -0.31827488 0.16688482 -0.009758111 0.02623765 -0.20585933 -0.34691576
## SR 1.00000000 -0.03438782 0.018894905 0.06241821 -0.22296350 0.04476040
## BP -0.03438782 1.00000000 -0.030328019 0.13451596 0.01322396 -0.23727079
## R 0.01889491 -0.03032802 1.000000000 0.03217615 -0.02332532 0.03162175
## NOX 0.06241821 0.13451596 0.032176145 1.00000000 -0.11447602 0.01880567
## SO2 -0.22296350 0.01322396 -0.023325318 -0.11447602 1.00000000 0.23627800
## O3 0.04476040 -0.23727079 0.031621747 0.01880567 0.23627800 1.00000000
## CO -0.04328178 -0.01646454 -0.029261298 0.23987194 0.10180341 0.05454941
## CO
## PM2.5 0.04802255
## WS -0.09658786
## WD 0.04161272
## T 0.02267939
## RH 0.07005124
## SR -0.04328178
## BP -0.01646454
## R -0.02926130
## NOX 0.23987194
## SO2 0.10180341
## O3 0.05454941
## CO 1.00000000#Plotting
# Convert the correlation matrices into long format for ggplot
long_spring <- melt(correlation_spring)
long_summer <- melt(correlation_summer)
long_autumn <- melt(correlation_autumn)
long_winter <- melt(correlation_winter)
# Function to create a heatmap
create_heatmap <- function(data, title) {
ggplot(data, aes(Var1, Var2, fill = value)) +
geom_tile() +
scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0, limit = c(-1, 1)) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
labs(fill = "Correlation", title = title)
}
# Create heatmaps for each season
heatmap_spring <- create_heatmap(long_spring, "Spring")
heatmap_summer <- create_heatmap(long_summer, "Summer")
heatmap_autumn <- create_heatmap(long_autumn, "Autumn")
heatmap_winter <- create_heatmap(long_winter, "Winter")
# Combine the heatmaps into one figure
combined_heatmap <- grid.arrange(heatmap_spring, heatmap_summer, heatmap_autumn, heatmap_winter, ncol = 2)
#Trend Test
#First checking the distribution
pollutants <- c("PM2.5", "NOX", "CO", "O3", "SO2")
# Function to check distribution
check_distribution <- function(data, pollutant) {
# Histogram
ggplot(data, aes_string(x = pollutant)) +
geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
labs(title = paste("Histogram of", pollutant), x = pollutant, y = "Frequency") +
theme_minimal()
# Q-Q Plot
qqnorm(data[[pollutant]], main = paste("Q-Q Plot of", pollutant))
qqline(data[[pollutant]], col = "red")
# Shapiro-Wilk Test
shapiro_test <- shapiro.test(data[[pollutant]])
print(paste("Shapiro-Wilk test for", pollutant, "p-value:", shapiro_test$p.value))
}
# Apply the function to each pollutant
for (pollutant in pollutants) {
check_distribution(data_cleaned, pollutant)
}
## [1] "Shapiro-Wilk test for PM2.5 p-value: 1.23489710160743e-38"
## [1] "Shapiro-Wilk test for NOX p-value: 2.21465209181895e-33"
## [1] "Shapiro-Wilk test for CO p-value: 1.37355309278081e-22"
## [1] "Shapiro-Wilk test for O3 p-value: 1.50415453463477e-27"
## [1] "Shapiro-Wilk test for SO2 p-value: 6.00675005349089e-35"#Combining all Q-Q plots for each pollutant
create_qq_plot <- function(data, pollutant) {
qq_plot <- ggplot(data, aes_string(sample = pollutant)) +
stat_qq() +
stat_qq_line() +
labs(title = paste("Q-Q Plot of", pollutant)) +
theme_minimal()
return(qq_plot)
}
qq_plots <- lapply(pollutants, function(p) create_qq_plot(data_cleaned, p))
combined_qq_plot <- do.call(grid.arrange, c(qq_plots, ncol = 2))
print(combined_qq_plot)## TableGrob (3 x 2) "arrange": 5 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
## 5 5 (3-3,1-1) arrange gtable[layout]#Pollutants are not normally distributed, according to the result.
#When dealing with non-normally distributed data like your PM2.5 measurements, it's essential to choose statistical methods and models that do not rely on the assumption of normality.
#Now, I am using Mann-Kendall trend test first as it is appropriate for this data.
# Function to perform Mann-Kendall test
perform_mann_kendall <- function(data, pollutant) {
# Extract the column for the pollutant
pollutant_data <- na.omit(data[[pollutant]]) # Remove NA values
# Conduct the Mann-Kendall test
mk_test <- MannKendall(pollutant_data)
return(mk_test)
}
# Apply the function to each pollutant
pollutants <- c("PM2.5", "NOX", "CO", "O3", "SO2")
mk_results <- lapply(pollutants, function(p) perform_mann_kendall(data_cleaned, p))
# Print the results
names(mk_results) <- pollutants
mk_results## $PM2.5
## tau = -0.0145, 2-sided pvalue =0.35736
##
## $NOX
## tau = -0.0488, 2-sided pvalue =0.0019756
##
## $CO
## tau = -0.0394, 2-sided pvalue =0.012659
##
## $O3
## tau = -0.0274, 2-sided pvalue =0.083273
##
## $SO2
## tau = 0.0386, 2-sided pvalue =0.014395#Plot the result
# Convert the Mann-Kendall test results into a data frame
mk_results_df <- data.frame(
Pollutant = c("PM2.5", "NOX", "CO", "O3", "SO2"),
Tau = c(-0.0145, -0.0488, -0.0394, -0.0274, 0.0386),
PValue = c(0.35736, 0.0019756, 0.012659, 0.083273, 0.014395)
)
# Add a column to indicate significance
mk_results_df$Significant <- mk_results_df$PValue < 0.05
# Create the plot
ggplot(mk_results_df, aes(x = Pollutant, y = Tau, fill = Significant)) +
geom_bar(stat = "identity", position = position_dodge()) +
scale_fill_manual(values = c("red", "blue"),
labels = c("Not Significant", "Significant"),
guide = guide_legend(title = "Trend Significance")) +
geom_text(aes(label = sprintf("tau=%.4f\np=%.4f", Tau, PValue), y = Tau + 0.01*sign(Tau)),
position = position_dodge(width = 0.9), size = 3, vjust = -0.5) +
labs(title = "Mann-Kendall Test Results for Pollutants",
x = "Pollutant",
y = "Tau Value") +
theme_minimal()
#VIF TEST
data_selected <- data_cleaned[c("PM2.5", "WS", "WD", "T", "RH", "SR", "BP", "R", "NOX", "SO2", "O3", "CO")]
# Fit a linear model with all selected variables
lm_model <- lm(PM2.5 ~ ., data = data_selected)
vif_results <- vif(lm_model)
print(vif_results)## WS WD T RH SR BP R NOX
## 1.538636 1.858846 2.763453 2.313589 1.921844 3.240407 1.083229 1.626086
## SO2 O3 CO
## 1.622984 1.207117 1.456148# plotting VIF Results
vif_df <- data.frame(Variable = names(vif_results), VIF = vif_results)
ggplot(vif_df, aes(x = Variable, y = VIF, group = 1)) +
geom_line(color = "blue") +
geom_point(color = "blue") +
geom_text(aes(label = sprintf("%.2f", VIF)), vjust = -0.5, hjust = 1.1) +
geom_hline(yintercept = 5, linetype = "dashed", color = "red") +
labs(title = "Variance Inflation Factor (VIF) for Each Variable",
x = "Variable",
y = "VIF Value") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
#MODEL IMPLEMENTATION
#I AM GOING TO USE RANDOM FOREST REGRESSION AND GENERALIZED LINEAR MODEL
# Set seed for reproducibility
set.seed(123)
# Splitting the data into training (80%) and testing (20%) sets
splitIndex <- createDataPartition(data_selected$PM2.5, p = 0.9, list = FALSE)
train_data <- data_selected[splitIndex, ]
test_data <- data_selected[-splitIndex, ]
# Random Forest
rf_model_train <- randomForest(PM2.5 ~ ., data = train_data, ntree = 500)
rf_predictions_test <- predict(rf_model_train, test_data)
# GLM
glm_model_train <- glm(PM2.5 ~ ., family = gaussian(), data = train_data)
glm_predictions_test <- predict(glm_model_train, test_data, type = "response")
# Calculate R-squared and RMSE for testing data
rf_r2_test <- cor(test_data$PM2.5, rf_predictions_test)^2
rf_rmse_test <- rmse(test_data$PM2.5, rf_predictions_test)
glm_r2_test <- cor(test_data$PM2.5, glm_predictions_test)^2
glm_rmse_test <- rmse(test_data$PM2.5, glm_predictions_test)
# Compare the results
comparison <- data.frame(Model = c("Random Forest", "GLM"),
R_Squared = c(rf_r2_test, glm_r2_test),
RMSE = c(rf_rmse_test, glm_rmse_test))
print(comparison)## Model R_Squared RMSE
## 1 Random Forest 0.8366860 26.38025
## 2 GLM 0.7232476 34.10853# Scatter plot for Random Forest results on test data
rf_test_plot <- ggplot(test_data, aes(x = PM2.5, y = rf_predictions_test)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
labs(title = "RFR: Observed vs Predicted (Test Data)",
x = "Observed PM2.5", y = "Predicted PM2.5") +
annotate("text", x = quantile(test_data$PM2.5, 0.1),
y = quantile(rf_predictions_test, 0.9),
label = sprintf("R² = %.2f\nRMSE = %.2f", rf_r2_test, rf_rmse_test),
hjust = 0, vjust = 1, size = 4, color = "blue") +
theme_minimal()
# Print the Random Forest plot
print(rf_test_plot)
# Scatter plot for GLM results on test data
glm_test_plot <- ggplot(test_data, aes(x = PM2.5, y = glm_predictions_test)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
labs(title = "GLM: Observed vs Predicted (Test Data)",
x = "Observed PM2.5", y = "Predicted PM2.5") +
annotate("text", x = quantile(test_data$PM2.5, 0.1),
y = quantile(glm_predictions_test, 0.9),
label = sprintf("R² = %.2f\nRMSE = %.2f", glm_r2_test, glm_rmse_test),
hjust = 0, vjust = 1, size = 4, color = "blue") +
theme_minimal()
print(glm_test_plot)
# Combine the two plots into one figure
combined_plot <- grid.arrange(rf_test_plot, glm_test_plot, ncol = 2)
# Print the combined plot
print(combined_plot)## TableGrob (1 x 2) "arrange": 2 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]