HANDS ON 2 MG 3 GRACILLIA N.V

library(ggplot2)
library(lmtest)

library(outliers)

library(car)

# Input data
df <- read.csv("C:/Users/User/Downloads/TravelInsurancePrediction.csv",
               sep=";",
               header=TRUE)
data <- df[c("Annual_Income", "Travel_Insurance")]

# Cek Missing Value
cat("\n=== Missing Value ===\n")

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

# Outlier
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 = "Price", 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)

# ANALISIS DESKRIPTIF DAN VISUALISASI
print("Statistik Deskriptif:")

summary(data)

# Scatter plot
ggplot(data, aes(x = Annual_Income, y = Travel_Insurance)) +
  geom_point(color = "blue", size = 3) +
  labs(title = "Hubungan Annual Income dan Travel Insurance",
       x = "Annual_Income", y = "Travel_Insurance") +
  theme_minimal()

# Korelasi
cor_test <- cor.test(data$Annual_Income, data$Travel_Insurance)
print(paste("Korelasi Pearson:", round(cor_test$estimate, 4)))

print(paste("p-value korelasi:", round(cor_test$p.value, 4)))

# MEMBANGUN MODEL REGRESI
model <- lm(Travel_Insurance ~ Annual_Income, data = data)
print("Ringkasan Model Regresi:")

summary(model)

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

# Normalitas Residual
shapiro_test <- shapiro.test(residuals(model))
cat("1. UJI NORMALITAS (Shapiro-Wilk):\n")

cat("   Statistik 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\n")
} else {
  cat("   Keputusan: Residual tidak normal\n")
}

# Q-Q Plot
qqnorm(residuals(model), main = "Q-Q Plot Residual")
qqline(residuals(model), col = "red")

# Homoskedastisitas
bp_test <- bptest(model)
cat("\n2. UJI HOMOSKEDASTISITAS (Breusch-Pagan):\n")

cat("   Statistik 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\n")
} else {
  cat("   Keputusan: Ada heteroskedastisitas\n")
}

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

# Tidak ada Autokorelasi
dw_test <- dwtest(model)
cat("\n3. UJI AUTOKORELASI (Durbin-Watson):\n")

cat("   Statistik DW =", round(dw_test$statistic, 4), "\n")

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

if(dw_test$p.value > 0.05) {
  cat("   Keputusan: Tidak ada autokorelasi\n")
} else {
  cat("   Keputusan: Ada autokorelasi\n")
}

# INTERPRETASI KOEFISIEN
cat("\n=== INTERPRETASI KOEFISIEN ===\n")

intercept <- coef(model)[1]
slope <- coef(model)[2]

cat("Persamaan Regresi: Travel_Insurance =", round(intercept, 2), "+", round(slope, 2), "* Annual_Income\n")

cat("\nInterpretasi:\n")

cat("1. Intercept (β0 =", round(intercept, 2), "):\n")

cat("  Travel_Insurance ketika Annual_Income = 0 adalah", round(intercept, 2))

cat("2. Slope (β1 =", round(slope, 2), "):\n")

cat("   Setiap kenaikan 1 satuan Annual_Income, nilai prediksi Travel_Insurance meningkat sebesar", round(slope, 6))

# ESTIMASI PARAMETER DAN INFERENSI
cat("\n=== ESTIMASI PARAMETER ===\n")

conf_int <- confint(model, level = 0.95)
cat("Interval Kepercayaan 95%:\n")

cat("   Intercept: [", round(conf_int[1,1], 3), ", ", round(conf_int[1,2], 3), "]\n", sep = "")

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

# Uji hipotesis untuk slope
cat("\nUji Hipotesis untuk Slope (β1):\n")

cat("   H0: β1 = 0 (tidak ada hubungan linear)\n")

cat("   H1: β1 ≠ 0 (ada hubungan linear)\n")

summary_model <- summary(model)
slope_pvalue <- summary_model$coefficients[2, 4]
cat("   p-value =", round(slope_pvalue, 6), "\n")

if(slope_pvalue < 0.05) {
  cat("   Keputusan: Tolak H0, ada hubungan linear signifikan\n")
} else {
  cat("   Keputusan: Gagal tolak H0, tidak ada hubungan linear signifikan\n")
}

# KOEFISIEN DETERMINASI
r_squared <- summary_model$r.squared
cat("\nKoefisien Determinasi (R²):\n")

cat("   R² =", round(r_squared, 4), "\n")

cat("   Artinya:", round(r_squared * 100, 2), "% variasi Travel Insurance dapat dijelaskan oleh Annual Income\n")

# VISUALISASI MODEL
ggplot(data, aes(x = Annual_Income, y = Travel_Insurance)) +
  geom_point(color = "blue", size = 3) +
  geom_smooth(method = "lm", se = TRUE, color = "red", fill = "pink") +
  labs(title = "Garis Regresi Linear",
       subtitle = paste("Y =", round(intercept, 2), "+", round(slope, 2), "X"),
       x = "Annual_Income", y = "Travel_Insurance") +
  theme_minimal()

# PREDIKSI
new_data <- data.frame(Annual_Income = c(400000, 700000))
prediction <- predict(model, newdata = new_data, interval = "confidence")
cat("\n=== PREDIKSI ===\n")

cat("Untuk Annual Income 400000, prediksi Travel Insurance =", round(prediction[1, "fit"], 2), "\n")

cat("Untuk Annual Income 700000, prediksi Travel Insurance =", round(prediction[2, "fit"], 2), "\n")

# DIAGNOSTIC PLOTS
par(mfrow = c(2, 2))
plot(model, which = 1:4)

par(mfrow = c(1, 1))

# RINGKASAN LENGKAP
cat("\n=== RINGKASAN ANALISIS ===\n")

cat("1. Model: Travel_Insurance = β0 + β1*Annual_Income + ε\n")

cat("2. Estimasi: Y =", round(intercept, 3), "+", round(slope, 3), "* X\n")

cat("3. R² =", round(r_squared, 4), "(", round(r_squared*100, 1), "%)\n")

f_stat <- summary_model$fstatistic
f_pvalue <- pf(f_stat[1], f_stat[2], f_stat[3], lower.tail = FALSE)
cat("4. Uji F (model): p-value =", round(f_pvalue, 6), "\n")

cat("5. Asumsi:\n")

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

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

cat("   - Autokorelasi: p =", round(dw_test$p.value, 4), "\n")

# Simpan hasil
hasil <- list(
  model = model,
  coefficients = coef(model),
  r_squared = r_squared,
  assumptions = list(
    normality = shapiro_test$p.value,
    homoscedasticity = bp_test$p.value,
    autocorrelation = dw_test$p.value
  ),
  confidence_intervals = conf_int
)