#input data

library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(outliers)
## Warning: package 'outliers' was built under R version 4.5.2
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
df <- read.csv("C:/Users/Msi user/Downloads/datania1k.csv")
data <- df[c("Age", "Annual_Premium")]

#cek missing velue

total_missing <- sum(is.na(data))
total_missing
## [1] 0

#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
  outlier_iqr <- data[[var]][data[[var]] < lower_bound | data [[var]] > upper_bound]
  
  #deteksi outlier dengan metode Z-score
  z_scores <- scale(data[[var]])
  outlier_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(outlier_iqr),
    n_outliers_z = length(outlier_z),
    grubbs_result = grubbs_outlier,
    lower_bound = lower_bound,
    upper_bound = upper_bound,
    outlier_values = unique(round(outlier_iqr, 2))
  )
  
  cat(" - Outlier (IQR method):", length(outlier_iqr), "\n")
  cat(" - Oulier (Z-score > 3):", length(outlier_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(outlier_iqr) > 0) {
    cat(" - Nilai outlier:", paste(head(unique(round(outlier_iqr, 2)), 5), collapse = ", "), "\n")
  }
  cat("\n")
  }
  
  return(outlier_report)
} 

#Deteksi Outlier
outlier_analysis <- detect_outliers(data)
## 
## === DETEKSI OUTLIER ===
## 
## Analisis outlier untuk variabel: Age 
##  - Outlier (IQR method): 0 
##  - Oulier (Z-score > 3): 0 
##  - Uji Grubbs: Tidak Terdeteksi 
##  - Batas Bawah: -12.88 
##  - Batas Atas: 88.12 
## 
## Analisis outlier untuk variabel: Annual_Premium 
##  - Outlier (IQR method): 26 
##  - Oulier (Z-score > 3): 5 
##  - Uji Grubbs: Terdeteksi 
##  - Batas Bawah: 1704.5 
##  - Batas Atas: 62266.5 
##  - Nilai outlier: 81192, 100278, 63273, 70452, 71918
#visual
visualize_outliers <- function(data) {
  cat("\n=== VISUALISASI OUTLIER ===\n")
  
  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)
## 
## === VISUALISASI OUTLIER ===

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#Statistik Deskriptif

cat("\n=== Statistik Deskriptif ===\n")
## 
## === Statistik Deskriptif ===
summary(data)
##       Age        Annual_Premium  
##  Min.   :20.00   Min.   :  2630  
##  1st Qu.:25.00   1st Qu.: 24415  
##  Median :38.00   Median : 31512  
##  Mean   :39.65   Mean   : 30364  
##  3rd Qu.:50.25   3rd Qu.: 39556  
##  Max.   :85.00   Max.   :100278

#Analisis Regresi

#Korelasi 
cat("\nKorelasi:\n")
## 
## Korelasi:
correlation <- cor(data$Age, data$Annual_Premium)
cat("Koefisien Korelasi (r) =", round(correlation, 4), "\n")
## Koefisien Korelasi (r) = 0.1388
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")
}
## Interpretasi: Korelasi Sangat Lemah/Tidak Ada
#Visual Data
cat("\n\n=== Visualisasi Data ===\n")
## 
## 
## === Visualisasi Data ===
#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 = "Äge (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")
## 
## PERSAMAAN REGRESI LINEAR:
beta0 <- coef(model)[1]
beta1 <- coef(model)[2]

#Ringkasan Model
cat("\nRINGKASAN MODEL:\n")
## 
## RINGKASAN MODEL:
summary_model <- summary(model)
print(summary_model)
## 
## Call:
## lm(formula = Age ~ Annual_Premium, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.016 -14.305  -0.939  10.210  45.079 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    3.559e+01  1.043e+00  34.114  < 2e-16 ***
## Annual_Premium 1.339e-04  3.025e-05   4.427 1.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.63 on 998 degrees of freedom
## Multiple R-squared:  0.01926,    Adjusted R-squared:  0.01827 
## F-statistic: 19.59 on 1 and 998 DF,  p-value: 1.063e-05
#Ekstrak Komponen Penting
cat("\nPERSAMAAN REGRESI:\n")
## 
## PERSAMAAN REGRESI:
intercept <- coef(model)[1]
slope <- coef(model)[2]
cat("Persamaan Regresi: Y =", round(intercept, 4), "+", round(slope, 4), "* X\n")
## Persamaan Regresi: Y = 35.5867 + 1e-04 * X
#Koefisien Determinasi
cat("\nKOEFISIEN DETERMINASI (R-squared):\n")
## 
## KOEFISIEN DETERMINASI (R-squared):
r_squared <- summary_model$r.squared
cat("R-squared =", round(r_squared, 4), "\n")
## R-squared = 0.0193
#ASUMSI REGRESI LINEAR
#-----------------------
cat("\n\n=== UJI ASUMSI REGRESI LINEAR ===\n")
## 
## 
## === UJI ASUMSI REGRESI LINEAR ===
#Normalitas Residual
cat("\n1. UJI NORMALITAS RESIDUAL (Shapiro-Wilk):\n")
## 
## 1. UJI NORMALITAS RESIDUAL (Shapiro-Wilk):
residuals <- resid(model)
shapiro_test <- shapiro.test(residuals)
cat("Statistik Uji W =", round(shapiro_test$statistic, 4), "\n")
## Statistik Uji W = 0.9339
cat("p-value =", round(shapiro_test$p.value, 4), "\n")
## p-value = 0
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")
}
## Keputusan: Residual tidak berdistribusi normal (p < 0.05)
# 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")
## 
## 2. UJI HOMOSKEDASTISITAS (Breusch-Pagan):
if (!require(lmtest)) {
  install.packages("lmtest")
  library(lmtest)
}
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 4.5.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.5.2
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bp_test <- bptest(model)
cat("Statistik Uji LM =", round(bp_test$statistic, 4), "\n")
## Statistik Uji LM = 9.2047
cat("p-value =", round(bp_test$p.value, 4), "\n")
## p-value = 0.0024
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")
}
## Keputusan: Varian residual tidak homogen (ada heteroskedastisitas)
# 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")
## 
## 
## === PENGUJIAN HIPOTESIS ===
# Hipotesis untuk Koefisien Regresi (Slope)
cat("\n1. UJI HIPOTESIS UNTUK SLOPE (β1):\n")
## 
## 1. UJI HIPOTESIS UNTUK SLOPE (β1):
cat("   Hipotesis Nol (H0): β1 = 0 (Tidak ada hubungan linear)\n")
##    Hipotesis Nol (H0): β1 = 0 (Tidak ada hubungan linear)
cat("   Hipotesis Alternatif (H1): β1 ≠ 0 (Ada hubungan linear)\n")
##    Hipotesis Alternatif (H1): β1 ≠ 0 (Ada hubungan linear)
# Dari output summary model
coef_table <- summary_model$coefficients
slope_p_value <- coef_table[2, 4]

cat("\n   Hasil Uji:\n")
## 
##    Hasil Uji:
cat("   t-statistik =", round(coef_table[2, 3], 4), "\n")
##    t-statistik = 4.4266
cat("   p-value =", round(slope_p_value, 6), "\n")
##    p-value = 1.1e-05
alpha <- 0.05
cat("   Tingkat signifikansi (α) =", alpha, "\n")
##    Tingkat signifikansi (α) = 0.05
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")
}
## 
##    KEPUTUSAN: Tolak H0
##    KESIMPULAN: Ada bukti statistik yang cukup untuk menyatakan bahwa
##                terdapat hubungan linear yang signifikan antara Age dan Annual Premium
# Interval Kepercayaan
cat("\n2. INTERVAL KEPERCAYAAN 95% UNTUK KOEFISIEN:\n")
## 
## 2. INTERVAL KEPERCAYAAN 95% UNTUK KOEFISIEN:
conf_int <- confint(model, level = 0.95)
cat("   Intercept (β0): [", round(conf_int[1, 1], 4), ", ", round(conf_int[1, 2], 4), "]\n", sep = "")
##    Intercept (β0): [33.5397, 37.6338]
cat("   Slope (β1):     [", round(conf_int[2, 1], 4), ", ", round(conf_int[2, 2], 4), "]\n", sep = "")
##    Slope (β1):     [1e-04, 2e-04]