# ============================================================================
# DATA SCREENING
# ============================================================================
# Results saved to: C:/Users/AWAN/Documents/Data_Screening_Results.xlsx
# ============================================================================

# STEP 0: Install packages (run ONCE only)
# Remove # symbols, highlight, press Ctrl+Enter, then put # back
# install.packages("tidyverse")
# install.packages("psych")
# install.packages("car")
# install.packages("mice")
# install.packages("openxlsx")

# STEP 1: Load packages
library(tidyverse)
library(psych)
library(car)
library(mice)
library(openxlsx)

cat("All packages loaded!\n\n")

# STEP 2: Set working directory
setwd("C:/Users/AWAN/Documents")
cat("Working directory:", getwd(), "\n\n")

# STEP 3: Load your data
# CHANGE "YOUR_FILE_NAME.csv" TO YOUR ACTUAL FILENAME
df <- read_csv("Survey-dataset.csv")

cat("Data loaded! Rows:", nrow(df), "| Columns:", ncol(df), "\n\n")

# STEP 4: Define variables
continuous_vars <- c("Age", "Education", "Income", "Length_Stay", 
                     "CEA_mean", "SCT_mean", "EC_mean", "WTP_Premium")

# ============================================================================
# STEP 5: RUN ALL SCREENING CHECKS
# ============================================================================

# --- 5a. Missing Data ---
missing_counts <- colSums(is.na(df))
missing_pct <- round(colSums(is.na(df)) / nrow(df) * 100, 2)

missing_table <- data.frame(
  Variable = names(missing_counts),
  Missing_Count = as.numeric(missing_counts),
  Missing_Percent = as.numeric(missing_pct)
)

original_n <- nrow(df)
df_clean <- df %>% drop_na()
final_n <- nrow(df_clean)
cases_removed <- original_n - final_n
attrition_pct <- round(cases_removed / original_n * 100, 2)

# --- 5b. Outliers ---
z_scores <- as.data.frame(lapply(df_clean[, continuous_vars], scale))
univariate_outliers <- sum(apply(abs(z_scores) > 3.29, 1, any))

pred_matrix <- as.matrix(df_clean[, continuous_vars[-length(continuous_vars)]])
maha_dist <- mahalanobis(pred_matrix, colMeans(pred_matrix), cov(pred_matrix))
critical_chisq <- qchisq(0.999, df = ncol(pred_matrix))
multivariate_outliers <- sum(maha_dist > critical_chisq)

# --- 5c. Normality ---
desc <- describe(df_clean[, continuous_vars], fast = TRUE)

normality_table <- data.frame(
  Variable = continuous_vars,
  N = desc$n,
  Mean = round(desc$mean, 2),
  SD = round(desc$sd, 2),
  Min = round(desc$min, 2),
  Max = round(desc$max, 2),
  Skewness = round(desc$skew, 3),
  Kurtosis = round(desc$kurtosis, 3),
  Skew_OK = abs(desc$skew) < 2.0,
  Kurt_OK = abs(desc$kurtosis) < 7.0
)

# --- 5d. Correlations ---
corr_matrix <- cor(df_clean[, continuous_vars], method = "pearson")
corr_df <- as.data.frame(round(corr_matrix, 3))
corr_df$Variable <- rownames(corr_df)
corr_df <- corr_df[, c("Variable", setdiff(names(corr_df), "Variable"))]

predictor_cors <- corr_matrix[1:7, 1:7]
max_cor <- max(abs(predictor_cors[upper.tri(predictor_cors)]))

# --- 5e. VIF ---
df_clean$CEA_c <- df_clean$CEA_mean - mean(df_clean$CEA_mean)
df_clean$SCT_c <- df_clean$SCT_mean - mean(df_clean$SCT_mean)

model3 <- lm(WTP_Premium ~ Age + Gender + Education + Income +
               Prior_Green_Exp + Trip_Purpose + Length_Stay +
               CEA_c + SCT_c + EC_mean + CEA_c:SCT_c, data = df_clean)

vif_values <- vif(model3)
vif_table <- data.frame(
  Variable = names(vif_values),
  VIF = round(as.numeric(vif_values), 4),
  Tolerance = round(1 / as.numeric(vif_values), 4)
)

# ============================================================================
# STEP 6: CREATE EXCEL WORKBOOK
# ============================================================================

wb <- createWorkbook()

# Sheet 1: Summary
summary_text <- data.frame(
  Check = c("Original Sample Size", 
            "Final Sample Size (after listwise deletion)",
            "Cases Removed",
            "Attrition Rate (%)",
            "",
            "Univariate Outliers (|z| > 3.29)",
            "Multivariate Outliers (Mahalanobis)",
            "",
            "All |Skewness| < 2.0",
            "All |Kurtosis| < 7.0",
            "",
            "Highest |r| among predictors",
            "",
            "Highest VIF",
            "All VIF < 5.0"),
  Result = c(as.character(original_n),
             as.character(final_n),
             as.character(cases_removed),
             paste0(attrition_pct, "%"),
             "",
             as.character(univariate_outliers),
             as.character(multivariate_outliers),
             "",
             as.character(all(abs(desc$skew) < 2.0)),
             as.character(all(abs(desc$kurtosis) < 7.0)),
             "",
             as.character(round(max_cor, 3)),
             "",
             as.character(round(max(vif_values), 4)),
             as.character(all(vif_values < 5.0)))
)

addWorksheet(wb, "Summary")
writeData(wb, "Summary", summary_text, startRow = 1, startCol = 1)

# Sheet 2: Missing Data
addWorksheet(wb, "Missing Data")
writeData(wb, "Missing Data", missing_table, startRow = 1, startCol = 1)

# Sheet 3: Descriptives & Normality
addWorksheet(wb, "Descriptives_Normality")
writeData(wb, "Descriptives_Normality", normality_table, startRow = 1, startCol = 1)

# Sheet 4: Correlations
addWorksheet(wb, "Correlations")
writeData(wb, "Correlations", corr_df, startRow = 1, startCol = 1)

# Sheet 5: VIF
addWorksheet(wb, "VIF")
writeData(wb, "VIF", vif_table, startRow = 1, startCol = 1)

# Save
output_file <- "Data_Screening_Results.xlsx"
saveWorkbook(wb, output_file, overwrite = TRUE)

cat("\n=================================================================\n")
cat("EXCEL FILE CREATED SUCCESSFULLY!\n")
cat("=================================================================\n")
cat("File name:", output_file, "\n")
cat("Saved to:", getwd(), "\n")
cat("\nSheets included:\n")
cat("  1. Summary - Overview of all checks\n")
cat("  2. Missing Data - Missing counts and percentages\n")
cat("  3. Descriptives_Normality - Means, SDs, skewness, kurtosis\n")
cat("  4. Correlations - Full correlation matrix\n")
cat("  5. VIF - Variance Inflation Factors\n")
cat("=================================================================\n")