# Load library
library(car) # Untuk VIF, Durbin-Watson, dll.
library(lmtest) # Untuk uji Breusch-Pagan, Durbin-Watson
library(nortest) # Untuk uji normalitas (Anderson-Darling)
library(MASS) # Untuk Box-Cox transformation
library(olsrr) # Untuk uji asumsi OLS lengkap
library(ggplot2) # Untuk visualisasi
library(corrplot) # Untuk matriks korelasi
# Set seed untuk reproduktifitas
set.seed(0803)
# ------------------------------------------------------------
# 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 = 10, sd = 2)
X3 <- rnorm(n, mean = 20, sd = 3)
X4 <- rnorm(n, mean = 10, sd = 2)
X5<- rexp(n, rate =1/10)
# Model: Y = 5 + 0.5*X1 + 0.3*X2 - 0.2*X3 + 0.1*X4 + error
Y <- 5 + 0.5*X1 + 0.1*X2 - 0.2*X3 + 0.1*X4 + 0.1*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")
cat("DATA YANG DIGUNAKAN\n")
## DATA YANG DIGUNAKAN
cat("=" %+% rep("=", 60) %+% "\n")
print(head(data))
## Y X1 X2 X3 X4 X5
## 1 29.94086 58.45793 9.570550 18.04178 10.206889 8.709971
## 2 45.41296 71.51275 9.576479 21.20520 8.259070 9.165373
## 3 26.31530 46.19398 12.512782 18.73754 10.398182 3.286498
## 4 30.96492 56.51039 11.249559 20.44780 8.673733 8.683955
## 5 14.15649 32.28167 12.197739 20.85138 10.317956 5.634298
## 6 28.75348 45.89046 11.573951 16.31628 11.426252 8.301876
# ------------------------------------------------------------
# 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. :11.08 Min. :22.93 Min. : 3.963 Min. :13.02
## 1st Qu.:25.69 1st Qu.:43.66 1st Qu.: 8.881 1st Qu.:17.80
## Median :29.77 Median :49.59 Median :10.088 Median :20.35
## Mean :29.59 Mean :49.80 Mean :10.032 Mean :20.03
## 3rd Qu.:34.03 3rd Qu.:57.70 3rd Qu.:11.410 3rd Qu.:22.21
## Max. :46.50 Max. :84.15 Max. :15.900 Max. :27.84
## X4 X5
## Min. : 3.422 Min. : 0.1068
## 1st Qu.: 8.354 1st Qu.: 2.9182
## Median : 9.656 Median : 6.9527
## Mean : 9.595 Mean : 9.4937
## 3rd Qu.:10.825 3rd Qu.:12.7905
## Max. :15.090 Max. :48.1652
# 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.896 -0.112 -0.085 -0.029 0.134
## X1 0.896 1.000 -0.123 -0.125 -0.067 0.047
## X2 -0.112 -0.123 1.000 0.172 0.183 -0.121
## X3 -0.085 -0.125 0.172 1.000 -0.072 0.021
## X4 -0.029 -0.067 0.183 -0.072 1.000 0.053
## X5 0.134 0.047 -0.121 0.021 0.053 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
## -8.3135 -2.0982 0.1183 2.2203 5.4204
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.220288 3.166545 0.070 0.9447
## X1 0.536426 0.027161 19.750 <2e-16 ***
## X2 -0.004402 0.142105 -0.031 0.9754
## X3 0.057208 0.097560 0.586 0.5590
## X4 0.096419 0.156060 0.618 0.5382
## X5 0.065836 0.033441 1.969 0.0519 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.963 on 94 degrees of freedom
## Multiple R-squared: 0.8122, Adjusted R-squared: 0.8023
## F-statistic: 81.33 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.2203 3.1665 0.0696 0.9447
## X1 X1 0.5364 0.0272 19.7497 0.0000
## X2 X2 -0.0044 0.1421 -0.0310 0.9754
## X3 X3 0.0572 0.0976 0.5864 0.5590
## X4 X4 0.0964 0.1561 0.6178 0.5382
## X5 X5 0.0658 0.0334 1.9687 0.0519
# 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) -6.0670 6.5075
## X1 0.4825 0.5904
## X2 -0.2866 0.2778
## X3 -0.1365 0.2509
## X4 -0.2134 0.4063
## X5 -0.0006 0.1322
# ------------------------------------------------------------
# 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: βj = 0 (Model tidak signifikan)\n")
## H0: βj = 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: 81.3323
cat("df1:", df1, "\n")
## df1: 5
cat("df2:", df2, "\n")
## df2: 94
cat("p-value:", format(f_pvalue, scientific = TRUE), "\n")
## p-value: 1.341839e-32
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 = 19.7497, p = 0.0000 ***
## X2 : t = -0.0310, p = 0.9754 ns
## X3 : t = 0.5864, p = 0.5590 ns
## X4 : t = 0.6178, p = 0.5382 ns
## X5 : t = 1.9687, p = 0.0519 .
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.8122
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.8023
cat("Residual Standard Error:", round(resid_se, 4), "\n")
## Residual Standard Error: 2.9628
cat("Interpretasi: Model mampu menjelaskan",
round(r_squared * 100, 2), "% variasi pada Y\n")
## Interpretasi: Model mampu menjelaskan 81.22 % variasi pada Y
# ------------------------------------------------------------
# 6. KRITERIA PEMILIHAN MODEL
# AIC, BIC, dan AICc
# ------------------------------------------------------------
cat("\n============================================================\n")
##
## ============================================================
cat("KRITERIA PEMILIHAN MODEL\n")
## KRITERIA PEMILIHAN MODEL
cat("============================================================\n")
## ============================================================
# AIC
aic_value <- AIC(model)
# BIC
bic_value <- BIC(model)
# AICc
n_obs <- nobs(model)
k <- length(coef(model))
aicc_value <- AIC(model) +
(2 * k * (k + 1)) / (n_obs - k - 1)
cat("\nAIC :", round(aic_value, 4), "\n")
##
## AIC : 508.8247
cat("BIC :", round(bic_value, 4), "\n")
## BIC : 527.0609
cat("AICc :", round(aicc_value, 4), "\n")
## AICc : 509.7279
# ------------------------------------------------------------
# 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.9796
cat("p-value =", round(shapiro_test$p.value, 4), "\n")
## p-value = 0.1245
cat("Kesimpulan:", ifelse(shapiro_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.8667
cat("Glejser p-value:", round(glejser_test$p.value, 4), "\n")
## Glejser p-value: 0.8037
# 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.8092
cat("p-value =", round(dw_test$p.value, 4), "\n")
## p-value = 0.1629
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
# --- 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.0321 1.1019 1.0588 1.0574 1.0252
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.0321 → Tidak ada multikolinearitas
## X2 : VIF = 1.1019 → Tidak ada multikolinearitas
## X3 : VIF = 1.0588 → Tidak ada multikolinearitas
## X4 : VIF = 1.0574 → Tidak ada multikolinearitas
## X5 : VIF = 1.0252 → Tidak ada multikolinearitas
# Tolerance
tolerance <- 1 / vif_values
cat("\nTolerance Values:\n")
##
## Tolerance Values:
print(round(tolerance, 4))
## X1 X2 X3 X4 X5
## 0.9689 0.9076 0.9445 0.9457 0.9754
# Condition Number
cat("\nCondition Number:", round(kappa(model), 4), "\n")
##
## Condition Number: 323.6919
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 = 7e-04
cat("p-value =", round(reset_test$p.value, 4), "\n")
## p-value = 0.9993
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.0602
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.2218
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: 6
# 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.965
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.1111
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=223.04
## Y ~ X1 + X2 + X3 + X4 + X5
##
## Df Sum of Sq RSS AIC
## - X2 1 0.0 825.1 221.04
## - X3 1 3.0 828.1 221.40
## - X4 1 3.4 828.5 221.44
## <none> 825.1 223.04
## - X5 1 34.0 859.2 225.08
## - X1 1 3423.8 4249.0 384.93
##
## Step: AIC=221.04
## Y ~ X1 + X3 + X4 + X5
##
## Df Sum of Sq RSS AIC
## - X3 1 3.1 828.2 219.41
## - X4 1 3.4 828.6 219.45
## <none> 825.1 221.04
## + X2 1 0.0 825.1 223.04
## - X5 1 34.8 859.9 223.17
## - X1 1 3449.0 4274.1 383.52
##
## Step: AIC=219.41
## Y ~ X1 + X4 + X5
##
## Df Sum of Sq RSS AIC
## - X4 1 2.9 831.1 217.76
## <none> 828.2 219.41
## + X3 1 3.1 825.1 221.04
## + X2 1 0.1 828.1 221.40
## - X5 1 35.5 863.7 221.61
## - X1 1 3482.3 4310.5 382.36
##
## Step: AIC=217.76
## Y ~ X1 + X5
##
## Df Sum of Sq RSS AIC
## <none> 831.1 217.76
## + X4 1 2.9 828.2 219.41
## + X3 1 2.6 828.6 219.45
## + X2 1 0.3 830.8 219.73
## - X5 1 36.8 867.9 220.09
## - X1 1 3485.1 4316.2 380.50
cat("\nModel Terbaik (Stepwise):\n")
##
## Model Terbaik (Stepwise):
print(summary(step_model))
##
## Call:
## lm(formula = Y ~ X1 + X5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.1540 -2.1037 0.1229 2.3296 5.8352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.38445 1.37030 1.740 0.085 .
## X1 0.53332 0.02644 20.168 <2e-16 ***
## X5 0.06766 0.03267 2.071 0.041 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.927 on 97 degrees of freedom
## Multiple R-squared: 0.8109, Adjusted R-squared: 0.807
## F-statistic: 208 on 2 and 97 DF, p-value: < 2.2e-16
# ------------------------------------------------------------
# STEPWISE BACKWARD
# ------------------------------------------------------------
cat("\n============================================================\n")
##
## ============================================================
cat("STEPWISE BACKWARD\n")
## STEPWISE BACKWARD
cat("============================================================\n")
## ============================================================
model_backward <- step(
model,
direction = "backward",
trace = 1
)
## Start: AIC=223.04
## Y ~ X1 + X2 + X3 + X4 + X5
##
## Df Sum of Sq RSS AIC
## - X2 1 0.0 825.1 221.04
## - X3 1 3.0 828.1 221.40
## - X4 1 3.4 828.5 221.44
## <none> 825.1 223.04
## - X5 1 34.0 859.2 225.08
## - X1 1 3423.8 4249.0 384.93
##
## Step: AIC=221.04
## Y ~ X1 + X3 + X4 + X5
##
## Df Sum of Sq RSS AIC
## - X3 1 3.1 828.2 219.41
## - X4 1 3.4 828.6 219.45
## <none> 825.1 221.04
## - X5 1 34.8 859.9 223.17
## - X1 1 3449.0 4274.1 383.52
##
## Step: AIC=219.41
## Y ~ X1 + X4 + X5
##
## Df Sum of Sq RSS AIC
## - X4 1 2.9 831.1 217.76
## <none> 828.2 219.41
## - X5 1 35.5 863.7 221.61
## - X1 1 3482.3 4310.5 382.36
##
## Step: AIC=217.76
## Y ~ X1 + X5
##
## Df Sum of Sq RSS AIC
## <none> 831.1 217.76
## - X5 1 36.8 867.9 220.09
## - X1 1 3485.1 4316.2 380.50
cat("\nModel Akhir Backward:\n")
##
## Model Akhir Backward:
print(formula(model_backward))
## Y ~ X1 + X5
cat("\nRingkasan Model Backward:\n")
##
## Ringkasan Model Backward:
print(summary(model_backward))
##
## Call:
## lm(formula = Y ~ X1 + X5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.1540 -2.1037 0.1229 2.3296 5.8352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.38445 1.37030 1.740 0.085 .
## X1 0.53332 0.02644 20.168 <2e-16 ***
## X5 0.06766 0.03267 2.071 0.041 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.927 on 97 degrees of freedom
## Multiple R-squared: 0.8109, Adjusted R-squared: 0.807
## F-statistic: 208 on 2 and 97 DF, p-value: < 2.2e-16
# ------------------------------------------------------------
# STEPWISE FORWARD
# ------------------------------------------------------------
cat("\n============================================================\n")
##
## ============================================================
cat("STEPWISE FORWARD\n")
## STEPWISE FORWARD
cat("============================================================\n")
## ============================================================
model_null <- lm(Y ~ 1, data = data)
model_forward <- step(
model_null,
scope = formula(model),
direction = "forward",
trace = 1
)
## Start: AIC=380.3
## Y ~ 1
##
## Df Sum of Sq RSS AIC
## + X1 1 3526.9 867.9 220.09
## <none> 4394.8 380.30
## + X5 1 78.6 4316.2 380.50
## + X2 1 55.6 4339.2 381.03
## + X3 1 32.0 4362.8 381.57
## + X4 1 3.7 4391.1 382.22
##
## Step: AIC=220.09
## Y ~ X1
##
## Df Sum of Sq RSS AIC
## + X5 1 36.760 831.13 217.76
## <none> 867.89 220.09
## + X4 1 4.190 863.70 221.60
## + X3 1 3.109 864.78 221.73
## + X2 1 0.027 867.86 222.09
##
## Step: AIC=217.76
## Y ~ X1 + X5
##
## Df Sum of Sq RSS AIC
## <none> 831.13 217.76
## + X4 1 2.92505 828.20 219.41
## + X3 1 2.56550 828.56 219.45
## + X2 1 0.29148 830.83 219.73
cat("\nModel Akhir Forward:\n")
##
## Model Akhir Forward:
print(formula(model_forward))
## Y ~ X1 + X5
cat("\nRingkasan Model Forward:\n")
##
## Ringkasan Model Forward:
print(summary(model_forward))
##
## Call:
## lm(formula = Y ~ X1 + X5, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.1540 -2.1037 0.1229 2.3296 5.8352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.38445 1.37030 1.740 0.085 .
## X1 0.53332 0.02644 20.168 <2e-16 ***
## X5 0.06766 0.03267 2.071 0.041 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.927 on 97 degrees of freedom
## Multiple R-squared: 0.8109, Adjusted R-squared: 0.807
## F-statistic: 208 on 2 and 97 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.2203 + 0.5364*X1 + -0.0044*X2 + 0.0572*X3 + 0.0964*X4 + 0.0658*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.341839e-32
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.9754
## X3: p-value = 0.5590
## X4: p-value = 0.5382
## X5: p-value = 0.0519
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.1245
cat("Heteroskedastisitas (Breusch-Pagan): p-value =",
round(bp_test$p.value, 4), "\n")
## Heteroskedastisitas (Breusch-Pagan): p-value = 0.8667
cat("Autokorelasi (Durbin-Watson): DW =",
round(dw_test$statistic, 4), "\n")
## Autokorelasi (Durbin-Watson): DW = 1.8092
cat("Multikolinearitas (VIF max):",
round(max(vif_values), 4), "\n")
## Multikolinearitas (VIF max): 1.1019
cat("Linearitas (Ramsey RESET): p-value =",
round(reset_test$p.value, 4), "\n")
## Linearitas (Ramsey RESET): p-value = 0.9993
cat("\n--- KESIMPULAN ---\n")
##
## --- KESIMPULAN ---
cat("R-squared:", round(r_squared, 4), "\n")
## R-squared: 0.8122
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.8023
cat("\n" %+% rep("=", 60) %+% "\n")
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
cat(rep("=", 60) %+% "\n")