Author:
- G1401221027 Nabil Ibni Nawawi
1. LIBRARY
if(!require(pacman)) install.packages("pacman")## Loading required package: pacman## Warning: package 'pacman' was built under R version 4.3.3pacman::p_load(readxl, sf, spdep, plm, splm, car, tidyverse, lmtest, tseries, GGally, gridExtra, ggplot2, dplyr, tidyr, DT, corrplot, dlookr, ggrepel, imputeTS, psych, ggtext, MASS, ggpubr, stargazer, caret, reshape2, sp, kableExtra, nortest)2. INPUT DATA & PETA
Import Data Excel
data_panel <- read_excel("C:/Users/Nabil Ibni Nawawi/Documents/KULIAH/Semester 8/Data/Data Skripsi tambah peubah operasional.xlsx")Rename nama peubah
data_panel <- data_panel %>%
rename(
UMK = Upah_minimum,
AMA = Persen_aksara,
PPM = Persen_miskin,
AngKer = Angka_ketergantungan,
LPP = Ljprtmbhn_penduduk
)
names(data_panel)## [1] "id" "Kab_Kota" "Tahun" "TPT" "TPAK"
## [6] "IPM" "PPM" "UMK" "AngKer" "AMA"
## [11] "GR" "LPP" "PDRB" "JumlahSMP" "JumlahSMA"
## [16] "JumlahFasKes"Import Peta (Shapefile)
peta_indonesia <- st_read("C:/Users/Nabil Ibni Nawawi/Documents/KULIAH/Semester 8/Batas kabkot/file_1758080022479_4916116/Batas Kabupaten.shp", quiet = TRUE)Mendefinisikan 27 Kab/Kota di Provinsi Jawa Barat
target_jabar <- c(
"Bogor", "Sukabumi", "Cianjur", "Bandung", "Garut", "Tasikmalaya",
"Ciamis", "Kuningan", "Cirebon", "Majalengka", "Sumedang", "Indramayu",
"Subang", "Purwakarta", "Karawang", "Bekasi", "Bandung Barat", "Pangandaran",
"Kota Bogor", "Kota Sukabumi", "Kota Bandung", "Kota Cirebon",
"Kota Bekasi", "Kota Depok", "Kota Cimahi", "Kota Tasikmalaya", "Kota Banjar"
)Cek Format Nama pada shp (DEBUGGING)
untuk melihat format penulisannya: Apakah “Kabupaten Bandung”, “KAB. BANDUNG”, atau “Bandung”
contoh_nama <- grep("Bandung|Bogor", unique(peta_indonesia$WADMKK), value = TRUE, ignore.case = TRUE)
cat("=== FORMAT NAMA DI FILE SHP ===\n")## === FORMAT NAMA DI FILE SHP ===print(contoh_nama)## [1] "Bandung" "Bandung Barat" "Bogor" "Kota Bandung"
## [5] "Kota Bogor"Filter hanya Kab/Kota di Provinsi Jawa Barat
peta_jabar <- peta_indonesia %>%
filter(WADMKK %in% target_jabar)Menyamakan urutan baris di Excel dan poligon di Peta
Rename WADMKK menjadi Kab_Kota
colnames(peta_jabar)[which(names(peta_jabar) == "WADMKK")] <- "Kab_Kota"
names(peta_jabar)## [1] "Kab_Kota" "geometry"data_panel <- data_panel %>% arrange(Kab_Kota, Tahun)
peta_jabar <- peta_jabar %>% arrange(Kab_Kota)print(peta_jabar, n = 27)## Simple feature collection with 27 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: 106.3703 ymin: -7.82099 xmax: 108.8468 ymax: -5.806538
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
## Kab_Kota geometry
## 1 Bandung MULTIPOLYGON Z (((107.75 -6...
## 2 Bandung Barat MULTIPOLYGON Z (((107.4121 ...
## 3 Bekasi MULTIPOLYGON Z (((107.0159 ...
## 4 Bogor MULTIPOLYGON Z (((106.9709 ...
## 5 Ciamis MULTIPOLYGON Z (((108.358 -...
## 6 Cianjur MULTIPOLYGON Z (((107.2302 ...
## 7 Cirebon MULTIPOLYGON Z (((108.5391 ...
## 8 Garut MULTIPOLYGON (((107.8994 -7...
## 9 Indramayu MULTIPOLYGON (((108.1988 -6...
## 10 Karawang MULTIPOLYGON Z (((107.0983 ...
## 11 Kota Bandung MULTIPOLYGON Z (((107.5979 ...
## 12 Kota Banjar MULTIPOLYGON Z (((108.5758 ...
## 13 Kota Bekasi MULTIPOLYGON Z (((107.0052 ...
## 14 Kota Bogor MULTIPOLYGON Z (((106.783 -...
## 15 Kota Cimahi MULTIPOLYGON Z (((107.5479 ...
## 16 Kota Cirebon MULTIPOLYGON Z (((108.5621 ...
## 17 Kota Depok MULTIPOLYGON (((106.8149 -6...
## 18 Kota Sukabumi MULTIPOLYGON Z (((106.9151 ...
## 19 Kota Tasikmalaya MULTIPOLYGON Z (((108.194 -...
## 20 Kuningan MULTIPOLYGON Z (((108.4456 ...
## 21 Majalengka MULTIPOLYGON Z (((108.1225 ...
## 22 Pangandaran MULTIPOLYGON (((108.4971 -7...
## 23 Purwakarta MULTIPOLYGON Z (((107.5028 ...
## 24 Subang MULTIPOLYGON Z (((107.8876 ...
## 25 Sukabumi MULTIPOLYGON (((106.581 -7....
## 26 Sumedang MULTIPOLYGON Z (((107.8566 ...
## 27 Tasikmalaya MULTIPOLYGON (((108.3231 -7...Mengecek Apakah Urutan Sudah Sama
urutan_excel <- unique(data_panel$Kab_Kota)
urutan_peta <- peta_jabar$Kab_Kota
# Cek Logika: Apakah Urutan Excel == Urutan Peta?
cek_validasi <- all(urutan_excel == urutan_peta)
if(cek_validasi == TRUE) {
cat("\nUrutan Peta dan Excel sudah SAMA PERSIS.\n")
} else {
cat("\nUrutan masih berantakan.\n")
cat("Cek 10 urutan pertama Excel:\n")
print(head(urutan_excel, 10))
cat("Cek 10 urutan pertama Peta:\n")
print(head(urutan_peta, 10))
}##
## Urutan Peta dan Excel sudah SAMA PERSIS.3. Eksplorasi Data
Deskriptif
cat("=== RINGKASAN DATA SEBELUM TRANSFORMASI ===\n")## === RINGKASAN DATA SEBELUM TRANSFORMASI ===summary(data_panel)## id Kab_Kota Tahun TPT
## Min. : 1.00 Length:162 Min. :2019 Min. : 1.520
## 1st Qu.: 41.25 Class :character 1st Qu.:2020 1st Qu.: 6.567
## Median : 81.50 Mode :character Median :2022 Median : 8.155
## Mean : 81.50 Mean :2022 Mean : 8.124
## 3rd Qu.:121.75 3rd Qu.:2023 3rd Qu.: 9.717
## Max. :162.00 Max. :2024 Max. :14.290
## TPAK IPM PPM UMK
## Min. :55.74 Min. :65.36 Min. : 2.070 Min. :1688218
## 1st Qu.:63.76 1st Qu.:69.06 1st Qu.: 6.652 1st Qu.:2202184
## Median :65.62 Median :71.60 Median : 8.410 Median :2893229
## Mean :66.12 Mean :72.53 Mean : 8.274 Mean :3079173
## 3rd Qu.:68.17 3rd Qu.:75.14 3rd Qu.:10.355 3rd Qu.:3865110
## Max. :80.15 Max. :83.29 Max. :13.130 Max. :5343430
## AngKer AMA GR LPP
## Min. :36.74 Min. :92.34 Min. :0.2840 Min. :0.2000
## 1st Qu.:41.48 1st Qu.:98.17 1st Qu.:0.3420 1st Qu.:0.9125
## Median :44.24 Median :99.18 Median :0.3630 Median :1.1750
## Mean :44.03 Mean :98.59 Mean :0.3714 Mean :1.1988
## 3rd Qu.:45.87 3rd Qu.:99.55 3rd Qu.:0.3997 3rd Qu.:1.4300
## Max. :57.59 Max. :99.95 Max. :0.4890 Max. :3.9500
## PDRB JumlahSMP JumlahSMA JumlahFasKes
## Min. : 3221 Min. : 25.0 Min. : 4.00 Min. : 212
## 1st Qu.: 22450 1st Qu.:113.2 1st Qu.: 28.00 1st Qu.:1037
## Median : 33338 Median :186.0 Median : 53.00 Median :1742
## Mean : 59150 Mean :216.5 Mean : 64.41 Mean :2009
## 3rd Qu.: 59459 3rd Qu.:302.8 3rd Qu.: 96.75 3rd Qu.:2662
## Max. :293665 Max. :768.0 Max. :227.00 Max. :5284aa
ggplot(data = data_panel, aes(x = Tahun, y = IPM, group = Kab_Kota, color = Kab_Kota)) +
geom_line() +
geom_point() +
labs(title = "IPM dari tahun ke tahun",
x = "Year",
y = "IPM") +
theme_minimal()
Hubungan IPM dengan peubah penjelas
vars <- c("IPM", "TPT", "TPAK", "PPM", "UMK", "AngKer", "AMA", "GR", "LPP", "PDRB", "JumlahSMP", "JumlahSMA", "JumlahFasKes")
data_subset <- data_panel[, vars]
# Membuat data dalam format long
data_long <- gather(data_subset, key = "Variable", value = "Value", -IPM)
# Plotting menggunakan facet_wrap untuk 2 baris dan 3 kolom
ggplot(data_long, aes(x = IPM, y = Value)) +
geom_point(color = "goldenrod3") + # Mengatur warna titik
facet_wrap(~ Variable, scales = "free", nrow = 4, ncol = 3) +
labs(title = "Hubungan antara IPM dan peubah penjelas",
x = "Dependent Variables",
y = "Values") +
theme_light()
Tes normalitas
pairs.panels(data_panel[, c("IPM", "TPT", "TPAK", "PPM", "UMK", "AngKer", "GR", "PDRB", "JumlahSMA", "JumlahFasKes")],
method = "pearson", hist.col = "darkgoldenrod", col = "orange", density=TRUE)
vars_of_interest <- c("IPM", "TPT", "TPAK", "PPM", "UMK", "AngKer", "GR", "PDRB", "JumlahSMA", "JumlahFasKes")
# Create histograms for each variable
histograms <- lapply(vars_of_interest, function(var) {
ggplot(data_panel, aes_string(x = var)) +
geom_histogram(fill = "darkgoldenrod", color = "orange") +
labs(title = var) +
theme_minimal()
})## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.gridExtra::grid.arrange(grobs = histograms, ncol = 3)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ks_test <- ks.test(data_panel$IPM, "pnorm", mean(data_panel$IPM), sd(data_panel$IPM))## Warning in ks.test.default(data_panel$IPM, "pnorm", mean(data_panel$IPM), :
## ties should not be present for the Kolmogorov-Smirnov testks_test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: data_panel$IPM
## D = 0.10443, p-value = 0.05841
## alternative hypothesis: two-sidedTak tolak Ho, IPM mengikuti sebaran normal
# Create QQ plot
qqplot <- ggqqplot(data_panel$IPM, color = "goldenrod2") +
theme_light()
print(qqplot)
# Create histogram for variable Y
histogram_IPM <- ggplot(data_panel, aes(x = IPM)) +
geom_histogram(fill = "darkgoldenrod", color = "darkgoldenrod") +
labs(title = "Histogram of IPM") +
theme_light()
print(histogram_IPM)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Distribusi PDRB dan UMK
p1 <- ggplot(data_panel, aes(x = PDRB)) +
geom_histogram(fill = "red", bins = 30) +
labs(title = "Distribusi PDRB (Miliar Rp)", subtitle = "Miring ke Kanan (Skewed)")
p2 <- ggplot(data_panel, aes(x = UMK)) +
geom_histogram(fill = "blue", bins = 30) +
labs(title = "Distribusi UMK (Juta Rp)", subtitle = "Miring ke Kanan (Skewed)")
grid.arrange(p1, p2, ncol = 2)
Distribusi AMA, LPP
p10 <- ggplot(data_panel, aes(x = AMA)) +
geom_histogram(fill = "red", bins = 30) +
labs(title = "Angka Melek Aksara", subtitle = "-")
p11 <- ggplot(data_panel, aes(x = LPP)) +
geom_histogram(fill = "blue", bins = 30) +
labs(title = "Laju Pertumbuhan Penduduk", subtitle = "-")
p12 <- ggplot(data_panel, aes(x = GR)) +
geom_histogram(fill = "green", bins = 30) +
labs(title = "Gini Ratio", subtitle = "-")
p13 <- ggplot(data_panel, aes(x = TPAK)) +
geom_histogram(fill = "yellow", bins = 30) +
labs(title = "Tingkat Partisipasi Angkatan Kerja", subtitle = "-")
grid.arrange(p10, p11, p12, p13, ncol = 2)
Distribusi TPT, PPM, AngKer
p14 <- ggplot(data_panel, aes(x = TPT)) +
geom_histogram(fill = "red", bins = 30) +
labs(title = "Tingkat Pengangguran Terbuka", subtitle = "-")
p15 <- ggplot(data_panel, aes(x = PPM)) +
geom_histogram(fill = "blue", bins = 30) +
labs(title = "Persentase Penduduk Miskin", subtitle = "-")
p16 <- ggplot(data_panel, aes(x = AngKer)) +
geom_histogram(fill = "green", bins = 30) +
labs(title = "Angka Ketergantungan", subtitle = "-")
p17 <- ggplot(data_panel, aes(x = IPM)) +
geom_histogram(fill = "yellow", bins = 30) +
labs(title = "IPM", subtitle = "-")
grid.arrange(p14, p15, p16, p17, ncol = 2)
Distribusi JumlahSMP, JumlahSMA, JumlahFasKes
p3 <- ggplot(data_panel, aes(x = JumlahSMP)) +
geom_histogram(fill = "red", bins = 30) +
labs(title = "Jumlah SMP (Unit)", subtitle = "-")
p4 <- ggplot(data_panel, aes(x = JumlahSMA)) +
geom_histogram(fill = "blue", bins = 30) +
labs(title = "Jumlah SMA (Unit)", subtitle = "-")
p5 <- ggplot(data_panel, aes(x = JumlahFasKes)) +
geom_histogram(fill = "green", bins = 30) +
labs(title = "Jumlah Fasilitas Kesehatan (Unit)", subtitle = "-")
grid.arrange(p3, p4, p5, ncol = 2)
4. Transformasi Data (Log Natural)
Transformasi Logaritma Natural pada Peubah PDRB dan UMK
data_panel_tr <- data_panel %>%
mutate(
Ln_PDRB = log(PDRB),
Ln_UMK = log(UMK),
)Cek hasil transformasi
p6 <- ggplot(data_panel_tr, aes(x = Ln_PDRB)) +
geom_histogram(fill = "green", bins = 30) +
labs(title = "Distribusi Ln_PDRB", subtitle = "Lebih Normal (Bagus untuk Regresi)")
print(p6)
p7 <- ggplot(data_panel_tr, aes(x = Ln_UMK)) +
geom_histogram(fill = "green", bins = 30) +
labs(title = "Distribusi Ln_UMK", subtitle = "Lebih Normal (Bagus untuk Regresi)")
print(p7)
RLB
cat("REGRESI LINIER BERGANDA (MODEL PENUH)\n")## REGRESI LINIER BERGANDA (MODEL PENUH)model_penuh <- lm(IPM ~ PPM + TPT + TPAK + AngKer + AMA + GR + LPP + Ln_PDRB + Ln_UMK + JumlahSMP + JumlahSMA + JumlahFasKes,
data = data_panel_tr)
summary(model_penuh)##
## Call:
## lm(formula = IPM ~ PPM + TPT + TPAK + AngKer + AMA + GR + LPP +
## Ln_PDRB + Ln_UMK + JumlahSMP + JumlahSMA + JumlahFasKes,
## data = data_panel_tr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2526 -1.2069 -0.2362 1.2088 4.1438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.7352116 14.7699426 4.654 7.15e-06 ***
## PPM -0.5475529 0.0765912 -7.149 3.64e-11 ***
## TPT -0.2826642 0.0820085 -3.447 0.000737 ***
## TPAK -0.1553391 0.0544253 -2.854 0.004930 **
## AngKer -0.3624693 0.0645180 -5.618 9.21e-08 ***
## AMA 0.0117406 0.1233406 0.095 0.924293
## GR 15.7089848 4.6164496 3.403 0.000857 ***
## LPP -0.4064039 0.2986723 -1.361 0.175663
## Ln_PDRB -0.1395488 0.3566309 -0.391 0.696136
## Ln_UMK 2.2579813 0.9052871 2.494 0.013716 *
## JumlahSMP -0.0213812 0.0036645 -5.835 3.24e-08 ***
## JumlahSMA 0.0899061 0.0123229 7.296 1.64e-11 ***
## JumlahFasKes -0.0014731 0.0003145 -4.684 6.30e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.643 on 149 degrees of freedom
## Multiple R-squared: 0.8863, Adjusted R-squared: 0.8771
## F-statistic: 96.77 on 12 and 149 DF, p-value: < 2.2e-16model_stepwise <- step(model_penuh, direction = "both", trace = 0)cat("HASIL MODEL TERBAIK DARI STEPWISE\n")## HASIL MODEL TERBAIK DARI STEPWISEsummary(model_stepwise)##
## Call:
## lm(formula = IPM ~ PPM + TPT + TPAK + AngKer + GR + Ln_UMK +
## JumlahSMP + JumlahSMA + JumlahFasKes, data = data_panel_tr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3400 -1.1954 -0.2191 1.3534 4.0897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.785627 11.783030 6.007 1.34e-08 ***
## PPM -0.557094 0.071758 -7.764 1.13e-12 ***
## TPT -0.264722 0.080745 -3.278 0.001294 **
## TPAK -0.133026 0.051497 -2.583 0.010732 *
## AngKer -0.340387 0.056402 -6.035 1.16e-08 ***
## GR 16.733456 4.427845 3.779 0.000225 ***
## Ln_UMK 1.880423 0.671848 2.799 0.005793 **
## JumlahSMP -0.021866 0.003276 -6.675 4.36e-10 ***
## JumlahSMA 0.090189 0.010698 8.431 2.46e-14 ***
## JumlahFasKes -0.001490 0.000285 -5.229 5.54e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.638 on 152 degrees of freedom
## Multiple R-squared: 0.8847, Adjusted R-squared: 0.8779
## F-statistic: 129.6 on 9 and 152 DF, p-value: < 2.2e-16vif(model_penuh)## PPM TPT TPAK AngKer AMA GR
## 2.619294 2.324597 2.427857 2.398060 2.277196 2.339865
## LPP Ln_PDRB Ln_UMK JumlahSMP JumlahSMA JumlahFasKes
## 1.331222 7.267250 5.204990 20.038963 21.267133 9.708547vif(model_stepwise)## PPM TPT TPAK AngKer GR Ln_UMK
## 2.314241 2.268328 2.187909 1.844697 2.166712 2.885563
## JumlahSMP JumlahSMA JumlahFasKes
## 16.120528 16.132680 8.022788formtanpasmp <- IPM ~ PPM + TPT + TPAK + AngKer + GR + Ln_UMK + JumlahSMA + JumlahFasKes
tanpasmp <- lm(IPM ~ PPM + TPT + TPAK + AngKer + GR + Ln_UMK + JumlahSMA + JumlahFasKes,
data = data_panel_tr)
vif(tanpasmp)## PPM TPT TPAK AngKer GR Ln_UMK
## 2.231444 2.259480 2.181118 1.838582 1.948282 2.768955
## JumlahSMA JumlahFasKes
## 7.115610 7.138493RLB hanya 2024
cat("SELEKSI STEPWISE: CROSS-SECTION (TAHUN 2024)\n")## SELEKSI STEPWISE: CROSS-SECTION (TAHUN 2024)# Filter data hanya untuk tahun 2024
data_2024 <- data_panel_tr %>% filter(Tahun == 2024)
model_2024 <- lm(IPM ~ PPM + TPT + TPAK + AngKer + AMA + GR + LPP + Ln_PDRB + Ln_UMK + JumlahSMP + JumlahSMA + JumlahFasKes,
data = data_2024)
stepwise_2024 <- step(model_2024, direction = "both", trace = 0)
summary(stepwise_2024)##
## Call:
## lm(formula = IPM ~ PPM + TPAK + AngKer + JumlahSMP + JumlahSMA +
## JumlahFasKes, data = data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0185 -0.6952 -0.2772 0.5627 2.6776
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.098e+02 6.801e+00 16.142 6.17e-13 ***
## PPM -5.535e-01 1.442e-01 -3.837 0.001028 **
## TPAK -2.060e-01 8.945e-02 -2.303 0.032177 *
## AngKer -3.410e-01 1.478e-01 -2.307 0.031853 *
## JumlahSMP -2.425e-02 6.093e-03 -3.980 0.000738 ***
## JumlahSMA 1.017e-01 1.823e-02 5.577 1.85e-05 ***
## JumlahFasKes -1.828e-03 5.232e-04 -3.493 0.002290 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.346 on 20 degrees of freedom
## Multiple R-squared: 0.9262, Adjusted R-squared: 0.9041
## F-statistic: 41.83 on 6 and 20 DF, p-value: 2.753e-10vif(stepwise_2024)## PPM TPAK AngKer JumlahSMP JumlahSMA JumlahFasKes
## 2.095204 1.360960 1.802832 14.727173 13.008643 6.941608formtanpasmp_2024 <- IPM ~ PPM + TPAK + AngKer + JumlahSMA + JumlahFasKes
tanpasmp_2024 <- lm(IPM ~ PPM + TPAK + AngKer + JumlahSMA + JumlahFasKes,
data = data_panel_tr)
vif(tanpasmp_2024)## PPM TPAK AngKer JumlahSMA JumlahFasKes
## 1.709011 1.039010 1.310925 5.592638 5.071455Model RLB dari semua tahun, peubah penjelas yang tepilih 8 peubah penjelas yaitu PPM, TPT, TPAK, AngKer, GR, Ln_UMK, JumlahSMA, JumlahFasKes
Sedangkan jika RLB hanya tahun terbatu yaitu 2024, peubah penjelas yang terpilih 5 peubah penjelas yaitu PPM, TPAK, AngKer, JumlahSMA, JumlahFasKes
5. Ekplorasi Spasial
# Menggabungkan data dengan Peta untuk tahun terbaru (2024)
tahun_terbaru <- max(data_panel$Tahun)
data_map <- data_panel_tr %>% filter(Tahun == tahun_terbaru)
# Join dengan Shapefile
peta_visual <- left_join(peta_jabar, data_map, by = "Kab_Kota")
# Memperbaiki poligon yang rusak/tidak tertutup
peta_visual <- st_make_valid(peta_visual)
# Plot Peta IPM (Y)
plot_ipm <- ggplot(peta_visual) +
geom_sf(aes(fill = IPM)) +
scale_fill_viridis_c(option = "plasma") +
labs(title = paste("Peta Sebaran IPM Jawa Barat", tahun_terbaru)) +
theme_minimal()
# Plot Peta Ln_PDRB (X)
plot_pdrb <- ggplot(peta_visual) +
geom_sf(aes(fill = Ln_PDRB)) +
scale_fill_viridis_c(option = "viridis") +
labs(title = paste("Peta Sebaran Ekonomi (Ln PDRB)", tahun_terbaru)) +
theme_minimal()
# Plot Peta PPM (X)
plot_ppm <- ggplot(peta_visual) +
geom_sf(aes(fill = PPM)) +
scale_fill_viridis_c(option = "plasma") +
labs(title = paste("Peta Sebaran PPM Jawa Barat", tahun_terbaru)) +
theme_minimal()
# Plot Peta TPT (X)
plot_tpt <- ggplot(peta_visual) +
geom_sf(aes(fill = TPT)) +
scale_fill_viridis_c(option = "viridis") +
labs(title = paste("Peta Sebaran TPT Jawa Barat", tahun_terbaru)) +
theme_minimal()
grid.arrange(plot_ipm, plot_pdrb, plot_ppm, plot_tpt, nrow = 2, ncol = 2)
# TPAK
plot_tpak <- ggplot(peta_visual) +
geom_sf(aes(fill = TPAK)) +
scale_fill_viridis_c(option = "plasma") +
labs(title = paste("Peta Sebaran TPAK", tahun_terbaru)) +
theme_minimal()
# UMK
plot_umk <- ggplot(peta_visual) +
geom_sf(aes(fill = Ln_UMK)) +
scale_fill_viridis_c(option = "viridis") +
labs(title = paste("Peta Sebaran (Ln UMK)", tahun_terbaru)) +
theme_minimal()
# Jumlah SMA
plot_sma <- ggplot(peta_visual) +
geom_sf(aes(fill = JumlahSMA)) +
scale_fill_viridis_c(option = "plasma") +
labs(title = paste("Peta Sebaran Jumlah SMA", tahun_terbaru)) +
theme_minimal()
# FasKes
plot_faskes <- ggplot(peta_visual) +
geom_sf(aes(fill = JumlahFasKes)) +
scale_fill_viridis_c(option = "viridis") +
labs(title = paste("Peta Sebaran Jumlah FasKes", tahun_terbaru)) +
theme_minimal()
grid.arrange(plot_tpak, plot_umk, plot_sma, plot_faskes, nrow = 2, ncol = 2)
6. Membuat Matriks Pembobot Spasial
Butuh sisaan (residual) dari model OLS biasa untuk diuji spasialnya, coba pakai model tanpasmp dan tanpasmp_2024
Ambil koordinat untuk perhitungan jarak
coords <- st_coordinates(st_centroid(peta_jabar))## Warning: st_centroid assumes attributes are constant over geometries## st_as_s2(): dropping Z and/or M coordinateTAHAP 1: MENCARI ‘K’ OPTIMAL (K-NEAREST NEIGHBOR)
Coba k dari 1 sampai 15 hingga dapat Moran’s I-nya paling tinggi
# LOOPING MENCARI K OPTIMAL (DENGAN DATA CROSS-SECTION)
#data_map berisi hanya tahun 2024
model_ols_2024 <- lm(formtanpasmp, data = data_map)
k_range <- 1:10
moran_k <- numeric(length(k_range))
cat("Mencari K optimal menggunakan data tahun", tahun_terbaru, "...\n")## Mencari K optimal menggunakan data tahun 2024 ...for (i in k_range) {
# a. Buat tetangga k
knn <- knearneigh(coords, k = i)
nb_k <- knn2nb(knn)
# b. Buat matriks (Row Standardized)
listw_k <- nb2listw(nb_k, style = "W")
# c. Uji Moran pada model Cross-Section
moran_test <- lm.morantest(model_ols_2024, listw = listw_k)
moran_k[i] <- moran_test$estimate[1]
}## Warning in knn2nb(knn): neighbour object has 8 sub-graphs## Warning in knearneigh(coords, k = i): k greater than one-third of the number of
## data points# Visualisasi Hasil
df_k <- data.frame(k = k_range, Moran_I = moran_k)
ggplot(df_k, aes(x=k, y=Moran_I)) +
geom_line(color="blue", linewidth=1) +
geom_point(color="red", size=3) +
labs(title = paste("Optimasi Matriks k-NN (Basis Data:", tahun_terbaru, ")"),
x = "Jumlah Tetangga (k)",
y = "Moran's I Index") +
scale_x_continuous(breaks = k_range) +
theme_minimal()
# LOOPING MENCARI K OPTIMAL (DENGAN DATA CROSS-SECTION)
#data_map berisi hanya tahun 2024
model_ols_2024_2 <- lm(formtanpasmp_2024, data = data_map)
k_range <- 1:10
moran_k <- numeric(length(k_range))
cat("Mencari K optimal menggunakan data tahun", tahun_terbaru, "...\n")## Mencari K optimal menggunakan data tahun 2024 ...for (i in k_range) {
# a. Buat tetangga k
knn <- knearneigh(coords, k = i)
nb_k <- knn2nb(knn)
# b. Buat matriks (Row Standardized)
listw_k <- nb2listw(nb_k, style = "W")
# c. Uji Moran pada model Cross-Section (27 baris vs 27x27 matriks sehingga pas)
moran_test <- lm.morantest(model_ols_2024_2, listw = listw_k)
moran_k[i] <- moran_test$estimate[1]
}## Warning in knn2nb(knn): neighbour object has 8 sub-graphs## Warning in knearneigh(coords, k = i): k greater than one-third of the number of
## data points# Visualisasi Hasil
df_k <- data.frame(k = k_range, Moran_I = moran_k)
ggplot(df_k, aes(x=k, y=Moran_I)) +
geom_line(color="blue", linewidth=1) +
geom_point(color="red", size=3) +
labs(title = paste("Optimasi Matriks k-NN (Basis Data:", tahun_terbaru, ")"),
x = "Jumlah Tetangga (k)",
y = "Moran's I Index") +
scale_x_continuous(breaks = k_range) +
theme_minimal()
# BUAT MATRIKS TETANGGA (k=2)
knn_final <- knearneigh(coords, k = 2)
nb_final <- knn2nb(knn_final)
# BUAT MATRIKS PEMBOBOT
# Ini adalah matriks final yang akan dipakai di spml()
listw_k2 <- nb2listw(nb_final, style = "W")
# CEK STRUKTUR TETANGGA (VISUALISASI)
# Melihat siapa tetangga siapa (2 terdekat)
plot(st_geometry(peta_jabar), border="grey")
plot(nb_final, coords, add=TRUE, col="red", length=0.08)
title(main="Peta Konektivitas 2 Tetangga Terdekat (k=2)")
#COBA COBA
# BUAT MATRIKS TETANGGA (k=6)
# Menggunakan koordinat centroid yang sudah dibuat sebelumnya
knn_final2 <- knearneigh(coords, k = 6)
nb_final2 <- knn2nb(knn_final2)
# BUAT MATRIKS PEMBOBOT
# Ini adalah matriks final yang akan dipakai di spml()
listw_k6 <- nb2listw(nb_final2, style = "W")
# CEK STRUKTUR TETANGGA (VISUALISASI)
# Melihat siapa tetangga siapa (6 terdekat)
plot(st_geometry(peta_jabar), border="grey")
plot(nb_final2, coords, add=TRUE, col="red", length=0.08)
title(main="Peta Konektivitas 6 Tetangga Terdekat (k=6)")
TAHAP 2: MENCARI ‘ALPHA’ OPTIMAL (INVERS JARAK)
Rumus: 1 / Jarak^alpha
Biasanya alpha antara 0.5 sampai 3.0 (Hukum Gravitasi Newton pakai alpha=2)
Menggunakan struktur tetangga dari K=2 tadi sebagai dasar
knn_final <- knearneigh(coords, k = 2)
nb_final <- knn2nb(knn_final)
dists <- nbdists(nb_final, coords)
alpha_range <- seq(0.5, 4.0, by = 0.25)
moran_idw <- numeric(length(alpha_range))
for (i in 1:length(alpha_range)) {
a <- alpha_range[i]
# Hitung bobot invers jarak pangkat alpha
weights <- lapply(dists, function(x) 1 / (x^a))
# Buat listw
listw_idw <- nb2listw(nb_final, glist = weights, style = "W")
# Uji Moran
moran_idw[i] <- lm.morantest(model_ols_2024, listw = listw_idw)$estimate[1]
}
# Plot Hasil Alpha Optimal
df_alpha <- data.frame(alpha = alpha_range, Moran_I = moran_idw)
ggplot(df_alpha, aes(x=alpha, y=Moran_I)) +
geom_line(color="green", size=1) +
geom_point(color="darkgreen", size=3) +
labs(title="Optimasi Parameter Invers Jarak (1/d^alpha)", x="Pangkat Alpha", y="Moran's I Index") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
PILIH ALPHA 0.5 atau 1
# HITUNG ULANG BOBOT DENGAN alpha = 0.5
weights_idw_final <- lapply(dists, function(x) 1 / (x^0.5))
# BUAT MATRIKS FINAL
listw_idw_opt <- nb2listw(nb_final, glist = weights_idw_final, style = "W")PAKAI YANG listw_idw_opt saat membandingkan MORAN I
TAHAP 3: MENCARI ‘BETA’ OPTIMAL (EKSPONENSIAL JARAK)
Rumus: exp(-beta * Jarak)
Beta biasanya kecil (misal 0.01 sampai 2) tergantung satuan jarak (derajat/meter)
beta_range <- seq(0.1, 5.0, by = 0.2)
moran_exp <- numeric(length(beta_range))
for (i in 1:length(beta_range)) {
b <- beta_range[i]
# Hitung bobot eksponensial
weights_exp <- lapply(dists, function(x) exp(-b * x))
# Buat listw
listw_exp <- nb2listw(nb_final, glist = weights_exp, style = "W")
# Uji Moran
moran_exp[i] <- lm.morantest(model_ols_2024, listw = listw_exp)$estimate[1]
}
# Plot Hasil Beta Optimal
df_beta <- data.frame(beta = beta_range, Moran_I = moran_exp)
ggplot(df_beta, aes(x=beta, y=Moran_I)) +
geom_line(color="purple", size=1) +
geom_point(color="black", size=3) +
labs(title="Optimasi Parameter Eksponensial (exp(-beta*d))", x="Beta", y="Moran's I Index") +
theme_minimal()
Terpilih 0.1
# HITUNG ULANG BOBOT DENGAN beta_terbaik
weights_exp_final <- lapply(dists, function(x) exp(-0.1 * x))
# BUAT MATRIKS FINAL
listw_exp_opt <- nb2listw(nb_final, glist = weights_exp_final, style = "W")Menggunakan Queen Contiguity
# MEMBUAT MATRIKS QUEEN CONTIGUITY
peta_jabar <- st_make_valid(peta_jabar)
# 1. Definisikan Tetangga (Berdasarkan Poligon yang Menempel)
nb_queen <- poly2nb(peta_jabar, queen = TRUE)
# 2. Cek Ringkasan Konektivitas
cat("RINGKASAN STRUKTUR TETANGGA QUEEN\n")## RINGKASAN STRUKTUR TETANGGA QUEENsummary(nb_queen)## Neighbour list object:
## Number of regions: 27
## Number of nonzero links: 106
## Percentage nonzero weights: 14.54047
## Average number of links: 3.925926
## Link number distribution:
##
## 1 2 3 4 5 6 7 8
## 4 3 6 4 1 7 1 1
## 4 least connected regions:
## 12 14 16 18 with 1 link
## 1 most connected region:
## 4 with 8 links# 3. Buat Matriks Pembobot (Row Standardized)
# zero.policy = TRUE penting jika ada wilayah kepulauan (tidak punya tetangga)
listw_queen <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)
# VISUALISASI KONEKSI QUEEN (UNTUK MEMBANDINGKAN DENGAN K-NN)
# Plot Peta Dasar
plot(st_geometry(peta_jabar), border = "grey", main = "Konektivitas Queen Contiguity")
# Plot Garis Hubungan Antar Tetangga
plot(nb_queen, coords, add = TRUE, col = "blue", length = 0.08)
# Legend
legend("bottomright", legend = "Tetangga (Menempel)", col = "blue", lty = 1, cex = 0.8)
“Average number of links” memiliki rata-rata sekitar 3-5
Penetuan matriks pembobot terbaik
# Ubah data ke format Panel Data Frame
pdata <- pdata.frame(data_panel_tr, index = c("Kab_Kota", "Tahun"))
# Memasukkan semua listw yang sudah dibuat
daftar_w <- list(
"Queen Contiguity" = listw_queen,
"k-NN (k=2)" = listw_k2, # k optimal 2
"k-NN (k=6)" = listw_k6, # k = 4 coba
"Invers Jarak" = listw_idw_opt, # alpha optimal 0.5
"Eksponensial" = listw_exp_opt # beta optimal 0.1
)
# Menyimpan Skor
skor_matriks <- data.frame(
Matriks = character(),
LM_Stat = numeric(),
P_Value = numeric(),
stringsAsFactors = FALSE
)
# Menggunakan uji LM (Lagrange Multiplier) standar untuk membandingkan
for (nama in names(daftar_w)) {
bobot <- daftar_w[[nama]]
# Kita ambil LM Error (SEM) atau LM Lag (SAR) sebagai indikator
# Biasanya jika matriksnya bagus, kedua nilai ini akan tinggi.
# Di sini kita pakai RLM (Robust LM) gabungan atau LM Error standar sebagai patokan
uji <- slmtest(formtanpasmp, data = pdata, listw = bobot, model = "within", test = "lme")
# Simpan hasilnya
skor_matriks <- rbind(skor_matriks, data.frame(
Matriks = nama,
LM_Stat = as.numeric(uji$statistic),
P_Value = as.numeric(uji$p.value)
))
}
# Urutkan dari LM Statistic Terbesar (Makin besar = Makin Kuat)
skor_final <- skor_matriks %>%
arrange(desc(LM_Stat))
print(skor_final)## Matriks LM_Stat P_Value
## 1 Queen Contiguity 43.72367 3.781698e-11
## 2 k-NN (k=2) 27.54711 1.533141e-07
## 3 Eksponensial 27.47360 1.592535e-07
## 4 Invers Jarak 24.47119 7.542934e-07
## 5 k-NN (k=6) 20.31433 6.570651e-06cat("Matriks pembobot terbaik:", skor_final$Matriks[1], "\n")## Matriks pembobot terbaik: Queen Contiguitycat("Nilai LM Statistic:", round(skor_final$LM_Stat[1], 2), "\n")## Nilai LM Statistic: 43.72cat("Gunakan matriks", skor_final$Matriks[1], "untuk analisis selanjutnya (SAR/SEM).")## Gunakan matriks Queen Contiguity untuk analisis selanjutnya (SAR/SEM).# Simpan yang tertinggi ke final
nama_pemenang <- skor_final$Matriks[1]
listw_terpilih <- daftar_w[[nama_pemenang]]Menguji autokorelasi pada peubah respon menggunakan indeks moran
# UJI MORAN'S I GLOBAL PADA 'IPM'
moran_ipm <- moran.test(data_map$IPM, listw = listw_terpilih)
print(moran_ipm)##
## Moran I test under randomisation
##
## data: data_map$IPM
## weights: listw_terpilih
##
## Moran I statistic standard deviate = 2.2533, p-value = 0.01212
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.27667928 -0.03846154 0.01956095# Interpretasi:
# Jika p-value < 0.05, berarti IPM memang terklaster secara spasial.
# (Wilayah IPM tinggi bertetangga dengan IPM tinggi).
# MEMBUAT MORAN SCATTERPLOT (VISUALISASI)
# Grafik (Kuadran I, II, III, IV)
moran.plot(data_map$IPM, listw = listw_terpilih,
labels = as.character(data_map$Kab_Kota), # Label nama kota
pch = 19, col = "blue",
main = paste("Moran Scatterplot IPM Tahun", tahun_terbaru),
xlab = "IPM (Standar)",
ylab = "Spatially Lagged IPM (Rata-rata Tetangga)")
# Garis referensi
abline(h=0, v=0, lty=2)
Menguji autokorelasi pada peubah respon yang di mean kan semua tahunnya menggunakan indeks moran
# Ambil rata-rata IPM per kota selama periode penelitian
data_rata2 <- data_panel %>%
group_by(Kab_Kota) %>%
summarise(IPM_Mean = mean(IPM))
# Uji Moran pada rata-rata tersebut
moran.test(data_rata2$IPM_Mean, listw = listw_terpilih)##
## Moran I test under randomisation
##
## data: data_rata2$IPM_Mean
## weights: listw_terpilih
##
## Moran I statistic standard deviate = 2.4218, p-value = 0.007721
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.29983406 -0.03846154 0.01951222Menguji autokorelasi pada peubah respon setiap tahunnya menggunakan indeks moran
list_tahun <- unique(data_panel$Tahun)
cat("HASIL UJI MORAN'S I DARI TAHUN KE TAHUN\n")## HASIL UJI MORAN'S I DARI TAHUN KE TAHUNfor (thn in list_tahun) {
# 1. Filter data tahun tersebut
data_subset <- data_panel %>% filter(Tahun == thn)
# 2. Uji Moran
# Pastikan urutan data subset SAMA dengan urutan di listw (biasanya aman)
uji <- moran.test(data_subset$IPM, listw = listw_terpilih)
# 3. Tampilkan Hasil
cat("Tahun:", thn,
"| Moran's I:", round(uji$estimate[1], 4),
"| p-value:", format.pval(uji$p.value, digits=4), "\n")
}## Tahun: 2019 | Moran's I: 0.317 | p-value: 0.00547
## Tahun: 2020 | Moran's I: 0.3038 | p-value: 0.007114
## Tahun: 2021 | Moran's I: 0.3015 | p-value: 0.007455
## Tahun: 2022 | Moran's I: 0.299 | p-value: 0.00784
## Tahun: 2023 | Moran's I: 0.2968 | p-value: 0.008197
## Tahun: 2024 | Moran's I: 0.2767 | p-value: 0.012127. Menentukan Model Panel Terbaik (UJI CHOW, HAUSMAN, LM)
Pemilihan Model Regresi Data Panel
# A. Common Effect Model (CEM) / Pooled OLS
# Menganggap tidak ada perbedaan karakteristik antar wilayah/waktu.
model_cem <- plm(formtanpasmp, data = pdata, model = "pooling")
# B. Fixed Effect Model (FEM) / LSDV
# Menganggap setiap wilayah punya karakteristik unik (intercept beda-beda).
model_fem <- plm(formtanpasmp, data = pdata, model = "within")
# C. Random Effect Model (REM) / GLS
# Menganggap perbedaan antar wilayah sebagai error random.
model_rem <- plm(formtanpasmp, data = pdata, model = "random")
# UJI PEMILIHAN MODEL (CEK MANA YANG TERBAIK)
# A. Uji Chow (CEM vs FEM)
# H0: CEM lebih baik
# H1: FEM lebih baik
uji_chow <- pooltest(model_cem, model_fem)
print(uji_chow)##
## F statistic
##
## data: formtanpasmp
## F = 163.19, df1 = 26, df2 = 127, p-value < 2.2e-16
## alternative hypothesis: unstabilitycat("Jika p-value < 0.05, pilih FEM.\n")## Jika p-value < 0.05, pilih FEM.# B. Uji Hausman (FEM vs REM)
# H0: REM lebih baik (konsisten)
# H1: FEM lebih baik (konsisten)
uji_hausman <- phtest(model_fem, model_rem)
print(uji_hausman)##
## Hausman Test
##
## data: formtanpasmp
## chisq = 12.078, df = 8, p-value = 0.1477
## alternative hypothesis: one model is inconsistentcat("Jika p-value < 0.05, pilih FEM. Jika > 0.05, pilih REM.\n")## Jika p-value < 0.05, pilih FEM. Jika > 0.05, pilih REM.# C. Uji Lagrange Multiplier (LM) (CEM vs REM) - Opsional jika Hausman pilih REM
# H0: CEM lebih baik
# H1: REM lebih baik
uji_lm <- plmtest(model_cem, type="bp") # Breusch-Pagan
print(uji_lm)##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: formtanpasmp
## chisq = 208.85, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectssummary(model_rem)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = formtanpasmp, data = pdata, model = "random")
##
## Balanced Panel: n = 27, T = 6, N = 162
##
## Effects:
## var std.dev share
## idiosyncratic 0.1207 0.3474 0.035
## individual 3.3238 1.8231 0.965
## theta: 0.9224
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.01036 -0.29676 -0.02587 0.32807 1.07958
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -6.7960e+01 9.6699e+00 -7.0280 2.095e-12 ***
## PPM -3.6191e-01 7.4024e-02 -4.8891 1.013e-06 ***
## TPT -2.1015e-01 3.2070e-02 -6.5528 5.647e-11 ***
## TPAK 1.5744e-02 2.0137e-02 0.7818 0.43432
## AngKer -3.2121e-02 2.6214e-02 -1.2254 0.22043
## GR 1.1002e+00 1.6479e+00 0.6676 0.50437
## Ln_UMK 9.8703e+00 6.3517e-01 15.5397 < 2.2e-16 ***
## JumlahSMA -5.8895e-03 8.4624e-03 -0.6960 0.48646
## JumlahFasKes -6.9758e-04 2.9373e-04 -2.3749 0.01755 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 162.56
## Residual Sum of Squares: 26.644
## R-Squared: 0.8361
## Adj. R-Squared: 0.82753
## Chisq: 780.487 on 8 DF, p-value: < 2.22e-16# A. Common Effect Model (CEM) / Pooled OLS
# Menganggap tidak ada perbedaan karakteristik antar wilayah/waktu.
model_cem2 <- plm(formtanpasmp_2024, data = pdata, model = "pooling")
# B. Fixed Effect Model (FEM) / LSDV
# Menganggap setiap wilayah punya karakteristik unik (intercept beda-beda).
model_fem2 <- plm(formtanpasmp_2024, data = pdata, model = "within")
# C. Random Effect Model (REM) / GLS
# Menganggap perbedaan antar wilayah sebagai error random.
model_rem2 <- plm(formtanpasmp_2024, data = pdata, model = "random")
# UJI PEMILIHAN MODEL (CEK MANA YANG TERBAIK)
# A. Uji Chow (CEM vs FEM)
# H0: CEM lebih baik
# H1: FEM lebih baik
uji_chow2 <- pooltest(model_cem2, model_fem2)
print(uji_chow2)##
## F statistic
##
## data: formtanpasmp_2024
## F = 38.83, df1 = 26, df2 = 130, p-value < 2.2e-16
## alternative hypothesis: unstabilitycat("Jika p-value < 0.05, pilih FEM.\n")## Jika p-value < 0.05, pilih FEM.# B. Uji Hausman (FEM vs REM)
# H0: REM lebih baik (konsisten)
# H1: FEM lebih baik (konsisten)
uji_hausman2 <- phtest(model_fem2, model_rem2)
print(uji_hausman2)##
## Hausman Test
##
## data: formtanpasmp_2024
## chisq = 31.728, df = 5, p-value = 6.725e-06
## alternative hypothesis: one model is inconsistentcat("Jika p-value < 0.05, pilih FEM. Jika > 0.05, pilih REM.\n")## Jika p-value < 0.05, pilih FEM. Jika > 0.05, pilih REM.# C. Uji Lagrange Multiplier (LM) (CEM vs REM) - Opsional jika Hausman pilih REM
# H0: CEM lebih baik
# H1: REM lebih baik
uji_lm2 <- plmtest(model_cem2, type="bp") # Breusch-Pagan
print(uji_lm2)##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: formtanpasmp_2024
## chisq = 196.66, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectssummary(model_fem2)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = formtanpasmp_2024, data = pdata, model = "within")
##
## Balanced Panel: n = 27, T = 6, N = 162
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.481234 -0.454017 -0.052999 0.424701 2.093978
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## PPM -0.21318369 0.15231890 -1.3996 0.1640
## TPAK 0.18052827 0.03227515 5.5934 1.259e-07 ***
## AngKer -0.04225065 0.04973521 -0.8495 0.3972
## JumlahSMA 0.08746932 0.01728552 5.0603 1.397e-06 ***
## JumlahFasKes 0.00044484 0.00069276 0.6421 0.5219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 142.13
## Residual Sum of Squares: 76.757
## R-Squared: 0.45996
## Adj. R-Squared: 0.33119
## F-statistic: 22.1449 on 5 and 130 DF, p-value: 5.2631e-16Melihat pengaruh efek spesifik
# FEM model tanpasmp_2024
# 1. UJI LM INDIVIDU
# H0: Tidak ada efek wilayah
uji_lm_ind2 <- plmtest(model_fem2, effect = "individual", type = "bp")
print(uji_lm_ind2)##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: formtanpasmp_2024
## chisq = 196.66, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects# 2. UJI LM WAKTU
# H0: Tidak ada efek tahun
uji_lm_waktu2 <- plmtest(model_fem2, effect = "time", type = "bp")
print(uji_lm_waktu2)##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: formtanpasmp_2024
## chisq = 40.486, df = 1, p-value = 1.981e-10
## alternative hypothesis: significant effects# 3. UJI LM DUA ARAH
# H0: Tidak ada efek wilayah ATAU waktu
uji_lm_dua_arah2 <- plmtest(model_fem2, effect = "twoways", type = "bp")
print(uji_lm_dua_arah2)##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: formtanpasmp_2024
## chisq = 237.15, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effectsEFEK two ways
# REM model tanpasmp
# 1. UJI LM INDIVIDU
# H0: Tidak ada efek wilayah
uji_lm_ind3 <- plmtest(model_rem, effect = "individual", type = "bp")
print(uji_lm_ind3)##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: formtanpasmp
## chisq = 208.85, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects# 2. UJI LM WAKTU
# H0: Tidak ada efek tahun
uji_lm_waktu3 <- plmtest(model_rem, effect = "time", type = "bp")
print(uji_lm_waktu3)##
## Lagrange Multiplier Test - time effects (Breusch-Pagan)
##
## data: formtanpasmp
## chisq = 2.2189, df = 1, p-value = 0.1363
## alternative hypothesis: significant effects# 3. UJI LM DUA ARAH
# H0: Tidak ada efek wilayah ATAU waktu
uji_lm_dua_arah3 <- plmtest(model_rem, effect = "twoways", type = "bp")
print(uji_lm_dua_arah3)##
## Lagrange Multiplier Test - two-ways effects (Breusch-Pagan)
##
## data: formtanpasmp
## chisq = 211.07, df = 2, p-value < 2.2e-16
## alternative hypothesis: significant effectsEFEK INDIVIDU
Uji Lagrange Multiplier (LM Test) untuk Efek Spasial
Menentukan apakah pakai SAR (Lag) atau SEM (Error)
cat("UJI LAGRANGE MULTIPLIER (SPASIAL)\n")## UJI LAGRANGE MULTIPLIER (SPASIAL)lm_lag <- slmtest(formtanpasmp, data = pdata, listw = listw_terpilih, model = "random", test = "lml")
lm_err <- slmtest(formtanpasmp, data = pdata, listw = listw_terpilih, model = "random", test = "lme")
lm_lag##
## LM test for spatial lag dependence
##
## data: formula (random transformation)
## LM = 43.977, df = 1, p-value = 3.322e-11
## alternative hypothesis: spatial lag dependencelm_err##
## LM test for spatial error dependence
##
## data: formula (random transformation)
## LM = 31.373, df = 1, p-value = 2.129e-08
## alternative hypothesis: spatial error dependencecat("UJI ROBUST LAGRANGE MULTIPLIER\n")## UJI ROBUST LAGRANGE MULTIPLIER# Robust LM Lag (Untuk SAR)
rlm_lag <- slmtest(formtanpasmp, data = pdata, listw = listw_terpilih,
model = "random", test = "rlml")
# Robust LM Error (Untuk SEM)
rlm_err <- slmtest(formtanpasmp, data = pdata, listw = listw_terpilih,
model = "random", test = "rlme")
print(rlm_lag)##
## Locally robust LM test for spatial lag dependence sub spatial error
##
## data: formula (random transformation)
## LM = 12.604, df = 1, p-value = 0.0003849
## alternative hypothesis: spatial lag dependenceprint(rlm_err)##
## Locally robust LM test for spatial error dependence sub spatial lag
##
## data: formula (random transformation)
## LM = 1.3729e-05, df = 1, p-value = 0.997
## alternative hypothesis: spatial error dependencecat("UJI LAGRANGE MULTIPLIER (SPASIAL) model tanpasmp_2024\n")## UJI LAGRANGE MULTIPLIER (SPASIAL) model tanpasmp_2024lm_lag2 <- slmtest(formtanpasmp_2024, data = pdata, listw = listw_terpilih, model = "within", test = "lml")
lm_err2 <- slmtest(formtanpasmp_2024, data = pdata, listw = listw_terpilih, model = "within", test = "lme")
lm_lag2##
## LM test for spatial lag dependence
##
## data: formula (within transformation)
## LM = 133.14, df = 1, p-value < 2.2e-16
## alternative hypothesis: spatial lag dependencelm_err2##
## LM test for spatial error dependence
##
## data: formula (within transformation)
## LM = 75.253, df = 1, p-value < 2.2e-16
## alternative hypothesis: spatial error dependencecat("UJI ROBUST LAGRANGE MULTIPLIER model tanpasmp_2024\n")## UJI ROBUST LAGRANGE MULTIPLIER model tanpasmp_2024# Robust LM Lag (Untuk SAR)
rlm_lag2 <- slmtest(formtanpasmp_2024, data = pdata, listw = listw_terpilih,
model = "within", test = "rlml")
# Robust LM Error (Untuk SEM)
rlm_err2 <- slmtest(formtanpasmp_2024, data = pdata, listw = listw_terpilih,
model = "within", test = "rlme")
print(rlm_lag2)##
## Locally robust LM test for spatial lag dependence sub spatial error
##
## data: formula (within transformation)
## LM = 57.922, df = 1, p-value = 2.727e-14
## alternative hypothesis: spatial lag dependenceprint(rlm_err2)##
## Locally robust LM test for spatial error dependence sub spatial lag
##
## data: formula (within transformation)
## LM = 0.035708, df = 1, p-value = 0.8501
## alternative hypothesis: spatial error dependence8. ESTIMASI MODEL SPASIAL FINAL
Kasus 1 Random Effect (model=“random”) untuk tanpasmp
cat("ESTIMASI MODEL SAR-REM (Spatial Autoregressive)\n")## ESTIMASI MODEL SAR-REM (Spatial Autoregressive)model_sar_rem <- spml(formula = formtanpasmp,
data = pdata,
listw = listw_terpilih,
model = "random",
effect = "individual",
lag = TRUE, # SAR hidup
spatial.error = "none") # Error mati
summary(model_sar_rem)## ML panel with spatial lag, random effects
##
## Call:
## spreml(formula = formula, data = data, index = index, w = listw2mat(listw),
## w2 = listw2mat(listw2), lag = lag, errors = errors, cl = cl)
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 41.1 44.5 46.7 47.6 50.4 57.0
##
## Error variance parameters:
## Estimate Std. Error t-value Pr(>|t|)
## phi 240.710 80.822 2.9783 0.002899 **
##
## Spatial autoregressive coefficient:
## Estimate Std. Error t-value Pr(>|t|)
## lambda 0.666223 0.055821 11.935 < 2.2e-16 ***
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -2.0945e+01 5.7687e+00 -3.6309 0.0002825 ***
## PPM -5.6963e-02 4.9038e-02 -1.1616 0.2453883
## TPT -6.7490e-02 1.9300e-02 -3.4969 0.0004707 ***
## TPAK 1.0564e-02 1.1341e-02 0.9316 0.3515683
## AngKer 3.3220e-05 1.5445e-02 0.0022 0.9982839
## GR 1.1005e+00 9.3240e-01 1.1803 0.2378736
## Ln_UMK 3.0710e+00 3.8513e-01 7.9739 1.537e-15 ***
## JumlahSMA -1.9753e-03 5.7467e-03 -0.3437 0.7310543
## JumlahFasKes 1.1325e-04 2.0097e-04 0.5635 0.5730640
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1rsquare_sar_rem<-1-var(model_sar_rem$resid)/var(pdata$IPM)
rsquare_sar_rem## [1] 0.203698aic_sar_rem <- -2 * model_sar_rem$logLik + 2 * length(coefficients(model_sar_rem))
bic_sar_rem <- -2 * model_sar_rem$logLik + length(coefficients(model_sar_rem)) * log(length(pdata$IPM))
# Tampilkan nilai AIC
print(aic_sar_rem)## [1] 223.9103# Tampilkan nilai BIC
print(bic_sar_rem)## [1] 251.6987Kasus 2 Fixed Effect (model=“within”) untuk tanpasmp_2024
cat("ESTIMASI MODEL SAR-FEM (Spatial Autoregressive) ===\n")## ESTIMASI MODEL SAR-FEM (Spatial Autoregressive) ===model_sar_fem <- spml(formula = formtanpasmp_2024,
data = pdata,
listw = listw_terpilih,
model = "within",
effect = "individual",
lag = TRUE, # SAR hidup
spatial.error = "none") # Error mati
summary(model_sar_fem)## Spatial panel fixed effects lag model
##
##
## Call:
## spml(formula = formtanpasmp_2024, data = pdata, listw = listw_terpilih,
## model = "within", effect = "individual", lag = TRUE, spatial.error = "none")
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.81292742 -0.13950341 0.00075076 0.11708068 0.89052716
##
## Spatial autoregressive coefficient:
## Estimate Std. Error t-value Pr(>|t|)
## lambda 0.885970 0.023044 38.447 < 2.2e-16 ***
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## PPM 0.05719274 0.04189694 1.3651 0.17223
## TPAK 0.02083577 0.00906530 2.2984 0.02154 *
## AngKer 0.01139597 0.01354395 0.8414 0.40012
## JumlahSMA 0.00991970 0.00495437 2.0022 0.04526 *
## JumlahFasKes 0.00035988 0.00018872 1.9069 0.05653 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1rsquare_sar_fem<-1-var(model_sar_fem$resid)/var(pdata$IPM)
rsquare_sar_fem## [1] 0.9979953aic_sar_fem <- -2 * model_sar_fem$logLik + 2 * length(coefficients(model_sar_fem))
bic_sar_fem <- -2 * model_sar_fem$logLik + length(coefficients(model_sar_fem)) * log(length(pdata$IPM))
# Tampilkan nilai AIC
print(aic_sar_fem)## [1] 16.76558# Tampilkan nilai BIC
print(bic_sar_fem)## [1] 35.29115Kasus 1 Random Effect (model=“random”) untuk tanpasmp
cat("ESTIMASI MODEL SEM-REM (Spatial Error Model)\n")## ESTIMASI MODEL SEM-REM (Spatial Error Model)model_sem_rem <- spml(formula = formtanpasmp,
data = pdata,
listw = listw_terpilih,
model = "random",
effect = "individual",
lag = FALSE, # SAR mati
spatial.error = "b") # 'b' = Baltagi (standar untuk REM)
summary(model_sem_rem)## ML panel with , random effects, spatial error correlation
##
## Call:
## spreml(formula = formula, data = data, index = index, w = listw2mat(listw),
## w2 = listw2mat(listw2), lag = lag, errors = errors, cl = cl)
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -6.6244 -3.4243 -0.7443 -0.0039 3.0316 8.4067
##
## Error variance parameters:
## Estimate Std. Error t-value Pr(>|t|)
## phi 319.12403 110.22640 2.8952 0.00379 **
## rho 0.85038 0.04448 19.1182 < 2e-16 ***
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -1.4262e+01 1.7072e+01 -0.8354 0.403475
## PPM 8.7860e-02 8.0191e-02 1.0956 0.273242
## TPT -6.0716e-02 3.0381e-02 -1.9985 0.045660 *
## TPAK 1.3970e-02 1.0939e-02 1.2771 0.201555
## AngKer -1.0416e-03 1.4711e-02 -0.0708 0.943552
## GR 2.5092e+00 8.9738e-01 2.7961 0.005172 **
## Ln_UMK 5.6602e+00 1.1312e+00 5.0039 5.619e-07 ***
## JumlahSMA 5.2195e-03 5.1506e-03 1.0134 0.310882
## JumlahFasKes 7.7454e-05 1.9033e-04 0.4069 0.684051
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1rsquare_sem_rem<-1-var(model_sem_rem$resid)/var(pdata$IPM)
rsquare_sem_rem## [1] 0.2121558aic_sem_rem <- -2 * model_sem_rem$logLik + 2 * length(coefficients(model_sem_rem))
bic_sem_rem <- -2 * model_sem_rem$logLik + length(coefficients(model_sem_rem)) * log(length(pdata$IPM))
# Tampilkan nilai AIC
print(aic_sem_rem)## [1] 243.3608# Tampilkan nilai BIC
print(bic_sem_rem)## [1] 271.1492Kasus 2 Fixed Effect (model=“within”) untuk tanpasmp_2024
cat("ESTIMASI MODEL SEM-FEM (Spatial Error Model)\n")## ESTIMASI MODEL SEM-FEM (Spatial Error Model)model_sem_fem <- spml(formula = formtanpasmp_2024,
data = pdata,
listw = listw_terpilih,
model = "within",
effect = "individual",
lag = FALSE, # SAR mati
spatial.error = "b")
summary(model_sem_fem)## Spatial panel fixed effects error model
##
##
## Call:
## spml(formula = formtanpasmp_2024, data = pdata, listw = listw_terpilih,
## model = "within", effect = "individual", lag = FALSE, spatial.error = "b")
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.18763 -0.79760 -0.29418 0.74853 2.53057
##
## Spatial error parameter:
## Estimate Std. Error t-value Pr(>|t|)
## rho 0.924876 0.018467 50.081 < 2.2e-16 ***
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## PPM 0.20849270 0.07443354 2.8011 0.005094 **
## TPAK 0.00556744 0.00944442 0.5895 0.555529
## AngKer 0.01773614 0.01293992 1.3707 0.170483
## JumlahSMA 0.00080713 0.00408231 0.1977 0.843269
## JumlahFasKes 0.00024551 0.00017432 1.4084 0.159000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1rsquare_sem_fem<-1-var(model_sem_fem$resid)/var(pdata$IPM)
rsquare_sem_fem## [1] 0.9591772aic_sem_fem <- -2 * model_sem_fem$logLik + 2 * length(coefficients(model_sem_fem))
bic_sem_fem <- -2 * model_sem_fem$logLik + length(coefficients(model_sem_fem)) * log(length(pdata$IPM))
# Tampilkan nilai AIC
print(aic_sem_fem)## [,1]
## [1,] 387.0717# Tampilkan nilai BIC
print(bic_sem_fem)## [,1]
## [1,] 405.5972Tabel Perbandingan R-squared, AIC, dan BIC
tabel_banding <- data.frame( Model = c("SAR-REM", "SAR-FEM", "SEM-REM", "SEM-FEM"), R_Squared = c(rsquare_sar_rem, rsquare_sar_fem, rsquare_sem_rem, rsquare_sem_fem), AIC = c(aic_sar_rem, aic_sar_fem, aic_sem_rem, aic_sem_fem), BIC = c(bic_sar_rem, bic_sar_fem, bic_sem_rem, bic_sem_fem),
stringsAsFactors = FALSE
)
print(tabel_banding)## Model R_Squared AIC BIC
## 1 SAR-REM 0.2036980 223.91030 251.69867
## 2 SAR-FEM 0.9979953 16.76558 35.29115
## 3 SEM-REM 0.2121558 243.36080 271.14917
## 4 SEM-FEM 0.9591772 387.07166 405.597239. PERHITUNGAN DAMPAK (DIRECT & INDIRECT)
cat("\nDAMPAK LANGSUNG & TIDAK LANGSUNG (SPILLOVER)\n")##
## DAMPAK LANGSUNG & TIDAK LANGSUNG (SPILLOVER)imp <- impacts(model_sar_fem, listw = listw_terpilih, time = length(unique(data_panel$Tahun)))
summary(imp, zstats = TRUE, R = 500)## Impact measures (lag, trace):
## Direct Indirect Total
## PPM 0.0901012737 0.39818651 0.488287782
## TPAK 0.0328246142 0.14506253 0.177887142
## AngKer 0.0179531776 0.07934087 0.097294044
## JumlahSMA 0.0156274621 0.06906278 0.084690244
## JumlahFasKes 0.0005669559 0.00250556 0.003072516
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:200
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 200
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## PPM 0.0936331 0.070021 4.951e-03 4.951e-03
## TPAK 0.0329396 0.014020 9.914e-04 9.914e-04
## AngKer 0.0191433 0.023027 1.628e-03 1.628e-03
## JumlahSMA 0.0162326 0.007939 5.613e-04 4.629e-04
## JumlahFasKes 0.0006172 0.000313 2.213e-05 2.213e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## PPM -5.311e-02 0.0515822 0.0930650 0.136763 0.238739
## TPAK 4.636e-03 0.0246877 0.0320825 0.042287 0.061641
## AngKer -2.482e-02 0.0036878 0.0172479 0.033955 0.067024
## JumlahSMA 1.503e-03 0.0101728 0.0164320 0.022211 0.031554
## JumlahFasKes -1.847e-05 0.0004268 0.0006213 0.000802 0.001297
##
## ========================================================
## Indirect:
##
## Iterations = 1:200
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 200
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## PPM 0.429443 0.352567 0.0249303 0.0202327
## TPAK 0.149865 0.073351 0.0051867 0.0051867
## AngKer 0.087943 0.113295 0.0080111 0.0080111
## JumlahSMA 0.073391 0.038841 0.0027465 0.0027465
## JumlahFasKes 0.002816 0.001609 0.0001138 0.0001138
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## PPM -2.004e-01 0.22033 0.389730 0.596443 1.175742
## TPAK 2.170e-02 0.10115 0.144182 0.194614 0.311542
## AngKer -9.958e-02 0.01529 0.075022 0.144130 0.341984
## JumlahSMA 6.408e-03 0.04647 0.071339 0.098820 0.162212
## JumlahFasKes -7.006e-05 0.00184 0.002638 0.003692 0.006757
##
## ========================================================
## Total:
##
## Iterations = 1:200
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 200
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## PPM 0.523076 0.420534 0.0297363 0.0297363
## TPAK 0.182805 0.086560 0.0061207 0.0061207
## AngKer 0.107087 0.135841 0.0096054 0.0096054
## JumlahSMA 0.089624 0.046407 0.0032815 0.0032815
## JumlahFasKes 0.003433 0.001908 0.0001349 0.0001349
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## PPM -2.461e-01 0.270532 0.488623 0.745636 1.445807
## TPAK 2.761e-02 0.126119 0.173699 0.235097 0.363233
## AngKer -1.245e-01 0.019355 0.091279 0.178763 0.393898
## JumlahSMA 7.911e-03 0.057386 0.088774 0.119229 0.193869
## JumlahFasKes -8.634e-05 0.002256 0.003304 0.004482 0.008175
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## PPM 0.0700205659 0.352567043 0.420534297
## TPAK 0.0140198927 0.073351096 0.086559611
## AngKer 0.0230270238 0.113294690 0.135841149
## JumlahSMA 0.0079386233 0.038841256 0.046407016
## JumlahFasKes 0.0003129854 0.001608794 0.001908347
##
## Simulated z-values:
## Direct Indirect Total
## PPM 1.3372225 1.2180466 1.2438371
## TPAK 2.3494880 2.0431205 2.1118937
## AngKer 0.8313388 0.7762367 0.7883234
## JumlahSMA 2.0447577 1.8895206 1.9312578
## JumlahFasKes 1.9720772 1.7502935 1.7989881
##
## Simulated p-values:
## Direct Indirect Total
## PPM 0.181150 0.223206 0.213560
## TPAK 0.018799 0.041041 0.034696
## AngKer 0.405782 0.437609 0.430508
## JumlahSMA 0.040879 0.058822 0.053451
## JumlahFasKes 0.048601 0.080068 0.0720219. UJI ASUMSI KLASIK (SISAAN)
#Uji Asumsi kenormalan
galat1 <- model_sar_fem$residuals
ks.test(galat1, "pnorm", mean=mean(galat1), sd=sd(galat1)) ##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: galat1
## D = 0.064188, p-value = 0.5168
## alternative hypothesis: two-sidedqqnorm(galat1)
qqline(galat1)
ad.test(galat1)##
## Anderson-Darling normality test
##
## data: galat1
## A = 0.63584, p-value = 0.09577residuals <- model_sar_fem$residuals
residuals_squared <- residuals^2
breusch_pagan_model <- lm(residuals_squared ~ PPM + TPAK + AngKer + JumlahSMA + JumlahFasKes, data=pdata)
summary(breusch_pagan_model)##
## Call:
## lm(formula = residuals_squared ~ PPM + TPAK + AngKer + JumlahSMA +
## JumlahFasKes, data = pdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06968 -0.03845 -0.02229 0.00534 0.73685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.583e-01 1.754e-01 -0.902 0.368
## PPM 2.073e-03 3.469e-03 0.597 0.551
## TPAK 2.230e-03 1.997e-03 1.117 0.266
## AngKer 7.871e-04 2.675e-03 0.294 0.769
## JumlahSMA -1.254e-04 3.544e-04 -0.354 0.724
## JumlahFasKes 5.400e-06 1.275e-05 0.424 0.672
##
## Residual standard error: 0.09215 on 156 degrees of freedom
## Multiple R-squared: 0.0233, Adjusted R-squared: -0.008001
## F-statistic: 0.7444 on 5 and 156 DF, p-value: 0.5913