library(olsrr)
## Warning: package 'olsrr' was built under R version 4.5.3
##
## Attaching package: 'olsrr'
## The following object is masked from 'package:datasets':
##
## rivers
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
set.seed(123)
#MEMBACA DATA PENELITIAN
file_data <- file.choose()
data <- read.csv2(
file_data,
header = TRUE,
stringsAsFactors = FALSE,
check.names = FALSE,
na.strings = c("", "NA")
)
names(data) <- trimws(names(data))
cat("\nNama kolom dari file:\n")
##
## Nama kolom dari file:
print(names(data))
## [1] "Provinsi" "Total fasilitas kesehatan (Y)"
## [3] "Penduduk (X1)" "Desa Kelurahan (X2)"
## [5] "Tenaga Kesehatan (X3)" "Luas wilayah (X4)"
names(data) <- c(
"Provinsi",
"Y",
"X1",
"X2",
"X3",
"X4"
)
data$Y <- as.numeric(
gsub(",", ".", data$Y, fixed = TRUE)
)
data$X1 <- as.numeric(
gsub(",", ".", data$X1, fixed = TRUE)
)
data$X2 <- as.numeric(
gsub(",", ".", data$X2, fixed = TRUE)
)
data$X3 <- as.numeric(
gsub(",", ".", data$X3, fixed = TRUE)
)
data$X4 <- as.numeric(
gsub(",", ".", data$X4, fixed = TRUE)
)
data <- data[
!is.na(data$Provinsi) &
trimws(data$Provinsi) != "",
]
data <- data[
complete.cases(
data[, c("Y", "X1", "X2", "X3", "X4")]
),
]
rownames(data) <- NULL
cat(
paste(rep("=", 60), collapse = ""),
"\n"
)
## ============================================================
cat("DATA YANG DIGUNAKAN\n")
## DATA YANG DIGUNAKAN
cat(
paste(rep("=", 60), collapse = ""),
"\n"
)
## ============================================================
print(data)
## Provinsi Y X1 X2 X3 X4
## 1 Aceh 452 5626.0 6513 623 56835.02
## 2 Sumatera Utara 827 15785.8 6113 1625 72437.76
## 3 Sumatera Barat 363 5914.3 1287 1088 42107.67
## 4 Riau 323 6811.2 1870 1821 89900.78
## 5 Jambi 253 3768.5 1586 2122 49023.04
## 6 Sumatera Selatan 443 8928.5 3285 519 86771.92
## 7 Bengkulu 207 2138.0 1513 2321 20122.21
## 8 Lampung 405 9522.9 2654 1951 33570.76
## 9 Kepulauan Bangka Belitung 93 1550.8 393 2011 16670.23
## 10 Kepulauan Riau 135 2213.5 430 2656 8170.38
## 11 DKI Jakarta 238 10678.0 267 1180 661.53
## 12 Jawa Barat 1542 50759.0 5957 461 37053.33
## 13 Jawa Tengah 1250 38233.9 8563 412 34347.43
## 14 DI Yogyakarta 1418 42089.3 8494 415 48055.88
## 15 Jawa Timur 389 12537.4 1552 2262 9355.76
## 16 Banten 202 4461.3 718 1859 5582.83
## 17 Bali 223 5731.1 1180 968 19631.99
## 18 Nusa Tenggara Barat 508 5742.6 3538 245 46378.11
## 19 Nusa Tenggara Timur 311 5766.0 2157 1765 147018.06
## 20 Kalimantan Barat 239 2845.0 1577 2449 153430.36
## 21 Kalimantan Tengah 294 4323.3 2015 1007 37125.43
## 22 Kalimantan Selatan 254 4267.6 1055 2310 126951.76
## 23 Kalimantan Timur 75 749.4 484 2719 69900.89
## 24 Kalimantan Utara 263 2721.4 1838 3072 14488.43
## 25 Sulawesi Utara 260 3156.1 2022 2163 61496.98
## 26 Sulawesi Tengah 598 9563.1 3060 711 45323.98
## 27 Sulawesi Selatan 351 2836.7 2292 971 36139.30
## 28 Sulawesi Tenggara 117 1242.2 732 3096 12024.98
## 29 Gorontalo 114 1525.3 650 3523 16590.67
## 30 Sulawesi Barat 270 1970.6 1262 3513 46133.83
## 31 Maluku 173 1373.8 1209 3568 31465.98
## 32 Maluku Utara 95 587.6 970 1968 60308.59
## 33 Papua Barat 131 636.4 1056 1980 39103.06
## 34 Papua Barat Daya 142 1073.6 1029 4122 81383.32
## 35 Papua 93 549.7 690 2394 117858.97
## 36 Papua Selatan 141 1492.3 1208 2392 61079.59
## 37 Papua Tengah 191 1484.9 2634 1671 52508.66
cat("\n6 BARIS PERTAMA:\n")
##
## 6 BARIS PERTAMA:
print(head(data))
## Provinsi Y X1 X2 X3 X4
## 1 Aceh 452 5626.0 6513 623 56835.02
## 2 Sumatera Utara 827 15785.8 6113 1625 72437.76
## 3 Sumatera Barat 363 5914.3 1287 1088 42107.67
## 4 Riau 323 6811.2 1870 1821 89900.78
## 5 Jambi 253 3768.5 1586 2122 49023.04
## 6 Sumatera Selatan 443 8928.5 3285 519 86771.92
cat("\nSTRUKTUR DATA:\n")
##
## STRUKTUR DATA:
str(data)
## 'data.frame': 37 obs. of 6 variables:
## $ Provinsi: chr "Aceh" "Sumatera Utara" "Sumatera Barat" "Riau" ...
## $ Y : num 452 827 363 323 253 443 207 405 93 135 ...
## $ X1 : num 5626 15786 5914 6811 3768 ...
## $ X2 : num 6513 6113 1287 1870 1586 ...
## $ X3 : num 623 1625 1088 1821 2122 ...
## $ X4 : num 56835 72438 42108 89901 49023 ...
n <- nrow(data)
cat("\nJumlah observasi =", n, "\n")
##
## Jumlah observasi = 37
cat("Jumlah variabel =", ncol(data), "\n")
## Jumlah variabel = 6
#DEFINISI VARIABEL
cat("\nKeterangan Variabel:\n")
##
## Keterangan Variabel:
cat("Y = Total Fasilitas Kesehatan\n")
## Y = Total Fasilitas Kesehatan
cat("X1 = Jumlah Penduduk (Ribu)\n")
## X1 = Jumlah Penduduk (Ribu)
cat("X2 = Jumlah Desa/Kelurahan\n")
## X2 = Jumlah Desa/Kelurahan
cat("X3 = Jumlah Tenaga Kesehatan\n")
## X3 = Jumlah Tenaga Kesehatan
cat("X4 = Luas Wilayah (km2)\n")
## X4 = Luas Wilayah (km2)
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
#STATISTIK DESKRIPTIF DAN EKSPLORASI DATA
cat("1. STATISTIK DESKRIPTIF\n")
## 1. STATISTIK DESKRIPTIF
# Statistik deskriptif
summary(data[, c("Y", "X1", "X2", "X3", "X4")])
## Y X1 X2 X3
## Min. : 75.0 Min. : 549.7 Min. : 267 Min. : 245
## 1st Qu.: 142.0 1st Qu.: 1525.3 1st Qu.:1029 1st Qu.:1007
## Median : 254.0 Median : 3768.5 Median :1552 Median :1968
## Mean : 361.7 Mean : 7585.3 Mean :2266 Mean :1890
## 3rd Qu.: 389.0 3rd Qu.: 6811.2 3rd Qu.:2634 3rd Qu.:2394
## Max. :1542.0 Max. :50759.0 Max. :8563 Max. :4122
## X4
## Min. : 661.5
## 1st Qu.: 20122.2
## Median : 45324.0
## Mean : 51000.3
## 3rd Qu.: 61497.0
## Max. :153430.4
# Matriks korelasi
cat("\nMatriks Korelasi:\n")
##
## Matriks Korelasi:
cor_matrix <- cor(
data[, c("Y", "X1", "X2", "X3", "X4")],
use = "complete.obs"
)
print(
round(cor_matrix, 3)
)
## Y X1 X2 X3 X4
## Y 1.000 0.966 0.893 -0.613 -0.035
## X1 0.966 1.000 0.799 -0.552 -0.096
## X2 0.893 0.799 1.000 -0.625 0.050
## X3 -0.613 -0.552 -0.625 1.000 0.033
## X4 -0.035 -0.096 0.050 0.033 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[, c("Y", "X1", "X2", "X3", "X4")],
main = "Scatter Plot Matrix",
pch = 19,
col = "steelblue"
)
#ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
cat("2. ESTIMASI MODEL MLR\n")
## 2. ESTIMASI MODEL MLR
model <- lm(
Y ~ X1 + X2 + X3 + X4,
data = data
)
cat("\nRingkasan Model:\n")
##
## Ringkasan Model:
print(
summary(model)
)
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -100.267 -32.771 -9.315 16.780 141.307
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 97.4711551 36.8127795 2.648 0.0125 *
## X1 0.0213641 0.0014718 14.516 1.25e-15 ***
## X2 0.0507886 0.0084492 6.011 1.05e-06 ***
## X3 -0.0113428 0.0127231 -0.892 0.3793
## X4 0.0001670 0.0002685 0.622 0.5385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 59.48 on 32 degrees of freedom
## Multiple R-squared: 0.9746, Adjusted R-squared: 0.9714
## F-statistic: 306.7 on 4 and 32 DF, p-value: < 2.2e-16
#KOEFISIEN REGRESI
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) 97.4712 36.8128 2.6478 0.0125
## X1 X1 0.0214 0.0015 14.5159 0.0000
## X2 X2 0.0508 0.0084 6.0110 0.0000
## X3 X3 -0.0113 0.0127 -0.8915 0.3793
## X4 X4 0.0002 0.0003 0.6217 0.5385
cat("\nInterval Kepercayaan 95%:\n")
##
## Interval Kepercayaan 95%:
ci <- confint(
model,
level = 0.95
)
print(
round(ci, 4)
)
## 2.5 % 97.5 %
## (Intercept) 22.4860 172.4563
## X1 0.0184 0.0244
## X2 0.0336 0.0680
## X3 -0.0373 0.0146
## X4 -0.0004 0.0007
#UJI SIGNIFIKANSI SIMULTAN (UJI F)
cat("3. UJI SIGNIFIKANSI SIMULTAN (UJI F)\n")
## 3. UJI SIGNIFIKANSI SIMULTAN (UJI F)
# 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 = 0 (Model tidak signifikan)\n"
)
## H0: β1 = β2 = β3 = β4 = 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: 306.711
cat(
"df1:",
df1,
"\n"
)
## df1: 4
cat(
"df2:",
df2,
"\n"
)
## df2: 32
cat(
"p-value:",
format(
f_pvalue,
scientific = TRUE
),
"\n"
)
## p-value: 5.044065e-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%
#UJI SIGNIFIKANSI PARSIAL (UJI T)
cat("4. UJI SIGNIFIKANSI PARSIAL (UJI T)\n")
## 4. UJI SIGNIFIKANSI PARSIAL (UJI T)
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 = 14.5159, p = 0.0000 ***
## X2 : t = 6.0110, p = 0.0000 ***
## X3 : t = -0.8915, p = 0.3793 ns
## X4 : t = 0.6217, p = 0.5385 ns
cat(
"\nKeterangan: *** p<0.001, ** p<0.01, * p<0.05, . p<0.1, ns tidak signifikan\n"
)
##
## Keterangan: *** p<0.001, ** p<0.01, * p<0.05, . p<0.1, ns tidak signifikan
#KOEFISIEN DETERMINASI (R² DAN ADJUSTED R²)
cat("5. KOEFISIEN DETERMINASI\n")
## 5. KOEFISIEN DETERMINASI
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.9746
cat(
"Adjusted R-squared:",
round(
adj_r_squared,
4
),
"\n"
)
## Adjusted R-squared: 0.9714
cat(
"Residual Standard Error:",
round(
resid_se,
4
),
"\n"
)
## Residual Standard Error: 59.4831
cat(
"Interpretasi: Model mampu menjelaskan",
round(
r_squared * 100,
2
),
"% variasi pada Y\n"
)
## Interpretasi: Model mampu menjelaskan 97.46 % variasi pada Y
library(nortest)
## Warning: package 'nortest' was built under R version 4.5.2
library(moments)
## Warning: package 'moments' was built under R version 4.5.2
library(lmtest)
## 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.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(car)
## 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.3
#UJI ASUMSI RESIDUAL
cat("6. UJI ASUMSI RESIDUAL\n")
## 6. UJI ASUMSI RESIDUAL
residuals <- residuals(model)
fitted_values <- fitted(model)
standardized_resid <- rstandard(model)
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.9523
cat(
"p-value =",
round(
shapiro_test$p.value,
4
),
"\n"
)
## p-value = 0.1142
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.1565
cat(
"p-value =",
round(
ks_test$p.value,
4
),
"\n"
)
## p-value = 0.2935
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.7614
cat(
"p-value =",
round(
ad_test$p.value,
4
),
"\n"
)
## p-value = 0.0434
cat(
"Kesimpulan:",
ifelse(
ad_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"
),
"\n"
)
## Kesimpulan: Residual TIDAK 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 = 1.8463
cat(
"p-value =",
round(
jb_test$p.value,
4
),
"\n"
)
## p-value = 0.3973
cat(
"Kesimpulan:",
ifelse(
jb_test$p.value > 0.05,
"Residual berdistribusi normal",
"Residual TIDAK berdistribusi normal"
),
"\n"
)
## Kesimpulan: Residual berdistribusi normal
cat(
"\n--- 6.2. UJI HETEROSKEDASTISITAS ---\n"
)
##
## --- 6.2. UJI HETEROSKEDASTISITAS ---
# Breusch-Pagan Test
bp_test <- bptest(model)
# Glejser Test
glejser_test <- bptest(
model,
~ fitted(model)
)
cat(
"Breusch-Pagan p-value:",
round(
bp_test$p.value,
4
),
"\n"
)
## Breusch-Pagan p-value: 0.0303
cat(
"Glejser p-value:",
round(
glejser_test$p.value,
4
),
"\n"
)
## Glejser p-value: 0.1136
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
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.9649
cat(
"p-value =",
round(
dw_test$p.value,
4
),
"\n"
)
## p-value = 0.3324
cat(
"Kesimpulan:",
ifelse(
dw_test$p.value > 0.05,
"Tidak ada autokorelasi",
"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.0788
cat(
"p-value =",
round(
bg_test$p.value,
4
),
"\n"
)
## p-value = 0.7789
cat(
"Kesimpulan:",
ifelse(
bg_test$p.value > 0.05,
"Tidak ada autokorelasi",
"Ada autokorelasi"
),
"\n"
)
## Kesimpulan: Tidak ada autokorelasi
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
## 2.9438 3.3570 1.6675 1.0601
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 = 2.9438 → Tidak ada multikolinearitas
## X2 : VIF = 3.3570 → Tidak ada multikolinearitas
## X3 : VIF = 1.6675 → Tidak ada multikolinearitas
## X4 : VIF = 1.0601 → Tidak ada multikolinearitas
tolerance <- 1 / vif_values
cat("\nTolerance Values:\n")
##
## Tolerance Values:
print(
round(
tolerance,
4
)
)
## X1 X2 X3 X4
## 0.3397 0.2979 0.5997 0.9433
condition_number <- kappa(
model.matrix(model)
)
cat(
"\nCondition Number:",
round(
condition_number,
4
),
"\n"
)
##
## Condition Number: 206951
cat(
"(Condition Number > 30 mengindikasikan multikolinearitas)\n"
)
## (Condition Number > 30 mengindikasikan multikolinearitas)
cat(
"\n--- 6.5. UJI LINEARITAS ---\n"
)
##
## --- 6.5. UJI LINEARITAS ---
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 = 1.5475
cat(
"p-value =",
round(
reset_test$p.value,
4
),
"\n"
)
## p-value = 0.2293
cat(
"Kesimpulan:",
ifelse(
reset_test$p.value > 0.05,
"Model linear (spesifikasi benar)",
"Model TIDAK linear (spesifikasi salah)"
),
"\n"
)
## Kesimpulan: Model linear (spesifikasi benar)
#DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
cat(
"7. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS\n"
)
## 7. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
# 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.5589
cat(
"Jumlah observasi dengan Cook's D > 1:",
sum(
cooks_d > 1
),
"\n"
)
## Jumlah observasi dengan Cook's D > 1: 0
# Leverage
leverage <- hatvalues(model)
cat("\nLeverage (Hat Values):\n")
##
## Leverage (Hat Values):
cat(
"Nilai maksimum:",
round(
max(leverage),
4
),
"\n"
)
## Nilai maksimum: 0.544
# Threshold 2(k+1)/n
leverage_threshold <-
2 * length(coef(model)) / n
cat(
"Threshold 2(k+1)/n:",
round(
leverage_threshold,
4
),
"\n"
)
## Threshold 2(k+1)/n: 0.2703
cat(
"Jumlah observasi dengan leverage > threshold:",
sum(
leverage > leverage_threshold
),
"\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.7024
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: 3
#VISUALISASI DIAGNOSTIK
cat("8. VISUALISASI DIAGNOSTIK\n")
## 8. VISUALISASI DIAGNOSTIK
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 = leverage_threshold,
col = "red",
lty = 2
)
# Reset par
par(
mfrow = c(1, 1)
)
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:olsrr':
##
## cement
#UJI BOX-COX UNTUK TRANSFORMASI
cat("9. UJI BOX-COX\n")
## 9. UJI BOX-COX
# Box-Cox transformation
bc <- boxcox(
model,
lambda = seq(
-2,
2,
0.1
)
)
# Lambda optimal
lambda_opt <-
bc$x[
which.max(bc$y)
]
cat(
"\nOptimal Lambda:",
round(
lambda_opt,
4
),
"\n"
)
##
## Optimal Lambda: 0.9495
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
#PEMILIHAN MODEL (STEPWISE)
cat("10. PEMILIHAN MODEL (STEPWISE)\n")
## 10. PEMILIHAN MODEL (STEPWISE)
step_model <- step(
model,
direction = "both",
trace = 1
)
## Start: AIC=306.97
## Y ~ X1 + X2 + X3 + X4
##
## Df Sum of Sq RSS AIC
## - X4 1 1368 114591 305.41
## - X3 1 2812 116036 305.88
## <none> 113224 306.97
## - X2 1 127846 241069 332.93
## - X1 1 745545 858769 379.94
##
## Step: AIC=305.41
## Y ~ X1 + X2 + X3
##
## Df Sum of Sq RSS AIC
## - X3 1 2597 117188 304.24
## <none> 114591 305.41
## + X4 1 1368 113224 306.97
## - X2 1 140315 254906 333.00
## - X1 1 768465 883056 378.97
##
## Step: AIC=304.24
## Y ~ X1 + X2
##
## Df Sum of Sq RSS AIC
## <none> 117188 304.24
## + X3 1 2597 114591 305.41
## + X4 1 1152 116036 305.88
## - X2 1 178719 295907 336.51
## - X1 1 788364 905552 377.90
cat("\nModel Terbaik (Stepwise):\n")
##
## Model Terbaik (Stepwise):
print(
summary(step_model)
)
##
## Call:
## lm(formula = Y ~ X1 + X2, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -107.879 -33.211 -8.738 15.856 150.892
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.557904 14.447656 5.299 7.03e-06 ***
## X1 0.021301 0.001408 15.124 < 2e-16 ***
## X2 0.054524 0.007572 7.201 2.49e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 58.71 on 34 degrees of freedom
## Multiple R-squared: 0.9737, Adjusted R-squared: 0.9721
## F-statistic: 629.1 on 2 and 34 DF, p-value: < 2.2e-16
#RINGKASAN HASIL
cat("11. RINGKASAN HASIL ANALISIS\n")
## 11. RINGKASAN HASIL ANALISIS
cat("\n--- MODEL AWAL ---\n")
##
## --- MODEL AWAL ---
cat("Persamaan Regresi:\n")
## Persamaan Regresi:
cat(
sprintf(
"Y = %.4f + %.4f*X1 + %.4f*X2 + %.4f*X3 + %.4f*X4\n",
coef(model)[1],
coef(model)[2],
coef(model)[3],
coef(model)[4],
coef(model)[5]
)
)
## Y = 97.4712 + 0.0214*X1 + 0.0508*X2 + -0.0113*X3 + 0.0002*X4
cat("\n--- UJI SIGNIFIKANSI ---\n")
##
## --- UJI SIGNIFIKANSI ---
cat(
"Uji F (Simultan): p-value =",
format(
f_pvalue,
scientific = TRUE
),
"\n"
)
## Uji F (Simultan): p-value = 5.044065e-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.3793
## X4: p-value = 0.5385
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.1142
cat(
"Heteroskedastisitas (Breusch-Pagan): p-value =",
round(
bp_test$p.value,
4
),
"\n"
)
## Heteroskedastisitas (Breusch-Pagan): p-value = 0.0303
cat(
"Autokorelasi (Durbin-Watson): DW =",
round(
dw_test$statistic,
4
),
"\n"
)
## Autokorelasi (Durbin-Watson): DW = 1.9649
cat(
"Multikolinearitas (VIF max):",
round(
max(vif_values),
4
),
"\n"
)
## Multikolinearitas (VIF max): 3.357
cat(
"Linearitas (Ramsey RESET): p-value =",
round(
reset_test$p.value,
4
),
"\n"
)
## Linearitas (Ramsey RESET): p-value = 0.2293
cat("\n--- MODEL HASIL STEPWISE ---\n")
##
## --- MODEL HASIL STEPWISE ---
cat(
"Formula model stepwise:\n"
)
## Formula model stepwise:
print(
formula(step_model)
)
## Y ~ X1 + X2
cat(
"\nR-squared model awal:",
round(
r_squared,
4
),
"\n"
)
##
## R-squared model awal: 0.9746
cat(
"Adjusted R-squared model awal:",
round(
adj_r_squared,
4
),
"\n"
)
## Adjusted R-squared model awal: 0.9714
cat(
"\nR-squared model stepwise:",
round(
summary(step_model)$r.squared,
4
),
"\n"
)
##
## R-squared model stepwise: 0.9737
cat(
"Adjusted R-squared model stepwise:",
round(
summary(step_model)$adj.r.squared,
4
),
"\n"
)
## Adjusted R-squared model stepwise: 0.9721
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
#PEMILIHAN MODEL (STEPWISE AIC DAN BIC)
# Stepwise berdasarkan AIC
model_aic <- step(
model,
direction = "both",
trace = 1
)
## Start: AIC=306.97
## Y ~ X1 + X2 + X3 + X4
##
## Df Sum of Sq RSS AIC
## - X4 1 1368 114591 305.41
## - X3 1 2812 116036 305.88
## <none> 113224 306.97
## - X2 1 127846 241069 332.93
## - X1 1 745545 858769 379.94
##
## Step: AIC=305.41
## Y ~ X1 + X2 + X3
##
## Df Sum of Sq RSS AIC
## - X3 1 2597 117188 304.24
## <none> 114591 305.41
## + X4 1 1368 113224 306.97
## - X2 1 140315 254906 333.00
## - X1 1 768465 883056 378.97
##
## Step: AIC=304.24
## Y ~ X1 + X2
##
## Df Sum of Sq RSS AIC
## <none> 117188 304.24
## + X3 1 2597 114591 305.41
## + X4 1 1152 116036 305.88
## - X2 1 178719 295907 336.51
## - X1 1 788364 905552 377.90
# Ringkasan model AIC
summary(model_aic)
##
## Call:
## lm(formula = Y ~ X1 + X2, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -107.879 -33.211 -8.738 15.856 150.892
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.557904 14.447656 5.299 7.03e-06 ***
## X1 0.021301 0.001408 15.124 < 2e-16 ***
## X2 0.054524 0.007572 7.201 2.49e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 58.71 on 34 degrees of freedom
## Multiple R-squared: 0.9737, Adjusted R-squared: 0.9721
## F-statistic: 629.1 on 2 and 34 DF, p-value: < 2.2e-16
# Nilai AIC model AIC
AIC(model_aic)
## [1] 411.2443
# Stepwise berdasarkan BIC
model_bic <- step(
model,
direction = "both",
k = log(nrow(data)),
trace = 1
)
## Start: AIC=315.02
## Y ~ X1 + X2 + X3 + X4
##
## Df Sum of Sq RSS AIC
## - X4 1 1368 114591 311.86
## - X3 1 2812 116036 312.32
## <none> 113224 315.02
## - X2 1 127846 241069 339.37
## - X1 1 745545 858769 386.38
##
## Step: AIC=311.86
## Y ~ X1 + X2 + X3
##
## Df Sum of Sq RSS AIC
## - X3 1 2597 117188 309.08
## <none> 114591 311.86
## + X4 1 1368 113224 315.02
## - X2 1 140315 254906 337.83
## - X1 1 768465 883056 383.80
##
## Step: AIC=309.08
## Y ~ X1 + X2
##
## Df Sum of Sq RSS AIC
## <none> 117188 309.08
## + X3 1 2597 114591 311.86
## + X4 1 1152 116036 312.32
## - X2 1 178719 295907 339.74
## - X1 1 788364 905552 381.12
# Ringkasan model BIC
summary(model_bic)
##
## Call:
## lm(formula = Y ~ X1 + X2, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -107.879 -33.211 -8.738 15.856 150.892
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.557904 14.447656 5.299 7.03e-06 ***
## X1 0.021301 0.001408 15.124 < 2e-16 ***
## X2 0.054524 0.007572 7.201 2.49e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 58.71 on 34 degrees of freedom
## Multiple R-squared: 0.9737, Adjusted R-squared: 0.9721
## F-statistic: 629.1 on 2 and 34 DF, p-value: < 2.2e-16
# Nilai BIC model BIC
BIC(model_bic)
## [1] 417.688
# PERBANDINGAN MODEL
AIC(model)
## [1] 413.9709
AIC(model_aic)
## [1] 411.2443
BIC(model)
## [1] 423.6364
BIC(model_bic)
## [1] 417.688
# Formula model
formula(model_aic)
## Y ~ X1 + X2
formula(model_bic)
## Y ~ X1 + X2