0. Persiapan: Library dan Data
# ============================================================
# ANALISIS MULTIPLE LINEAR REGRESSION (MLR)
# Data: Food Delivery Times (waktu pengantaran makanan)
# Alur: EDA -> Estimasi -> Uji F & t -> R2 -> Uji Asumsi ->
# Diagnostik Outlier -> Pemilihan Model -> Ringkasan
# ============================================================
# ------------------------------------------------------------
# 0. PERSIAPAN: LOAD LIBRARY
# ------------------------------------------------------------
# Install package jika belum ada (uncomment jika perlu)
# install.packages(c("car", "lmtest", "nortest", "moments", "corrplot"))
library(car) # Untuk VIF
## Warning: package 'car' was built under R version 4.5.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.5.2
library(lmtest) # Untuk uji Breusch-Pagan dan Ramsey RESET
## Warning: package 'lmtest' was built under R version 4.5.3
## 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
library(nortest) # Untuk uji normalitas Anderson-Darling
## Warning: package 'nortest' was built under R version 4.5.2
library(moments) # Untuk skewness dan kurtosis residual
## Warning: package 'moments' was built under R version 4.5.2
library(corrplot) # Untuk matriks korelasi
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
1. Data dan Preprocessing
# ------------------------------------------------------------
# 1. MEMBACA DATA
# ------------------------------------------------------------
# Nilai kosong pada file CSV berupa sel kosong, sehingga perlu
# didefinisikan sebagai NA saat membaca data
data_asli <- read.csv("Food_Delivery_Times.csv",
na.strings = c("", "NA"),
stringsAsFactors = FALSE)
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("DATA YANG DIGUNAKAN\n")
## DATA YANG DIGUNAKAN
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
str(data_asli)
## 'data.frame': 1000 obs. of 9 variables:
## $ Order_ID : int 522 738 741 661 412 679 627 514 860 137 ...
## $ Distance_km : num 7.93 16.42 9.52 7.44 19.03 ...
## $ Weather : chr "Windy" "Clear" "Foggy" "Rainy" ...
## $ Traffic_Level : chr "Low" "Medium" "Low" "Medium" ...
## $ Time_of_Day : chr "Afternoon" "Evening" "Night" "Afternoon" ...
## $ Vehicle_Type : chr "Scooter" "Bike" "Scooter" "Scooter" ...
## $ Preparation_Time_min : int 12 20 28 5 16 8 12 5 20 29 ...
## $ Courier_Experience_yrs: num 1 2 1 1 5 9 1 6 6 1 ...
## $ Delivery_Time_min : int 43 84 59 37 68 57 49 46 35 73 ...
print(head(data_asli))
## Order_ID Distance_km Weather Traffic_Level Time_of_Day Vehicle_Type
## 1 522 7.93 Windy Low Afternoon Scooter
## 2 738 16.42 Clear Medium Evening Bike
## 3 741 9.52 Foggy Low Night Scooter
## 4 661 7.44 Rainy Medium Afternoon Scooter
## 5 412 19.03 Clear Low Morning Bike
## 6 679 19.40 Clear Low Evening Scooter
## Preparation_Time_min Courier_Experience_yrs Delivery_Time_min
## 1 12 1 43
## 2 20 2 84
## 3 28 1 59
## 4 5 1 37
## 5 16 5 68
## 6 8 9 57
# ------------------------------------------------------------
# 1.1 MISSING VALUE
# ------------------------------------------------------------
cat("\nJumlah missing value per variabel:\n")
##
## Jumlah missing value per variabel:
print(colSums(is.na(data_asli)))
## Order_ID Distance_km Weather
## 0 0 30
## Traffic_Level Time_of_Day Vehicle_Type
## 30 30 0
## Preparation_Time_min Courier_Experience_yrs Delivery_Time_min
## 0 30 0
cat("\nBaris dengan minimal 1 missing value:",
sum(!complete.cases(data_asli)), "dari", nrow(data_asli), "baris\n")
##
## Baris dengan minimal 1 missing value: 117 dari 1000 baris
# ------------------------------------------------------------
# 1.2 PREPROCESSING
# ------------------------------------------------------------
# Order_ID hanya penanda pesanan (bukan variabel penjelas) -> dibuang
data <- data_asli[, setdiff(names(data_asli), "Order_ID")]
# Missing value hanya sekitar 3% per variabel, dan jumlah data
# setelah penghapusan masih besar -> baris yang memuat missing dihapus
data <- na.omit(data)
# Variabel kategorik diubah menjadi factor. Urutan level menentukan
# kategori referensi (level pertama) pada koefisien regresi
data$Weather <- factor(data$Weather,
levels = c("Clear", "Foggy", "Rainy", "Snowy", "Windy"))
data$Traffic_Level <- factor(data$Traffic_Level,
levels = c("Low", "Medium", "High"))
data$Time_of_Day <- factor(data$Time_of_Day,
levels = c("Morning", "Afternoon", "Evening", "Night"))
data$Vehicle_Type <- factor(data$Vehicle_Type,
levels = c("Bike", "Scooter", "Car"))
n <- nrow(data)
cat("\nJumlah observasi setelah preprocessing (n):", n, "\n")
##
## Jumlah observasi setelah preprocessing (n): 883
cat("Jumlah variabel prediktor:", ncol(data) - 1, "\n")
## Jumlah variabel prediktor: 7
2. Statistik Deskriptif dan Eksplorasi Data (EDA)
# ------------------------------------------------------------
# 2. STATISTIK DESKRIPTIF DAN EKSPLORASI DATA
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("1. STATISTIK DESKRIPTIF\n")
## 1. STATISTIK DESKRIPTIF
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Statistik deskriptif
summary(data)
## Distance_km Weather Traffic_Level Time_of_Day Vehicle_Type
## Min. : 0.59 Clear:425 Low :347 Morning :281 Bike :451
## 1st Qu.: 5.13 Foggy: 98 Medium:355 Afternoon:257 Scooter:260
## Median :10.28 Rainy:188 High :181 Evening :269 Car :172
## Mean :10.05 Snowy: 86 Night : 76
## 3rd Qu.:15.03 Windy: 86
## Max. :19.99
## Preparation_Time_min Courier_Experience_yrs Delivery_Time_min
## Min. : 5.00 Min. :0.00 Min. : 8.00
## 1st Qu.:11.00 1st Qu.:2.00 1st Qu.: 41.00
## Median :17.00 Median :5.00 Median : 55.00
## Mean :17.02 Mean :4.64 Mean : 56.43
## 3rd Qu.:23.00 3rd Qu.:7.00 3rd Qu.: 71.00
## Max. :29.00 Max. :9.00 Max. :141.00
# Matriks korelasi (hanya untuk variabel numerik)
var_numerik <- c("Delivery_Time_min", "Distance_km",
"Preparation_Time_min", "Courier_Experience_yrs")
cat("\nMatriks Korelasi:\n")
##
## Matriks Korelasi:
cor_matrix <- cor(data[, var_numerik])
print(round(cor_matrix, 3))
## Delivery_Time_min Distance_km Preparation_Time_min
## Delivery_Time_min 1.000 0.783 0.297
## Distance_km 0.783 1.000 -0.026
## Preparation_Time_min 0.297 -0.026 1.000
## Courier_Experience_yrs -0.077 -0.003 -0.031
## Courier_Experience_yrs
## Delivery_Time_min -0.077
## Distance_km -0.003
## Preparation_Time_min -0.031
## Courier_Experience_yrs 1.000
# Visualisasi matriks korelasi
corrplot(cor_matrix, method = "color", type = "upper",
addCoef.col = "black", tl.col = "black", tl.srt = 45,
title = "Matriks Korelasi Antar Variabel Numerik", mar = c(0, 0, 2, 0))
# Scatter plot matrix
pairs(data[, var_numerik], main = "Scatter Plot Matrix", pch = 19, col = "steelblue")
# Distribusi dan outlier variabel respon
par(mfrow = c(1, 2))
hist(data$Delivery_Time_min, main = "Histogram Waktu Pengantaran",
xlab = "Waktu pengantaran (menit)", col = "steelblue", border = "white")
boxplot(data$Delivery_Time_min, main = "Boxplot Waktu Pengantaran",
ylab = "Waktu pengantaran (menit)", col = "steelblue")
par(mfrow = c(1, 1))
outlier_y <- boxplot.stats(data$Delivery_Time_min)$out
cat("\nJumlah outlier Delivery_Time_min (aturan boxplot):", length(outlier_y), "\n")
##
## Jumlah outlier Delivery_Time_min (aturan boxplot): 4
cat("(Outlier tidak dihapus karena merupakan nilai yang masuk akal, bukan salah input;\n")
## (Outlier tidak dihapus karena merupakan nilai yang masuk akal, bukan salah input;
cat(" pengaruhnya terhadap model diperiksa pada bagian diagnostik.)\n")
## pengaruhnya terhadap model diperiksa pada bagian diagnostik.)
# Waktu pengantaran menurut variabel kategorik
par(mfrow = c(2, 2))
boxplot(Delivery_Time_min ~ Weather, data = data, main = "Menurut Cuaca",
xlab = "", ylab = "Menit", col = "lightblue")
boxplot(Delivery_Time_min ~ Traffic_Level, data = data, main = "Menurut Lalu Lintas",
xlab = "", ylab = "Menit", col = "lightblue")
boxplot(Delivery_Time_min ~ Time_of_Day, data = data, main = "Menurut Waktu",
xlab = "", ylab = "Menit", col = "lightblue")
boxplot(Delivery_Time_min ~ Vehicle_Type, data = data, main = "Menurut Kendaraan",
xlab = "", ylab = "Menit", col = "lightblue")
par(mfrow = c(1, 1))
3. Estimasi Model MLR
# ------------------------------------------------------------
# 3. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("2. ESTIMASI MODEL MLR\n")
## 2. ESTIMASI MODEL MLR
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Fit model MLR (model penuh, semua prediktor)
model <- lm(Delivery_Time_min ~ Distance_km + Preparation_Time_min +
Courier_Experience_yrs + Weather + Traffic_Level +
Time_of_Day + Vehicle_Type,
data = data)
# Ringkasan model
cat("\nRingkasan Model:\n")
##
## Ringkasan Model:
print(summary(model))
##
## Call:
## lm(formula = Delivery_Time_min ~ Distance_km + Preparation_Time_min +
## Courier_Experience_yrs + Weather + Traffic_Level + Time_of_Day +
## Vehicle_Type, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.488 -4.684 -1.189 3.012 61.446
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.80052 1.45345 3.303 0.000996 ***
## Distance_km 2.95778 0.05959 49.632 < 2e-16 ***
## Preparation_Time_min 0.95199 0.04664 20.412 < 2e-16 ***
## Courier_Experience_yrs -0.55123 0.11623 -4.743 2.46e-06 ***
## WeatherFoggy 7.78608 1.12570 6.917 8.97e-12 ***
## WeatherRainy 4.90641 0.88229 5.561 3.57e-08 ***
## WeatherSnowy 9.84553 1.19586 8.233 6.67e-16 ***
## WeatherWindy 2.13038 1.19369 1.785 0.074659 .
## Traffic_LevelMedium 6.04437 0.76392 7.912 7.68e-15 ***
## Traffic_LevelHigh 12.09852 0.92239 13.116 < 2e-16 ***
## Time_of_DayAfternoon 0.70893 0.86777 0.817 0.414174
## Time_of_DayEvening 1.35702 0.85929 1.579 0.114647
## Time_of_DayNight -0.60854 1.30464 -0.466 0.641014
## Vehicle_TypeScooter -1.20253 0.78460 -1.533 0.125723
## Vehicle_TypeCar 0.26682 0.90421 0.295 0.768001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.01 on 868 degrees of freedom
## Multiple R-squared: 0.7879, Adjusted R-squared: 0.7845
## F-statistic: 230.3 on 14 and 868 DF, p-value: < 2.2e-16
# Ekstrak koefisien
coef_tab <- summary(model)$coefficients
cat("\nKoefisien Regresi:\n")
##
## Koefisien Regresi:
coef_df <- data.frame(
Variabel = rownames(coef_tab),
Koefisien = round(coef_tab[, 1], 4),
Std_Error = round(coef_tab[, 2], 4),
t_value = round(coef_tab[, 3], 4),
p_value = round(coef_tab[, 4], 4),
row.names = NULL
)
print(coef_df)
## Variabel Koefisien Std_Error t_value p_value
## 1 (Intercept) 4.8005 1.4534 3.3028 0.0010
## 2 Distance_km 2.9578 0.0596 49.6320 0.0000
## 3 Preparation_Time_min 0.9520 0.0466 20.4120 0.0000
## 4 Courier_Experience_yrs -0.5512 0.1162 -4.7428 0.0000
## 5 WeatherFoggy 7.7861 1.1257 6.9167 0.0000
## 6 WeatherRainy 4.9064 0.8823 5.5610 0.0000
## 7 WeatherSnowy 9.8455 1.1959 8.2330 0.0000
## 8 WeatherWindy 2.1304 1.1937 1.7847 0.0747
## 9 Traffic_LevelMedium 6.0444 0.7639 7.9123 0.0000
## 10 Traffic_LevelHigh 12.0985 0.9224 13.1165 0.0000
## 11 Time_of_DayAfternoon 0.7089 0.8678 0.8170 0.4142
## 12 Time_of_DayEvening 1.3570 0.8593 1.5792 0.1146
## 13 Time_of_DayNight -0.6085 1.3046 -0.4664 0.6410
## 14 Vehicle_TypeScooter -1.2025 0.7846 -1.5327 0.1257
## 15 Vehicle_TypeCar 0.2668 0.9042 0.2951 0.7680
# Interval kepercayaan 95%
cat("\nInterval Kepercayaan 95%:\n")
##
## Interval Kepercayaan 95%:
ci <- confint(model, level = 0.95)
print(round(ci, 4))
## 2.5 % 97.5 %
## (Intercept) 1.9478 7.6532
## Distance_km 2.8408 3.0747
## Preparation_Time_min 0.8604 1.0435
## Courier_Experience_yrs -0.7794 -0.3231
## WeatherFoggy 5.5767 9.9955
## WeatherRainy 3.1747 6.6381
## WeatherSnowy 7.4984 12.1926
## WeatherWindy -0.2125 4.4732
## Traffic_LevelMedium 4.5450 7.5437
## Traffic_LevelHigh 10.2881 13.9089
## Time_of_DayAfternoon -0.9942 2.4121
## Time_of_DayEvening -0.3295 3.0435
## Time_of_DayNight -3.1692 1.9521
## Vehicle_TypeScooter -2.7425 0.3374
## Vehicle_TypeCar -1.5079 2.0415
4. Uji Signifikansi Simultan (Uji F)
# ------------------------------------------------------------
# 4. UJI SIGNIFIKANSI SIMULTAN (UJI F)
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("3. UJI SIGNIFIKANSI SIMULTAN (UJI F)\n")
## 3. UJI SIGNIFIKANSI SIMULTAN (UJI F)
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Ekstrak informasi F-test
f_stat <- summary(model)$fstatistic
f_value <- unname(f_stat[1])
df1 <- unname(f_stat[2])
df2 <- unname(f_stat[3])
f_pvalue <- pf(f_value, df1, df2, lower.tail = FALSE)
cat("\nHipotesis:\n")
##
## Hipotesis:
cat("H0: β1 = β2 = ... = βp = 0 (Model tidak signifikan)\n")
## H0: β1 = β2 = ... = βp = 0 (Model tidak signifikan)
cat("H1: Minimal ada satu βj ≠ 0 (Model signifikan)\n")
## H1: Minimal ada satu βj ≠ 0 (Model signifikan)
cat("\nHasil Uji F:\n")
##
## Hasil Uji F:
cat("F-statistic:", round(f_value, 4), "\n")
## F-statistic: 230.3132
cat("df1:", df1, "\n")
## df1: 14
cat("df2:", df2, "\n")
## df2: 868
cat("p-value:", format(f_pvalue, scientific = TRUE), "\n")
## p-value: 1.226723e-280
if (f_pvalue < 0.05) {
cat("Kesimpulan: Tolak H0 → Model signifikan pada α = 5%\n")
} else {
cat("Kesimpulan: Gagal Tolak H0 → Model tidak signifikan\n")
}
## Kesimpulan: Tolak H0 → Model signifikan pada α = 5%
5. Uji Signifikansi Parsial (Uji t)
# ------------------------------------------------------------
# 5. UJI SIGNIFIKANSI PARSIAL (UJI T)
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("4. UJI SIGNIFIKANSI PARSIAL (UJI T)\n")
## 4. UJI SIGNIFIKANSI PARSIAL (UJI T)
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("\nHipotesis untuk setiap βj:\n")
##
## Hipotesis untuk setiap βj:
cat("H0: βj = 0 (Variabel tidak signifikan)\n")
## H0: βj = 0 (Variabel tidak signifikan)
cat("H1: βj ≠ 0 (Variabel signifikan)\n")
## H1: βj ≠ 0 (Variabel signifikan)
cat("(Untuk variabel kategorik, setiap koefisien membandingkan satu kategori\n")
## (Untuk variabel kategorik, setiap koefisien membandingkan satu kategori
cat(" dengan kategori referensinya: Clear, Low, Morning, dan Bike)\n")
## dengan kategori referensinya: Clear, Low, Morning, dan Bike)
cat("\nHasil Uji T:\n")
##
## Hasil Uji T:
for (i in 2:nrow(coef_tab)) {
var_name <- rownames(coef_tab)[i]
t_val <- coef_tab[i, 3]
p_val <- coef_tab[i, 4]
sig <- ifelse(p_val < 0.001, "***",
ifelse(p_val < 0.01, "**",
ifelse(p_val < 0.05, "*",
ifelse(p_val < 0.1, ".", "ns"))))
cat(sprintf("%-28s: t = %8.4f, p = %s %s\n",
var_name, t_val, format.pval(p_val, digits = 3, eps = 0.0001), sig))
}
## Distance_km : t = 49.6320, p = <1e-04 ***
## Preparation_Time_min : t = 20.4120, p = <1e-04 ***
## Courier_Experience_yrs : t = -4.7428, p = <1e-04 ***
## WeatherFoggy : t = 6.9167, p = <1e-04 ***
## WeatherRainy : t = 5.5610, p = <1e-04 ***
## WeatherSnowy : t = 8.2330, p = <1e-04 ***
## WeatherWindy : t = 1.7847, p = 0.0747 .
## Traffic_LevelMedium : t = 7.9123, p = <1e-04 ***
## Traffic_LevelHigh : t = 13.1165, p = <1e-04 ***
## Time_of_DayAfternoon : t = 0.8170, p = 0.414 ns
## Time_of_DayEvening : t = 1.5792, p = 0.115 ns
## Time_of_DayNight : t = -0.4664, p = 0.641 ns
## Vehicle_TypeScooter : t = -1.5327, p = 0.126 ns
## Vehicle_TypeCar : t = 0.2951, p = 0.768 ns
cat("\nKeterangan: *** p<0.001, ** p<0.01, * p<0.05, . p<0.1, ns tidak signifikan\n")
##
## Keterangan: *** p<0.001, ** p<0.01, * p<0.05, . p<0.1, ns tidak signifikan
# Uji parsial per VARIABEL (bukan per kategori). Diperlukan karena
# variabel kategorik memiliki lebih dari satu koefisien dummy
cat("\nUji F parsial per variabel (setiap variabel diuji utuh):\n")
##
## Uji F parsial per variabel (setiap variabel diuji utuh):
print(drop1(model, test = "F"))
## Single term deletions
##
## Model:
## Delivery_Time_min ~ Distance_km + Preparation_Time_min + Courier_Experience_yrs +
## Weather + Traffic_Level + Time_of_Day + Vehicle_Type
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 87026 4083.5
## Distance_km 1 246976 334002 5269.1 2463.3352 < 2.2e-16 ***
## Preparation_Time_min 1 41774 128800 4427.7 416.6505 < 2.2e-16 ***
## Courier_Experience_yrs 1 2255 89282 4104.1 22.4942 2.463e-06 ***
## Weather 4 10655 97682 4177.5 26.5688 < 2.2e-16 ***
## Traffic_Level 2 17953 104980 4245.1 89.5330 < 2.2e-16 ***
## Time_of_Day 3 366 87392 4081.2 1.2157 0.3029
## Vehicle_Type 2 303 87329 4082.6 1.5095 0.2216
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
6. Koefisien Determinasi
# ------------------------------------------------------------
# 6. KOEFISIEN DETERMINASI (R² dan Adjusted R²)
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("5. KOEFISIEN DETERMINASI\n")
## 5. KOEFISIEN DETERMINASI
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
r_squared <- summary(model)$r.squared
adj_r_squared <- summary(model)$adj.r.squared
resid_se <- summary(model)$sigma
cat("\nR-squared:", round(r_squared, 4), "\n")
##
## R-squared: 0.7879
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.7845
cat("Residual Standard Error:", round(resid_se, 4), "menit\n")
## Residual Standard Error: 10.013 menit
cat("Interpretasi: Model mampu menjelaskan",
round(r_squared * 100, 2), "% variasi waktu pengantaran\n")
## Interpretasi: Model mampu menjelaskan 78.79 % variasi waktu pengantaran
7. Uji Asumsi Residual
# ------------------------------------------------------------
# 7. UJI ASUMSI RESIDUAL
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("6. UJI ASUMSI RESIDUAL\n")
## 6. UJI ASUMSI RESIDUAL
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Ekstrak residual
resid_model <- residuals(model)
fitted_values <- fitted(model)
standardized_resid <- rstandard(model)
# --- 7.1. UJI NORMALITAS RESIDUAL ---
cat("\n--- 6.1. UJI NORMALITAS RESIDUAL ---\n")
##
## --- 6.1. UJI NORMALITAS RESIDUAL ---
# Pemilihan uji: n = 883 (sampel besar), sehingga dipakai Anderson-Darling.
# Uji ini sensitif pada bagian ekor sebaran, dan kinerjanya membaik pada
# sampel yang lebih besar
# Tidak dipakai: Shapiro-Wilk (umumnya dipakai untuk sampel kecil-sedang),
# Kolmogorov-Smirnov (parameter mean & sd diestimasi dari data sehingga
# p-value tidak valid), dan Jarque-Bera (hasilnya sejalan, tidak diulang).
ad_test <- ad.test(resid_model)
cat("\nAnderson-Darling Test:\n")
##
## Anderson-Darling Test:
cat("H0: Residual berdistribusi normal\n")
## H0: Residual berdistribusi normal
cat("A =", round(ad_test$statistic, 4), "\n")
## A = 33.3302
cat("p-value =", format(ad_test$p.value, scientific = TRUE, digits = 4), "\n")
## p-value = 3.7e-24
cat("Kesimpulan:", ifelse(ad_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"), "\n")
## Kesimpulan: Residual TIDAK berdistribusi normal
# Skewness dan kurtosis sebagai bukti pendukung bentuk sebaran residual
skew_resid <- skewness(resid_model)
kurt_resid <- kurtosis(resid_model)
cat("\nSkewness residual:", round(skew_resid, 4), "(normal = 0)\n")
##
## Skewness residual: 2.5496 (normal = 0)
cat("Kurtosis residual:", round(kurt_resid, 4), "(normal = 3)\n")
## Kurtosis residual: 14.3563 (normal = 3)
# --- 7.2. UJI HETEROSKEDASTISITAS ---
cat("\n--- 6.2. UJI HETEROSKEDASTISITAS ---\n")
##
## --- 6.2. UJI HETEROSKEDASTISITAS ---
# Breusch-Pagan (versi studentized/Koenker pada bptest), cocok untuk
# n besar dan tidak mensyaratkan residual normal
bp_test <- bptest(model)
cat("\nBreusch-Pagan Test:\n")
##
## Breusch-Pagan Test:
cat("H0: Varians residual konstan (homoskedastis)\n")
## H0: Varians residual konstan (homoskedastis)
cat("BP =", round(bp_test$statistic, 4), ", df =", bp_test$parameter, "\n")
## BP = 20.4758 , df = 14
cat("p-value =", round(bp_test$p.value, 4), "\n")
## p-value = 0.1158
cat("Kesimpulan:", ifelse(bp_test$p.value > 0.05,
"Tidak ada heteroskedastisitas",
"Ada heteroskedastisitas"), "\n")
## Kesimpulan: Tidak ada heteroskedastisitas
# --- 7.3. UJI AUTOKORELASI ---
# Tidak dilakukan. Uji Durbin-Watson/Breusch-Godfrey untuk data yang
# berurutan (deret waktu). Data ini cross-section: setiap baris adalah
# pesanan berbeda tanpa urutan waktu, sehingga uji tersebut tidak bermakna.
# --- 7.4. UJI MULTIKOLINEARITAS ---
cat("\n--- 6.4. UJI MULTIKOLINEARITAS ---\n")
##
## --- 6.4. UJI MULTIKOLINEARITAS ---
# Karena ada variabel kategorik (lebih dari 1 df), vif() dari package car
# menghasilkan GVIF. Nilai GVIF^(1/Df) dikuadratkan agar sebanding dengan
# VIF biasa, sehingga batas 5 dan 10 tetap dapat dipakai.
vif_raw <- vif(model)
cat("\nGVIF (keluaran vif):\n")
##
## GVIF (keluaran vif):
print(round(vif_raw, 4))
## GVIF Df GVIF^(1/(2*Df))
## Distance_km 1.0110 1 1.0055
## Preparation_Time_min 1.0086 1 1.0043
## Courier_Experience_yrs 1.0147 1 1.0073
## Weather 1.0538 4 1.0066
## Traffic_Level 1.0274 2 1.0068
## Time_of_Day 1.0307 3 1.0051
## Vehicle_Type 1.0306 2 1.0076
if (is.matrix(vif_raw)) {
vif_values <- vif_raw[, 3]^2
} else {
vif_values <- vif_raw
}
cat("\nInterpretasi VIF:\n")
##
## Interpretasi VIF:
for (i in 1:length(vif_values)) {
vif_val <- vif_values[i]
var_name <- names(vif_values)[i]
status <- ifelse(vif_val < 5, "Tidak ada multikolinearitas",
ifelse(vif_val < 10, "Multikolinearitas moderat",
"Multikolinearitas serius"))
cat(sprintf("%-25s: VIF = %7.4f → %s\n", var_name, vif_val, status))
}
## Distance_km : VIF = 1.0110 → Tidak ada multikolinearitas
## Preparation_Time_min : VIF = 1.0086 → Tidak ada multikolinearitas
## Courier_Experience_yrs : VIF = 1.0147 → Tidak ada multikolinearitas
## Weather : VIF = 1.0132 → Tidak ada multikolinearitas
## Traffic_Level : VIF = 1.0136 → Tidak ada multikolinearitas
## Time_of_Day : VIF = 1.0101 → Tidak ada multikolinearitas
## Vehicle_Type : VIF = 1.0152 → Tidak ada multikolinearitas
# Tolerance
tolerance <- 1 / vif_values
cat("\nTolerance Values:\n")
##
## Tolerance Values:
print(round(tolerance, 4))
## Distance_km Preparation_Time_min Courier_Experience_yrs
## 0.9891 0.9915 0.9855
## Weather Traffic_Level Time_of_Day
## 0.9870 0.9866 0.9900
## Vehicle_Type
## 0.9850
# --- 7.5. UJI LINEARITAS ---
cat("\n--- 6.5. UJI LINEARITAS ---\n")
##
## --- 6.5. UJI LINEARITAS ---
# Ramsey RESET Test
reset_test <- resettest(model, power = 2:3, type = "fitted")
cat("\nRamsey RESET Test:\n")
##
## Ramsey RESET Test:
cat("H0: Model linear (spesifikasi benar)\n")
## H0: Model linear (spesifikasi benar)
cat("F =", round(reset_test$statistic, 4), "\n")
## F = 1.0352
cat("p-value =", round(reset_test$p.value, 4), "\n")
## p-value = 0.3556
cat("Kesimpulan:", ifelse(reset_test$p.value > 0.05,
"Model linear (spesifikasi benar)",
"Model TIDAK linear (spesifikasi salah)"), "\n")
## Kesimpulan: Model linear (spesifikasi benar)
8. Diagnostik Outlier dan Influential Points
# ------------------------------------------------------------
# 8. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("7. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS\n")
## 7. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
p_param <- length(coef(model)) # jumlah parameter (k+1)
# Cook's Distance
cooks_d <- cooks.distance(model)
cat("\nCook's Distance:\n")
##
## Cook's Distance:
cat("Nilai maksimum:", round(max(cooks_d), 4), "\n")
## Nilai maksimum: 0.0501
cat("Jumlah observasi dengan Cook's D > 1:", sum(cooks_d > 1), "\n")
## Jumlah observasi dengan Cook's D > 1: 0
# Leverage (Hat Values)
leverage <- hatvalues(model)
cat("\nLeverage (Hat Values):\n")
##
## Leverage (Hat Values):
cat("Nilai maksimum:", round(max(leverage), 4), "\n")
## Nilai maksimum: 0.0338
cat("Threshold 2(k+1)/n:", round(2 * p_param / n, 4), "\n")
## Threshold 2(k+1)/n: 0.034
cat("Jumlah observasi dengan leverage > threshold:",
sum(leverage > 2 * p_param / n), "\n")
## Jumlah observasi dengan leverage > threshold: 0
# Studentized Residuals
std_resid <- rstudent(model)
cat("\nStudentized Residuals:\n")
##
## Studentized Residuals:
cat("Nilai maksimum absolut:", round(max(abs(std_resid)), 4), "\n")
## Nilai maksimum absolut: 6.3344
cat("Jumlah observasi dengan |std_resid| > 3:",
sum(abs(std_resid) > 3), "\n")
## Jumlah observasi dengan |std_resid| > 3: 21
# DFBETAS
dfbetas_val <- dfbetas(model)
cat("\nDFBETAS:\n")
##
## DFBETAS:
cat("Jumlah nilai dengan |DFBETAS| > 1:", sum(abs(dfbetas_val) > 1), "\n")
## Jumlah nilai dengan |DFBETAS| > 1: 0
9. Visualisasi Diagnostik
# ------------------------------------------------------------
# 9. VISUALISASI DIAGNOSTIK
# ------------------------------------------------------------
# Set par untuk 2x2 plot
par(mfrow = c(2, 2))
# Plot 1: Residual vs Fitted
plot(fitted_values, resid_model,
xlab = "Fitted Values", ylab = "Residuals",
main = "Residual vs Fitted",
pch = 19, col = "steelblue")
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted_values, resid_model), col = "green", lwd = 2)
# Plot 2: QQ-Plot
qqnorm(resid_model, pch = 19, col = "steelblue",
main = "Normal Q-Q Plot")
qqline(resid_model, col = "red", lwd = 2)
# Plot 3: Scale-Location
plot(fitted_values, sqrt(abs(standardized_resid)),
xlab = "Fitted Values", ylab = "√|Standardized Residuals|",
main = "Scale-Location",
pch = 19, col = "steelblue")
lines(lowess(fitted_values, sqrt(abs(standardized_resid))),
col = "green", lwd = 2)
# Plot 4: Residual vs Leverage
plot(leverage, standardized_resid,
xlab = "Leverage", ylab = "Standardized Residuals",
main = "Residual vs Leverage",
pch = 19, col = "steelblue")
abline(h = 0, col = "red", lty = 2)
abline(v = 2 * p_param / n, col = "red", lty = 2)
# Reset par
par(mfrow = c(1, 1))
10. Pemilihan Model (AIC, BIC, AICc dengan Backward dan
Forward)
# ------------------------------------------------------------
# 10. PEMILIHAN MODEL (STEPWISE)
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("8. PEMILIHAN MODEL (STEPWISE)\n")
## 8. PEMILIHAN MODEL (STEPWISE)
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
# Kriteria: AIC -> step(..., k = 2)
# BIC -> step(..., k = log(n))
# AICc tidak tersedia di step(), sehingga dihitung manual dengan rumus
# AICc = AIC + 2K(K+1) / (n - K - 1)
# (K = jumlah parameter yang diestimasi, termasuk varians galat)
aicc <- function(m) {
K <- attr(logLik(m), "df")
AIC(m) + (2 * K * (K + 1)) / (nobs(m) - K - 1)
}
# Model awal untuk forward selection (hanya intersep) dan batas modelnya
model_null <- lm(Delivery_Time_min ~ 1, data = data)
batas <- list(lower = ~ 1, upper = formula(model))
# --- Berdasarkan AIC ---
aic_backward <- step(model, direction = "backward", k = 2, trace = 0)
aic_forward <- step(model_null, scope = batas, direction = "forward", k = 2, trace = 0)
aic_both <- step(model, scope = batas, direction = "both", k = 2, trace = 0)
# --- Berdasarkan BIC ---
bic_backward <- step(model, direction = "backward", k = log(n), trace = 0)
bic_forward <- step(model_null, scope = batas, direction = "forward", k = log(n), trace = 0)
bic_both <- step(model, scope = batas, direction = "both", k = log(n), trace = 0)
# Langkah-langkah seleksi
cat("\nLangkah Backward (AIC):\n"); print(aic_backward$anova)
##
## Langkah Backward (AIC):
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 NA NA 868 87026.41 4083.537
## 2 - Time_of_Day 3 365.6656 871 87392.07 4081.239
## 3 - Vehicle_Type 2 275.0251 873 87667.10 4080.013
cat("\nLangkah Forward (AIC):\n"); print(aic_forward$anova)
##
## Langkah Forward (AIC):
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 NA NA 882 410305.89 5424.797
## 2 + Distance_km -1 251699.413 881 158606.48 4587.526
## 3 + Preparation_Time_min -1 41132.045 880 117474.43 4324.445
## 4 + Traffic_Level -2 16640.469 878 100833.96 4193.570
## 5 + Weather -4 10938.713 874 89895.25 4100.175
## 6 + Courier_Experience_yrs -1 2228.152 873 87667.10 4080.013
cat("\nLangkah Backward (BIC):\n"); print(bic_backward$anova)
##
## Langkah Backward (BIC):
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 NA NA 868 87026.41 4155.286
## 2 - Time_of_Day 3 365.6656 871 87392.07 4138.639
## 3 - Vehicle_Type 2 275.0251 873 87667.10 4127.847
cat("\nLangkah Forward (BIC):\n"); print(bic_forward$anova)
##
## Langkah Forward (BIC):
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 NA NA 882 410305.89 5429.580
## 2 + Distance_km -1 251699.413 881 158606.48 4597.093
## 3 + Preparation_Time_min -1 41132.045 880 117474.43 4338.795
## 4 + Traffic_Level -2 16640.469 878 100833.96 4217.487
## 5 + Weather -4 10938.713 874 89895.25 4143.225
## 6 + Courier_Experience_yrs -1 2228.152 873 87667.10 4127.847
# Perbandingan semua model
daftar_model <- list(
"Model penuh" = model,
"AIC - Backward" = aic_backward,
"AIC - Forward" = aic_forward,
"AIC - Both" = aic_both,
"BIC - Backward" = bic_backward,
"BIC - Forward" = bic_forward,
"BIC - Both" = bic_both
)
tabel_model <- data.frame(
Model = names(daftar_model),
Variabel = sapply(daftar_model, function(m) paste(attr(terms(m), "term.labels"), collapse = " + ")),
K = sapply(daftar_model, function(m) attr(logLik(m), "df")),
AIC = sapply(daftar_model, AIC),
BIC = sapply(daftar_model, BIC),
AICc = sapply(daftar_model, aicc),
Adj_R2 = sapply(daftar_model, function(m) summary(m)$adj.r.squared),
row.names = NULL
)
tabel_model[, c("AIC", "BIC", "AICc", "Adj_R2")] <-
round(tabel_model[, c("AIC", "BIC", "AICc", "Adj_R2")], 4)
cat("\nPerbandingan Model:\n")
##
## Perbandingan Model:
print(tabel_model, right = FALSE)
## Model
## 1 Model penuh
## 2 AIC - Backward
## 3 AIC - Forward
## 4 AIC - Both
## 5 BIC - Backward
## 6 BIC - Forward
## 7 BIC - Both
## Variabel
## 1 Distance_km + Preparation_Time_min + Courier_Experience_yrs + Weather + Traffic_Level + Time_of_Day + Vehicle_Type
## 2 Distance_km + Preparation_Time_min + Courier_Experience_yrs + Weather + Traffic_Level
## 3 Distance_km + Preparation_Time_min + Traffic_Level + Weather + Courier_Experience_yrs
## 4 Distance_km + Preparation_Time_min + Courier_Experience_yrs + Weather + Traffic_Level
## 5 Distance_km + Preparation_Time_min + Courier_Experience_yrs + Weather + Traffic_Level
## 6 Distance_km + Preparation_Time_min + Traffic_Level + Weather + Courier_Experience_yrs
## 7 Distance_km + Preparation_Time_min + Courier_Experience_yrs + Weather + Traffic_Level
## K AIC BIC AICc Adj_R2
## 1 16 6591.382 6667.915 6592.010 0.7845
## 2 11 6587.859 6640.475 6588.162 0.7841
## 3 11 6587.859 6640.475 6588.162 0.7841
## 4 11 6587.859 6640.475 6588.162 0.7841
## 5 11 6587.859 6640.475 6588.162 0.7841
## 6 11 6587.859 6640.475 6588.162 0.7841
## 7 11 6587.859 6640.475 6588.162 0.7841
# Pemilihan AIC vs AICc: AICc disarankan bila n/K < 40 (Burnham & Anderson, 2002)
K_penuh <- attr(logLik(model), "df")
cat("\nn / K (model penuh) =", round(n / K_penuh, 1), "\n")
##
## n / K (model penuh) = 55.2
cat(ifelse(n / K_penuh > 40,
"n/K > 40 → koreksi AICc sangat kecil, AIC dan AICc memberi pilihan yang sama\n",
"n/K < 40 → gunakan AICc sebagai kriteria utama\n"))
## n/K > 40 → koreksi AICc sangat kecil, AIC dan AICc memberi pilihan yang sama
# Model terbaik menurut masing-masing kriteria
idx_aic <- which.min(tabel_model$AIC)
idx_bic <- which.min(tabel_model$BIC)
idx_aicc <- which.min(tabel_model$AICc)
cat("\nAIC terkecil :", tabel_model$Model[idx_aic], "\n")
##
## AIC terkecil : AIC - Backward
cat("BIC terkecil :", tabel_model$Model[idx_bic], "\n")
## BIC terkecil : AIC - Backward
cat("AICc terkecil:", tabel_model$Model[idx_aicc], "\n")
## AICc terkecil: AIC - Backward
# Model akhir (kriteria utama: AICc, sesuai n/K di atas; ganti idx_aic
# jika memilih AIC). Jika AIC, BIC, dan AICc memilih himpunan variabel
# yang sama, maka pilihan model sudah konsisten.
final_model <- daftar_model[[idx_aicc]]
cat("\nModel Terbaik (Stepwise):\n")
##
## Model Terbaik (Stepwise):
print(summary(final_model))
##
## Call:
## lm(formula = Delivery_Time_min ~ Distance_km + Preparation_Time_min +
## Courier_Experience_yrs + Weather + Traffic_Level, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.855 -4.925 -1.191 3.101 62.311
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.03563 1.32542 3.799 0.000155 ***
## Distance_km 2.95856 0.05957 49.667 < 2e-16 ***
## Preparation_Time_min 0.95311 0.04664 20.435 < 2e-16 ***
## Courier_Experience_yrs -0.54740 0.11621 -4.710 2.87e-06 ***
## WeatherFoggy 7.77951 1.12595 6.909 9.39e-12 ***
## WeatherRainy 4.99198 0.88143 5.663 2.01e-08 ***
## WeatherSnowy 10.07569 1.19060 8.463 < 2e-16 ***
## WeatherWindy 2.17772 1.19110 1.828 0.067841 .
## Traffic_LevelMedium 5.92850 0.76073 7.793 1.85e-14 ***
## Traffic_LevelHigh 12.03682 0.92188 13.057 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.02 on 873 degrees of freedom
## Multiple R-squared: 0.7863, Adjusted R-squared: 0.7841
## F-statistic: 357 on 9 and 873 DF, p-value: < 2.2e-16
11. Ringkasan Hasil
# ------------------------------------------------------------
# 11. RINGKASAN HASIL
# ------------------------------------------------------------
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("9. RINGKASAN HASIL ANALISIS\n")
## 9. RINGKASAN HASIL ANALISIS
cat(strrep("=", 60), "\n", sep = "")
## ============================================================
cat("\n--- MODEL AKHIR ---\n")
##
## --- MODEL AKHIR ---
koef_final <- coef(final_model)
persamaan <- paste0(
"Y = ", round(koef_final[1], 4),
paste0(sprintf(" %s %.4f*%s",
ifelse(koef_final[-1] >= 0, "+", "-"),
abs(koef_final[-1]), names(koef_final)[-1]),
collapse = "")
)
cat("Persamaan Regresi:\n")
## Persamaan Regresi:
cat(persamaan, "\n")
## Y = 5.0356 + 2.9586*Distance_km + 0.9531*Preparation_Time_min - 0.5474*Courier_Experience_yrs + 7.7795*WeatherFoggy + 4.9920*WeatherRainy + 10.0757*WeatherSnowy + 2.1777*WeatherWindy + 5.9285*Traffic_LevelMedium + 12.0368*Traffic_LevelHigh
cat("\n--- UJI SIGNIFIKANSI (MODEL PENUH) ---\n")
##
## --- UJI SIGNIFIKANSI (MODEL PENUH) ---
cat("Uji F (Simultan): p-value =", format(f_pvalue, scientific = TRUE), "\n")
## Uji F (Simultan): p-value = 1.226723e-280
cat("\n--- UJI ASUMSI RESIDUAL (MODEL PENUH) ---\n")
##
## --- UJI ASUMSI RESIDUAL (MODEL PENUH) ---
cat("Normalitas (Anderson-Darling): p-value =",
format(ad_test$p.value, scientific = TRUE, digits = 4), "\n")
## Normalitas (Anderson-Darling): p-value = 3.7e-24
cat("Heteroskedastisitas (Breusch-Pagan): p-value =",
round(bp_test$p.value, 4), "\n")
## Heteroskedastisitas (Breusch-Pagan): p-value = 0.1158
cat("Multikolinearitas (VIF max):", round(max(vif_values), 4), "\n")
## Multikolinearitas (VIF max): 1.0152
cat("Linearitas (Ramsey RESET): p-value =",
round(reset_test$p.value, 4), "\n")
## Linearitas (Ramsey RESET): p-value = 0.3556
cat("\n--- KESIMPULAN ---\n")
##
## --- KESIMPULAN ---
cat("Model penuh : R-squared =", round(r_squared, 4),
", Adjusted R-squared =", round(adj_r_squared, 4), "\n")
## Model penuh : R-squared = 0.7879 , Adjusted R-squared = 0.7845
cat("Model akhir : R-squared =", round(summary(final_model)$r.squared, 4),
", Adjusted R-squared =", round(summary(final_model)$adj.r.squared, 4), "\n")
## Model akhir : R-squared = 0.7863 , Adjusted R-squared = 0.7841
cat("\n", strrep("=", 60), "\n", sep = "")
##
## ============================================================
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
cat(strrep("=", 60), "\n", sep = "")
## ============================================================