# LOAD LIBRARY
library(readxl) # Membaca file Excel
## Warning: package 'readxl' was built under R version 4.5.2
library(car) # VIF
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.5.2
library(lmtest) # Breusch-Pagan, Durbin-Watson, RESET
## Warning: package 'lmtest' was built under R version 4.5.2
## 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(nortest) # Uji normalitas tambahan
## Warning: package 'nortest' was built under R version 4.5.2
library(MASS) # Box-Cox
library(olsrr) # Diagnostik OLS
## Warning: package 'olsrr' was built under R version 4.5.3
##
## Attaching package: 'olsrr'
## The following object is masked from 'package:MASS':
##
## cement
## The following object is masked from 'package:datasets':
##
## rivers
library(ggplot2) # Visualisasi
## Warning: package 'ggplot2' was built under R version 4.5.3
library(corrplot) # Matriks korelasi
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
# Set seed untuk reproduktifitas
set.seed(0803)
# 1. MEMBACA DATA
data_raw <- read_excel("C:/Users/ASUS/Downloads/Analisis Multivariat/Dataset_MLR_Kalimantan_Tengah_2024.xlsx")
# Menyamakan nama kolom agar mudah digunakan dalam R
data <- data.frame(
Wilayah = data_raw[["Kabupaten_Kota"]],
Y = data_raw[["LPE"]],
X1 = data_raw[["IPM"]],
X2 = data_raw[["Kemiskinan"]],
X3 = data_raw[["TPT"]]
)
# Jumlah observasi
n <- nrow(data)
print(data)
## Wilayah Y X1 X2 X3
## 1 Kotawaringin Barat 0.98 73.95 3.59 4.76
## 2 Kotawaringin Timur -3.06 73.18 5.62 5.25
## 3 Kapuas -1.04 71.18 5.04 4.98
## 4 Barito Selatan -2.90 72.73 4.45 4.21
## 5 Barito Utara -2.24 71.44 5.17 5.29
## 6 Sukamara 1.98 69.04 3.23 4.70
## 7 Lamandau 1.85 72.21 3.09 2.83
## 8 Seruyan -2.23 69.22 6.85 4.30
## 9 Katingan -3.18 72.45 4.79 5.69
## 10 Pulang Pisau 2.68 70.57 4.09 2.63
## 11 Gunung Mas 3.39 72.00 4.75 2.49
## 12 Barito Timur -2.73 73.09 6.09 2.91
## 13 Murung Raya -2.45 69.54 5.85 3.10
## 14 Palangka Raya -2.85 81.17 3.44 5.95
## 15 Kotawaringin Barat 5.61 74.15 3.95 4.70
## 16 Kotawaringin Timur 2.10 73.25 5.91 5.15
## 17 Kapuas 4.71 71.34 5.35 4.91
## 18 Barito Selatan 2.13 73.05 4.62 4.16
## 19 Barito Utara 2.82 71.64 5.61 5.14
## 20 Sukamara 4.74 69.28 3.66 4.65
## 21 Lamandau 4.01 72.28 3.56 2.30
## 22 Seruyan 2.12 69.31 7.22 4.25
## 23 Katingan 2.90 72.66 5.25 5.50
## 24 Pulang Pisau 3.24 70.65 4.24 2.60
## 25 Gunung Mas 5.09 72.22 5.35 3.11
## 26 Barito Timur 2.97 73.17 6.38 3.22
## 27 Murung Raya 4.38 69.67 6.15 3.03
## 28 Palangka Raya 4.32 81.22 3.75 5.86
## 29 Kotawaringin Barat 6.01 74.39 3.93 4.51
## 30 Kotawaringin Timur 7.41 73.45 5.95 5.00
## 31 Kapuas 7.04 71.72 5.52 3.91
## 32 Barito Selatan 6.28 73.45 4.88 3.53
## 33 Barito Utara 6.24 72.16 5.80 4.82
## 34 Sukamara 5.62 69.86 3.72 6.46
## 35 Lamandau 6.05 72.81 3.34 3.41
## 36 Seruyan 4.01 69.81 7.43 3.96
## 37 Katingan 5.58 73.43 5.50 5.33
## 38 Pulang Pisau 4.68 71.05 4.70 1.96
## 39 Gunung Mas 6.47 72.50 5.64 2.96
## 40 Barito Timur 6.06 73.69 6.59 2.95
## 41 Murung Raya 7.03 70.13 6.40 2.77
## 42 Palangka Raya 6.25 81.47 3.61 5.64
## 43 Kotawaringin Barat 6.10 74.92 4.18 4.45
## 44 Kotawaringin Timur 1.81 73.99 5.69 4.77
## 45 Kapuas 5.71 72.40 5.21 3.66
## 46 Barito Selatan 3.27 74.01 4.72 4.33
## 47 Barito Utara 5.49 72.71 5.35 4.85
## 48 Sukamara 5.64 70.35 3.96 5.23
## 49 Lamandau 1.59 73.44 3.12 3.32
## 50 Seruyan 4.55 70.24 7.12 3.61
## 51 Katingan 5.98 73.90 4.99 4.96
## 52 Pulang Pisau 4.84 71.62 4.58 2.07
## 53 Gunung Mas 4.25 73.18 5.47 3.24
## 54 Barito Timur 3.47 74.21 6.63 3.37
## 55 Murung Raya 5.46 70.91 6.44 2.75
## 56 Palangka Raya 6.57 81.95 3.44 5.13
## 57 Kotawaringin Barat 4.10 75.35 4.11 4.42
## 58 Kotawaringin Timur 4.00 74.47 5.66 4.63
## 59 Kapuas 4.95 72.98 5.25 3.61
## 60 Barito Selatan 4.70 74.76 4.83 4.12
## 61 Barito Utara 5.08 73.17 5.67 4.71
## 62 Sukamara 3.89 70.83 4.14 4.95
## 63 Lamandau 3.64 73.95 3.25 3.17
## 64 Seruyan 3.04 70.66 7.08 3.47
## 65 Katingan 4.67 74.37 5.26 4.88
## 66 Pulang Pisau 4.41 72.36 4.56 1.99
## 67 Gunung Mas 4.48 73.88 5.68 3.12
## 68 Barito Timur 4.29 74.81 6.66 3.26
## 69 Murung Raya 5.05 71.58 6.58 2.90
## 70 Palangka Raya 6.62 82.53 3.52 5.02
# 2. STATISTIK DESKRIPTIF
# Statistik deskriptif dasar
summary(data[, c("Y", "X1", "X2", "X3")])
## Y X1 X2 X3
## Min. :-3.180 Min. :69.04 Min. :3.090 Min. :1.960
## 1st Qu.: 2.715 1st Qu.:71.22 1st Qu.:4.095 1st Qu.:3.132
## Median : 4.350 Median :72.72 Median :5.190 Median :4.230
## Mean : 3.596 Mean :72.99 Mean :5.046 Mean :4.070
## 3rd Qu.: 5.603 3rd Qu.:73.94 3rd Qu.:5.772 3rd Qu.:4.940
## Max. : 7.410 Max. :82.53 Max. :7.430 Max. :6.460
# ------------------------------------------------------------
# Membuat tabel statistik deskriptif
# ------------------------------------------------------------
deskriptif <- data.frame(
Variabel = c("LPE", "IPM", "Kemiskinan", "TPT"),
Mean = c(
mean(data$Y, na.rm = TRUE),
mean(data$X1, na.rm = TRUE),
mean(data$X2, na.rm = TRUE),
mean(data$X3, na.rm = TRUE)
),
SD = c(
sd(data$Y, na.rm = TRUE),
sd(data$X1, na.rm = TRUE),
sd(data$X2, na.rm = TRUE),
sd(data$X3, na.rm = TRUE)
),
Minimum = c(
min(data$Y, na.rm = TRUE),
min(data$X1, na.rm = TRUE),
min(data$X2, na.rm = TRUE),
min(data$X3, na.rm = TRUE)
),
Maximum = c(
max(data$Y, na.rm = TRUE),
max(data$X1, na.rm = TRUE),
max(data$X2, na.rm = TRUE),
max(data$X3, na.rm = TRUE)
)
)
# Membulatkan HANYA kolom numerik
deskriptif[, 2:5] <- round(deskriptif[, 2:5], 3)
# Menampilkan tabel
print(deskriptif)
## Variabel Mean SD Minimum Maximum
## 1 LPE 3.596 2.782 -3.18 7.41
## 2 IPM 72.987 2.905 69.04 82.53
## 3 Kemiskinan 5.046 1.152 3.09 7.43
## 4 TPT 4.070 1.085 1.96 6.46
# 3. MATRIKS KORELASI
# Membuat matriks korelasi
cor_matrix <- cor(
data[, c("Y", "X1", "X2", "X3")],
use = "complete.obs"
)
print(round(cor_matrix, 3))
## Y X1 X2 X3
## Y 1.000 0.138 -0.028 -0.127
## X1 0.138 1.000 -0.380 0.352
## X2 -0.028 -0.380 1.000 -0.203
## X3 -0.127 0.352 -0.203 1.000
# 4. VISUALISASI MATRIKS KORELASI
# Pastikan package corrplot sudah terpasang
if (!requireNamespace("corrplot", quietly = TRUE)) {
install.packages("corrplot")
}
library(corrplot)
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)
)
# 5. SCATTER PLOT MATRIX
pairs(
data[, c("Y", "X1", "X2", "X3")],
main = "Scatter Plot Matrix",
pch = 19
)
# 6. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
# Model regresi
#
# Y = LPE 2024
# X1 = IPM 2024
# X2 = Persentase Penduduk Miskin 2024
# X3 = TPT 2024
model <- lm(
Y ~ X1 + X2 + X3,
data = data
)
# Ringkasan model
summary(model)
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1062 -0.7963 0.3911 1.8445 4.1686
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.33293 9.86029 -0.947 0.347
## X1 0.20360 0.13004 1.566 0.122
## X2 0.02879 0.31355 0.092 0.927
## X3 -0.51031 0.32877 -1.552 0.125
##
## Residual standard error: 2.766 on 66 degrees of freedom
## Multiple R-squared: 0.05425, Adjusted R-squared: 0.01126
## F-statistic: 1.262 on 3 and 66 DF, p-value: 0.2947
# 7. KOEFISIEN REGRESI
coef_model <- summary(model)$coefficients
coef_df <- data.frame(
Variabel = rownames(coef_model),
Koefisien = coef_model[, 1],
Std_Error = coef_model[, 2],
t_value = coef_model[, 3],
p_value = coef_model[, 4],
row.names = NULL
)
# Membulatkan kolom numerik saja
coef_df[, 2:5] <- round(
coef_df[, 2:5],
4
)
print(coef_df)
## Variabel Koefisien Std_Error t_value p_value
## 1 (Intercept) -9.3329 9.8603 -0.9465 0.3473
## 2 X1 0.2036 0.1300 1.5657 0.1222
## 3 X2 0.0288 0.3135 0.0918 0.9271
## 4 X3 -0.5103 0.3288 -1.5522 0.1254
cat("\nInterval Kepercayaan 95%:\n")
##
## Interval Kepercayaan 95%:
print(
round(
confint(model, level = 0.95),
4
)
)
## 2.5 % 97.5 %
## (Intercept) -29.0196 10.3538
## X1 -0.0560 0.4632
## X2 -0.5972 0.6548
## X3 -1.1667 0.1461
# 8. UJI SIGNIFIKANSI SIMULTAN (UJI F)
f_stat <- summary(model)$fstatistic
f_value <- unname(f_stat[1])
df1 <- unname(f_stat[2])
df2 <- unname(f_stat[3])
f_pvalue <- pf(
f_value,
df1,
df2,
lower.tail = FALSE
)
cat("H0 : β1 = β2 = β3 = 0\n")
## H0 : β1 = β2 = β3 = 0
cat("H1 : Minimal terdapat satu βj ≠0\n\n")
## H1 : Minimal terdapat satu βj ≠0
cat(
"F-statistic =",
round(f_value, 4),
"\n"
)
## F-statistic = 1.262
cat(
"df1 =",
df1,
"\n"
)
## df1 = 3
cat(
"df2 =",
df2,
"\n"
)
## df2 = 66
cat(
"p-value =",
format.pval(f_pvalue, digits = 4),
"\n\n"
)
## p-value = 0.2947
if (f_pvalue < 0.05) {
cat(
"Kesimpulan: H0 ditolak. Model signifikan pada taraf 5%.\n"
)
} else {
cat(
"Kesimpulan: H0 tidak ditolak. Model tidak signifikan pada taraf 5%.\n"
)
}
## Kesimpulan: H0 tidak ditolak. Model tidak signifikan pada taraf 5%.
# 9. UJI SIGNIFIKANSI PARSIAL (UJI T)
cat("H0 : βj = 0\n")
## H0 : βj = 0
cat("H1 : βj ≠0\n\n")
## H1 : βj ≠0
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]
cat(
var_name,
": t =",
round(t_val, 4),
", p-value =",
format.pval(
p_val,
digits = 4
),
"\n"
)
}
## X1 : t = 1.5657 , p-value = 0.1222
## X2 : t = 0.0918 , p-value = 0.9271
## X3 : t = -1.5522 , p-value = 0.1254
# 10. KOEFISIEN DETERMINASI
r_squared <- summary(model)$r.squared
adj_r_squared <- summary(model)$adj.r.squared
resid_se <- summary(model)$sigma
cat(
"R-squared =",
round(r_squared, 4),
"\n"
)
## R-squared = 0.0543
cat(
"Adjusted R-squared =",
round(adj_r_squared, 4),
"\n"
)
## Adjusted R-squared = 0.0113
cat(
"Residual Standard Error =",
round(resid_se, 4),
"\n\n"
)
## Residual Standard Error = 2.7659
cat(
"Model menjelaskan sekitar",
round(r_squared * 100, 2),
"% variasi LPE 2024.\n"
)
## Model menjelaskan sekitar 5.43 % variasi LPE 2024.
# 11. UJI NORMALITAS RESIDUAL
residuals_model <- residuals(model)
# Shapiro-Wilk
shapiro_test <- shapiro.test(
residuals_model
)
cat("Shapiro-Wilk Test\n")
## Shapiro-Wilk Test
cat(
"W =",
round(
unname(shapiro_test$statistic),
4
),
"\n"
)
## W = 0.898
cat(
"p-value =",
format.pval(
shapiro_test$p.value,
digits = 4
),
"\n\n"
)
## p-value = 3.225e-05
if (shapiro_test$p.value > 0.05) {
cat(
"Kesimpulan: Tidak terdapat bukti yang cukup bahwa residual tidak berdistribusi normal.\n"
)
} else {
cat(
"Kesimpulan: Terdapat indikasi residual tidak berdistribusi normal.\n"
)
}
## Kesimpulan: Terdapat indikasi residual tidak berdistribusi normal.
# Anderson-Darling
ad_test <- ad.test(
residuals_model
)
cat("\nAnderson-Darling Test\n")
##
## Anderson-Darling Test
cat(
"A =",
round(
unname(ad_test$statistic),
4
),
"\n"
)
## A = 2.3081
cat(
"p-value =",
format.pval(
ad_test$p.value,
digits = 4
),
"\n"
)
## p-value = 6.596e-06
# 12. UJI HETEROSKEDASTISITAS
bp_test <- bptest(model)
cat("Breusch-Pagan Test\n")
## Breusch-Pagan Test
cat(
"BP =",
round(
unname(bp_test$statistic),
4
),
"\n"
)
## BP = 5.0148
cat(
"df =",
unname(bp_test$parameter),
"\n"
)
## df = 3
cat(
"p-value =",
format.pval(
bp_test$p.value,
digits = 4
),
"\n\n"
)
## p-value = 0.1707
if (bp_test$p.value > 0.05) {
cat(
"Kesimpulan: Tidak terdapat bukti yang cukup adanya heteroskedastisitas.\n"
)
} else {
cat(
"Kesimpulan: Terdapat indikasi heteroskedastisitas.\n"
)
}
## Kesimpulan: Tidak terdapat bukti yang cukup adanya heteroskedastisitas.
# 13. UJI MULTIKOLINEARITAS
vif_values <- vif(model)
vif_table <- data.frame(
Variabel = names(vif_values),
VIF = as.numeric(vif_values)
)
vif_table$VIF <- round(
vif_table$VIF,
4
)
print(vif_table)
## Variabel VIF
## 1 X1 1.2872
## 2 X2 1.1762
## 3 X3 1.1487
cat("\nInterpretasi umum:\n")
##
## Interpretasi umum:
cat(
"VIF < 5 : tidak terdapat indikasi multikolinearitas yang kuat\n"
)
## VIF < 5 : tidak terdapat indikasi multikolinearitas yang kuat
cat(
"VIF 5-10 : indikasi multikolinearitas sedang\n"
)
## VIF 5-10 : indikasi multikolinearitas sedang
cat(
"VIF > 10 : indikasi multikolinearitas kuat\n"
)
## VIF > 10 : indikasi multikolinearitas kuat
# 14. UJI AUTOKORELASI
dw_test <- dwtest(model)
cat("Durbin-Watson Test\n")
## Durbin-Watson Test
cat(
"DW =",
round(
unname(dw_test$statistic),
4
),
"\n"
)
## DW = 0.8184
cat(
"p-value =",
format.pval(
dw_test$p.value,
digits = 4
),
"\n\n"
)
## p-value = 1.877e-08
cat(
"Catatan: data merupakan cross-section kabupaten/kota, sehingga autokorelasi bukan fokus utama seperti pada data time series.\n"
)
## Catatan: data merupakan cross-section kabupaten/kota, sehingga autokorelasi bukan fokus utama seperti pada data time series.
# 15. UJI LINEARITAS / SPESIFIKASI MODEL
reset_test <- resettest(
model,
power = 2:3,
type = "fitted"
)
cat("Ramsey RESET Test\n")
## Ramsey RESET Test
cat(
"F =",
round(
unname(reset_test$statistic),
4
),
"\n"
)
## F = 0.6614
cat(
"p-value =",
format.pval(
reset_test$p.value,
digits = 4
),
"\n\n"
)
## p-value = 0.5196
if (reset_test$p.value > 0.05) {
cat(
"Kesimpulan: Tidak terdapat bukti yang cukup adanya kesalahan spesifikasi model.\n"
)
} else {
cat(
"Kesimpulan: Terdapat indikasi kesalahan spesifikasi model.\n"
)
}
## Kesimpulan: Tidak terdapat bukti yang cukup adanya kesalahan spesifikasi model.
# 16. DIAGNOSTIK OUTLIER DAN OBSERVASI INFLUENSIAL
# Cook's Distance
cooks_d <- cooks.distance(model)
cat(
"Cook's Distance maksimum =",
round(max(cooks_d), 4),
"\n"
)
## Cook's Distance maksimum = 0.3111
cat(
"Jumlah observasi dengan Cook's D > 1 =",
sum(cooks_d > 1),
"\n\n"
)
## Jumlah observasi dengan Cook's D > 1 = 0
# Leverage
leverage <- hatvalues(model)
n <- nrow(data)
k <- length(coef(model))
leverage_threshold <- 2 * k / n
cat(
"Leverage maksimum =",
round(max(leverage), 4),
"\n"
)
## Leverage maksimum = 0.1726
cat(
"Threshold leverage =",
round(leverage_threshold, 4),
"\n"
)
## Threshold leverage = 0.1143
cat(
"Jumlah observasi di atas threshold =",
sum(leverage > leverage_threshold),
"\n\n"
)
## Jumlah observasi di atas threshold = 7
# Studentized residual
student_resid <- rstudent(model)
cat(
"Nilai maksimum absolut studentized residual =",
round(
max(abs(student_resid)),
4
),
"\n"
)
## Nilai maksimum absolut studentized residual = 2.9239
cat(
"Jumlah |studentized residual| > 3 =",
sum(abs(student_resid) > 3),
"\n"
)
## Jumlah |studentized residual| > 3 = 0
# 17. VISUALISASI DIAGNOSTIK
par(mfrow = c(2, 2))
# 1. Residual vs Fitted
plot(
fitted(model),
residuals(model),
xlab = "Fitted Values",
ylab = "Residuals",
main = "Residual vs Fitted",
pch = 19
)
abline(
h = 0,
lty = 2
)
lines(
lowess(
fitted(model),
residuals(model)
),
lwd = 2
)
# 2. Normal Q-Q
qqnorm(
residuals(model),
pch = 19,
main = "Normal Q-Q Plot"
)
qqline(
residuals(model),
lwd = 2
)
# 3. Scale-Location
plot(
fitted(model),
sqrt(abs(rstandard(model))),
xlab = "Fitted Values",
ylab = "√|Standardized Residuals|",
main = "Scale-Location",
pch = 19
)
lines(
lowess(
fitted(model),
sqrt(abs(rstandard(model)))
),
lwd = 2
)
# 4. Residual vs Leverage
plot(
hatvalues(model),
rstandard(model),
xlab = "Leverage",
ylab = "Standardized Residuals",
main = "Residual vs Leverage",
pch = 19
)
abline(
h = 0,
lty = 2
)
abline(
v = leverage_threshold,
lty = 2
)
par(mfrow = c(1, 1))
# 18. UJI BOX-COX
library(readxl)
library(MASS)
data_mlr <- read_excel(
"Dataset_MLR_Kalimantan_Tengah_2024.xlsx"
)
# Model regresi untuk Box-Cox
model_bc <- lm(
I(LPE + abs(min(LPE)) + 1) ~ IPM + Kemiskinan + TPT,
data = data_mlr
)
# Uji Box-Cox
bc <- boxcox(
model_bc,
lambda = seq(-2, 2, by = 0.1),
plotit = TRUE
)
# Lambda optimal
lambda_opt <- bc$x[which.max(bc$y)]
cat(
"Lambda optimal =",
round(lambda_opt, 4),
"\n\n"
)
## Lambda optimal = 1.8384
cat("Interpretasi:\n")
## Interpretasi:
if (abs(lambda_opt - 1) < 0.2) {
cat("Lambda mendekati 1 -> tidak perlu transformasi.\n")
} else if (abs(lambda_opt) < 0.2) {
cat("Lambda mendekati 0 -> pertimbangkan transformasi log.\n")
} else if (abs(lambda_opt - 0.5) < 0.2) {
cat("Lambda mendekati 0.5 -> pertimbangkan transformasi akar.\n")
} else {
cat("Lambda tidak mendekati 0, 0.5, atau 1.\n")
cat("Pertimbangkan transformasi Box-Cox sesuai lambda optimal.\n")
}
## Lambda tidak mendekati 0, 0.5, atau 1.
## Pertimbangkan transformasi Box-Cox sesuai lambda optimal.
# 19. PEMILIHAN MODEL DENGAN STEPWISE
model_full <- lm(
Y ~ X1 + X2 + X3,
data = data
)
step_model <- step(
model_full,
direction = "both",
trace = 1
)
## Start: AIC=146.31
## Y ~ X1 + X2 + X3
##
## Df Sum of Sq RSS AIC
## - X2 1 0.0645 504.97 144.32
## <none> 504.91 146.31
## - X3 1 18.4308 523.34 146.82
## - X1 1 18.7543 523.66 146.86
##
## Step: AIC=144.32
## Y ~ X1 + X3
##
## Df Sum of Sq RSS AIC
## <none> 504.97 144.32
## - X3 1 18.7254 523.70 144.87
## - X1 1 20.3241 525.29 145.08
## + X2 1 0.0645 504.91 146.31
cat("\nFormula model hasil stepwise:\n")
##
## Formula model hasil stepwise:
print(
formula(step_model)
)
## Y ~ X1 + X3
cat("\nRingkasan model hasil stepwise:\n")
##
## Ringkasan model hasil stepwise:
print(
summary(step_model)
)
##
## Call:
## lm(formula = Y ~ X1 + X3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1150 -0.7963 0.3922 1.8492 4.1987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.8845 8.5021 -1.045 0.300
## X1 0.1996 0.1215 1.642 0.105
## X3 -0.5127 0.3253 -1.576 0.120
##
## Residual standard error: 2.745 on 67 degrees of freedom
## Multiple R-squared: 0.05413, Adjusted R-squared: 0.0259
## F-statistic: 1.917 on 2 and 67 DF, p-value: 0.155
# 20. STEPWISE BACKWARD
model_backward <- step(
model_full,
direction = "backward",
trace = 1
)
## Start: AIC=146.31
## Y ~ X1 + X2 + X3
##
## Df Sum of Sq RSS AIC
## - X2 1 0.0645 504.97 144.32
## <none> 504.91 146.31
## - X3 1 18.4308 523.34 146.82
## - X1 1 18.7543 523.66 146.86
##
## Step: AIC=144.32
## Y ~ X1 + X3
##
## Df Sum of Sq RSS AIC
## <none> 504.97 144.32
## - X3 1 18.725 523.70 144.87
## - X1 1 20.324 525.29 145.08
cat(
"Formula model backward:\n"
)
## Formula model backward:
print(
formula(model_backward)
)
## Y ~ X1 + X3
cat(
"\nAIC model backward =",
round(
AIC(model_backward),
4
),
"\n"
)
##
## AIC model backward = 344.9718
# 21. STEPWISE FORWARD
model_null <- lm(
Y ~ 1,
data = data
)
model_forward <- step(
model_null,
scope = list(
lower = formula(model_null),
upper = formula(model_full)
),
direction = "forward",
trace = 1
)
## Start: AIC=144.22
## Y ~ 1
##
## Df Sum of Sq RSS AIC
## <none> 533.87 144.22
## + X1 1 10.1730 523.70 144.87
## + X3 1 8.5743 525.29 145.08
## + X2 1 0.4323 533.44 146.16
cat(
"Formula model forward:\n"
)
## Formula model forward:
print(
formula(model_forward)
)
## Y ~ 1
cat(
"\nAIC model forward =",
round(
AIC(model_forward),
4
),
"\n"
)
##
## AIC model forward = 344.8673
# 22. PERBANDINGAN MODEL
# Fungsi menghitung AICc
hitung_aicc <- function(model_object) {
n_model <- nobs(model_object)
k_model <- length(
coef(model_object)
)
if (
(n_model - k_model - 1) <= 0
) {
return(NA)
}
AIC(model_object) +
(
2 * k_model * (k_model + 1)
) /
(
n_model - k_model - 1
)
}
# Membuat tabel perbandingan
perbandingan_model <- data.frame(
Model = c(
"Model Awal",
"Stepwise",
"Backward",
"Forward"
),
Formula = c(
paste(
deparse(formula(model_full)),
collapse = ""
),
paste(
deparse(formula(step_model)),
collapse = ""
),
paste(
deparse(formula(model_backward)),
collapse = ""
),
paste(
deparse(formula(model_forward)),
collapse = ""
)
),
AIC = c(
AIC(model_full),
AIC(step_model),
AIC(model_backward),
AIC(model_forward)
),
BIC = c(
BIC(model_full),
BIC(step_model),
BIC(model_backward),
BIC(model_forward)
),
AICc = c(
hitung_aicc(model_full),
hitung_aicc(step_model),
hitung_aicc(model_backward),
hitung_aicc(model_forward)
)
)
# Membulatkan hanya kolom numerik
perbandingan_model[, 3:5] <-
round(
perbandingan_model[, 3:5],
4
)
print(perbandingan_model)
## Model Formula AIC BIC AICc
## 1 Model Awal Y ~ X1 + X2 + X3 346.9628 358.2053 347.5782
## 2 Stepwise Y ~ X1 + X3 344.9718 353.9657 345.3354
## 3 Backward Y ~ X1 + X3 344.9718 353.9657 345.3354
## 4 Forward Y ~ 1 344.8673 349.3643 344.9261
# 23. PREDIKSI DAN RESIDUAL
data_hasil <- data
data_hasil$Prediksi <-
predict(model_full)
data_hasil$Residual <-
residuals(model_full)
hasil_prediksi <- data_hasil[
,
c(
"Wilayah",
"Y",
"Prediksi",
"Residual"
)
]
# Membulatkan nilai numerik
hasil_prediksi$Y <-
round(
hasil_prediksi$Y,
4
)
hasil_prediksi$Prediksi <-
round(
hasil_prediksi$Prediksi,
4
)
hasil_prediksi$Residual <-
round(
hasil_prediksi$Residual,
4
)
print(hasil_prediksi)
## Wilayah Y Prediksi Residual
## 1 Kotawaringin Barat 0.98 3.3978 -2.4178
## 2 Kotawaringin Timur -3.06 3.0494 -6.1094
## 3 Kapuas -1.04 2.7633 -3.8033
## 4 Barito Selatan -2.90 3.4548 -6.3548
## 5 Barito Utara -2.24 2.6618 -4.9018
## 6 Sukamara 1.98 2.4183 -0.4383
## 7 Lamandau 1.85 4.0140 -2.1640
## 8 Seruyan -2.23 2.7633 -4.9933
## 9 Katingan -3.18 2.6523 -5.8323
## 10 Pulang Pisau 2.68 3.8109 -1.1309
## 11 Gunung Mas 3.39 4.1925 -0.8025
## 12 Barito Timur -2.73 4.2387 -6.9687
## 13 Murung Raya -2.45 3.4121 -5.8621
## 14 Palangka Raya -2.85 4.2562 -7.1062
## 15 Kotawaringin Barat 5.61 3.4795 2.1305
## 16 Kotawaringin Timur 2.10 3.1230 -1.0230
## 17 Kapuas 4.71 2.8405 1.8695
## 18 Barito Selatan 2.13 3.5504 -1.4204
## 19 Barito Utara 2.82 2.7917 0.0283
## 20 Sukamara 4.74 2.5051 2.2349
## 21 Lamandau 4.01 4.3122 -0.3022
## 22 Seruyan 2.12 2.8178 -0.6978
## 23 Katingan 2.90 2.8053 0.0947
## 24 Pulang Pisau 3.24 3.8469 -0.6069
## 25 Gunung Mas 5.09 3.9382 1.1518
## 26 Barito Timur 2.97 4.1052 -1.1352
## 27 Murung Raya 4.38 3.4829 0.8971
## 28 Palangka Raya 4.32 4.3212 -0.0012
## 29 Kotawaringin Barat 6.01 3.6247 2.3853
## 30 Kotawaringin Timur 7.41 3.2414 4.1686
## 31 Kapuas 7.04 3.4331 3.6069
## 32 Barito Selatan 6.28 3.9608 2.3192
## 33 Barito Utara 6.24 3.0663 3.1737
## 34 Sukamara 5.62 1.7013 3.9187
## 35 Lamandau 6.05 3.8474 2.2026
## 36 Seruyan 4.01 3.0737 0.9363
## 37 Katingan 5.58 3.0560 2.5240
## 38 Pulang Pisau 4.68 4.2681 0.4119
## 39 Gunung Mas 6.47 4.0801 2.3899
## 40 Barito Timur 6.06 4.3549 1.7051
## 41 Murung Raya 7.03 3.7164 3.3136
## 42 Palangka Raya 6.25 4.4804 1.7696
## 43 Kotawaringin Barat 6.10 3.7704 2.3296
## 44 Kotawaringin Timur 1.81 3.4613 -1.6513
## 45 Kapuas 5.71 3.6902 2.0198
## 46 Barito Selatan 3.27 3.6619 -0.3919
## 47 Barito Utara 5.49 3.1500 2.3400
## 48 Sukamara 5.64 2.4356 3.2044
## 49 Lamandau 1.59 4.0152 -2.4252
## 50 Seruyan 4.55 3.3309 1.2191
## 51 Katingan 5.98 3.3258 2.6542
## 52 Pulang Pisau 4.84 4.3246 0.5154
## 53 Gunung Mas 4.25 4.0708 0.1792
## 54 Barito Timur 3.47 4.2476 -0.7776
## 55 Murung Raya 5.46 3.8866 1.5734
## 56 Palangka Raya 6.57 4.8335 1.7365
## 57 Kotawaringin Barat 4.10 3.8713 0.2287
## 58 Kotawaringin Timur 4.00 3.6296 0.3704
## 59 Kapuas 4.95 3.8349 1.1151
## 60 Barito Selatan 4.70 3.9250 0.7750
## 61 Barito Utara 5.08 3.3244 1.7556
## 62 Sukamara 3.89 2.6814 1.2086
## 63 Lamandau 3.64 4.1994 -0.5594
## 64 Seruyan 3.04 3.4867 -0.4467
## 65 Katingan 4.67 3.4701 1.1999
## 66 Pulang Pisau 4.41 4.5155 -0.1055
## 67 Gunung Mas 4.48 4.2806 0.1994
## 68 Barito Timur 4.29 4.4267 -0.1367
## 69 Murung Raya 5.05 3.9505 1.0995
## 70 Palangka Raya 6.62 5.0100 1.6100
# 24. PERSAMAAN REGRESI
b <- coef(model_full)
cat(
"Persamaan regresi:\n\n"
)
## Persamaan regresi:
cat(
"Y = ",
round(b[1], 4),
" + (",
round(b[2], 4),
")X1 + (",
round(b[3], 4),
")X2 + (",
round(b[4], 4),
")X3\n",
sep = ""
)
## Y = -9.3329 + (0.2036)X1 + (0.0288)X2 + (-0.5103)X3
cat("\nKeterangan:\n")
##
## Keterangan:
cat(
"Y = LPE 2024\n"
)
## Y = LPE 2024
cat(
"X1 = IPM 2024\n"
)
## X1 = IPM 2024
cat(
"X2 = Persentase Penduduk Miskin 2024\n"
)
## X2 = Persentase Penduduk Miskin 2024
cat(
"X3 = TPT 2024\n"
)
## X3 = TPT 2024
# 25. RINGKASAN HASIL ANALISIS
cat("\n--- MODEL ---\n")
##
## --- MODEL ---
cat(
"Y = ",
round(b[1], 4),
" + ",
round(b[2], 4),
"X1 + ",
round(b[3], 4),
"X2 + ",
round(b[4], 4),
"X3\n",
sep = ""
)
## Y = -9.3329 + 0.2036X1 + 0.0288X2 + -0.5103X3
cat("\n--- GOODNESS OF FIT ---\n")
##
## --- GOODNESS OF FIT ---
cat(
"R-squared =",
round(r_squared, 4),
"\n"
)
## R-squared = 0.0543
cat(
"Adjusted R-squared =",
round(adj_r_squared, 4),
"\n"
)
## Adjusted R-squared = 0.0113
cat(
"AIC =",
round(
AIC(model_full),
4
),
"\n"
)
## AIC = 346.9628
cat(
"BIC =",
round(
BIC(model_full),
4
),
"\n"
)
## BIC = 358.2053
cat(
"AICc =",
round(
hitung_aicc(model_full),
4
),
"\n"
)
## AICc = 347.5782
cat("\n--- UJI SIMULTAN ---\n")
##
## --- UJI SIMULTAN ---
cat(
"Uji F p-value =",
format.pval(
f_pvalue,
digits = 4
),
"\n"
)
## Uji F p-value = 0.2947
cat("\n--- UJI PARSIAL ---\n")
##
## --- UJI PARSIAL ---
for (i in 2:nrow(coef_df)) {
cat(
coef_df$Variabel[i],
": p-value =",
format.pval(
coef_df$p_value[i],
digits = 4
),
"\n"
)
}
## X1 : p-value = 0.1222
## X2 : p-value = 0.9271
## X3 : p-value = 0.1254
cat("\n--- UJI ASUMSI ---\n")
##
## --- UJI ASUMSI ---
cat(
"Normalitas Shapiro-Wilk p-value =",
format.pval(
shapiro_test$p.value,
digits = 4
),
"\n"
)
## Normalitas Shapiro-Wilk p-value = 3.225e-05
cat(
"Heteroskedastisitas Breusch-Pagan p-value =",
format.pval(
bp_test$p.value,
digits = 4
),
"\n"
)
## Heteroskedastisitas Breusch-Pagan p-value = 0.1707
cat(
"VIF maksimum =",
round(
max(vif_values),
4
),
"\n"
)
## VIF maksimum = 1.2872
cat(
"Ramsey RESET p-value =",
format.pval(
reset_test$p.value,
digits = 4
),
"\n"
)
## Ramsey RESET p-value = 0.5196
cat(
"\nAnalisis Multiple Linear Regression selesai.\n"
)
##
## Analisis Multiple Linear Regression selesai.