# Cronbach's Alpha Reliability Analysis
# For 4 domains (preferences, personality, values, life satisfaction) at 2 time points
# 5 items per domain per time point

# Load required libraries
library(psych)
library(dplyr)
library(tidyr)

# Read the data
data <- read.csv("exp1.csv")

# Define the scale variables for each domain and time point
# Past time point scales
past_pref_vars <- c("NPastDiff_pref_read", "NPastDiff_pref_music", "NPastDiff_pref_tv", 
                   "NPastDiff_pref_nap", "NPastDiff_pref_travel")

past_pers_vars <- c("NPastDiff_pers_extravert", "NPastDiff_pers_critical", "NPastDiff_pers_dependable", 
                   "NPastDiff_pers_anxious", "NPastDiff_pers_complex")

past_val_vars <- c("NPastDiff_val_obey", "NPastDiff_val_trad", "NPastDiff_val_opinion", 
                  "NPastDiff_val_performance", "NPastDiff_val_justice")

past_life_vars <- c("NPastDiff_life_ideal", "NPastDiff_life_excellent", "NPastDiff_life_satisfied", 
                   "NPastDiff_life_important", "NPastDiff_life_change")

# Future time point scales
fut_pref_vars <- c("NFutDiff_pref_read", "NFutDiff_pref_music", "NFutDiff_pref_tv", 
                  "NFutDiff_pref_nap", "NFutDiff_pref_travel")

fut_pers_vars <- c("NFutDiff_pers_extravert", "NFutDiff_pers_critical", "NFutDiff_pers_dependable", 
                  "NFutDiff_pers_anxious", "NFutDiff_pers_complex")

fut_val_vars <- c("NFutDiff_val_obey", "NFutDiff_val_trad", "NFutDiff_val_opinion", 
                 "NFutDiff_val_performance", "NFutDiff_val_justice")

fut_life_vars <- c("NFutDiff_life_ideal", "NFutDiff_life_excellent", "NFutDiff_life_satisfied", 
                  "NFutDiff_life_important", "NFutDiff_life_change")

# Function to calculate Cronbach's alpha and return detailed results
calc_cronbach_alpha <- function(data, var_names, scale_name) {
  # Check for missing values
  scale_data <- data[, var_names]
  missing_info <- data.frame(
    Variable = var_names,
    Missing_Count = colSums(is.na(scale_data)),
    Missing_Percent = round(colSums(is.na(scale_data)) / nrow(scale_data) * 100, 2)
  )
  
  # Remove rows with any missing values for reliability analysis
  complete_data <- scale_data[complete.cases(scale_data), ]
  
  cat("\n", "="*60, "\n")
  cat("SCALE:", scale_name, "\n")
  cat("="*60, "\n")
  
  cat("Sample size for reliability analysis:", nrow(complete_data), "\n")
  cat("Original sample size:", nrow(data), "\n")
  cat("Cases removed due to missing data:", nrow(data) - nrow(complete_data), "\n\n")
  
  cat("Missing data summary:\n")
  print(missing_info)
  
  if(nrow(complete_data) < 3) {
    cat("\nWARNING: Insufficient complete cases for reliability analysis\n")
    return(NULL)
  }
  
  # Calculate Cronbach's alpha
  alpha_result <- alpha(complete_data, check.keys = TRUE)
  
  cat("\nCronbach's Alpha Results:\n")
  cat("Raw alpha:", round(alpha_result$total$raw_alpha, 4), "\n")
  cat("Standardized alpha:", round(alpha_result$total$std.alpha, 4), "\n")
  cat("Average inter-item correlation:", round(alpha_result$total$average_r, 4), "\n")
  cat("Number of items:", alpha_result$total$nvar, "\n")
  
  # Item statistics
  cat("\nItem Statistics:\n")
  item_stats <- data.frame(
    Item = var_names,
    Alpha_if_deleted = round(alpha_result$alpha.drop$raw_alpha, 4),
    Item_total_correlation = round(alpha_result$item.stats$r.drop, 4),
    Mean = round(alpha_result$item.stats$mean, 4),
    SD = round(alpha_result$item.stats$sd, 4)
  )
  print(item_stats)
  
  # Check assumptions
  cat("\nAssumption Checks:\n")
  
  # 1. Check for sufficient sample size (minimum 30 recommended)
  sample_size_ok <- nrow(complete_data) >= 30
  cat("Sample size adequate (≥30):", sample_size_ok, "\n")
  
  # 2. Check for adequate inter-item correlations (should be > 0.30)
  inter_item_cors <- cor(complete_data)
  inter_item_cors[lower.tri(inter_item_cors)] <- NA
  diag(inter_item_cors) <- NA
  inter_item_cors_flat <- as.vector(inter_item_cors)
  inter_item_cors_flat <- inter_item_cors_flat[!is.na(inter_item_cors_flat)]
  adequate_cors <- sum(inter_item_cors_flat > 0.30) / length(inter_item_cors_flat)
  cat("Proportion of inter-item correlations > 0.30:", round(adequate_cors, 4), "\n")
  
  # 3. Check for negative correlations (concerning for unidimensionality)
  negative_cors <- sum(inter_item_cors_flat < 0, na.rm = TRUE)
  cat("Number of negative inter-item correlations:", negative_cors, "\n")
  
  # 4. Check item variances (should be roughly similar)
  item_vars <- apply(complete_data, 2, var)
  var_ratio <- max(item_vars) / min(item_vars)
  cat("Ratio of highest to lowest item variance:", round(var_ratio, 4), "\n")
  
  return(alpha_result)
}

# Calculate Cronbach's alpha for all scales
cat("CRONBACH'S ALPHA RELIABILITY ANALYSIS")
cat("\nData: exp1.csv")
cat("\nTotal sample size:", nrow(data))

# Past time point analyses
past_pref_alpha <- calc_cronbach_alpha(data, past_pref_vars, "Past Preferences")
past_pers_alpha <- calc_cronbach_alpha(data, past_pers_vars, "Past Personality")
past_val_alpha <- calc_cronbach_alpha(data, past_val_vars, "Past Values")
past_life_alpha <- calc_cronbach_alpha(data, past_life_vars, "Past Life Satisfaction")

# Future time point analyses
fut_pref_alpha <- calc_cronbach_alpha(data, fut_pref_vars, "Future Preferences")
fut_pers_alpha <- calc_cronbach_alpha(data, fut_pers_vars, "Future Personality")
fut_val_alpha <- calc_cronbach_alpha(data, fut_val_vars, "Future Values")
fut_life_alpha <- calc_cronbach_alpha(data, fut_life_vars, "Future Life Satisfaction")

# Summary table
cat("\n", "="*80, "\n")
cat("SUMMARY OF CRONBACH'S ALPHA COEFFICIENTS")
cat("\n", "="*80, "\n")

summary_results <- data.frame(
  Scale = c("Past Preferences", "Past Personality", "Past Values", "Past Life Satisfaction",
           "Future Preferences", "Future Personality", "Future Values", "Future Life Satisfaction"),
  Alpha = c(
    if(!is.null(past_pref_alpha)) round(past_pref_alpha$total$raw_alpha, 4) else NA,
    if(!is.null(past_pers_alpha)) round(past_pers_alpha$total$raw_alpha, 4) else NA,
    if(!is.null(past_val_alpha)) round(past_val_alpha$total$raw_alpha, 4) else NA,
    if(!is.null(past_life_alpha)) round(past_life_alpha$total$raw_alpha, 4) else NA,
    if(!is.null(fut_pref_alpha)) round(fut_pref_alpha$total$raw_alpha, 4) else NA,
    if(!is.null(fut_pers_alpha)) round(fut_pers_alpha$total$raw_alpha, 4) else NA,
    if(!is.null(fut_val_alpha)) round(fut_val_alpha$total$raw_alpha, 4) else NA,
    if(!is.null(fut_life_alpha)) round(fut_life_alpha$total$raw_alpha, 4) else NA
  ),
  Items = rep(5, 8),
  Interpretation = c(
    if(!is.null(past_pref_alpha)) {
      alpha_val <- past_pref_alpha$total$raw_alpha
      if(alpha_val >= 0.90) "Excellent" else if(alpha_val >= 0.80) "Good" else if(alpha_val >= 0.70) "Acceptable" else if(alpha_val >= 0.60) "Questionable" else "Poor"
    } else "Insufficient data",
    if(!is.null(past_pers_alpha)) {
      alpha_val <- past_pers_alpha$total$raw_alpha
      if(alpha_val >= 0.90) "Excellent" else if(alpha_val >= 0.80) "Good" else if(alpha_val >= 0.70) "Acceptable" else if(alpha_val >= 0.60) "Questionable" else "Poor"
    } else "Insufficient data",
    if(!is.null(past_val_alpha)) {
      alpha_val <- past_val_alpha$total$raw_alpha
      if(alpha_val >= 0.90) "Excellent" else if(alpha_val >= 0.80) "Good" else if(alpha_val >= 0.70) "Acceptable" else if(alpha_val >= 0.60) "Questionable" else "Poor"
    } else "Insufficient data",
    if(!is.null(past_life_alpha)) {
      alpha_val <- past_life_alpha$total$raw_alpha
      if(alpha_val >= 0.90) "Excellent" else if(alpha_val >= 0.80) "Good" else if(alpha_val >= 0.70) "Acceptable" else if(alpha_val >= 0.60) "Questionable" else "Poor"
    } else "Insufficient data",
    if(!is.null(fut_pref_alpha)) {
      alpha_val <- fut_pref_alpha$total$raw_alpha
      if(alpha_val >= 0.90) "Excellent" else if(alpha_val >= 0.80) "Good" else if(alpha_val >= 0.70) "Acceptable" else if(alpha_val >= 0.60) "Questionable" else "Poor"
    } else "Insufficient data",
    if(!is.null(fut_pers_alpha)) {
      alpha_val <- fut_pers_alpha$total$raw_alpha
      if(alpha_val >= 0.90) "Excellent" else if(alpha_val >= 0.80) "Good" else if(alpha_val >= 0.70) "Acceptable" else if(alpha_val >= 0.60) "Questionable" else "Poor"
    } else "Insufficient data",
    if(!is.null(fut_val_alpha)) {
      alpha_val <- fut_val_alpha$total$raw_alpha
      if(alpha_val >= 0.90) "Excellent" else if(alpha_val >= 0.80) "Good" else if(alpha_val >= 0.70) "Acceptable" else if(alpha_val >= 0.60) "Questionable" else "Poor"
    } else "Insufficient data",
    if(!is.null(fut_life_alpha)) {
      alpha_val <- fut_life_alpha$total$raw_alpha
      if(alpha_val >= 0.90) "Excellent" else if(alpha_val >= 0.80) "Good" else if(alpha_val >= 0.70) "Acceptable" else if(alpha_val >= 0.60) "Questionable" else "Poor"
    } else "Insufficient data"
  )
)

print(summary_results)

# Save results to CSV
write.csv(summary_results, "cronbach_alpha_summary.csv", row.names = FALSE)
cat("\nSummary results saved to: cronbach_alpha_summary.csv\n")

cat("\n", "="*80, "\n")
cat("INTERPRETATION GUIDE FOR CRONBACH'S ALPHA")
cat("\n", "="*80, "\n")
cat("α ≥ 0.90: Excellent reliability\n")
cat("α ≥ 0.80: Good reliability\n") 
cat("α ≥ 0.70: Acceptable reliability\n")
cat("α ≥ 0.60: Questionable reliability\n")
cat("α < 0.60: Poor reliability\n")