rm(list = ls())
# Set seed for reproducibility (important for random numbers)
set.seed(123)
# ============================================================================
# Task 1A: Generate Production Data
# ============================================================================
# Initialize an empty dataframe to store production data
production_data <- data.frame(Sol = integer(),
                              Station = integer(),
                              Team = integer(),
                              Output = numeric(),
                              stringsAsFactors = FALSE)
# Nested loops to populate the dataframe
# Loop through each Sol (day) 1 to 30
for(sol in 1:30) {
  # Loop through each Station 1 to 6
  for(station in 1:6) {
    # Loop through each Team 1 to 4
    for(team in 1:4) {
      # Calculate RandomNoise: random integer between -5 and +5
      random_noise <- sample(-5:5, 1)
      # Calculate Output using the formula: Output = BaseOutput + (Station × 3) + (Team × 2) + RandomNoise
      base_output <- 50
      output <- base_output + (station * 3) + (team * 2) + random_noise
      # Create a temporary dataframe for this row
      temp_row <- data.frame(Sol = sol,
                            Station = station,
                            Team = team,
                            Output = output)
      # Append the row to production_data
      production_data <- rbind(production_data, temp_row)
    }
  }
}
# Display the first 10 rows of the dataframe
cat("\n============================================================================\n")
cat("Task 1A: First 10 rows of production_data\n")
cat("============================================================================\n")
print(head(production_data, 10))
# Verify the structure and dimensions
cat("\nStructure of production_data:\n")
str(production_data)
cat("\nTotal rows (should be 30 × 6 × 4 = 720):", nrow(production_data), "\n")
# ============================================================================
# Task 1B: Total Production per Station
# ============================================================================
# Initialize an empty dataframe for station totals
station_totals <- data.frame(Station = integer(),
                            TotalOutput = numeric(),
                            stringsAsFactors = FALSE)
# Calculate total output for each station using loops
for(station in 1:6) {
  # Initialize sum for this station
  station_sum <- 0
  # Loop through all rows in production_data
  for(row in 1:nrow(production_data)) {
    # If the station matches, add its output to the sum
    if(production_data$Station[row] == station) {
      station_sum <- station_sum + production_data$Output[row]
    }
  }
  # Add the result to station_totals dataframe
  temp_row <- data.frame(Station = station, TotalOutput = station_sum)
  station_totals <- rbind(station_totals, temp_row)
}
# Display station_totals
cat("\n============================================================================\n")
cat("Task 1B: Total Production per Station\n")
cat("============================================================================\n")
print(station_totals)
# ============================================================================
# Task 1C: Detect Production Anomalies
# ============================================================================
# Initialize an empty dataframe for anomalies
anomalies <- data.frame(Sol = integer(),
                        Station = integer(),
                        Team = integer(),
                        Output = numeric(),
                        Anomaly = logical(),
                        stringsAsFactors = FALSE)
# Calculate daily average for each station and detect anomalies
# Loop through each Sol
for(sol in 1:30) {
  # Loop through each Station
  for(station in 1:6) {
    # Calculate sum of outputs for this station on this sol
    daily_station_sum <- 0
    team_count <- 0
    # Loop through all rows to find matching sol and station
    for(row in 1:nrow(production_data)) {
      if(production_data$Sol[row] == sol && production_data$Station[row] == station) {
        daily_station_sum <- daily_station_sum + production_data$Output[row]
        team_count <- team_count + 1
      }
    }
    # Calculate daily average (sum divided by number of teams, which should be 4)
    daily_avg <- daily_station_sum / team_count
    # Now check each team's output for this station and sol
    for(team in 1:4) {
      # Find the specific row for this team
      for(row in 1:nrow(production_data)) {
        if(production_data$Sol[row] == sol && 
           production_data$Station[row] == station && 
           production_data$Team[row] == team) {
          # Get the output
          team_output <- production_data$Output[row]
          # Calculate 10% threshold
          threshold <- daily_avg * 0.10
          # Check if deviation is more than 10%
          is_anomaly <- abs(team_output - daily_avg) > threshold
          # Add to anomalies dataframe
          temp_row <- data.frame(Sol = sol,
                                Station = station,
                                Team = team,
                                Output = team_output,
                                Anomaly = is_anomaly)
          anomalies <- rbind(anomalies, temp_row)
        }
      }
    }
  }
}
# Display the first 10 anomalies
cat("\n============================================================================\n")
cat("Task 1C: First 10 Detected Anomalies\n")
cat("============================================================================\n")
print(head(anomalies, 10))
# Summary of anomalies
cat("\nTotal anomalies detected:", sum(anomalies$Anomaly), "out of", nrow(anomalies), "records\n")
# ============================================================================
# Task 1D: Identify Peak Productivity Days
# ============================================================================
# Initialize an empty dataframe for peak days
peak_days <- data.frame(Station = integer(),
                        Sol = integer(),
                        TotalOutput = numeric(),
                        stringsAsFactors = FALSE)
# For each station, find the sol with maximum total output
for(station in 1:6) {
  # Initialize variables to track maximum
  max_output <- -Inf
  max_sol <- NA
  # Check each sol for this station
  for(sol in 1:30) {
    # Calculate total output for this station on this sol
    sol_total <- 0
    # Sum all teams' outputs for this station and sol
    for(row in 1:nrow(production_data)) {
      if(production_data$Station[row] == station && production_data$Sol[row] == sol) {
        sol_total <- sol_total + production_data$Output[row]
      }
    }
    # Update maximum if this sol has higher total
    if(sol_total > max_output) {
      max_output <- sol_total
      max_sol <- sol
    }
  }
  # Add the peak day for this station to the dataframe
  temp_row <- data.frame(Station = station,
                        Sol = max_sol,
                        TotalOutput = max_output)
  peak_days <- rbind(peak_days, temp_row)
}
# Display peak_days
cat("\n============================================================================\n")
cat("Task 1D: Peak Productivity Days per Station\n")
cat("============================================================================\n")
print(peak_days)
# ============================================================================
# Task 1E: Calculate Efficiency Score
# ============================================================================
# Add EfficiencyScore column to production_data
# First, create an empty vector for the new column
efficiency_scores <- numeric(nrow(production_data))
# Calculate EfficiencyScore for each row
for(row in 1:nrow(production_data)) {
  output <- production_data$Output[row]
  station <- production_data$Station[row]
  team <- production_data$Team[row]
  # Calculate EfficiencyScore = Output / (Station + Team)
  efficiency_scores[row] <- output / (station + team)
}
# Add the new column to production_data
production_data$EfficiencyScore <- efficiency_scores
# Display the first 10 rows of the updated dataframe
cat("\n============================================================================\n")
cat("Task 1E: First 10 rows with EfficiencyScore\n")
cat("============================================================================\n")
print(head(production_data, 10))
# ============================================================================
# Task 1F: Code Presentation
# ============================================================================
# The code above already includes:
# 1. Clear, informative comments explaining each step
# 2. Consistent indentation
# 3. Readable variable names
# 4. Proper formatting with section headers
cat("\n============================================================================\n")
cat("Task 1F: Code Presentation\n")
cat("============================================================================\n")
cat("
################################################################################
# Air Quality Analysis - UCI Dataset
# Author: Malik Haroon 
# Date: March 2026
# Description: Comprehensive analysis of air quality data including exploratory
#              data analysis, visualization, correlation analysis, and statistical tests
################################################################################
# FIRST, INSTALL THE REQUIRED PACKAGES (run this section once)
# Uncomment the lines below to install packages if you haven't already
install.packages("tidyverse")
install.packages("lubridate") 
install.packages("corrplot")
install.packages("gridExtra")
# NOW LOAD THE LIBRARIES
library(tidyverse)    # For data manipulation and visualization
library(lubridate)    # For date/time handling
library(corrplot)     # For correlation visualizations
library(gridExtra)    # For arranging multiple plots
################################################################################
# TASK 1: Data Import and Initial Exploration
################################################################################
# Import the dataset from your Downloads folder
AirQuality <- read.csv("C:/Users/Administrator/Downloads/air+quality/AirQualityUCI.csv", 
                       sep = ";",           # Dataset uses semicolon separator
                       dec = ",",            # Decimal separator is comma
                       na.strings = c("", "NA", "-200"))  # Handle missing values
# Check column names and structure
cat("\n=== TASK 1: DATA OVERVIEW ===\n")
cat("\nColumn Names:\n")
print(names(AirQuality))
cat("\nData Types:\n")
print(str(AirQuality))
# Display first few rows
cat("\nFirst 6 rows of the dataset:\n")
print(head(AirQuality, 6))
# Dataset dimensions
cat("\nDataset Dimensions:\n")
cat("Rows:", nrow(AirQuality), "\n")
cat("Columns:", ncol(AirQuality), "\n")
# Basic descriptive statistics for numeric columns
cat("\nDescriptive Statistics:\n")
# Select only numeric columns for summary
numeric_cols <- AirQuality[, sapply(AirQuality, is.numeric)]
print(summary(numeric_cols))
# Check for missing values
cat("\nMissing Values per Column:\n")
print(colSums(is.na(AirQuality)))
################################################################################
# TASK 2: Exploratory Visualization
################################################################################
cat("\n=== TASK 2: EXPLORATORY VISUALIZATION ===\n")
# Data preprocessing for time series analysis
# Remove the last two columns (X and X.1) which are all NA
AirQuality <- AirQuality[, 1:15]
# Create proper datetime column
AirQuality$DateTime <- as.POSIXct(paste(AirQuality$Date, AirQuality$Time), 
                                   format = "%d/%m/%Y %H.%M.%S")
# Remove rows with missing datetime or pollutant values
# First, convert -200 values (used as missing value indicator) to NA
AirQuality[AirQuality == -200] <- NA
# Create clean dataset
AirQuality_clean <- AirQuality %>%
  filter(!is.na(DateTime)) %>%
  filter(!is.na(CO.GT.), !is.na(C6H6.GT.), !is.na(T)) %>%
  filter(CO.GT. > 0, C6H6.GT. > 0, T > -100)  # Remove unrealistic values
# Time series plot for CO and Benzene (two pollutants)
# Justification: Line plots are ideal for showing trends over time
p1 <- ggplot(AirQuality_clean, aes(x = DateTime)) +
  geom_line(aes(y = CO.GT., color = "CO"), alpha = 0.7) +
  geom_line(aes(y = C6H6.GT., color = "Benzene"), alpha = 0.7) +
  labs(title = "Time Trends: CO and Benzene Levels",
       x = "Date/Time",
       y = "Concentration",
       color = "Pollutant") +
  theme_minimal() +
  theme(legend.position = "bottom")
# Scatterplot: Temperature vs CO
# Justification: Scatterplots reveal relationships between continuous variables
p2 <- ggplot(AirQuality_clean, aes(x = T, y = CO.GT.)) +
  geom_point(aes(color = RH), alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  labs(title = "Temperature vs CO Levels",
       x = "Temperature (°C)",
       y = "CO (mg/m³)",
       color = "Relative Humidity (%)") +
  theme_minimal()
# Boxplot: Wind direction categories vs CO sensor response
# Create wind direction categories (if WD column exists)
if("WD" %in% names(AirQuality_clean)) {
  AirQuality_clean$WindDir_Cat <- cut(AirQuality_clean$WD,
                                       breaks = c(0, 90, 180, 270, 360),
                                       labels = c("NE", "SE", "SW", "NW"),
                                       include.lowest = TRUE)
  # Justification: Boxplots show distribution differences across categories
  p3 <- ggplot(AirQuality_clean %>% filter(!is.na(WindDir_Cat)), 
               aes(x = WindDir_Cat, y = PT08.S1.CO.)) +
    geom_boxplot(fill = "skyblue", outlier.color = "red", outlier.alpha = 0.3) +
    labs(title = "CO Sensor Response by Wind Direction",
         x = "Wind Direction",
         y = "CO Sensor Response") +
    theme_minimal()
} else {
  # If WD column doesn't exist, create a different plot
  p3 <- ggplot(AirQuality_clean, aes(x = factor(floor(RH/10)*10), y = CO.GT.)) +
    geom_boxplot(fill = "skyblue", outlier.color = "red", outlier.alpha = 0.3) +
    labs(title = "CO Levels by Humidity Range",
         x = "Humidity Range (%)",
         y = "CO (mg/m³)") +
    theme_minimal()
}
# Display plots
grid.arrange(p1, p2, p3, ncol = 2, nrow = 2)
################################################################################
# TASK 3: Correlation Analysis
################################################################################
cat("\n=== TASK 3: CORRELATION ANALYSIS ===\n")
# Select pollutants for correlation analysis
pollutants <- AirQuality_clean %>%
  select(CO.GT., PT08.S1.CO., C6H6.GT., PT08.S2.NMHC., NOx.GT., PT08.S3.NOx.)
# Calculate correlation matrix
cor_matrix <- cor(pollutants, use = "complete.obs")
# Display correlation matrix
cat("\nCorrelation Matrix:\n")
print(round(cor_matrix, 3))
# Visualize correlations using heatmap
# Justification: Heatmaps provide an intuitive overview of correlation patterns
par(mfrow = c(1, 2))
# Heatmap
corrplot(cor_matrix, method = "color", type = "upper", 
         tl.col = "black", tl.srt = 45,
         title = "Pollutant Correlations",
         mar = c(0, 0, 2, 0))
# Scatterplot matrix for selected pollutants
# Justification: Shows detailed relationships between variables
pairs(pollutants[, 1:4], 
      pch = 19, 
      col = rgb(0, 0, 1, 0.3),
      main = "Scatterplot Matrix: Pollutant Relationships")
cat("\nInterpretation:\n")
cat("- Strong positive correlations (>0.7) indicate closely related pollutants\n")
cat("- Negative correlations suggest inverse relationships\n")
cat("- The sensor responses (PT08 series) show strong correlations with\n")
cat("  corresponding gas concentrations, validating measurement accuracy\n")
################################################################################
# TASK 4: Hypothesis Testing - Two-Group Comparison
################################################################################
cat("\n=== TASK 4: HYPOTHESIS TESTING (TWO-GROUP COMPARISON) ===\n")
# Compare CO levels based on temperature (high vs low)
# Create temperature categories (median split)
median_temp <- median(AirQuality_clean$T, na.rm = TRUE)
AirQuality_clean$Temp_Category <- ifelse(AirQuality_clean$T > median_temp, 
                                          "High Temperature", "Low Temperature")
# Extract CO levels for each group
co_low_temp <- AirQuality_clean$CO.GT.[AirQuality_clean$Temp_Category == "Low Temperature"]
co_high_temp <- AirQuality_clean$CO.GT.[AirQuality_clean$Temp_Category == "High Temperature"]
# Remove NA values
co_low_temp <- co_low_temp[!is.na(co_low_temp)]
co_high_temp <- co_high_temp[!is.na(co_high_temp)]
# Check normality assumption (using sample due to large dataset)
set.seed(123)  # for reproducibility
shapiro_low <- shapiro.test(sample(co_low_temp, min(500, length(co_low_temp))))
shapiro_high <- shapiro.test(sample(co_high_temp, min(500, length(co_high_temp))))
cat("\nNormality Test Results:\n")
cat("Low Temperature p-value:", shapiro_low$p.value, "\n")
cat("High Temperature p-value:", shapiro_high$p.value, "\n")
# Since data may not be normal, use Wilcoxon test
# Hypotheses:
# H0: There is no difference in CO levels between high and low temperature conditions
# H1: CO levels differ between high and low temperature conditions
wilcox_result <- wilcox.test(co_low_temp, co_high_temp)
cat("\nWilcoxon Rank Sum Test Results:\n")
cat("W statistic:", wilcox_result$statistic, "\n")
cat("p-value:", wilcox_result$p.value, "\n")
# Visualize the comparison
p4 <- ggplot(AirQuality_clean %>% filter(!is.na(Temp_Category)), 
             aes(x = Temp_Category, y = CO.GT., fill = Temp_Category)) +
  geom_boxplot() +
  labs(title = "CO Levels: High vs Low Temperature",
       subtitle = paste("Wilcoxon test p-value =", round(wilcox_result$p.value, 4)),
       x = "Temperature Condition",
       y = "CO (mg/m³)") +
  theme_minimal() +
  theme(legend.position = "none")
print(p4)
cat("\nInterpretation:\n")
if(wilcox_result$p.value < 0.05) {
  cat("Reject H0: Significant difference in CO levels between high and low temperature conditions.\n")
  cat("This suggests temperature affects pollutant levels.\n")
} else {
  cat("Fail to reject H0: No significant difference in CO levels between temperature conditions.\n")
}
################################################################################
# TASK 5: Group Comparison Across Seasons
################################################################################
cat("\n=== TASK 5: SEASONAL COMPARISON ===\n")
# Create season variable from dates
AirQuality_clean$Date <- dmy(AirQuality_clean$Date)
AirQuality_clean$Month <- month(AirQuality_clean$Date)
AirQuality_clean$Season <- case_when(
  AirQuality_clean$Month %in% c(12, 1, 2) ~ "Winter",
  AirQuality_clean$Month %in% c(3, 4, 5) ~ "Spring",
  AirQuality_clean$Month %in% c(6, 7, 8) ~ "Summer",
  AirQuality_clean$Month %in% c(9, 10, 11) ~ "Autumn"
)
# Select pollutant for analysis (using Benzene - C6H6.GT.)
# Remove NA values
seasonal_data <- AirQuality_clean %>%
  filter(!is.na(Season), !is.na(C6H6.GT.))
# Boxplot for seasonal comparison
# Justification: Boxplots effectively show distribution differences across groups
p5 <- ggplot(seasonal_data, aes(x = Season, y = C6H6.GT., fill = Season)) +
  geom_boxplot() +
  labs(title = "Benzene Levels Across Seasons",
       x = "Season",
       y = "Benzene (µg/m³)") +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_fill_manual(values = c("Winter" = "lightblue", 
                                "Spring" = "lightgreen",
                                "Summer" = "gold", 
                                "Autumn" = "orange"))
print(p5)
# Statistical test: ANOVA for seasonal differences
# Hypotheses:
# H0: Mean benzene levels are equal across all seasons
# H1: At least one season has different mean benzene levels
# Test for homogeneity of variances
bartlett_test <- bartlett.test(C6H6.GT. ~ Season, data = seasonal_data)
cat("\nBartlett Test for Homogeneity of Variances:\n")
cat("p-value:", bartlett_test$p.value, "\n")
# Perform ANOVA
anova_result <- aov(C6H6.GT. ~ Season, data = seasonal_data)
cat("\nANOVA Results:\n")
print(summary(anova_result))
# If ANOVA significant, perform post-hoc test
if(summary(anova_result)[[1]][["Pr(>F)"]][1] < 0.05) {
  cat("\nPost-hoc Test (Tukey HSD):\n")
  tukey_result <- TukeyHSD(anova_result)
  print(tukey_result)
}
# Check normality of residuals
residuals <- residuals(anova_result)
shapiro_resid <- shapiro.test(sample(residuals, min(500, length(residuals))))
cat("\nNormality of Residuals (Shapiro-Wilk test):\n")
cat("p-value:", shapiro_resid$p.value, "\n")
cat("\nInterpretation:\n")
if(summary(anova_result)[[1]][["Pr(>F)"]][1] < 0.05) {
  cat("Reject H0: Significant differences in benzene levels across seasons.\n")
  cat("This suggests seasonal factors (e.g., temperature, sunlight) affect pollutant levels.\n")
} else {
  cat("Fail to reject H0: No significant seasonal differences in benzene levels.\n")
}
################################################################################
# TASK 6: Summary and Code Quality
################################################################################
cat("\n=== TASK 6: ANALYSIS SUMMARY ===\n")
cat("\nAnalysis completed successfully!\n")
cat("Key findings:\n")
cat("1. Dataset contains", nrow(AirQuality_clean), "valid observations\n")
cat("2. Strong correlations observed between related pollutants\n")
cat("3. Temperature shows", ifelse(wilcox_result$p.value < 0.05, 
                                    "significant", "no significant"), 
    "effect on CO levels\n")
cat("4. Seasonal patterns in benzene levels are", 
    ifelse(summary(anova_result)[[1]][["Pr(>F)"]][1] < 0.05,
           "significant", "not significant"), "\n")
# Clean up (optional)
save.image("C:\\Ehtisham\\Haroon\\New\\2\\Code")