# ============================================================
# ANALISIS MULTIPLE LINEAR REGRESSION (MLR) LENGKAP
# Mulai dari Estimasi hingga Uji Hipotesis & Residual
# ============================================================
# ------------------------------------------------------------
# 0. PERSIAPAN: LOAD LIBRARY DAN DATA
# ------------------------------------------------------------
# Load library
library(car)
## Loading required package: carData
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(nortest)
library(MASS)
library(olsrr)
##
## Attaching package: 'olsrr'
## The following object is masked from 'package:MASS':
##
## cement
## The following object is masked from 'package:datasets':
##
## rivers
library(ggplot2)
library(corrplot)
## corrplot 0.95 loaded
library(moments)
# Set seed untuk reproduktifitas
set.seed(123)
# ------------------------------------------------------------
# 1. MEMBUAT DATA SIMULASI
# ------------------------------------------------------------
# Simulasi data: 100 observasi, 5 variabel independen
n <- 100
# X1 = Normal
X1 <- rnorm(n, mean = 50, sd = 10)
# X2 = Normal
X2 <- rnorm(n, mean = 30, sd = 5)
# X3 = Weibull
X3 <- rweibull(n, shape = 2, scale = 8)
# X4 = Beta
X4 <- rbeta(n, shape1 = 2, shape2 = 5) * 20
# X5 = Chi-Square
# X5 digunakan sebagai variabel dengan distribusi tidak normal
X5 <- rchisq(n, df = 5)
# Model:
# Y = 5 + 0.5*X1 + 0.3*X2 - 0.2*X3 + 0.1*X4
# + 0.2*X5 + error
Y <- 5 +
0.5*X1 +
0.3*X2 -
0.2*X3 +
0.1*X4 +
0.2*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 32.80633 44.39524 26.44797 0.9480543 4.315718 5.141125
## 2 38.62270 47.69823 31.28442 11.2776791 3.707507 1.772879
## 3 56.27093 65.58708 28.76654 2.5232146 4.424703 13.094574
## 4 40.51181 50.70508 28.26229 5.9390167 6.086879 19.656701
## 5 41.94787 51.29288 25.24191 7.7055151 14.002152 6.675153
## 6 45.38196 67.15065 29.77486 7.1507056 5.363228 5.572990
# ------------------------------------------------------------
# 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. :25.36 Min. :26.91 Min. :19.73 Min. : 0.4657
## 1st Qu.:34.92 1st Qu.:45.06 1st Qu.:25.99 1st Qu.: 4.7105
## Median :38.57 Median :50.62 Median :28.87 Median : 7.0309
## Mean :39.01 Mean :50.90 Mean :29.46 Mean : 7.2130
## 3rd Qu.:42.48 3rd Qu.:56.92 3rd Qu.:32.34 3rd Qu.: 9.6466
## Max. :56.27 Max. :71.87 Max. :46.21 Max. :22.1597
## X4 X5
## Min. : 0.3583 Min. : 0.1821
## 1st Qu.: 3.9110 1st Qu.: 2.4663
## Median : 5.8606 Median : 4.1222
## Mean : 6.2836 Mean : 4.6251
## 3rd Qu.: 8.6400 3rd Qu.: 6.3293
## Max. :16.2690 Max. :19.6567
# 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.760 0.263 -0.092 0.009 0.089
## X1 0.760 1.000 -0.050 0.011 -0.059 -0.052
## X2 0.263 -0.050 1.000 -0.101 -0.089 0.063
## X3 -0.092 0.011 -0.101 1.000 0.094 0.107
## X4 0.009 -0.059 -0.089 0.094 1.000 -0.085
## X5 0.089 -0.052 0.063 0.107 -0.085 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.3186 -2.4312 -0.0366 2.1377 8.0578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.41563 2.95303 1.157 0.2503
## X1 0.48151 0.03472 13.867 < 2e-16 ***
## X2 0.33859 0.06602 5.129 1.56e-06 ***
## X3 -0.14242 0.08645 -1.647 0.1028
## X4 0.16578 0.09343 1.774 0.0792 .
## X5 0.23691 0.10406 2.277 0.0251 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.138 on 94 degrees of freedom
## Multiple R-squared: 0.6988, Adjusted R-squared: 0.6827
## F-statistic: 43.61 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) 3.4156 2.9530 1.1567 0.2503
## X1 X1 0.4815 0.0347 13.8674 0.0000
## X2 X2 0.3386 0.0660 5.1289 0.0000
## X3 X3 -0.1424 0.0864 -1.6475 0.1028
## X4 X4 0.1658 0.0934 1.7745 0.0792
## X5 X5 0.2369 0.1041 2.2767 0.0251
# 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) -2.4477 9.2789
## X1 0.4126 0.5505
## X2 0.2075 0.4697
## X3 -0.3141 0.0292
## X4 -0.0197 0.3513
## X5 0.0303 0.4435
# ------------------------------------------------------------
# 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: 43.6093
cat(
"df1:",
df1,
"\n"
)
## df1: 5
cat(
"df2:",
df2,
"\n"
)
## df2: 94
cat(
"p-value:",
format(f_pvalue, scientific = TRUE),
"\n"
)
## p-value: 4.815451e-23
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 = 13.8674, p = 0.0000 ***
## X2 : t = 5.1289, p = 0.0000 ***
## X3 : t = -1.6475, p = 0.1028 ns
## X4 : t = 1.7745, p = 0.0792 .
## X5 : t = 2.2767, p = 0.0251 *
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.6988
cat(
"Adjusted R-squared:",
round(adj_r_squared, 4),
"\n"
)
## Adjusted R-squared: 0.6827
cat(
"Residual Standard Error:",
round(resid_se, 4),
"\n"
)
## Residual Standard Error: 3.1382
cat(
"Interpretasi: Model mampu menjelaskan",
round(r_squared * 100, 2),
"% variasi pada Y\n"
)
## Interpretasi: Model mampu menjelaskan 69.88 % 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.9936
cat(
"p-value =",
round(shapiro_test$p.value, 4),
"\n"
)
## p-value = 0.9237
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.0458
cat(
"p-value =",
round(ks_test$p.value, 4),
"\n"
)
## p-value = 0.9848
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.2166
cat(
"p-value =",
round(ad_test$p.value, 4),
"\n"
)
## p-value = 0.8401
cat(
"Kesimpulan:",
ifelse(
ad_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"
),
"\n"
)
## Kesimpulan: Residual berdistribusi normal
# Jarque-Bera Test
jb_test <- jarque.test(residuals)
cat("\nJarque-Bera Test:\n")
##
## Jarque-Bera Test:
cat(
"JB =",
round(jb_test$statistic, 4),
"\n"
)
## JB = 0.5471
cat(
"p-value =",
round(jb_test$p.value, 4),
"\n"
)
## p-value = 0.7607
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.0029
cat(
"Glejser p-value:",
round(glejser_test$p.value, 4),
"\n"
)
## Glejser p-value: 0.0237
# 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.2335
cat(
"p-value =",
round(dw_test$p.value, 4),
"\n"
)
## p-value = 0.8816
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 = 1.7107
cat(
"p-value =",
round(bg_test$p.value, 4),
"\n"
)
## p-value = 0.1909
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
vif_values <- vif(model)
cat("\nVIF Values:\n")
##
## VIF Values:
print(round(vif_values, 4))
## X1 X2 X3 X4 X5
## 1.0099 1.0241 1.0335 1.0288 1.0288
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.0099 → Tidak ada multikolinearitas
## X2 : VIF = 1.0241 → Tidak ada multikolinearitas
## X3 : VIF = 1.0335 → Tidak ada multikolinearitas
## X4 : VIF = 1.0288 → Tidak ada multikolinearitas
## X5 : VIF = 1.0288 → Tidak ada multikolinearitas
# Tolerance
tolerance <- 1 / vif_values
cat("\nTolerance Values:\n")
##
## Tolerance Values:
print(round(tolerance, 4))
## X1 X2 X3 X4 X5
## 0.9902 0.9764 0.9676 0.9720 0.9720
# Condition Number
cat(
"\nCondition Number:",
round(kappa(model), 4),
"\n"
)
##
## Condition Number: 279.6448
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.1979
cat(
"p-value =",
round(reset_test$p.value, 4),
"\n"
)
## p-value = 0.8208
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.2185
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.2657
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.8833
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: 0.101
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=234.54
## Y ~ X1 + X2 + X3 + X4 + X5
##
## Df Sum of Sq RSS AIC
## <none> 925.74 234.54
## - X3 1 26.73 952.47 235.39
## - X4 1 31.01 956.75 235.84
## - X5 1 51.05 976.79 237.91
## - X2 1 259.06 1184.80 257.22
## - X1 1 1893.87 2819.61 343.92
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.3186 -2.4312 -0.0366 2.1377 8.0578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.41563 2.95303 1.157 0.2503
## X1 0.48151 0.03472 13.867 < 2e-16 ***
## X2 0.33859 0.06602 5.129 1.56e-06 ***
## X3 -0.14242 0.08645 -1.647 0.1028
## X4 0.16578 0.09343 1.774 0.0792 .
## X5 0.23691 0.10406 2.277 0.0251 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.138 on 94 degrees of freedom
## Multiple R-squared: 0.6988, Adjusted R-squared: 0.6827
## F-statistic: 43.61 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 = 3.4156 + 0.4815*X1 + 0.3386*X2 + -0.1424*X3 + 0.1658*X4 + 0.2369*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 = 4.815451e-23
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.1028
## X4: p-value = 0.0792
## X5: p-value = 0.0251
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.9237
cat(
"Heteroskedastisitas (Breusch-Pagan): p-value =",
round(bp_test$p.value, 4),
"\n"
)
## Heteroskedastisitas (Breusch-Pagan): p-value = 0.0029
cat(
"Autokorelasi (Durbin-Watson): DW =",
round(dw_test$statistic, 4),
"\n"
)
## Autokorelasi (Durbin-Watson): DW = 2.2335
cat(
"Multikolinearitas (VIF max):",
round(max(vif_values), 4),
"\n"
)
## Multikolinearitas (VIF max): 1.0335
cat(
"Linearitas (Ramsey RESET): p-value =",
round(reset_test$p.value, 4),
"\n"
)
## Linearitas (Ramsey RESET): p-value = 0.8208
cat("\n--- KESIMPULAN ---\n")
##
## --- KESIMPULAN ---
cat(
"R-squared:",
round(r_squared, 4),
"\n"
)
## R-squared: 0.6988
cat(
"Adjusted R-squared:",
round(adj_r_squared, 4),
"\n"
)
## Adjusted R-squared: 0.6827
cat(
"Distribusi X1: Normal\n"
)
## Distribusi X1: Normal
cat(
"Distribusi X2: Normal\n"
)
## Distribusi X2: Normal
cat(
"Distribusi X3: Weibull\n"
)
## Distribusi X3: Weibull
cat(
"Distribusi X4: Beta\n"
)
## Distribusi X4: Beta
cat(
"Distribusi X5: Chi-Square (non-normal)\n"
)
## Distribusi X5: Chi-Square (non-normal)
cat("\n", rep("=", 60), "\n")
##
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =