# ============================================================
# ANALISIS MULTIPLE LINEAR REGRESSION (MLR) LENGKAP
# Mulai dari Estimasi hingga Uji Hipotesis & Residual
# ============================================================
# ------------------------------------------------------------
# 0. PERSIAPAN: LOAD LIBRARY DAN DATA
# ------------------------------------------------------------
# Install package jika belum ada (uncomment jika perlu)
# install.packages(c("car", "lmtest", "nortest", "MASS",
# "olsrr", "ggplot2", "corrplot", "moments"))
# Load library
library(car) # Untuk VIF, Durbin-Watson, dll.
## 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, Durbin-Watson
## 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(MASS) # Untuk Box-Cox transformation
## Warning: package 'MASS' was built under R version 4.5.3
library(olsrr) # Untuk uji asumsi OLS lengkap
## Warning: package 'olsrr' was built under R version 4.5.3
##
## Attaching package: 'olsrr'
## The following object is masked from 'package:MASS':
##
## cement
## The following object is masked from 'package:datasets':
##
## rivers
library(ggplot2) # Untuk visualisasi
## Warning: package 'ggplot2' was built under R version 4.5.3
library(corrplot) # Untuk matriks korelasi
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
# Set seed untuk reproduktifitas
set.seed(123)
# Set seed untuk reproduktifitas
set.seed(123)
# ------------------------------------------------------------
# 1. MEMBUAT DATA SIMULASI (Ganti dengan data Anda)
# ------------------------------------------------------------
# Simulasi data: 100 observasi, 4 variabel independen
n <- 100
X1 <- rnorm(n, mean = 50, sd = 10)
X2 <- rnorm(n, mean = 70, sd = 15)
X3 <- rnorm(n, mean = 100, sd = 20)
X4 <- rnorm(n, mean = 60, sd = 8)
X5 <- runif(n, min = 5, max = 100) # non-normal (distribusi rata/flat)
# Model:Y <- 5 + 0.5*X1 + 0.3*X2 - 0.2*X3 + 0.1*X4 + 0.4*X5 + rnorm(n, mean = 0, sd = 3)
Y <- 5 + 0.5*X1 + 0.3*X2 - 0.2*X3 + 0.1*X4 + 0.4*X5 + rnorm(n, mean = 0, sd = 3)
# Buat data frame
data <- data.frame(Y = Y, X1 = X1, X2 = X2, X3 = X3, X4 = X4, X5 =X5)
# Tampilkan 6 baris pertama
cat("=" %+% rep("=", 60) %+% "\n")
## Warning: <ggplot> %+% x was deprecated in ggplot2 4.0.0.
## ℹ Please use <ggplot> + x instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
cat("DATA YANG DIGUNAKAN\n")
## DATA YANG DIGUNAKAN
cat("=" %+% rep("=", 60) %+% "\n")
print(head(data))
## Y X1 X2 X3 X4 X5
## 1 45.81047 44.39524 59.34390 143.97621 54.27806 49.71477
## 2 50.19530 47.69823 73.85326 126.24826 53.97849 39.75532
## 3 51.90738 65.58708 66.29962 94.69710 52.49169 16.52085
## 4 38.70511 50.70508 64.78686 110.86388 51.57989 9.46440
## 5 49.40957 51.29288 55.72572 91.71320 56.50272 29.96565
## 6 78.36823 67.15065 69.32458 90.47506 62.64943 97.02091
# ------------------------------------------------------------
# 2. STATISTIK DESKRIPTIF DAN EKSPLORASI DATA
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("1. STATISTIK DESKRIPTIF\n")
## 1. STATISTIK DESKRIPTIF
cat(rep("=", 60) %+% "\n")
# Statistik deskriptif
summary(data)
## Y X1 X2 X3
## Min. :27.89 Min. :26.91 Min. : 39.20 Min. : 64.87
## 1st Qu.:45.75 1st Qu.:45.06 1st Qu.: 57.98 1st Qu.: 89.37
## Median :60.64 Median :50.62 Median : 66.61 Median :100.72
## Mean :57.36 Mean :50.90 Mean : 68.39 Mean :102.41
## 3rd Qu.:67.92 3rd Qu.:56.92 3rd Qu.: 77.02 3rd Qu.:115.27
## Max. :82.08 Max. :71.87 Max. :118.62 Max. :145.86
## X4 X5
## Min. :40.27 Min. : 5.39
## 1st Qu.:54.16 1st Qu.:27.13
## Median :59.97 Median :54.17
## Mean :59.71 Mean :51.12
## 3rd Qu.:65.51 3rd Qu.:73.61
## Max. :80.57 Max. :99.22
# Matriks korelasi
cat("\nMatriks Korelasi:\n")
##
## Matriks Korelasi:
cor_matrix <- cor(data)
print(round(cor_matrix, 3))
## Y X1 X2 X3 X4 X5
## Y 1.000 0.279 0.320 -0.334 0.225 0.804
## X1 0.279 1.000 -0.050 -0.129 -0.044 -0.081
## X2 0.320 -0.050 1.000 0.031 0.044 0.036
## X3 -0.334 -0.129 0.031 1.000 -0.045 0.007
## X4 0.225 -0.044 0.044 -0.045 1.000 0.162
## X5 0.804 -0.081 0.036 0.007 0.162 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", mar = c(0,0,2,0))
# Scatter plot matrix
pairs(data, main = "Scatter Plot Matrix", pch = 19, col = "steelblue")
# ------------------------------------------------------------
# 3. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("2. ESTIMASI MODEL MLR\n")
## 2. ESTIMASI MODEL MLR
cat(rep("=", 60) %+% "\n")
# Fit model MLR
model <- lm(Y ~ X1 + X2 + X3 + X4 + X5, data = data)
# Ringkasan model
cat("\nRingkasan Model:\n")
##
## Ringkasan Model:
print(summary(model))
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5771 -1.7043 -0.2487 1.9458 6.3534
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.07586 3.42130 2.360 0.020318 *
## X1 0.47059 0.03083 15.262 < 2e-16 ***
## X2 0.28564 0.01920 14.875 < 2e-16 ***
## X3 -0.21155 0.01477 -14.324 < 2e-16 ***
## X4 0.12850 0.03396 3.784 0.000271 ***
## X5 0.38710 0.01021 37.923 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.764 on 94 degrees of freedom
## Multiple R-squared: 0.9587, Adjusted R-squared: 0.9565
## F-statistic: 436.8 on 5 and 94 DF, p-value: < 2.2e-16
# Ekstrak koefisien
cat("\nKoefisien Regresi:\n")
##
## Koefisien Regresi:
coef_df <- data.frame(
Variabel = names(coef(model)),
Koefisien = round(coef(model), 4),
Std_Error = round(summary(model)$coefficients[, 2], 4),
t_value = round(summary(model)$coefficients[, 3], 4),
p_value = round(summary(model)$coefficients[, 4], 4)
)
print(coef_df)
## Variabel Koefisien Std_Error t_value p_value
## (Intercept) (Intercept) 8.0759 3.4213 2.3605 0.0203
## X1 X1 0.4706 0.0308 15.2618 0.0000
## X2 X2 0.2856 0.0192 14.8754 0.0000
## X3 X3 -0.2115 0.0148 -14.3242 0.0000
## X4 X4 0.1285 0.0340 3.7842 0.0003
## X5 X5 0.3871 0.0102 37.9228 0.0000
# 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.2828 14.8689
## X1 0.4094 0.5318
## X2 0.2475 0.3238
## X3 -0.2409 -0.1822
## X4 0.0611 0.1959
## X5 0.3668 0.4074
# ------------------------------------------------------------
# 4. UJI SIGNIFIKANSI SIMULTAN (UJI F)
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("3. UJI SIGNIFIKANSI SIMULTAN (UJI F)\n")
## 3. UJI SIGNIFIKANSI SIMULTAN (UJI F)
cat(rep("=", 60) %+% "\n")
# Ekstrak informasi F-test
f_stat <- summary(model)$fstatistic
f_value <- f_stat[1]
df1 <- f_stat[2]
df2 <- f_stat[3]
f_pvalue <- pf(f_value, df1, df2, lower.tail = FALSE)
cat("\nHipotesis:\n")
##
## Hipotesis:
cat("H0: β1 = β2 = β3 = β4 = β5 = 0 (Model tidak signifikan)\n")
## H0: β1 = β2 = β3 = β4 = β5 = 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: 436.7941
cat("df1:", df1, "\n")
## df1: 5
cat("df2:", df2, "\n")
## df2: 94
cat("p-value:", format(f_pvalue, scientific = TRUE), "\n")
## p-value: 2.027572e-63
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)
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("4. UJI SIGNIFIKANSI PARSIAL (UJI T)\n")
## 4. UJI SIGNIFIKANSI PARSIAL (UJI T)
cat(rep("=", 60) %+% "\n")
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("\nHasil Uji T:\n")
##
## Hasil Uji T:
for (i in 2:nrow(coef_df)) {
var_name <- coef_df$Variabel[i]
t_val <- coef_df$t_value[i]
p_val <- coef_df$p_value[i]
sig <- ifelse(p_val < 0.001, "***",
ifelse(p_val < 0.01, "**",
ifelse(p_val < 0.05, "*",
ifelse(p_val < 0.1, ".", "ns"))))
cat(sprintf("%-10s: t = %7.4f, p = %8.4f %s\n",
var_name, t_val, p_val, sig))
}
## X1 : t = 15.2618, p = 0.0000 ***
## X2 : t = 14.8754, p = 0.0000 ***
## X3 : t = -14.3242, p = 0.0000 ***
## X4 : t = 3.7842, p = 0.0003 ***
## X5 : t = 37.9228, p = 0.0000 ***
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
# ------------------------------------------------------------
# 6. KOEFISIEN DETERMINASI (R² dan Adjusted R²)
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("5. KOEFISIEN DETERMINASI\n")
## 5. KOEFISIEN DETERMINASI
cat(rep("=", 60) %+% "\n")
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.9587
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.9565
cat("Residual Standard Error:", round(resid_se, 4), "\n")
## Residual Standard Error: 2.7636
cat("Interpretasi: Model mampu menjelaskan",
round(r_squared * 100, 2), "% variasi pada Y\n")
## Interpretasi: Model mampu menjelaskan 95.87 % variasi pada Y
# ------------------------------------------------------------
# 7. UJI ASUMSI RESIDUAL
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("6. UJI ASUMSI RESIDUAL\n")
## 6. UJI ASUMSI RESIDUAL
cat(rep("=", 60) %+% "\n")
# Ekstrak residual
residuals <- 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 ---
# Shapiro-Wilk Test
shapiro_test <- shapiro.test(residuals)
cat("\nShapiro-Wilk Test:\n")
##
## Shapiro-Wilk Test:
cat("W =", round(shapiro_test$statistic, 4), "\n")
## W = 0.9928
cat("p-value =", round(shapiro_test$p.value, 4), "\n")
## p-value = 0.8753
cat("Kesimpulan:", ifelse(shapiro_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"), "\n")
## Kesimpulan: Residual berdistribusi normal
# Kolmogorov-Smirnov Test
ks_test <- ks.test(residuals, "pnorm", mean(residuals), sd(residuals))
cat("\nKolmogorov-Smirnov Test:\n")
##
## Kolmogorov-Smirnov Test:
cat("D =", round(ks_test$statistic, 4), "\n")
## D = 0.0595
cat("p-value =", round(ks_test$p.value, 4), "\n")
## p-value = 0.8707
cat("Kesimpulan:", ifelse(ks_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"), "\n")
## Kesimpulan: Residual berdistribusi normal
# Anderson-Darling Test
ad_test <- ad.test(residuals)
cat("\nAnderson-Darling Test:\n")
##
## Anderson-Darling Test:
cat("A =", round(ad_test$statistic, 4), "\n")
## A = 0.2243
cat("p-value =", round(ad_test$p.value, 4), "\n")
## p-value = 0.8187
cat("Kesimpulan:", ifelse(ad_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"), "\n")
## Kesimpulan: Residual berdistribusi normal
# Jarque-Bera Test (dari package moments)
library(moments)
## Warning: package 'moments' was built under R version 4.5.2
jb_test <- jarque.test(residuals)
cat("\nJarque-Bera Test:\n")
##
## Jarque-Bera Test:
cat("JB =", round(jb_test$statistic, 4), "\n")
## JB = 0.3151
cat("p-value =", round(jb_test$p.value, 4), "\n")
## p-value = 0.8542
cat("Kesimpulan:", ifelse(jb_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"), "\n")
## Kesimpulan: Residual berdistribusi normal
# --- 7.2. UJI HETEROSKEDASTISITAS ---
bp_test <- bptest(model)
glejser_test <- bptest(model, ~ fitted(model))
cat("Breusch-Pagan p-value:", round(bp_test$p.value, 4), "\n")
## Breusch-Pagan p-value: 0.8599
cat("Glejser p-value:", round(glejser_test$p.value, 4), "\n")
## Glejser p-value: 0.8408
# Kesimpulan
if (bp_test$p.value > 0.05 & glejser_test$p.value > 0.05) {
cat("Kesimpulan: Tidak ada heteroskedastisitas\n")
} else {
cat("Kesimpulan: Ada heteroskedastisitas\n")
}
## Kesimpulan: Tidak ada heteroskedastisitas
# --- 7.3. UJI AUTOKORELASI ---
cat("\n--- 6.3. UJI AUTOKORELASI ---\n")
##
## --- 6.3. UJI AUTOKORELASI ---
# Durbin-Watson Test
dw_test <- dwtest(model)
cat("\nDurbin-Watson Test:\n")
##
## Durbin-Watson Test:
cat("DW =", round(dw_test$statistic, 4), "\n")
## DW = 1.983
cat("p-value =", round(dw_test$p.value, 4), "\n")
## p-value = 0.4705
cat("Kesimpulan:",
ifelse(dw_test$statistic < 1.5, "Ada autokorelasi positif",
ifelse(dw_test$statistic > 2.5, "Ada autokorelasi negatif",
"Tidak ada autokorelasi")), "\n")
## Kesimpulan: Tidak ada autokorelasi
# Breusch-Godfrey Test
bg_test <- bgtest(model, order = 1)
cat("\nBreusch-Godfrey Test (Lag 1):\n")
##
## Breusch-Godfrey Test (Lag 1):
cat("LM =", round(bg_test$statistic, 4), "\n")
## LM = 0.0039
cat("p-value =", round(bg_test$p.value, 4), "\n")
## p-value = 0.9502
cat("Kesimpulan:", ifelse(bg_test$p.value > 0.05,
"Tidak ada autokorelasi",
"Ada autokorelasi"), "\n")
## Kesimpulan: Tidak ada autokorelasi
# --- 7.4. UJI MULTIKOLINEARITAS ---
cat("\n--- 6.4. UJI MULTIKOLINEARITAS ---\n")
##
## --- 6.4. UJI MULTIKOLINEARITAS ---
# VIF (Variance Inflation Factor)
vif_values <- vif(model)
cat("\nVIF Values:\n")
##
## VIF Values:
print(round(vif_values, 4))
## X1 X2 X3 X4 X5
## 1.0269 1.0056 1.0204 1.0322 1.0334
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("%-10s: VIF = %7.4f → %s\n", var_name, vif_val, status))
}
## X1 : VIF = 1.0269 → Tidak ada multikolinearitas
## X2 : VIF = 1.0056 → Tidak ada multikolinearitas
## X3 : VIF = 1.0204 → Tidak ada multikolinearitas
## X4 : VIF = 1.0322 → Tidak ada multikolinearitas
## X5 : VIF = 1.0334 → Tidak ada multikolinearitas
# Tolerance
tolerance <- 1 / vif_values
cat("\nTolerance Values:\n")
##
## Tolerance Values:
print(round(tolerance, 4))
## X1 X2 X3 X4 X5
## 0.9738 0.9945 0.9800 0.9688 0.9676
# Condition Number
cat("\nCondition Number:", round(kappa(model), 4), "\n")
##
## Condition Number: 734.3349
cat("(Condition Number > 30 mengindikasikan multikolinearitas)\n")
## (Condition Number > 30 mengindikasikan multikolinearitas)
# --- 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("F =", round(reset_test$statistic, 4), "\n")
## F = 3.056
cat("p-value =", round(reset_test$p.value, 4), "\n")
## p-value = 0.0519
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
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("7. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS\n")
## 7. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
cat(rep("=", 60) %+% "\n")
# 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.1131
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.162
cat("Threshold 2(k+1)/n:", round(2 * (length(coef(model))) / n, 4), "\n")
## Threshold 2(k+1)/n: 0.12
cat("Jumlah observasi dengan leverage > threshold:",
sum(leverage > 2 * length(coef(model)) / n), "\n")
## Jumlah observasi dengan leverage > threshold: 4
# 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: 2.97
cat("Jumlah observasi dengan |std_resid| > 3:",
sum(abs(std_resid) > 3), "\n")
## Jumlah observasi dengan |std_resid| > 3: 0
# DFBETAS
dfbetas_val <- dfbetas(model)
cat("\nDFBETAS:\n")
##
## DFBETAS:
cat("Jumlah observasi dengan |DFBETAS| > 1:",
sum(abs(dfbetas_val) > 1), "\n")
## Jumlah observasi dengan |DFBETAS| > 1: 0
# ------------------------------------------------------------
# 9. VISUALISASI DIAGNOSTIK
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("8. VISUALISASI DIAGNOSTIK\n")
## 8. VISUALISASI DIAGNOSTIK
cat(rep("=", 60) %+% "\n")
# Set par untuk 2x2 plot
par(mfrow = c(2, 2))
# Plot 1: Residual vs Fitted
plot(fitted_values, residuals,
xlab = "Fitted Values", ylab = "Residuals",
main = "Residual vs Fitted",
pch = 19, col = "steelblue")
abline(h = 0, col = "red", lty = 2)
lines(lowess(fitted_values, residuals), col = "green", lwd = 2)
# Plot 2: QQ-Plot
qqnorm(residuals, pch = 19, col = "steelblue",
main = "Normal Q-Q Plot")
qqline(residuals, 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 * length(coef(model)) / n, col = "red", lty = 2)
# Reset par
par(mfrow = c(1, 1))
# ------------------------------------------------------------
# 10. UJI BOX-COX UNTUK TRANSFORMASI
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("9. UJI BOX-COX\n")
## 9. UJI BOX-COX
cat(rep("=", 60) %+% "\n")
# Box-Cox transformation
bc <- boxcox(model, lambda = seq(-2, 2, 0.1))
lambda_opt <- bc$x[which.max(bc$y)]
cat("\nOptimal Lambda:", round(lambda_opt, 4), "\n")
##
## Optimal Lambda: 1.1919
cat("Interpretasi:\n")
## Interpretasi:
cat("- Jika λ ≈ 1: Tidak perlu transformasi\n")
## - Jika λ ≈ 1: Tidak perlu transformasi
cat("- Jika λ ≈ 0: Gunakan log transformation\n")
## - Jika λ ≈ 0: Gunakan log transformation
cat("- Jika λ ≈ 0.5: Gunakan square root transformation\n")
## - Jika λ ≈ 0.5: Gunakan square root transformation
# ------------------------------------------------------------
# 11. PEMILIHAN MODEL (STEPWISE)
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("10. PEMILIHAN MODEL (STEPWISE)\n")
## 10. PEMILIHAN MODEL (STEPWISE)
cat(rep("=", 60) %+% "\n")
# Stepwise selection
step_model <- step(model, direction = "both", trace = 1)
## Start: AIC=209.12
## Y ~ X1 + X2 + X3 + X4 + X5
##
## Df Sum of Sq RSS AIC
## <none> 717.9 209.12
## - X4 1 109.4 827.3 221.30
## - X3 1 1567.1 2285.0 322.89
## - X2 1 1690.0 2407.9 328.13
## - X1 1 1778.9 2496.8 331.76
## - X5 1 10983.6 11701.5 486.23
cat("\nModel Terbaik (Stepwise):\n")
##
## Model Terbaik (Stepwise):
print(summary(step_model))
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5771 -1.7043 -0.2487 1.9458 6.3534
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.07586 3.42130 2.360 0.020318 *
## X1 0.47059 0.03083 15.262 < 2e-16 ***
## X2 0.28564 0.01920 14.875 < 2e-16 ***
## X3 -0.21155 0.01477 -14.324 < 2e-16 ***
## X4 0.12850 0.03396 3.784 0.000271 ***
## X5 0.38710 0.01021 37.923 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.764 on 94 degrees of freedom
## Multiple R-squared: 0.9587, Adjusted R-squared: 0.9565
## F-statistic: 436.8 on 5 and 94 DF, p-value: < 2.2e-16
# ------------------------------------------------------------
# 12. RINGKASAN HASIL
# ------------------------------------------------------------
cat("\n" %+% rep("=", 60) %+% "\n")
cat("11. RINGKASAN HASIL ANALISIS\n")
## 11. RINGKASAN HASIL ANALISIS
cat(rep("=", 60) %+% "\n")
cat("\n--- MODEL AKHIR ---\n")
##
## --- MODEL AKHIR ---
cat("Persamaan Regresi:\n")
## Persamaan Regresi:
cat(sprintf("Y = %.4f + %.4f*X1 + %.4f*X2 + %.4f*X3 + %.4f*X4 + %.4f*X5\n",
coef(model)[1], coef(model)[2], coef(model)[3],
coef(model)[4], coef(model)[5], coef(model)[6]))
## Y = 8.0759 + 0.4706*X1 + 0.2856*X2 + -0.2115*X3 + 0.1285*X4 + 0.3871*X5
cat("\n--- UJI SIGNIFIKANSI ---\n")
##
## --- UJI SIGNIFIKANSI ---
cat("Uji F (Simultan): p-value =", format(f_pvalue, scientific = TRUE), "\n")
## Uji F (Simultan): p-value = 2.027572e-63
cat("Uji t (Parsial):\n")
## Uji t (Parsial):
for (i in 2:nrow(coef_df)) {
cat(sprintf(" %s: p-value = %.4f\n",
coef_df$Variabel[i], coef_df$p_value[i]))
}
## X1: p-value = 0.0000
## X2: p-value = 0.0000
## X3: p-value = 0.0000
## X4: p-value = 0.0003
## X5: p-value = 0.0000
cat("\n--- UJI ASUMSI RESIDUAL ---\n")
##
## --- UJI ASUMSI RESIDUAL ---
cat("Normalitas (Shapiro-Wilk): p-value =",
round(shapiro_test$p.value, 4), "\n")
## Normalitas (Shapiro-Wilk): p-value = 0.8753
cat("Heteroskedastisitas (Breusch-Pagan): p-value =",
round(bp_test$p.value, 4), "\n")
## Heteroskedastisitas (Breusch-Pagan): p-value = 0.8599
cat("Autokorelasi (Durbin-Watson): DW =",
round(dw_test$statistic, 4), "\n")
## Autokorelasi (Durbin-Watson): DW = 1.983
cat("Multikolinearitas (VIF max):",
round(max(vif_values), 4), "\n")
## Multikolinearitas (VIF max): 1.0334
cat("Linearitas (Ramsey RESET): p-value =",
round(reset_test$p.value, 4), "\n")
## Linearitas (Ramsey RESET): p-value = 0.0519
cat("\n--- KESIMPULAN ---\n")
##
## --- KESIMPULAN ---
cat("R-squared:", round(r_squared, 4), "\n")
## R-squared: 0.9587
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.9565
cat("\n" %+% rep("=", 60) %+% "\n")
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
cat(rep("=", 60) %+% "\n")