# ============================================================
# 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.
## Loading required package: carData
library(lmtest) # Untuk uji Breusch-Pagan, Durbin-Watson
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(nortest) # Untuk uji normalitas (Anderson-Darling)
library(MASS) # Untuk Box-Cox transformation
library(olsrr) # Untuk uji asumsi OLS lengkap
##
## 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
library(corrplot) # Untuk matriks korelasi
## corrplot 0.95 loaded
# Set seed untuk reproduktifitas
set.seed(0701)
# ------------------------------------------------------------
# 1. MEMBUAT DATA SIMULASI (Ganti dengan data Anda)
# ------------------------------------------------------------
# Simulasi data: 100 observasi, 5 variabel independen
n <- 100
X1 <- rnorm(n, mean = 50, sd = 10)
X2 <- rnorm(n, mean = 30, sd = 5)
X3 <- rnorm(n, mean = 20, sd = 3)
X4 <- rnorm(n, mean = 10, sd = 2)
X5 <- rchisq(n, df = 5)
# Model: Y = 5 + 0.5*X1 + 0.3*X2 - 0.2*X3 + 0.1*X4 + 0.4*X5 + error
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 35.33329 51.48466 23.22128 17.83111 11.524209 1.981234
## 2 41.82149 51.45661 37.92569 21.68509 8.637324 11.629715
## 3 37.69388 48.29904 31.13409 21.53040 10.298892 3.000536
## 4 29.42875 36.36672 30.07415 23.33477 11.603805 1.403770
## 5 39.02744 57.70044 32.69679 17.12079 8.879617 3.675208
## 6 32.14399 45.28521 36.69981 16.47896 7.177080 2.389776
# ------------------------------------------------------------
# 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. :22.95 Min. :29.56 Min. :17.05 Min. :13.57
## 1st Qu.:34.36 1st Qu.:42.27 1st Qu.:25.28 1st Qu.:17.50
## Median :38.31 Median :48.59 Median :29.83 Median :20.07
## Mean :37.78 Mean :49.58 Mean :29.74 Mean :19.83
## 3rd Qu.:41.00 3rd Qu.:55.55 3rd Qu.:33.46 3rd Qu.:21.96
## Max. :55.84 Max. :84.89 Max. :40.65 Max. :26.25
## X4 X5
## Min. : 6.063 Min. : 0.2871
## 1st Qu.: 8.494 1st Qu.: 2.7276
## Median : 9.722 Median : 3.6594
## Mean : 9.907 Mean : 4.9389
## 3rd Qu.:11.268 3rd Qu.: 5.8853
## Max. :17.314 Max. :19.4463
# 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.789 0.128 -0.063 0.065 0.170
## X1 0.789 1.000 -0.155 -0.086 0.020 -0.054
## X2 0.128 -0.155 1.000 0.057 -0.143 0.070
## X3 -0.063 -0.086 0.057 1.000 -0.027 0.049
## X4 0.065 0.020 -0.143 -0.027 1.000 -0.087
## X5 0.170 -0.054 0.070 0.049 -0.087 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
## -5.9637 -2.1441 -0.0045 2.4474 8.5030
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29507 3.93588 0.075 0.940399
## X1 0.49582 0.03191 15.536 < 2e-16 ***
## X2 0.29022 0.06100 4.757 7.07e-06 ***
## X3 -0.02721 0.10950 -0.248 0.804293
## X4 0.31590 0.16476 1.917 0.058235 .
## X5 0.34099 0.08812 3.870 0.000201 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.198 on 94 degrees of freedom
## Multiple R-squared: 0.7359, Adjusted R-squared: 0.7218
## F-statistic: 52.37 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) 0.2951 3.9359 0.0750 0.9404
## X1 X1 0.4958 0.0319 15.5363 0.0000
## X2 X2 0.2902 0.0610 4.7574 0.0000
## X3 X3 -0.0272 0.1095 -0.2485 0.8043
## X4 X4 0.3159 0.1648 1.9173 0.0582
## X5 X5 0.3410 0.0881 3.8696 0.0002
# 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) -7.5197 8.1098
## X1 0.4325 0.5592
## X2 0.1691 0.4113
## X3 -0.2446 0.1902
## X4 -0.0112 0.6430
## X5 0.1660 0.5159
# ------------------------------------------------------------
# 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: 52.3732
cat("df1:", df1, "\n")
## df1: 5
cat("df2:", df2, "\n")
## df2: 94
cat("p-value:", format(f_pvalue, scientific = TRUE), "\n")
## p-value: 1.079434e-25
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.5363, p = 0.0000 ***
## X2 : t = 4.7574, p = 0.0000 ***
## X3 : t = -0.2485, p = 0.8043 ns
## X4 : t = 1.9173, p = 0.0582 .
## X5 : t = 3.8696, p = 0.0002 ***
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.7359
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.7218
cat("Residual Standard Error:", round(resid_se, 4), "\n")
## Residual Standard Error: 3.198
cat("Interpretasi: Model mampu menjelaskan",
round(r_squared * 100, 2), "% variasi pada Y\n")
## Interpretasi: Model mampu menjelaskan 73.59 % 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.9851
cat("p-value =", round(shapiro_test$p.value, 4), "\n")
## p-value = 0.3213
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.063
cat("p-value =", round(ks_test$p.value, 4), "\n")
## p-value = 0.8223
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.3477
cat("p-value =", round(ad_test$p.value, 4), "\n")
## p-value = 0.472
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)
jb_test <- jarque.test(residuals)
cat("\nJarque-Bera Test:\n")
##
## Jarque-Bera Test:
cat("JB =", round(jb_test$statistic, 4), "\n")
## JB = 0.9572
cat("p-value =", round(jb_test$p.value, 4), "\n")
## p-value = 0.6197
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.1503
cat("Glejser p-value:", round(glejser_test$p.value, 4), "\n")
## Glejser p-value: 0.0362
# 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: 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 = 2.1382
cat("p-value =", round(dw_test$p.value, 4), "\n")
## p-value = 0.7481
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.5748
cat("p-value =", round(bg_test$p.value, 4), "\n")
## p-value = 0.4484
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.0328 1.0499 1.0116 1.0272 1.0147
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.0328 → Tidak ada multikolinearitas
## X2 : VIF = 1.0499 → Tidak ada multikolinearitas
## X3 : VIF = 1.0116 → Tidak ada multikolinearitas
## X4 : VIF = 1.0272 → Tidak ada multikolinearitas
## X5 : VIF = 1.0147 → Tidak ada multikolinearitas
# Tolerance
tolerance <- 1 / vif_values
cat("\nTolerance Values:\n")
##
## Tolerance Values:
print(round(tolerance, 4))
## X1 X2 X3 X4 X5
## 0.9683 0.9524 0.9886 0.9735 0.9855
# Condition Number
cat("\nCondition Number:", round(kappa(model), 4), "\n")
##
## Condition Number: 357.0637
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 = 0.254
cat("p-value =", round(reset_test$p.value, 4), "\n")
## p-value = 0.7763
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.0665
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.1994
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: 7
# 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.8123
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.3939
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=238.32
## Y ~ X1 + X2 + X3 + X4 + X5
##
## Df Sum of Sq RSS AIC
## - X3 1 0.63 962.0 236.38
## <none> 961.3 238.32
## - X4 1 37.60 998.9 240.15
## - X5 1 153.13 1114.5 251.10
## - X2 1 231.46 1192.8 257.89
## - X1 1 2468.56 3429.9 363.51
##
## Step: AIC=236.38
## Y ~ X1 + X2 + X4 + X5
##
## Df Sum of Sq RSS AIC
## <none> 962.0 236.38
## - X4 1 37.76 999.7 238.23
## + X3 1 0.63 961.3 238.32
## - X5 1 152.59 1114.6 249.10
## - X2 1 230.87 1192.8 255.89
## - X1 1 2489.30 3451.3 362.13
cat("\nModel Terbaik (Stepwise):\n")
##
## Model Terbaik (Stepwise):
print(summary(step_model))
##
## Call:
## lm(formula = Y ~ X1 + X2 + X4 + X5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0194 -2.1173 -0.0634 2.4878 8.5657
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.25879 3.22791 -0.080 0.936268
## X1 0.49643 0.03166 15.679 < 2e-16 ***
## X2 0.28962 0.06065 4.775 6.52e-06 ***
## X4 0.31654 0.16392 1.931 0.056464 .
## X5 0.34010 0.08761 3.882 0.000191 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.182 on 95 degrees of freedom
## Multiple R-squared: 0.7357, Adjusted R-squared: 0.7246
## F-statistic: 66.1 on 4 and 95 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 = 0.2951 + 0.4958*X1 + 0.2902*X2 + -0.0272*X3 + 0.3159*X4 + 0.3410*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 = 1.079434e-25
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.8043
## X4: p-value = 0.0582
## X5: p-value = 0.0002
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.3213
cat("Heteroskedastisitas (Breusch-Pagan): p-value =",
round(bp_test$p.value, 4), "\n")
## Heteroskedastisitas (Breusch-Pagan): p-value = 0.1503
cat("Autokorelasi (Durbin-Watson): DW =",
round(dw_test$statistic, 4), "\n")
## Autokorelasi (Durbin-Watson): DW = 2.1382
cat("Multikolinearitas (VIF max):",
round(max(vif_values), 4), "\n")
## Multikolinearitas (VIF max): 1.0499
cat("Linearitas (Ramsey RESET): p-value =",
round(reset_test$p.value, 4), "\n")
## Linearitas (Ramsey RESET): p-value = 0.7763
cat("\n--- KESIMPULAN ---\n")
##
## --- KESIMPULAN ---
cat("R-squared:", round(r_squared, 4), "\n")
## R-squared: 0.7359
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.7218
cat("\n" %+% rep("=", 60) %+% "\n")
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
cat(rep("=", 60) %+% "\n")