1. 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(corrplot)
## corrplot 0.95 loaded
2. Import Data
# Mencari file CSV secara otomatis
file_csv <- list.files(
pattern = "\\.csv$",
full.names = TRUE
)
# Memastikan hanya ada satu file CSV yang digunakan
if (length(file_csv) == 0) {
stop("File CSV tidak ditemukan. Pastikan file CSV berada di folder yang sama dengan file Rmd.")
}
# Jika terdapat lebih dari satu file CSV
if (length(file_csv) > 1) {
stop("Terdapat lebih dari satu file CSV. Pastikan hanya file dataset yang digunakan berada di folder yang sama dengan file Rmd.")
}
# Membaca data
data <- read.csv(
file_csv[1],
header = TRUE,
check.names = FALSE,
stringsAsFactors = FALSE
)
# Informasi data
cat("File yang digunakan:", basename(file_csv[1]), "\n")
## File yang digunakan: HDI Unemployment and Education - Data Indonesia Provinces (2017-2023).csv
cat("Jumlah baris:", nrow(data), "\n")
## Jumlah baris: 245
cat("Jumlah kolom:", ncol(data), "\n")
## Jumlah kolom: 8
3. Pemeriksaan Data
# Melihat 6 baris pertama
head(data)
## Provinsi Tingkat Pengangguran Terbuka (TPT) - Februari
## 1 ACEH 7.39
## 2 SUMATERA UTARA 6.41
## 3 SUMATERA BARAT 5.80
## 4 RIAU 5.76
## 5 JAMBI 3.67
## 6 SUMATERA SELATAN 3.80
## Tingkat Pengangguran Terbuka (TPT) - Agustus
## 1 6.57
## 2 5.60
## 3 5.58
## 4 6.22
## 5 3.87
## 6 4.39
## Tingkat Partisipasi Angkatan Kerja (TPAK) - Februari
## 1 65.59
## 2 69.13
## 3 70.42
## 4 68.42
## 5 70.84
## 6 72.12
## Tingkat Partisipasi Angkatan Kerja (TPAK) - Agustus Tahun
## 1 63.74 2017
## 2 68.88 2017
## 3 66.29 2017
## 4 64.00 2017
## 5 67.52 2017
## 6 69.50 2017
## Indeks Pembangunan Manusia Rata-rata Lama Sekolah
## 1 70.60 8.98
## 2 70.57 9.25
## 3 71.24 8.72
## 4 71.79 8.76
## 5 69.99 8.15
## 6 68.86 7.99
# Melihat struktur data
str(data)
## 'data.frame': 245 obs. of 8 variables:
## $ Provinsi : chr "ACEH" "SUMATERA UTARA" "SUMATERA BARAT" "RIAU" ...
## $ Tingkat Pengangguran Terbuka (TPT) - Februari : num 7.39 6.41 5.8 5.76 3.67 3.8 2.81 4.43 4.46 6.44 ...
## $ Tingkat Pengangguran Terbuka (TPT) - Agustus : num 6.57 5.6 5.58 6.22 3.87 4.39 3.74 4.33 3.78 7.16 ...
## $ Tingkat Partisipasi Angkatan Kerja (TPAK) - Februari: num 65.6 69.1 70.4 68.4 70.8 ...
## $ Tingkat Partisipasi Angkatan Kerja (TPAK) - Agustus : num 63.7 68.9 66.3 64 67.5 ...
## $ Tahun : int 2017 2017 2017 2017 2017 2017 2017 2017 2017 2017 ...
## $ Indeks Pembangunan Manusia : num 70.6 70.6 71.2 71.8 70 ...
## $ Rata-rata Lama Sekolah : num 8.98 9.25 8.72 8.76 8.15 7.99 8.47 7.79 7.78 9.79 ...
# Melihat nama variabel
names(data)
## [1] "Provinsi"
## [2] "Tingkat Pengangguran Terbuka (TPT) - Februari"
## [3] "Tingkat Pengangguran Terbuka (TPT) - Agustus"
## [4] "Tingkat Partisipasi Angkatan Kerja (TPAK) - Februari"
## [5] "Tingkat Partisipasi Angkatan Kerja (TPAK) - Agustus"
## [6] "Tahun"
## [7] "Indeks Pembangunan Manusia"
## [8] "Rata-rata Lama Sekolah"
# Melihat ukuran data
dim(data)
## [1] 245 8
# Melihat ringkasan data
summary(data)
## Provinsi Tingkat Pengangguran Terbuka (TPT) - Februari
## Length :245 Min. : 0.880
## N.unique : 35 1st Qu.: 3.590
## N.blank : 0 Median : 4.540
## Min.nchar: 4 Mean : 4.827
## Max.nchar: 20 3rd Qu.: 5.800
## Max. :10.120
## Tingkat Pengangguran Terbuka (TPT) - Agustus
## Min. : 1.400
## 1st Qu.: 3.870
## Median : 4.740
## Mean : 5.124
## 3rd Qu.: 6.280
## Max. :10.950
## Tingkat Partisipasi Angkatan Kerja (TPAK) - Februari
## Min. :61.85
## 1st Qu.:67.23
## Median :69.74
## Mean :69.60
## 3rd Qu.:72.01
## Max. :80.23
## Tingkat Partisipasi Angkatan Kerja (TPAK) - Agustus Tahun
## Min. :60.18 Min. :2017
## 1st Qu.:65.46 1st Qu.:2018
## Median :68.40 Median :2020
## Mean :68.11 Mean :2020
## 3rd Qu.:69.83 3rd Qu.:2022
## Max. :79.02 Max. :2023
## Indeks Pembangunan Manusia Rata-rata Lama Sekolah
## Min. :59.09 Min. : 6.270
## 1st Qu.:69.55 1st Qu.: 7.960
## Median :71.30 Median : 8.610
## Mean :71.28 Mean : 8.611
## 3rd Qu.:72.79 3rd Qu.: 9.220
## Max. :82.46 Max. :11.450
# Memeriksa data kosong
colSums(is.na(data))
## Provinsi
## 0
## Tingkat Pengangguran Terbuka (TPT) - Februari
## 0
## Tingkat Pengangguran Terbuka (TPT) - Agustus
## 0
## Tingkat Partisipasi Angkatan Kerja (TPAK) - Februari
## 0
## Tingkat Partisipasi Angkatan Kerja (TPAK) - Agustus
## 0
## Tahun
## 0
## Indeks Pembangunan Manusia
## 0
## Rata-rata Lama Sekolah
## 0
————————————————————
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)
## Provinsi Tingkat Pengangguran Terbuka (TPT) - Februari
## Length :245 Min. : 0.880
## N.unique : 35 1st Qu.: 3.590
## N.blank : 0 Median : 4.540
## Min.nchar: 4 Mean : 4.827
## Max.nchar: 20 3rd Qu.: 5.800
## Max. :10.120
## Tingkat Pengangguran Terbuka (TPT) - Agustus
## Min. : 1.400
## 1st Qu.: 3.870
## Median : 4.740
## Mean : 5.124
## 3rd Qu.: 6.280
## Max. :10.950
## Tingkat Partisipasi Angkatan Kerja (TPAK) - Februari
## Min. :61.85
## 1st Qu.:67.23
## Median :69.74
## Mean :69.60
## 3rd Qu.:72.01
## Max. :80.23
## Tingkat Partisipasi Angkatan Kerja (TPAK) - Agustus Tahun
## Min. :60.18 Min. :2017
## 1st Qu.:65.46 1st Qu.:2018
## Median :68.40 Median :2020
## Mean :68.11 Mean :2020
## 3rd Qu.:69.83 3rd Qu.:2022
## Max. :79.02 Max. :2023
## Indeks Pembangunan Manusia Rata-rata Lama Sekolah
## Min. :59.09 Min. : 6.270
## 1st Qu.:69.55 1st Qu.: 7.960
## Median :71.30 Median : 8.610
## Mean :71.28 Mean : 8.611
## 3rd Qu.:72.79 3rd Qu.: 9.220
## Max. :82.46 Max. :11.450
# Memilih variabel yang digunakan dalam model
data_korelasi <- data.frame(
TPT_Ags = data[[3]],
TPAK_Ags = data[[5]],
IPM = data[[7]],
RLS = data[[8]]
)
# Mengubah nama variabel agar lebih sederhana
names(data_korelasi) <- c(
"TPT_Ags",
"TPAK_Ags",
"IPM",
"RLS"
)
# Matriks korelasi
cat("\nMatriks Korelasi:\n")
##
## Matriks Korelasi:
cor_matrix <- cor(data_korelasi)
# Menampilkan matriks korelasi
print(round(cor_matrix, 3))
## TPT_Ags TPAK_Ags IPM RLS
## TPT_Ags 1.000 -0.619 0.326 0.474
## TPAK_Ags -0.619 1.000 -0.221 -0.424
## IPM 0.326 -0.221 1.000 0.779
## RLS 0.474 -0.424 0.779 1.000
# Visualisasi matriks korelasi
corrplot(
cor_matrix,
method = "color",
type = "upper",
col = colorRampPalette(c("blue", "white", "red"))(200),
addCoef.col = "darkblue",
tl.col = "darkblue",
tl.srt = 0,
tl.cex = 0.9,
title = "Matriks Korelasi Antar Variabel",
mar = c(0, 0, 2, 0)
)
# Scatter plot matrix
pairs(
data_korelasi,
main = "Scatter Plot Matrix",
pch = 19
)
============================================================
3. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
============================================================
cat("\n============================================================\n")
##
## ============================================================
cat("3. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION\n")
## 3. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
cat("============================================================\n")
## ============================================================
# Membuat data yang digunakan dalam model
data_model <- data.frame(
TPT_Ags = data[["Tingkat Pengangguran Terbuka (TPT) - Agustus"]],
TPAK_Ags = data[["Tingkat Partisipasi Angkatan Kerja (TPAK) - Agustus"]],
IPM = data[["Indeks Pembangunan Manusia"]],
RLS = data[["Rata-rata Lama Sekolah"]]
)
# Fit model MLR
model <- lm(
TPT_Ags ~ TPAK_Ags + IPM + RLS,
data = data_model
)
# Ringkasan 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, row.names = FALSE)
## Variabel Koefisien Std_Error t_value p_value
## (Intercept) 18.0189 2.5801 6.9837 0.0000
## TPAK_Ags -0.2588 0.0273 -9.4665 0.0000
## IPM 0.0145 0.0355 0.4100 0.6822
## RLS 0.4294 0.1570 2.7343 0.0067
# Interval kepercayaan 95%
cat("\nInterval Kepercayaan 95%:\n")
##
## Interval Kepercayaan 95%:
ci <- round(confint(model, level = 0.95), 4)
colnames(ci) <- c("Batas_Bawah", "Batas_Atas")
print(ci)
## Batas_Bawah Batas_Atas
## (Intercept) 12.9364 23.1014
## TPAK_Ags -0.3127 -0.2050
## IPM -0.0553 0.0844
## RLS 0.1200 0.7387
# Persamaan regresi
cat("\nPersamaan Regresi Linear Berganda:\n")
##
## Persamaan Regresi Linear Berganda:
cat(
"TPT_Ags = 18.0189 - 0.2588TPAK_Ags + 0.0145IPM + 0.4294RLS\n"
)
## TPT_Ags = 18.0189 - 0.2588TPAK_Ags + 0.0145IPM + 0.4294RLS
————————————————————
4. UJI SIGNIFIKANSI SIMULTAN (UJI F)
————————————————————
cat("\n============================================================\n")
##
## ============================================================
cat("4. UJI SIGNIFIKANSI SIMULTAN (UJI F)\n")
## 4. UJI SIGNIFIKANSI SIMULTAN (UJI F)
cat("============================================================\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 = 0 (Model tidak signifikan)\n")
## H0: β1 = β2 = β3 = 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: 62.6443
cat("df1:", df1, "\n")
## df1: 3
cat("df2:", df2, "\n")
## df2: 241
cat("p-value:", format(f_pvalue, scientific = TRUE), "\n")
## p-value: 5.591516e-30
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============================================================\n")
##
## ============================================================
cat("5. UJI SIGNIFIKANSI PARSIAL (UJI t)\n")
## 5. UJI SIGNIFIKANSI PARSIAL (UJI t)
cat("============================================================\n")
## ============================================================
# Ekstrak koefisien dan hasil uji t
t_test <- summary(model)$coefficients
cat("\nHipotesis:\n")
##
## Hipotesis:
cat("H0: βj = 0 (Variabel tidak berpengaruh signifikan)\n")
## H0: βj = 0 (Variabel tidak berpengaruh signifikan)
cat("H1: βj ≠0 (Variabel berpengaruh signifikan)\n")
## H1: βj ≠0 (Variabel berpengaruh signifikan)
cat("\nHasil Uji t:\n")
##
## Hasil Uji t:
t_table <- data.frame(
Variabel = rownames(t_test),
Koefisien = round(t_test[, 1], 4),
t_value = round(t_test[, 3], 4),
p_value = round(t_test[, 4], 4)
)
print(t_table, row.names = FALSE)
## Variabel Koefisien t_value p_value
## (Intercept) 18.0189 6.9837 0.0000
## TPAK_Ags -0.2588 -9.4665 0.0000
## IPM 0.0145 0.4100 0.6822
## RLS 0.4294 2.7343 0.0067
cat("\nKeputusan pada α = 5%:\n")
##
## Keputusan pada α = 5%:
for (i in 2:nrow(t_test)) {
if (t_test[i, 4] < 0.05) {
cat(rownames(t_test)[i], ": Tolak H0 → signifikan\n")
} else {
cat(rownames(t_test)[i], ": Gagal Tolak H0 → tidak signifikan\n")
}
}
## TPAK_Ags : Tolak H0 → signifikan
## IPM : Gagal Tolak H0 → tidak signifikan
## RLS : Tolak H0 → signifikan
————————————————————
6. KOEFISIEN DETERMINASI
————————————————————
cat("\n============================================================\n")
##
## ============================================================
cat("6. KOEFISIEN DETERMINASI\n")
## 6. KOEFISIEN DETERMINASI
cat("============================================================\n")
## ============================================================
# Ekstrak R-squared dan Adjusted R-squared
r_squared <- summary(model)$r.squared
adj_r_squared <- summary(model)$adj.r.squared
cat("\nR-squared:", round(r_squared, 4), "\n")
##
## R-squared: 0.4381
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.4311
cat("\nInterpretasi:\n")
##
## Interpretasi:
cat(
"R-squared menunjukkan proporsi variasi TPT yang dapat dijelaskan ",
"oleh variabel TPAK, IPM, dan RLS dalam model.\n"
)
## R-squared menunjukkan proporsi variasi TPT yang dapat dijelaskan oleh variabel TPAK, IPM, dan RLS dalam model.
————————————————————
7. UJI ASUMSI RESIDUAL
————————————————————
cat("\n============================================================\n")
##
## ============================================================
cat("7. UJI ASUMSI RESIDUAL\n")
## 7. UJI ASUMSI RESIDUAL
cat("============================================================\n")
## ============================================================
# 7.1 Uji Normalitas Residual - Shapiro-Wilk
cat("\n7.1 UJI NORMALITAS RESIDUAL (SHAPIRO-WILK)\n")
##
## 7.1 UJI NORMALITAS RESIDUAL (SHAPIRO-WILK)
shapiro_test <- shapiro.test(residuals(model))
cat("W =", round(shapiro_test$statistic, 4), "\n")
## W = 0.953
cat("p-value =", format(shapiro_test$p.value, scientific = TRUE), "\n")
## p-value = 3.934847e-07
if (shapiro_test$p.value > 0.05) {
cat("Kesimpulan: Gagal Tolak H0 → residual berdistribusi normal\n")
} else {
cat("Kesimpulan: Tolak H0 → residual tidak berdistribusi normal\n")
}
## Kesimpulan: Tolak H0 → residual tidak berdistribusi normal
# Q-Q Plot Residual
qqnorm(
residuals(model),
main = "Q-Q Plot Residual"
)
qqline(
residuals(model),
col = "red",
lwd = 2
)
# 7.2 UJI HETEROSKEDASTISITAS (BREUSCH-PAGAN)
cat("\n7.2 UJI HETEROSKEDASTISITAS (BREUSCH-PAGAN)\n")
##
## 7.2 UJI HETEROSKEDASTISITAS (BREUSCH-PAGAN)
bp_test <- bptest(model)
cat("BP =", round(bp_test$statistic, 4), "\n")
## BP = 2.7667
cat("df =", bp_test$parameter, "\n")
## df = 3
cat("p-value =", format(bp_test$p.value, scientific = TRUE), "\n")
## p-value = 4.29011e-01
if (bp_test$p.value > 0.05) {
cat("Kesimpulan: Gagal Tolak H0 → tidak terdapat heteroskedastisitas\n")
} else {
cat("Kesimpulan: Tolak H0 → terdapat heteroskedastisitas\n")
}
## Kesimpulan: Gagal Tolak H0 → tidak terdapat heteroskedastisitas
# 7.3 UJI AUTOKORELASI (DURBIN-WATSON)
cat("\n7.3 UJI AUTOKORELASI (DURBIN-WATSON)\n")
##
## 7.3 UJI AUTOKORELASI (DURBIN-WATSON)
dw_test <- dwtest(model)
cat("DW =", round(dw_test$statistic, 4), "\n")
## DW = 1.4364
cat("p-value =", format(dw_test$p.value, scientific = TRUE), "\n")
## p-value = 3.958082e-06
cat("\nKesimpulan:\n")
##
## Kesimpulan:
if (dw_test$p.value > 0.05) {
cat("Gagal Tolak H0 → tidak terdapat autokorelasi\n")
} else {
cat("Tolak H0 → terdapat autokorelasi\n")
}
## Tolak H0 → terdapat autokorelasi
# 7.4 UJI MULTIKOLINEARITAS (VIF)
cat("\n7.4 UJI MULTIKOLINEARITAS (VIF)\n")
##
## 7.4 UJI MULTIKOLINEARITAS (VIF)
vif_value <- vif(model)
cat("\nNilai VIF:\n")
##
## Nilai VIF:
print(round(vif_value, 4))
## TPAK_Ags IPM RLS
## 1.2653 2.6378 3.0579
cat("\nKesimpulan:\n")
##
## Kesimpulan:
if (all(vif_value < 10)) {
cat("Seluruh variabel memiliki VIF < 10 → tidak terdapat masalah multikolinearitas\n")
} else {
cat("Terdapat variabel dengan VIF >= 10 → terdapat indikasi multikolinearitas\n")
}
## Seluruh variabel memiliki VIF < 10 → tidak terdapat masalah multikolinearitas
# 7.5 UJI SPESIFIKASI MODEL (RAMSEY RESET)
cat("\n7.5 UJI SPESIFIKASI MODEL (RAMSEY RESET)\n")
##
## 7.5 UJI SPESIFIKASI MODEL (RAMSEY RESET)
reset_test <- resettest(model, power = 2:3, type = "fitted")
cat("\nHasil Uji Ramsey RESET:\n")
##
## Hasil Uji Ramsey RESET:
cat("F =", round(reset_test$statistic, 4), "\n")
## F = 7.2353
cat("df1 =", reset_test$parameter[1], "\n")
## df1 = 2
cat("df2 =", reset_test$parameter[2], "\n")
## df2 = 239
cat("p-value =", format(round(reset_test$p.value, 5), nsmall = 5), "\n")
## p-value = 0.00089
cat("\nKesimpulan:\n")
##
## Kesimpulan:
if (reset_test$p.value > 0.05) {
cat("Gagal Tolak H0 → tidak terdapat indikasi kesalahan spesifikasi model\n")
} else {
cat("Tolak H0 → terdapat indikasi kesalahan spesifikasi model\n")
}
## Tolak H0 → terdapat indikasi kesalahan spesifikasi model
============================================================
8. DIAGNOSTIK OUTLIER DAN OBSERVASI BERPENGARUH
============================================================
cat("\n============================================================\n")
##
## ============================================================
cat("8. DIAGNOSTIK OUTLIER DAN OBSERVASI BERPENGARUH\n")
## 8. DIAGNOSTIK OUTLIER DAN OBSERVASI BERPENGARUH
cat("============================================================\n")
## ============================================================
n <- nrow(data_model)
# Cook's Distance
cook_dist <- cooks.distance(model)
cat("\nCook's Distance:\n")
##
## Cook's Distance:
cat("Nilai maksimum:", round(max(cook_dist), 4), "\n")
## Nilai maksimum: 0.0673
cat("Threshold 4/n:", round(4 / n, 4), "\n")
## Threshold 4/n: 0.0163
cat("Jumlah observasi dengan Cook's Distance > threshold:",
sum(cook_dist > 4 / n), "\n")
## Jumlah observasi dengan Cook's Distance > threshold: 15
# 8.2 LEVERAGE
leverage <- hatvalues(model)
cat("\nLeverage (Hat Values):\n")
##
## Leverage (Hat Values):
cat("Nilai maksimum:", round(max(leverage), 4), "\n")
## Nilai maksimum: 0.0722
cat("Threshold 2(k+1)/n:", round(2 * (length(coef(model))) / n, 4), "\n")
## Threshold 2(k+1)/n: 0.0327
cat("Jumlah observasi dengan leverage > threshold:",
sum(leverage > 2 * length(coef(model)) / n), "\n")
## Jumlah observasi dengan leverage > threshold: 30
# 8.3 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: 3.473
cat("Jumlah observasi dengan |std_resid| > 3:",
sum(abs(std_resid) > 3), "\n")
## Jumlah observasi dengan |std_resid| > 3: 5
# 8.4 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
————————————————————
# Nilai yang digunakan untuk visualisasi
fitted_values <- fitted(model)
residuals <- residuals(model)
standardized_resid <- rstandard(model)
leverage <- hatvalues(model)
n <- nrow(data_model)
# 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. TRANSFORMASI BOX-COX
============================================================
library(MASS)
# Mencari nilai lambda optimal
boxcox_result <- boxcox(model,
lambda = seq(-2, 2, by = 0.1),
plotit = TRUE)
# Nilai lambda dengan log-likelihood maksimum
lambda_opt <- boxcox_result$x[
which.max(boxcox_result$y)
]
cat("Lambda optimal:", round(lambda_opt, 4), "\n")
## Lambda optimal: 0.2626
============================================================
11. STEPWISE SELECTION
============================================================
stepwise_model <- step(
model,
direction = "both",
trace = 0
)
cat("\nModel hasil Stepwise Selection:\n")
##
## Model hasil Stepwise Selection:
print(formula(stepwise_model))
## TPT_Ags ~ TPAK_Ags + RLS
cat("\nRingkasan Model:\n")
##
## Ringkasan Model:
summary(stepwise_model)
##
## Call:
## lm(formula = TPT_Ags ~ TPAK_Ags + RLS, data = data_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3292 -0.9741 -0.2335 0.7439 4.4515
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.47892 2.31938 7.967 6.36e-14 ***
## TPAK_Ags -0.25669 0.02679 -9.582 < 2e-16 ***
## RLS 0.47928 0.09896 4.843 2.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.311 on 242 degrees of freedom
## Multiple R-squared: 0.4377, Adjusted R-squared: 0.4331
## F-statistic: 94.21 on 2 and 242 DF, p-value: < 2.2e-16
============================================================
12. RINGKASAN
============================================================
cat("\n============================================================\n")
##
## ============================================================
cat("12. RINGKASAN HASIL\n")
## 12. RINGKASAN HASIL
cat("============================================================\n")
## ============================================================
cat("\nJumlah observasi:", nrow(data_model), "\n")
##
## Jumlah observasi: 245
cat("Jumlah variabel prediktor:", length(coef(model)) - 1, "\n")
## Jumlah variabel prediktor: 3
cat("\nModel regresi:\n")
##
## Model regresi:
cat("TPT_Ags = 18.0189 - 0.2588TPAK_Ags + 0.0145IPM + 0.4294RLS\n")
## TPT_Ags = 18.0189 - 0.2588TPAK_Ags + 0.0145IPM + 0.4294RLS
cat("\nKoefisien Determinasi:\n")
##
## Koefisien Determinasi:
cat("R-squared =", round(summary(model)$r.squared, 4), "\n")
## R-squared = 0.4381
cat("Adjusted R-squared =", round(summary(model)$adj.r.squared, 4), "\n")
## Adjusted R-squared = 0.4311
cat("\nUji F:\n")
##
## Uji F:
cat("F-statistic =", round(summary(model)$fstatistic[1], 4), "\n")
## F-statistic = 62.6443
cat("\nHasil Uji t:\n")
##
## Hasil Uji t:
cat("TPAK_Ags: signifikan\n")
## TPAK_Ags: signifikan
cat("IPM: tidak signifikan\n")
## IPM: tidak signifikan
cat("RLS: signifikan\n")
## RLS: signifikan
cat("\nHasil Uji Asumsi:\n")
##
## Hasil Uji Asumsi:
cat("Normalitas residual: tidak terpenuhi\n")
## Normalitas residual: tidak terpenuhi
cat("Heteroskedastisitas: tidak terdeteksi\n")
## Heteroskedastisitas: tidak terdeteksi
cat("Autokorelasi: terindikasi, dengan catatan data bersifat panel-like\n")
## Autokorelasi: terindikasi, dengan catatan data bersifat panel-like
cat("Multikolinearitas: tidak terdeteksi\n")
## Multikolinearitas: tidak terdeteksi
cat("Ramsey RESET: terdapat indikasi masalah spesifikasi\n")
## Ramsey RESET: terdapat indikasi masalah spesifikasi
cat("\nDiagnostik Outlier:\n")
##
## Diagnostik Outlier:
cat("Observasi dengan Cook's Distance > 4/n: 15\n")
## Observasi dengan Cook's Distance > 4/n: 15
cat("Observasi dengan leverage > threshold: 30\n")
## Observasi dengan leverage > threshold: 30
cat("Observasi dengan |Studentized Residual| > 3: 5\n")
## Observasi dengan |Studentized Residual| > 3: 5
cat("Observasi dengan |DFBETAS| > 1: 0\n")
## Observasi dengan |DFBETAS| > 1: 0
cat("\nBox-Cox:\n")
##
## Box-Cox:
cat("Lambda optimal =", round(lambda_opt, 4), "\n")
## Lambda optimal = 0.2626
cat("\nStepwise Selection:\n")
##
## Stepwise Selection:
cat("Model terpilih: TPT_Ags ~ TPAK_Ags + RLS\n")
## Model terpilih: TPT_Ags ~ TPAK_Ags + RLS
cat("Variabel IPM dikeluarkan dari model.\n")
## Variabel IPM dikeluarkan dari model.