library(dplyr)

library(outliers)
library(ggplot2)

df <- read.csv("C:/Users/User/Documents/SEMESTER 6/Analisis Regresi/datania1k.csv")
data <- df[c("Age", "Annual_Premium")]

total_missing <- sum(is.na(data))
total_missing

detect_outliers <- function(data) {
  cat("\n=== DETEKSI OUTLIER ===\n\n")
  
  outlier_report <- list()
  
  for(var in names(data)[sapply(data, is.numeric)]) {
    cat("Analisis outlier untuk variabel:", var, "\n")
    
    # Statistik deskriptif
    stats <- summary(data[[var]])
    iqr_val <- IQR(data[[var]], na.rm = TRUE)
    q1 <- quantile(data[[var]], 0.25, na.rm = TRUE)
    q3 <- quantile(data[[var]], 0.75, na.rm = TRUE)
    lower_bound <- q1 - 1.5 * iqr_val
    upper_bound <- q3 + 1.5 * iqr_val
    
    # Deteksi outlier dengan metode IQR
    outliers_iqr <- data[[var]][data[[var]] < lower_bound | data[[var]] > upper_bound]
    
    # Deteksi outlier dengan metode Z-score
    z_scores <- scale(data[[var]])
    outliers_z <- data[[var]][abs(z_scores) > 3]
    
    # Deteksi outlier dengan metode Grubbs (uji statistik)
    if(length(na.omit(data[[var]])) > 6) {
      tryCatch({
        grubbs_test <- grubbs.test(na.omit(data[[var]]))
        grubbs_outlier <- ifelse(grubbs_test$p.value < 0.05, "Terdeteksi", "Tidak terdeteksi")
      }, error = function(e) {
        grubbs_outlier <- "Tidak dapat dihitung"
      })
    } else {
      grubbs_outlier <- "Data tidak cukup"
    }
    
    # Ringkasan
    outlier_report[[var]] <- list(
      n_outliers_iqr = length(outliers_iqr),
      n_outliers_z = length(outliers_z),
      grubbs_result = grubbs_outlier,
      lower_bound = lower_bound,
      upper_bound = upper_bound,
      outlier_values = unique(round(outliers_iqr, 2))
    )
    
    cat("  - Outlier (IQR method):", length(outliers_iqr), "\n")
    cat("  - Outlier (Z-score > 3):", length(outliers_z), "\n")
    cat("  - Uji Grubbs:", grubbs_outlier, "\n")
    cat("  - Batas bawah:", round(lower_bound, 2), "\n")
    cat("  - Batas atas:", round(upper_bound, 2), "\n")
    
    if(length(outliers_iqr) > 0) {
      cat("  - Nilai outlier:", paste(head(unique(round(outliers_iqr, 2)), 5), collapse = ", "), "\n")
    }
    cat("\n")
  }
  
  return(outlier_report)
}

# Deteksi outlier
outlier_analysis <- detect_outliers(data)

# Visual
visualize_outliers <- function(data) {
  
  numeric_vars <- names(data)[sapply(data, is.numeric)]
  
  # Boxplot untuk setiap variabel numerik
  par(mfrow = c(2, 3))
  for(var in numeric_vars[1:min(6, length(numeric_vars))]) {
    boxplot(data[[var]], main = var, col = "lightblue", 
            ylab = "Nilai", outline = TRUE)
    grid()
  }
  par(mfrow = c(1, 1))
  
  # Histogram dengan overlay outlier
  for(var in numeric_vars[1:min(3, length(numeric_vars))]) {
    # Hitung batas outlier
    q1 <- quantile(data[[var]], 0.25, na.rm = TRUE)
    q3 <- quantile(data[[var]], 0.75, na.rm = TRUE)
    iqr_val <- IQR(data[[var]], na.rm = TRUE)
    lower_bound <- q1 - 1.5 * iqr_val
    upper_bound <- q3 + 1.5 * iqr_val
    
    # Identifikasi outlier
    is_outlier <- data[[var]] < lower_bound | data[[var]] > upper_bound
    
    # Plot histogram
    hist_data <- ggplot(data.frame(value = data[[var]]), aes(x = value)) +
      geom_histogram(aes(y = ..density..), bins = 30, fill = "lightblue", alpha = 0.7) +
      geom_density(color = "darkblue", linewidth = 1) +
      geom_vline(xintercept = c(lower_bound, upper_bound), 
                 color = "red", linetype = "dashed", linewidth = 1) +
      labs(title = paste("Distribusi dan Outlier:", var),
           x = var, y = "Density") +
      theme_minimal() +
      annotate("text", x = lower_bound, y = 0, 
               label = "Bawah", vjust = 2, color = "red") +
      annotate("text", x = upper_bound, y = 0, 
               label = "Atas", vjust = 2, color = "red")
    
    print(hist_data)
  }
  
  # Scatter plot matrix untuk melihat outlier multivariat
  if(length(numeric_vars) >= 3) {
    pairs(data[, numeric_vars[1:min(4, length(numeric_vars))]], 
          main = "Scatter Plot Matrix untuk Deteksi Outlier",
          pch = 19, col = alpha("blue", 0.6))
  }
}

# Jalankan visualisasi
visualize_outliers(data)

cat("\n==== STATISTIK DESKRIPTIF ====\n")

summary(data)

# Korelasi
cat("\nKORELASI:\n")

correlation <- cor(data$Age, data$Annual_Premium)
cat("Koefisien Korelasi (r) =", round(correlation, 4), "\n")

if (abs(correlation) >= 0.9) {
  cat("Interpretasi: Korelasi sangat kuat\n")
} else if (abs(correlation) >= 0.7) {
  cat("Interpretasi: Korelasi kuat\n")
} else if (abs(correlation) >= 0.5) {
  cat("Interpretasi: Korelasi sedang\n")
} else if (abs(correlation) >= 0.3) {
  cat("Interpretasi: Korelasi lemah\n")
} else {
  cat("Interpretasi: Korelasi sangat lemah/tidak ada\n")
}

# Visual Data
cat("\n\n=== VISUALISASI DATA ===\n")

# Scatter plot
par(mfrow = c(1, 2)) # Dua plot berdampingan

# Plot 1: Scatter plot tanpa garis regresi
plot(data$Age, data$Annual_Premium,
     main = "Scatter Plot: Age vs Annual Premium",
     xlab = "Age (X)",
     ylab = "Annual Premium (Y)",
     pch = 19, col = "blue", cex = 1.5)
grid()

# Plot 2: Scatter plot dengan garis regresi
plot(data$Age, data$Annual_Premium,
     main = "Data dengan Garis Regresi",
     xlab = "Age (X)",
     ylab = "Annual Premium (Y)",
     pch = 19, col = "blue", cex = 1.5)
grid()

# Kembalikan ke satu plot per window
par(mfrow = c(1, 1))

# Model regresi linear sederhana
model <- lm(Age ~ Annual_Premium, data = data)

# Persamaan regresi
cat("\nPERSAMAAN REGRESI LINEAR:\n")

beta0 <- coef(model)[1]
beta1 <- coef(model)[2]

# Ringkasan model
cat("\nRINGKASAN MODEL:\n")

summary_model <- summary(model)
print(summary_model)

# Ekstrak komponen penting
cat("\nPERSAMAAN REGRESI:\n")

intercept <- coef(model)[1]
slope <- coef(model)[2]
cat("Persamaan Regresi: Y =", round(intercept, 4), "+", round(slope, 4), "* X\n")

# Koefisien determinasi
cat("\nKOEFISIEN DETERMINASI (R-squared):\n")

r_squared <- summary_model$r.squared
cat("R-squared =", round(r_squared, 4), "\n")

# ASUMSI REGRESI LINEAR
# ------------------------
cat("\n\n=== UJI ASUMSI REGRESI LINEAR ===\n")

# Normalitas Residual
cat("\n1. UJI NORMALITAS RESIDUAL (Shapiro-Wilk):\n")

residuals <- resid(model)
shapiro_test <- shapiro.test(residuals)
cat("Statistik Uji W =", round(shapiro_test$statistic, 4), "\n")

cat("p-value =", round(shapiro_test$p.value, 4), "\n")

if (shapiro_test$p.value > 0.05) {
  cat("Keputusan: Residual berdistribusi normal (p > 0.05)\n")
} else {
  cat("Keputusan: Residual tidak berdistribusi normal (p < 0.05)\n")
}

# Plot normalitas residual
par(mfrow = c(1, 2))
qqnorm(residuals, main = "Q-Q Plot Residual")
qqline(residuals, col = "red")
hist(residuals, main = "Histogram Residual", 
     xlab = "Residual", col = "lightblue", border = "black")

# Homoskedastisitas (Uji Breusch-Pagan)
cat("\n2. UJI HOMOSKEDASTISITAS (Breusch-Pagan):\n")

if (!require(lmtest)) {
  install.packages("lmtest")
  library(lmtest)
}
bp_test <- bptest(model)
cat("Statistik Uji LM =", round(bp_test$statistic, 4), "\n")

cat("p-value =", round(bp_test$p.value, 4), "\n")

if (bp_test$p.value > 0.05) {
  cat("Keputusan: Varian residual homogen (homoskedastisitas terpenuhi)\n")
} else {
  cat("Keputusan: Varian residual tidak homogen (ada heteroskedastisitas)\n")
}

# Plot residual vs fitted
par(mfrow = c(1, 2))
plot(fitted(model), residuals,
     main = "Residual vs Fitted Values",
     xlab = "Fitted Values", ylab = "Residuals",
     pch = 19, col = "blue")
abline(h = 0, col = "red", lty = 2)
grid()

# Plot residual vs prediktor
plot(data$Age, residuals,
     main = "Residual vs Prediktor (X)",
     xlab = "Age", ylab = "Residuals",
     pch = 19, col = "blue")
abline(h = 0, col = "red", lty = 2)
grid()

par(mfrow = c(1, 1))

# PENGUJIAN HIPOTESIS
# ----------------------
cat("\n\n=== PENGUJIAN HIPOTESIS ===\n")

# Hipotesis untuk Koefisien Regresi (Slope)
cat("\n1. UJI HIPOTESIS UNTUK SLOPE (β1):\n")

cat("   Hipotesis Nol (H0): β1 = 0 (Tidak ada hubungan linear)\n")

cat("   Hipotesis Alternatif (H1): β1 ≠ 0 (Ada hubungan linear)\n")

# Dari output summary model
coef_table <- summary_model$coefficients
slope_p_value <- coef_table[2, 4]

cat("\n   Hasil Uji:\n")

cat("   t-statistik =", round(coef_table[2, 3], 4), "\n")

cat("   p-value =", round(slope_p_value, 6), "\n")

alpha <- 0.05
cat("   Tingkat signifikansi (α) =", alpha, "\n")

if (slope_p_value < alpha) {
  cat("\n   KEPUTUSAN: Tolak H0\n")
  cat("   KESIMPULAN: Ada bukti statistik yang cukup untuk menyatakan bahwa\n")
  cat("               terdapat hubungan linear yang signifikan antara Age dan Annual Premium\n")
} else {
  cat("\n   KEPUTUSAN: Gagal tolak H0\n")
  cat("   KESIMPULAN: Tidak ada bukti statistik yang cukup untuk menyatakan bahwa\n")
  cat("               terdapat hubungan linear yang signifikan antara Age dan Annual Premium\n")
}

# Interval Kepercayaan
cat("\n2. INTERVAL KEPERCAYAAN 95% UNTUK KOEFISIEN:\n")

conf_int <- confint(model, level = 0.95)
cat("   Intercept (β0): [", round(conf_int[1, 1], 4), ", ", round(conf_int[1, 2], 4), "]\n", sep = "")

cat("   Slope (β1):     [", round(conf_int[2, 1], 4), ", ", round(conf_int[2, 2], 4), "]\n", sep = "")