Analisis Regresi
Ghea Gabriella Karowe, Rusdianto
2026-02-04
#DATA
# Install dan load packages yang diperlukan
if (!require("tidyverse")) install.packages("tidyverse")## Loading required package: tidyverse## Warning: package 'tidyverse' was built under R version 4.5.2## Warning: package 'ggplot2' was built under R version 4.5.2## Warning: package 'tibble' was built under R version 4.5.2## Warning: package 'tidyr' was built under R version 4.5.2## Warning: package 'readr' was built under R version 4.5.2## Warning: package 'purrr' was built under R version 4.5.2## Warning: package 'dplyr' was built under R version 4.5.2## Warning: package 'stringr' was built under R version 4.5.2## Warning: package 'forcats' was built under R version 4.5.2## Warning: package 'lubridate' was built under R version 4.5.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.6
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.1 ✔ tibble 3.3.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsif (!require("VIM")) install.packages("VIM")## Loading required package: VIM## Warning: package 'VIM' was built under R version 4.5.2## Loading required package: colorspace## Warning: package 'colorspace' was built under R version 4.5.2## Loading required package: grid
## VIM is ready to use.
##
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
##
## Attaching package: 'VIM'
##
## The following object is masked from 'package:datasets':
##
## sleepif (!require("naniar")) install.packages("naniar")## Loading required package: naniar## Warning: package 'naniar' was built under R version 4.5.2if (!require("outliers")) install.packages("outliers")## Loading required package: outliers## Warning: package 'outliers' was built under R version 4.5.2if (!require("ggplot2")) install.packages("ggplot2")
library(tidyverse)
library(VIM)
library(naniar)
library(outliers)
library(ggplot2)
library(mice)## Warning: package 'mice' was built under R version 4.5.2##
## Attaching package: 'mice'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## cbind, rbindlibrary(dplyr)
data <- read.csv("C:/Users/Asus/Downloads//datania1k.csv")
Y <- data[c("Age","Annual_Premium")]
head(Y)## Age Annual_Premium
## 1 22 36513
## 2 24 2630
## 3 22 35832
## 4 72 36685
## 5 66 2630
## 6 42 31226#DATA HILANG
# 1. Ringkasan missing values
analyze_missing_values <- function(Y) {
cat("=== ANALISIS DATA HILANG ===\n\n")
# Total missing values
total_missing <- sum(is.na(Y))
cat("Total nilai hilang:", total_missing, "\n")
cat("Persentase nilai hilang:", round(total_missing / (nrow(Y) * ncol(Y)) * 100, 2), "%\n\n")
# Missing values per variabel
missing_per_var <- sapply(Y, function(x) sum(is.na(x)))
missing_df <- data.frame(
Variable = names(missing_per_var),
Missing_Count = missing_per_var,
Missing_Percentage = round(missing_per_var / nrow(Y) * 100, 2)
) %>%
arrange(desc(Missing_Count))
print("Missing values per variabel:")
print(missing_df)
# Pattern missing values
cat("\nPattern data hilang (5 observasi pertama):\n")
print(head(md.pattern(Y, plot = FALSE), 5))
# Visualisasi missing values
par(mfrow = c(1, 2))
# Plot 1: Aggregation plot
tryCatch({
aggr_plot <- aggr(Y, col = c('navyblue', 'red'),
numbers = TRUE, sortVars = TRUE,
labels = names(Y), cex.axis = 0.7,
gap = 3, ylab = c("Histogram data hilang", "Pattern"))
}, error = function(e) {
cat("Gagal membuat plot aggr, menggunakan alternatif...\n")
})
# Plot 2: Heatmap missing values
vis_miss(Y, warn_large_data = FALSE) +
theme_minimal() +
labs(title = "Heatmap Data Hilang")
par(mfrow = c(1, 1))
return(missing_df)
}
# Jalankan analisis missing values
missing_analysis <- analyze_missing_values(Y)## === ANALISIS DATA HILANG ===
##
## Total nilai hilang: 0
## Persentase nilai hilang: 0 %
##
## [1] "Missing values per variabel:"
## Variable Missing_Count Missing_Percentage
## Age Age 0 0
## Annual_Premium Annual_Premium 0 0
##
## Pattern data hilang (5 observasi pertama):
## /\ /\
## { `---' }
## { O O }
## ==> V <== No need for mice. This data set is completely observed.
## \ \|/ /
## `-----'
##
## Age Annual_Premium
## 1000 1 1 0
## 0 0 0
##
## Variables sorted by number of missings:
## Variable Count
## Age 0
## Annual_Premium 0#OUTLIER
# 3. Deteksi outlier
detect_outliers <- function(Y) {
cat("\n=== DETEKSI OUTLIER ===\n\n")
outlier_report <- list()
for(var in names(Y)[sapply(Y, is.numeric)]) {
cat("Analisis outlier untuk variabel:", var, "\n")
# Statistik deskriptif
stats <- summary(Y[[var]])
iqr_val <- IQR(Y[[var]], na.rm = TRUE)
q1 <- quantile(Y[[var]], 0.25, na.rm = TRUE)
q3 <- quantile(Y[[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 <- Y[[var]][Y[[var]] < lower_bound | Y[[var]] > upper_bound]
# Deteksi outlier dengan metode Z-score
z_scores <- scale(Y[[var]])
outliers_z <- Y[[var]][abs(z_scores) > 3]
# Deteksi outlier dengan metode Grubbs (uji statistik)
if(length(na.omit(Y[[var]])) > 6) {
tryCatch({
grubbs_test <- grubbs.test(na.omit(Y[[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(Y)##
## === DETEKSI OUTLIER ===
##
## Analisis outlier untuk variabel: Age
## - Outlier (IQR method): 0
## - Outlier (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
## - Outlier (Z-score > 3): 5
## - Uji Grubbs: Terdeteksi
## - Batas bawah: 1704.5
## - Batas atas: 62266.5
## - Nilai outlier: 81192, 100278, 63273, 70452, 71918# 4. Visualisasi outlier
visualize_outliers <- function(Y) {
cat("\n=== VISUALISASI OUTLIER ===\n")
numeric_vars <- names(Y)[sapply(Y, 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 <- Y[[var]] < lower_bound | Y[[var]] > upper_bound
# Plot histogram
hist_data <- ggplot(data.frame(value = Y[[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(Y[, 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(Y)##
## === 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.
# 5. Metode penanganan outlier
handle_outliers <- function(Y, method = "cap") {
cat("\n=== PENANGANAN OUTLIER ===\n")
cat("Metode yang digunakan:", method, "\n\n")
data_processed <- Y
for(var in names(data_processed)[sapply(data_processed, is.numeric)]) {
# Hitung batas outlier
q1 <- quantile(data_processed[[var]], 0.25, na.rm = TRUE)
q3 <- quantile(data_processed[[var]], 0.75, na.rm = TRUE)
iqr_val <- IQR(data_processed[[var]], na.rm = TRUE)
lower_bound <- q1 - 1.5 * iqr_val
upper_bound <- q3 + 1.5 * iqr_val
# Identifikasi outlier
is_outlier <- data_processed[[var]] < lower_bound | data_processed[[var]] > upper_bound
n_outliers <- sum(is_outlier, na.rm = TRUE)
if(n_outliers > 0) {
cat("Variabel", var, "memiliki", n_outliers, "outlier\n")
if(method == "cap") {
# Cap outliers (Winsorizing)
data_processed[[var]][data_processed[[var]] < lower_bound] <- lower_bound
data_processed[[var]][data_processed[[var]] > upper_bound] <- upper_bound
cat(" Outlier di-cap ke batas [", round(lower_bound, 2), ", ",
round(upper_bound, 2), "]\n", sep = "")
} else if(method == "remove") {
# Hapus outlier (set sebagai NA)
data_processed[[var]][is_outlier] <- NA
cat(" Outlier dihapus (dijadikan NA)\n")
} else if(method == "transform") {
# Transformasi logaritmik
if(all(data_processed[[var]] > 0, na.rm = TRUE)) {
data_processed[[var]] <- log(data_processed[[var]] + 1)
cat(" Dilakukan transformasi log\n")
} else {
cat(" Transformasi log tidak dapat dilakukan (nilai negatif)\n")
}
} else if(method == "median") {
# Ganti dengan median
median_val <- median(data_processed[[var]], na.rm = TRUE)
data_processed[[var]][is_outlier] <- median_val
cat(" Diganti dengan median:", round(median_val, 2), "\n")
}
}
}
return(data_processed)
}
# Contoh penggunaan metode penanganan outlier
data_capped <- handle_outliers(Y, method = "cap")##
## === PENANGANAN OUTLIER ===
## Metode yang digunakan: cap
##
## Variabel Annual_Premium memiliki 26 outlier
## Outlier di-cap ke batas [1704.5, 62266.5]data_no_outliers <- handle_outliers(Y, method = "remove")##
## === PENANGANAN OUTLIER ===
## Metode yang digunakan: remove
##
## Variabel Annual_Premium memiliki 26 outlier
## Outlier dihapus (dijadikan NA)#ANALISIS REGRESI
# ====================================================
# ANALISIS REGRESI LINEAR SEDERHANA
# Data: Age vs Annual_Premium
# ====================================================
# 1. MEMUAT DATA
# ----------------------
# Baca file CSV (sesuaikan path jika perlu)
data_studi <- read.csv("datania1k.csv")
# Ambil hanya kolom yang dibutuhkan
data_studi <- data_studi[, c("Age", "Annual_Premium")]
# Ubah nama variabel agar mudah dipanggil
colnames(data_studi) <- c("Age", "Annual_Premium")
cat("DATA YANG DIGUNAKAN:\n")## DATA YANG DIGUNAKAN:print(head(data_studi))## Age Annual_Premium
## 1 22 36513
## 2 24 2630
## 3 22 35832
## 4 72 36685
## 5 66 2630
## 6 42 31226# 2. ANALISIS DESKRIPTIF
# ----------------------
cat("\n\n=== ANALISIS DESKRIPTIF ===\n")##
##
## === ANALISIS DESKRIPTIF ===desc_stats <- data.frame(
Variabel = c("Age (X)", "Annual Premium (Y)"),
Mean = c(mean(data_studi$Age),
mean(data_studi$Annual_Premium)),
SD = c(sd(data_studi$Age),
sd(data_studi$Annual_Premium)),
Min = c(min(data_studi$Age),
min(data_studi$Annual_Premium)),
Max = c(max(data_studi$Age),
max(data_studi$Annual_Premium))
)
print(desc_stats)## Variabel Mean SD Min Max
## 1 Age (X) 39.653 15.77693 20 85
## 2 Annual Premium (Y) 30364.102 16348.15212 2630 100278# Korelasi
correlation <- cor(data_studi$Age,
data_studi$Annual_Premium)
cat("\nKoefisien Korelasi (r) =",
round(correlation, 4), "\n")##
## Koefisien Korelasi (r) = 0.1388# 3. VISUALISASI
# ----------------------
plot(data_studi$Age,
data_studi$Annual_Premium,
main = "Scatter Plot: Age vs Annual Premium",
xlab = "Age",
ylab = "Annual Premium",
pch = 19,
col = "blue")
# 4. REGRESI LINEAR
# ----------------------
model <- lm(Annual_Premium ~ Age,
data = data_studi)
summary_model <- summary(model)
cat("\n=== RINGKASAN MODEL ===\n")##
## === RINGKASAN MODEL ===print(summary_model)##
## Call:
## lm(formula = Annual_Premium ~ Age, data = data_studi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32960 -5971 1517 9265 72308
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24662.40 1386.17 17.792 < 2e-16 ***
## Age 143.79 32.48 4.427 1.06e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16200 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# Persamaan regresi
intercept <- coef(model)[1]
slope <- coef(model)[2]
cat("\nPersamaan Regresi:\n")##
## Persamaan Regresi:cat("Y =", round(intercept,4),
"+", round(slope,4), "X\n")## Y = 24662.4 + 143.7899 X# R-squared
r_squared <- summary_model$r.squared
cat("\nR-squared =",
round(r_squared,4), "\n")##
## R-squared = 0.0193# 5. UJI ASUMSI
# ----------------------
# Normalitas
residuals <- resid(model)
shapiro_test <- shapiro.test(residuals)
cat("\nUji Normalitas Shapiro-Wilk\n")##
## Uji Normalitas Shapiro-Wilkprint(shapiro_test)##
## Shapiro-Wilk normality test
##
## data: residuals
## W = 0.94079, p-value < 2.2e-16# Homoskedastisitas
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.numericbp_test <- bptest(model)
cat("\nUji Breusch-Pagan\n")##
## Uji Breusch-Paganprint(bp_test)##
## studentized Breusch-Pagan test
##
## data: model
## BP = 4.9487, df = 1, p-value = 0.02611# Plot residual
plot(fitted(model), residuals,
main = "Residual vs Fitted",
xlab = "Fitted",
ylab = "Residual",
pch = 19,
col = "blue")
abline(h = 0, col = "red")
# 6. UJI HIPOTESIS KOEFISIEN
# ----------------------
coef_table <- summary_model$coefficients
p_value <- coef_table[2,4]
cat("\n=== UJI HIPOTESIS SLOPE ===\n")##
## === UJI HIPOTESIS SLOPE ===cat("p-value =", p_value, "\n")## p-value = 1.062773e-05if (p_value < 0.05) {
cat("Kesimpulan: Signifikan\n")
} else {
cat("Kesimpulan: Tidak signifikan\n")
}## Kesimpulan: Signifikan# 7. GARIS REGRESI
# ----------------------
plot(data_studi$Age,
data_studi$Annual_Premium,
main = "Regresi Linear Age vs Premium",
xlab = "Age",
ylab = "Annual Premium",
pch = 19,
col = "blue")
abline(model, col = "red", lwd = 2)
# 8. INTERPRETASI SINGKAT
# ----------------------
cat("\n=== INTERPRETASI ===\n")##
## === INTERPRETASI ===cat("Setiap kenaikan 1 tahun usia,\n")## Setiap kenaikan 1 tahun usia,cat("premium berubah sebesar",
round(slope,2), "\n")## premium berubah sebesar 143.79cat("Model menjelaskan",
round(r_squared*100,2),
"% variasi premium\n")## Model menjelaskan 1.93 % variasi premium