1. Persiapan Data
1.1 Import data
panel
data_nsc <- readr::read_csv(
"C:/Users/nyayu/Statistika/Semester 6/NSC/Data/data nsc.csv",
show_col_types = FALSE
)
# Pastikan seluruh kolom selain kabkot berbentuk numerik
data_nsc <- data_nsc %>%
mutate(across(-kabkot, ~ readr::parse_number(as.character(.))))
str(data_nsc)
## tibble [514 × 26] (S3: tbl_df/tbl/data.frame)
## $ kabkot : chr [1:514] "Simeulue" "Aceh Singkil" "Aceh Selatan" "Aceh Tenggara" ...
## $ ipm_2020 : num [1:514] 67.9 69.9 69.8 71.2 69 ...
## $ ipm_2021 : num [1:514] 68.3 70.2 70.1 71.3 69.2 ...
## $ ipm_2022 : num [1:514] 69.2 70.6 70.6 72.2 70.1 ...
## $ ipm_2023 : num [1:514] 70 71.1 71.1 72.9 70.7 ...
## $ ipm_2024 : num [1:514] 71 71.8 71.8 73.6 71.3 ...
## $ pdrb_2020: num [1:514] 1602 1714 4241 3436 8273 ...
## $ pdrb_2021: num [1:514] 1648 1780 4346 3487 8434 ...
## $ pdrb_2022: num [1:514] 1708 1845 4481 3584 8748 ...
## $ pdrb_2023: num [1:514] 1807 1910 4639 3693 8931 ...
## $ pdrb_2024: num [1:514] 1887 1976 4788 3826 9345 ...
## $ tpt_2020 : num [1:514] 5.47 8.24 6.54 5.72 7.26 3.05 7.3 7.62 6.45 4.12 ...
## $ tpt_2021 : num [1:514] 5.71 8.36 6.46 6.43 7.13 2.61 7.09 7.7 7.28 4.32 ...
## ..- attr(*, "problems")= tibble [2 × 4] (S3: tbl_df/tbl/data.frame)
## .. ..$ row : int [1:2] 490 492
## .. ..$ col : int [1:2] NA NA
## .. ..$ expected: chr [1:2] "a number" "a number"
## .. ..$ actual : chr [1:2] "-" "-"
## $ tpt_2022 : num [1:514] 6 6.88 4.82 5.09 8.07 4.44 6.09 8.28 5.94 4.2 ...
## ..- attr(*, "problems")= tibble [6 × 4] (S3: tbl_df/tbl/data.frame)
## .. ..$ row : int [1:6] 489 492 495 503 504 505
## .. ..$ col : int [1:6] NA NA NA NA NA NA
## .. ..$ expected: chr [1:6] "a number" "a number" "a number" "a number" ...
## .. ..$ actual : chr [1:6] "-" "-" "-" "-" ...
## $ tpt_2023 : num [1:514] 5.85 6.84 4.73 5 8.03 4.42 6.07 8.17 5.92 4.14 ...
## ..- attr(*, "problems")= tibble [12 × 4] (S3: tbl_df/tbl/data.frame)
## .. ..$ row : int [1:12] 489 490 491 492 493 495 499 500 502 503 ...
## .. ..$ col : int [1:12] NA NA NA NA NA NA NA NA NA NA ...
## .. ..$ expected: chr [1:12] "a number" "a number" "a number" "a number" ...
## .. ..$ actual : chr [1:12] "-" "-" "-" "-" ...
## $ tpt_2024 : num [1:514] 5.5 6.44 4.56 4.79 7.75 4.21 5.58 7.93 5.74 3.93 ...
## ..- attr(*, "problems")= tibble [1 × 4] (S3: tbl_df/tbl/data.frame)
## .. ..$ row : int 504
## .. ..$ col : int NA
## .. ..$ expected: chr "a number"
## .. ..$ actual : chr "-"
## $ tpak_2020: num [1:514] 70.4 62 61.4 71.3 61.9 ...
## $ tpak_2021: num [1:514] 71.1 62.8 60.8 69.6 59.5 ...
## $ tpak_2022: num [1:514] 64.4 57.3 60.7 67.8 58.4 ...
## $ tpak_2023: num [1:514] 70.6 57.8 58.9 70.4 60.5 ...
## ..- attr(*, "problems")= tibble [13 × 4] (S3: tbl_df/tbl/data.frame)
## .. ..$ row : int [1:13] 490 491 492 493 495 498 499 500 501 502 ...
## .. ..$ col : int [1:13] NA NA NA NA NA NA NA NA NA NA ...
## .. ..$ expected: chr [1:13] "a number" "a number" "a number" "a number" ...
## .. ..$ actual : chr [1:13] "-" "-" "-" "-" ...
## $ tpak_2024: num [1:514] 71.1 57.8 59.2 70.5 60.8 ...
## $ san_2020 : num [1:514] 66.5 64.2 68.8 54 74.6 ...
## $ san_2021 : num [1:514] 71.6 69.6 62.5 62.7 66.8 ...
## ..- attr(*, "problems")= tibble [1 × 4] (S3: tbl_df/tbl/data.frame)
## .. ..$ row : int 500
## .. ..$ col : int NA
## .. ..$ expected: chr "a number"
## .. ..$ actual : chr "-"
## $ san_2022 : num [1:514] 70.2 78.9 70.8 62.9 69.1 ...
## $ san_2023 : num [1:514] 72.7 67.9 74.5 65.1 70 ...
## $ san_2024 : num [1:514] 76.5 63.2 78 64 73.6 ...
## ..- attr(*, "problems")= tibble [1 × 4] (S3: tbl_df/tbl/data.frame)
## .. ..$ row : int 504
## .. ..$ col : int NA
## .. ..$ expected: chr "a number"
## .. ..$ actual : chr "-"
1.3 Imputasi missing
value
data_long <- data_long %>%
group_by(kabkot) %>%
arrange(tahun, .by_group = TRUE) %>%
mutate(
tpt = na.approx(tpt, x = tahun, na.rm = FALSE),
tpak = na.approx(tpak, x = tahun, na.rm = FALSE),
san = na.approx(san, x = tahun, na.rm = FALSE)
) %>%
fill(tpt, tpak, san, .direction = "downup") %>%
ungroup()
colSums(is.na(data_long))
## kabkot tahun ipm pdrb tpt tpak san
## 0 0 0 0 0 0 0
1.4 Cleaning nilai
TPAK yang tidak wajar
data_long <- data_long %>%
mutate(tpak = ifelse(tpak < 10, NA_real_, tpak)) %>%
group_by(kabkot) %>%
arrange(tahun, .by_group = TRUE) %>%
mutate(tpak = na.approx(tpak, x = tahun, na.rm = FALSE)) %>%
fill(tpak, .direction = "downup") %>%
ungroup()
colSums(is.na(data_long))
## kabkot tahun ipm pdrb tpt tpak san
## 0 0 0 0 0 0 0
1.6 Import data
spasial
shp_kabkota <- sf::st_read(
"C:/Users/nyayu/Statistika/Semester 6/NSC/Data/shp indonesia/Kab_Kota.shp",
quiet = TRUE
)
names(shp_kabkota)
## [1] "KODE_KK" "KODE_PROV" "KAB_KOTA" "PROVINSI" "FID" "geometry"
## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.27241 ymin: -9.829316 xmax: 119.9263 ymax: 5.61174
## Geodetic CRS: WGS 84
## KODE_KK KODE_PROV KAB_KOTA PROVINSI FID
## 1 11.71 11 Kota Banda Aceh Aceh 11.71
## 2 12.72 12 Kota Pematangsiantar Sumatera Utara 12.72
## 3 35.28 35 Pamekasan Jawa Timur 35.28
## 4 32.16 32 Bekasi Jawa Barat 32.16
## 5 53.17 53 Sumba Tengah Nusa Tenggara Timur 53.17
## 6 61.07 61 Bengkayang Kalimantan Barat 61.07
## geometry
## 1 MULTIPOLYGON (((95.30342 5....
## 2 MULTIPOLYGON (((99.03411 2....
## 3 MULTIPOLYGON (((113.4331 -7...
## 4 MULTIPOLYGON (((107.0517 -6...
## 5 MULTIPOLYGON (((119.5395 -9...
## 6 MULTIPOLYGON (((109.0315 0....
1.7 Harmonisasi nama
kabupaten/kota
data_long <- data_long %>%
mutate(
kabkot = dplyr::recode(
kabkot,
"Banyu Asin" = "Banyuasin",
"Batang Hari" = "Batanghari",
"Fakfak" = "Fak Fak",
"Gunung Kidul" = "Gunungkidul",
"Karang Asem" = "Karangasem",
"Kepulauan Seribu" = "Administrasi Kepulauan Seribu",
"Kepulauan Sitaro" = "Kepulauan Siau Tagulandang Biaro",
"Kepulauan Tanimbar (Maluku Tenggara Barat)" = "Kepulauan Tanimbar",
"Kota Banjar Baru" = "Kota Banjarbaru",
"Kota Baru" = "Kotabaru",
"Kota Baubau" = "Kota Bau Bau",
"Kota Jakarta Barat" = "Kota Administrasi Jakarta Barat",
"Kota Jakarta Pusat" = "Kota Administrasi Jakarta Pusat",
"Kota Jakarta Selatan" = "Kota Administrasi Jakarta Selatan",
"Kota Jakarta Timur" = "Kota Administrasi Jakarta Timur",
"Kota Jakarta Utara" = "Kota Administrasi Jakarta Utara",
"Kota Lubuklinggau" = "Kota Lubuk Linggau",
"Kota Padangsidimpuan" = "Kota Padang Sidempuan",
"Kota Palangka Raya" = "Kota Palangkaraya",
"Kota Parepare" = "Kota Pare Pare",
"Kota Pematang Siantar" = "Kota Pematangsiantar",
"Kota Sawah Lunto" = "Kota Sawahlunto",
"Labuhan Batu" = "Labuhanbatu",
"Labuhan Batu Selatan" = "Labuhanbatu Selatan",
"Labuhan Batu Utara" = "Labuhanbatu Utara",
"Mahakam Hulu" = "Mahakam Ulu",
"Mukomuko" = "Muko Muko",
"Pangkajene Dan Kepulauan" = "Pangkajene Kepulauan",
"Pasangkayu (Mamuju Utara)" = "Pasangkayu",
"Toba Samosir" = "Toba",
"Tojo Una-una" = "Tojo Una Una",
"Tolitoli" = "Toli Toli",
"Tulangbawang" = "Tulang Bawang"
)
)
setdiff(sort(unique(data_long$kabkot)), sort(unique(shp_kabkota$KAB_KOTA)))
## character(0)
setdiff(sort(unique(shp_kabkota$KAB_KOTA)), sort(unique(data_long$kabkot)))
## character(0)
1.8 Merge data panel
dan data spasial
data_spasial <- shp_kabkota %>%
left_join(data_long, by = c("KAB_KOTA" = "kabkot"))
sum(is.na(data_spasial$ipm))
## [1] 0
3. Analisis Kawasan
Indonesia Timur
timur <- c(
"Nusa Tenggara Timur",
"Maluku",
"Maluku Utara",
"Papua",
"Papua Barat",
"Papua Selatan",
"Papua Tengah",
"Papua Pegunungan",
"Papua Barat Daya"
)
data_timur <- data_spasial %>%
filter(PROVINSI %in% timur) %>%
arrange(KAB_KOTA, tahun)
panel_timur <- st_drop_geometry(data_timur)
pdata_timur <- pdata.frame(
panel_timur,
index = c("KAB_KOTA", "tahun")
)
table(panel_timur$tahun)
##
## 2020 2021 2022 2023 2024
## 85 85 85 85 85
n_distinct(panel_timur$KAB_KOTA)
## [1] 85
3.1 Choropleth IPM
kawasan timur tahun 2024
data_timur_2024 <- data_timur %>%
filter(tahun == 2024) %>%
arrange(KAB_KOTA)
tmap_mode("plot")
peta_ipm_timur <- tm_shape(data_timur_2024) +
tm_polygons(
fill = "ipm",
fill.scale = tm_scale_intervals(
style = "jenks",
n = 5,
values = "Blues"
),
fill.legend = tm_legend(title = "IPM")
) +
tm_title("Peta IPM Kawasan Indonesia Timur Tahun 2024") +
tm_layout(frame = FALSE, legend.outside = TRUE)
peta_ipm_timur
3.2 Spatial weights
dengan KNN
coords_timur <- st_coordinates(
st_centroid(st_geometry(data_timur_2024))
)
k <- 8
repeat {
knn_timur <- suppressWarnings(knearneigh(coords_timur, k = k))
neighbors_timur <- suppressWarnings(knn2nb(knn_timur))
if (spdep::n.comp.nb(neighbors_timur)$nc == 1) break
k <- k + 1
if (k >= nrow(coords_timur)) stop("KNN tidak dapat membentuk graph yang terhubung.")
}
listw_timur <- nb2listw(neighbors_timur, style = "W")
k
## [1] 8
spdep::n.comp.nb(neighbors_timur)$nc
## [1] 1
3.3 Uji Moran’s
I
moran_ipm_timur <- moran.test(
data_timur_2024$ipm,
listw_timur
)
moran_ipm_timur
##
## Moran I test under randomisation
##
## data: data_timur_2024$ipm
## weights: listw_timur
##
## Moran I statistic standard deviate = 8.6522, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.399751197 -0.011904762 0.002263662
4. Regresi Panel
Global
4.1 Common Effect
Model
model_cem_timur <- plm(
ipm ~ log_pdrb + tpt + tpak + san,
data = pdata_timur,
model = "pooling"
)
summary(model_cem_timur)
## Pooling Model
##
## Call:
## plm(formula = ipm ~ log_pdrb + tpt + tpak + san, data = pdata_timur,
## model = "pooling")
##
## Balanced Panel: n = 85, T = 5, N = 425
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -13.55976 -2.06294 0.18557 2.44127 9.24152
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 37.9460196 2.5205812 15.0545 < 2.2e-16 ***
## log_pdrb 2.6906244 0.2084365 12.9086 < 2.2e-16 ***
## tpt 0.6073960 0.0945485 6.4242 3.598e-10 ***
## tpak -0.0956089 0.0232195 -4.1176 4.611e-05 ***
## san 0.1733031 0.0082514 21.0029 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 31584
## Residual Sum of Squares: 5254.8
## R-Squared: 0.83363
## Adj. R-Squared: 0.83204
## F-statistic: 526.115 on 4 and 420 DF, p-value: < 2.22e-16
4.2 Fixed Effect
Model
model_fem_timur <- plm(
ipm ~ log_pdrb + tpt + tpak + san,
data = pdata_timur,
model = "within"
)
summary(model_fem_timur)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = ipm ~ log_pdrb + tpt + tpak + san, data = pdata_timur,
## model = "within")
##
## Balanced Panel: n = 85, T = 5, N = 425
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.133609 -0.554622 -0.082448 0.578238 4.026226
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## log_pdrb 0.9961034 0.1548847 6.4313 4.35e-10 ***
## tpt -0.2855049 0.0495665 -5.7600 1.90e-08 ***
## tpak 0.0031148 0.0098506 0.3162 0.752
## san 0.0407611 0.0069914 5.8302 1.30e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 399.47
## Residual Sum of Squares: 289.79
## R-Squared: 0.27457
## Adj. R-Squared: 0.084573
## F-statistic: 31.793 on 4 and 336 DF, p-value: < 2.22e-16
4.3 Random Effect
Model
model_rem_timur <- plm(
ipm ~ log_pdrb + tpt + tpak + san,
data = pdata_timur,
model = "random"
)
summary(model_rem_timur)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = ipm ~ log_pdrb + tpt + tpak + san, data = pdata_timur,
## model = "random")
##
## Balanced Panel: n = 85, T = 5, N = 425
##
## Effects:
## var std.dev share
## idiosyncratic 0.8625 0.9287 0.078
## individual 10.1404 3.1844 0.922
## theta: 0.8707
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -4.52790 -0.62807 0.12740 0.80180 3.65973
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 48.6181723 1.8388213 26.4399 <2e-16 ***
## log_pdrb 1.6207514 0.1844172 8.7885 <2e-16 ***
## tpt -0.0639767 0.0606084 -1.0556 0.2912
## tpak -0.0170710 0.0124380 -1.3725 0.1699
## san 0.0851405 0.0080488 10.5781 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 921.06
## Residual Sum of Squares: 599.54
## R-Squared: 0.34907
## Adj. R-Squared: 0.34287
## Chisq: 225.234 on 4 DF, p-value: < 2.22e-16
4.4 Uji Chow dan
Hausman
chow_timur <- pFtest(model_fem_timur, model_cem_timur)
hausman_timur <- phtest(model_fem_timur, model_rem_timur)
chow_timur
##
## F test for individual effects
##
## data: ipm ~ log_pdrb + tpt + tpak + san
## F = 68.533, df1 = 84, df2 = 336, p-value < 2.2e-16
## alternative hypothesis: significant effects
##
## Hausman Test
##
## data: ipm ~ log_pdrb + tpt + tpak + san
## chisq = 484.71, df = 4, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent
4.5 Ringkasan
pemilihan model
tabel_pemilihan <- data.frame(
Uji = c("Chow", "Hausman"),
Statistik = c(
round(chow_timur$statistic, 3),
round(hausman_timur$statistic, 3)
),
p_value = c(
format(chow_timur$p.value, scientific = TRUE, digits = 3),
format(hausman_timur$p.value, scientific = TRUE, digits = 3)
)
)
knitr::kable(
tabel_pemilihan,
caption = "Hasil Uji Pemilihan Model Panel"
)
Hasil Uji Pemilihan Model Panel
| F |
Chow |
68.533 |
2.23e-169 |
| chisq |
Hausman |
484.712 |
1.36e-103 |
5. Diagnostik
Panel
5.1
Heteroskedastisitas
bp_timur <- bptest(model_fem_timur)
bp_timur
##
## studentized Breusch-Pagan test
##
## data: model_fem_timur
## BP = 43.874, df = 4, p-value = 6.816e-09
5.2 Cross-sectional
dependence
pcd_timur <- pcdtest(model_fem_timur, test = "cd")
pcd_timur
##
## Pesaran CD test for cross-sectional dependence in panels
##
## data: ipm ~ log_pdrb + tpt + tpak + san
## z = 79.087, p-value < 2.2e-16
## alternative hypothesis: cross-sectional dependence
5.3 Moran residual
FEM
panel_timur$res_fem <- residuals(model_fem_timur)
res_fem_2024 <- panel_timur %>%
filter(tahun == 2024) %>%
arrange(KAB_KOTA)
moran_res_fem <- moran.test(
res_fem_2024$res_fem,
listw_timur
)
moran_res_fem
##
## Moran I test under randomisation
##
## data: res_fem_2024$res_fem
## weights: listw_timur
##
## Moran I statistic standard deviate = -0.58219, p-value = 0.7198
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.039063581 -0.011904762 0.002176164
5.4 Normalitas
residual
shapiro_timur <- shapiro.test(residuals(model_fem_timur))
shapiro_timur
##
## Shapiro-Wilk normality test
##
## data: residuals(model_fem_timur)
## W = 0.98729, p-value = 0.0009076
5.5 Driscoll-Kraay
robust standard errors
dk_timur <- coeftest(
model_fem_timur,
vcov. = vcovSCC(model_fem_timur, type = "HC1", maxlag = 1)
)
dk_timur
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## log_pdrb 0.9961034 0.2410890 4.1317 4.551e-05 ***
## tpt -0.2855049 0.0485624 -5.8791 9.949e-09 ***
## tpak 0.0031148 0.0120935 0.2576 0.7969047
## san 0.0407611 0.0114223 3.5686 0.0004112 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
6. Spatial Panel
Regression
6.1 SAR-FEM
sar_fem_timur <- spml(
ipm ~ log_pdrb + tpt + tpak + san,
data = pdata_timur,
listw = listw_timur,
model = "within",
effect = "individual",
lag = TRUE,
spatial.error = "none"
)
summary(sar_fem_timur)
## Spatial panel fixed effects lag model
##
##
## Call:
## spml(formula = ipm ~ log_pdrb + tpt + tpak + san, data = pdata_timur,
## listw = listw_timur, model = "within", effect = "individual",
## lag = TRUE, spatial.error = "none")
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.4750114 -0.1853183 0.0024392 0.1593558 4.0683960
##
## Spatial autoregressive coefficient:
## Estimate Std. Error t-value Pr(>|t|)
## lambda 0.794645 0.027696 28.692 < 2.2e-16 ***
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## log_pdrb 0.5079701 0.0790677 6.4245 1.323e-10 ***
## tpt -0.0423404 0.0253878 -1.6677 0.09537 .
## tpak -0.0046953 0.0050066 -0.9378 0.34834
## san 0.0103915 0.0035924 2.8926 0.00382 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
6.2 SEM-FEM
sem_fem_timur <- spml(
ipm ~ log_pdrb + tpt + tpak + san,
data = pdata_timur,
listw = listw_timur,
model = "within",
effect = "individual",
lag = FALSE,
spatial.error = "b"
)
summary(sem_fem_timur)
## Spatial panel fixed effects error model
##
##
## Call:
## spml(formula = ipm ~ log_pdrb + tpt + tpak + san, data = pdata_timur,
## listw = listw_timur, model = "within", effect = "individual",
## lag = FALSE, spatial.error = "b")
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.18225 -0.78169 -0.10675 0.66069 4.80708
##
## Spatial error parameter:
## Estimate Std. Error t-value Pr(>|t|)
## rho 0.840116 0.024634 34.103 < 2.2e-16 ***
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## log_pdrb 0.4085869 0.0749813 5.4492 5.06e-08 ***
## tpt -0.0065347 0.0267345 -0.2444 0.80690
## tpak -0.0076356 0.0050021 -1.5265 0.12690
## san 0.0074774 0.0037888 1.9736 0.04843 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
6.3 Perbandingan
model spasial
aic_sar <- tryCatch(AIC(sar_fem_timur), error = function(e) NA_real_)
aic_sem <- tryCatch(AIC(sem_fem_timur), error = function(e) NA_real_)
ll_sar <- tryCatch(as.numeric(logLik(sar_fem_timur)), error = function(e) NA_real_)
ll_sem <- tryCatch(as.numeric(logLik(sem_fem_timur)), error = function(e) NA_real_)
tabel_spatial <- data.frame(
Model = c("SAR-FEM", "SEM-FEM"),
AIC = c(aic_sar, aic_sem),
LogLik = c(ll_sar, ll_sem)
)
knitr::kable(
tabel_spatial,
caption = "Perbandingan Model Spatial Panel"
)
Perbandingan Model Spatial Panel
| SAR-FEM |
NA |
NA |
| SEM-FEM |
NA |
NA |
6.4 Moran residual
SAR-FEM
panel_timur$res_sar <- residuals(sar_fem_timur)
res_timur_2024 <- panel_timur %>%
filter(tahun == 2024) %>%
arrange(KAB_KOTA)
data_timur_2024$res_sar <- res_timur_2024$res_sar
moran_res_sar <- moran.test(
data_timur_2024$res_sar,
listw_timur
)
moran_res_sar
##
## Moran I test under randomisation
##
## data: data_timur_2024$res_sar
## weights: listw_timur
##
## Moran I statistic standard deviate = 0.14722, p-value = 0.4415
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.005452596 -0.011904762 0.001920779
## [1] -308.2943
## [,1]
## [1,] -999.6365
loglik_sar <- as.numeric(summary(sar_fem_timur)$logLik)
loglik_sem <- as.numeric(summary(sem_fem_timur)$logLik)
k_sar <- length(coef(sar_fem_timur)) + 1
k_sem <- length(coef(sem_fem_timur)) + 1
aic_sar <- -2 * loglik_sar + 2 * k_sar
aic_sem <- -2 * loglik_sem + 2 * k_sem
data.frame(
Model = c("SAR-FEM", "SEM-FEM"),
AIC = c(aic_sar, aic_sem),
LogLik = c(loglik_sar, loglik_sem)
)
## Model AIC LogLik
## 1 SAR-FEM 628.5886 -308.2943
## 2 SEM-FEM 2011.2729 -999.6365
7. Ringkasan hasil
utama
cat("========================================\n")
## ========================================
cat(" RINGKASAN HASIL PENELITIAN\n")
## RINGKASAN HASIL PENELITIAN
cat("========================================\n\n")
## ========================================
cat("1. Pemilihan model panel:\n")
## 1. Pemilihan model panel:
cat(" Chow : F = ", round(chow_timur$statistic, 3),
", p-value = ", format(chow_timur$p.value, scientific = TRUE, digits = 3), "\n", sep = "")
## Chow : F = 68.533, p-value = 2.23e-169
cat(" Hausman : Chi-square = ", round(hausman_timur$statistic, 3),
", p-value = ", format(hausman_timur$p.value, scientific = TRUE, digits = 3), "\n", sep = "")
## Hausman : Chi-square = 484.712, p-value = 1.36e-103
cat(" Model panel terbaik: Fixed Effect Model (FEM)\n\n")
## Model panel terbaik: Fixed Effect Model (FEM)
cat("2. Diagnostik panel:\n")
## 2. Diagnostik panel:
cat(" Breusch-Pagan : BP = ", round(bp_timur$statistic, 3),
", p-value = ", format(bp_timur$p.value, scientific = TRUE, digits = 3), "\n", sep = "")
## Breusch-Pagan : BP = 43.874, p-value = 6.82e-09
cat(" Pesaran CD : z = ", round(pcd_timur$statistic, 3),
", p-value = ", format(pcd_timur$p.value, scientific = TRUE, digits = 3), "\n", sep = "")
## Pesaran CD : z = 79.087, p-value = 0e+00
cat(" Driscoll-Kraay SE sudah dihitung untuk FEM.\n\n")
## Driscoll-Kraay SE sudah dihitung untuk FEM.
cat("3. Spatial panel:\n")
## 3. Spatial panel:
cat(" SAR-FEM dan SEM-FEM telah diestimasi.\n")
## SAR-FEM dan SEM-FEM telah diestimasi.
cat(" Pilih model dengan AIC lebih kecil dan LogLik lebih besar.\n")
## Pilih model dengan AIC lebih kecil dan LogLik lebih besar.
