Library
library(spdep)
## Warning: package 'spdep' was built under R version 4.4.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.4.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Warning: package 'sf' was built under R version 4.4.3
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.4.3
library(sf)
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(spData)
library(geodata)
## Warning: package 'geodata' was built under R version 4.4.3
## Loading required package: terra
## Warning: package 'terra' was built under R version 4.4.3
## terra 1.8.60
library(tmap)
## Warning: package 'tmap' was built under R version 4.4.3
library(skimr)
## Warning: package 'skimr' was built under R version 4.4.3
library(psych)
## Warning: package 'psych' was built under R version 4.4.3
##
## Attaching package: 'psych'
## The following objects are masked from 'package:terra':
##
## describe, distance, rescale
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.4.3
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(grid)
##
## Attaching package: 'grid'
## The following object is masked from 'package:terra':
##
## depth
library(viridis)
## Warning: package 'viridis' was built under R version 4.4.3
## Loading required package: viridisLite
library(terra)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:psych':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
library(Matrix)
library(spatial)
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.4.3
##
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption
Upload Data
setwd("C:\\Users\\Ghozzy\\OneDrive\\Documents\\SMESTER 5\\Sparsial Stat")
indo_adm2 <- st_read("IDN_AdminBoundaries_candidate.gdb",
layer = "idn_admbnda_adm2_bps_20200401")
## Reading layer `idn_admbnda_adm2_bps_20200401' from data source
## `C:\Users\Ghozzy\OneDrive\Documents\SMESTER 5\Sparsial Stat\IDN_AdminBoundaries_candidate.gdb'
## using driver `OpenFileGDB'
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received a polygon with more than 100 parts. The
## processing may be really slow. You can skip the processing by setting
## METHOD=SKIP.
## Simple feature collection with 522 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS: WGS 84
jabar <- indo_adm2[indo_adm2$admin1Name_en == "Jawa Barat", ]
jabar <- jabar[jabar$admin2Name_en != "Waduk Cirata", ]
plot(st_geometry(jabar))
data <- na.omit(read.csv("spatial_data_uts.csv", sep=";", dec="."))
jabar$id <- c(1:27)
jabar_sf <- st_as_sf(jabar)
data$id <- c(1:27)
jabar_merged <- jabar_sf %>%
left_join(data, by = "id")
Peta IPM
ggplot(data = jabar_merged) +
geom_sf(aes(fill = Y..IPM.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(
option = "plasma",
direction = -1,
name = "Indeks Pembangunan Manusia (IPM)"
) +
labs(
title = "Peta IPM Pulau Jawa 2024",
subtitle = "Nilai Indeks Pembangunan Manusia per Kabupaten/Kota",
caption = "Sumber: BPS (2024)"
) +
theme_minimal(base_size = 14) +
theme(
panel.background = element_rect(fill = "white", color = NA),
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
Peta PDRB
ggplot(data = jabar_merged) +
geom_sf(aes(fill = X1..PDRB.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "Produk Domestik Regional Bruto (PDRB)") +
labs(
title = "Sebaran PDRB di Jawa Barat",
subtitle = "Nilai Produk Domestik Regional Bruto per kabupaten/kota",
caption = "Sumber: BPS (2024)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
Peta Indeks Pengetahuan
ggplot(data = jabar_merged) +
geom_sf(aes(fill = X2..Indeks.Pengetahuan.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "inferno", direction = -1,
name = "Indeks Pengetahuan") +
labs(
title = "Sebaran Indeks Pengetahuan di Jawa Barat",
subtitle = "Nilai indeks pengetahuan per kabupaten/kota",
caption = "Sumber: BPS (2024)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
Peta UMR
ggplot(data = jabar_merged) +
geom_sf(aes(fill = X3..UMR.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "cividis", direction = -1,
name = "Upah Minimum Rakyat (UMR)") +
labs(
title = "Sebaran UMR di Jawa Barat",
subtitle = "Nilai upah minimum per kabupaten/kota",
caption = "Sumber: BPS (2024)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
Peta Pengeluaran per Kapita
ggplot(data = jabar_merged) +
geom_sf(aes(fill = X4..Pengeluaran.per.kapita.),
color = "white", linewidth = 0.3) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "Pengeluaran per Kapita") +
labs(
title = "Sebaran Pengeluaran per Kapita di Jawa Barat",
subtitle = "Nilai pengeluaran per kapita per kabupaten/kota",
caption = "Sumber: BPS (2024)"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.major = element_line(color = "gray90"),
panel.grid.minor = element_line(color = "gray95"),
plot.title = element_text(face = "bold", hjust = 0.5)
)
Deskriptif Data
data_desc <- jabar_merged %>%
st_drop_geometry() %>%
select(
Y..IPM.,
X1..PDRB.,
X2..Indeks.Pengetahuan.,
X3..UMR.,
X4..Pengeluaran.per.kapita.
)
describe(data_desc)
## vars n mean sd median trimmed
## Y..IPM. 1 27 74.68 4.33 73.82 74.36
## X1..PDRB. 2 27 105053.23 112708.74 61236.34 87923.66
## X2..Indeks.Pengetahuan. 3 27 65.60 6.78 63.38 65.20
## X3..UMR. 4 27 3370534.44 1106474.44 3294485.00 3315578.61
## X4..Pengeluaran.per.kapita. 5 27 11928.85 2389.91 11563.00 11674.70
## mad min max range skew
## Y..IPM. 4.42 68.89 83.75 14.86 0.73
## X1..PDRB. 40057.20 5549.26 421323.50 415774.24 1.64
## X2..Indeks.Pengetahuan. 6.23 57.39 78.58 21.19 0.59
## X3..UMR. 1356318.06 2070192.00 5343430.00 3273238.00 0.45
## X4..Pengeluaran.per.kapita. 2044.51 8965.00 18795.00 9830.00 1.17
## kurtosis se
## Y..IPM. -0.59 0.83
## X1..PDRB. 1.37 21690.81
## X2..Indeks.Pengetahuan. -1.12 1.30
## X3..UMR. -1.27 212941.11
## X4..Pengeluaran.per.kapita. 1.06 459.94
# Rebuild queen contiguity weights
# --- 1. spatial weights (queen) ---
nb <- poly2nb(as_Spatial(jabar_merged), queen = TRUE)
lwW <- nb2listw(nb, style = "W", zero.policy = TRUE)
lwB <- nb2listw(nb, style = "B", zero.policy = TRUE)
summary(nb)
## 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
which(card(nb) == 0)
## integer(0)
FUNCTION
analyze_spatial <- function(data_sf, var_col, label) {
cat("\n\n==============================\n", label, "\n==============================\n")
x_raw <- data_sf[[var_col]]
# ---------- Global tests ----------
moran_res <- moran.test(x_raw, lwW, randomisation = TRUE, alternative = "two.sided")
geary_res <- geary.test(x_raw, lwW, randomisation = TRUE, alternative = "two.sided")
print(moran_res)
print(geary_res)
# ---------- Local Moran ----------
x <- scale(x_raw)[, 1]
lagx <- lag.listw(lwW, x)
lisa <- localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
lisa_df <- as.data.frame(lisa)
names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
alpha <- 0.05
quad <- dplyr::case_when(
x >= 0 & lagx >= 0 ~ "High-High",
x < 0 & lagx < 0 ~ "Low-Low",
x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
)
LISA_sf <- bind_cols(data_sf, lisa_df) |>
mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))
# ---------- Local Geary’s C ----------
localC_vals <- localC(x_raw, lwW)
Geary_sf <- mutate(data_sf,
LocalC = as.numeric(localC_vals))
# ---------- Getis–Ord Gi* ----------
sum_x <- sum(x_raw)
Wb <- listw2mat(lwB)
Wb_star <- Wb; diag(Wb_star) <- 1
G_star_raw <- as.numeric(Wb_star %*% x_raw) / sum_x
Gz <- localG(x_raw, listw = lwW)
G_sf <- mutate(data_sf,
z_Gistar = as.numeric(Gz),
hotcold = case_when(
z_Gistar >= 1.96 ~ "Hot spot (p<0.05)",
z_Gistar <= -1.96 ~ "Cold spot (p<0.05)",
TRUE ~ "Not significant"
))
# ---------- MAPS ----------
p1 <- ggplot(LISA_sf) +
geom_sf(aes(fill = quad), color = "white", size = 0.2) +
scale_fill_manual(values = c(
"High-High" = "#d73027",
"Low-Low" = "#4575b4",
"High-Low (Outlier)" = "#fdae61",
"Low-High (Outlier)" = "#74add1",
"Not significant" = "grey85"
)) +
labs(title = paste("Local Moran's I (LISA) —", label), fill = "Kategori") +
theme_minimal()
p2 <- ggplot(Geary_sf) +
geom_sf(aes(fill = LocalC), color = "white", size = 0.2) +
scale_fill_viridis_c(option = "magma", direction = -1,
name = "Local Geary’s C") +
labs(title = paste("Local Geary’s C —", label)) +
theme_minimal()
p3 <- ggplot(G_sf) +
geom_sf(aes(fill = hotcold), color = "white", size = 0.2) +
scale_fill_manual(values = c(
"Hot spot (p<0.05)" = "#b2182b",
"Cold spot (p<0.05)" = "#2166ac",
"Not significant" = "grey85"
)) +
labs(title = paste("Getis–Ord Gi* —", label), fill = NULL) +
theme_minimal()
print(p1); print(p2); print(p3)
}
RUN for each variable (West Java)
analyze_spatial(jabar_merged, "Y..IPM.", "Indeks Pembangunan Manusia (IPM)")
##
##
## ==============================
## Indeks Pembangunan Manusia (IPM)
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
##
## Moran I statistic standard deviate = 2.236, p-value = 0.02535
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.27459101 -0.03846154 0.01960115
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
##
## Geary C statistic standard deviate = 2.415, p-value = 0.01573
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic Expectation Variance
## 0.6069834 1.0000000 0.0264836
analyze_spatial(jabar_merged, "X1..PDRB.", "Produk Domestik Regional Bruto (PDRB)")
##
##
## ==============================
## Produk Domestik Regional Bruto (PDRB)
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
##
## Moran I statistic standard deviate = 2.4963, p-value = 0.01255
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.29479733 -0.03846154 0.01782296
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
##
## Geary C statistic standard deviate = 1.6091, p-value = 0.1076
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic Expectation Variance
## 0.69552243 1.00000000 0.03580473
analyze_spatial(jabar_merged, "X2..Indeks.Pengetahuan.", "Indeks Pengetahuan")
##
##
## ==============================
## Indeks Pengetahuan
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
##
## Moran I statistic standard deviate = -0.090652, p-value = 0.9278
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## -0.05130804 -0.03846154 0.02008253
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
##
## Geary C statistic standard deviate = 0.21152, p-value = 0.8325
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic Expectation Variance
## 0.96725881 1.00000000 0.02396022
analyze_spatial(jabar_merged, "X3..UMR.", "Upah Minimum Rakyat (UMR)")
##
##
## ==============================
## Upah Minimum Rakyat (UMR)
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
##
## Moran I statistic standard deviate = 0.037084, p-value = 0.9704
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## -0.03318861 -0.03846154 0.02021745
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
##
## Geary C statistic standard deviate = -0.051232, p-value = 0.9591
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic Expectation Variance
## 1.00781240 1.00000000 0.02325299
analyze_spatial(jabar_merged, "X4..Pengeluaran.per.kapita.", "Pengeluaran per Kapita")
##
##
## ==============================
## Pengeluaran per Kapita
## ==============================
##
## Moran I test under randomisation
##
## data: x_raw
## weights: lwW
##
## Moran I statistic standard deviate = -1.4182, p-value = 0.1561
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## -0.22929156 -0.03846154 0.01810562
##
##
## Geary C test under randomisation
##
## data: x_raw
## weights: lwW
##
## Geary C statistic standard deviate = -0.57649, p-value = 0.5643
## alternative hypothesis: two.sided
## sample estimates:
## Geary C statistic Expectation Variance
## 1.10680275 1.00000000 0.03432302
Standar
# Extract relevant variables from your data
data_subset <- jabar_merged %>%
st_drop_geometry() %>% # Remove geometry for correlation
select(Y..IPM., X4..Pengeluaran.per.kapita.,X2..Indeks.Pengetahuan., X1..PDRB., X3..UMR.)
# Check correlation matrix
cor_matrix <- cor(data_subset, use = "complete.obs") # Use "complete.obs" to handle NA if any
print(cor_matrix)
## Y..IPM. X4..Pengeluaran.per.kapita.
## Y..IPM. 1.0000000 -0.13022832
## X4..Pengeluaran.per.kapita. -0.1302283 1.00000000
## X2..Indeks.Pengetahuan. -0.2517940 0.77408165
## X1..PDRB. 0.2780267 -0.04548669
## X3..UMR. -0.1053600 0.64015150
## X2..Indeks.Pengetahuan. X1..PDRB. X3..UMR.
## Y..IPM. -0.2517940 0.27802673 -0.1053600
## X4..Pengeluaran.per.kapita. 0.7740816 -0.04548669 0.6401515
## X2..Indeks.Pengetahuan. 1.0000000 -0.12609677 0.4428079
## X1..PDRB. -0.1260968 1.00000000 0.3119602
## X3..UMR. 0.4428079 0.31196017 1.0000000
# Install and use VIF if not already (from package 'car')
# Jika belum terinstal
library(car)
# Run a simple OLS model first to check VIF
ols_model_temp <- lm(Y..IPM. ~ X4..Pengeluaran.per.kapita. + X2..Indeks.Pengetahuan. + X1..PDRB. + X3..UMR., data = data_subset)
vif_results <- vif(ols_model_temp)
print(vif_results)
## X4..Pengeluaran.per.kapita. X2..Indeks.Pengetahuan.
## 3.544125 2.554084
## X1..PDRB. X3..UMR.
## 1.263254 2.116692
# Interpret: Jika VIF > 5 atau 10, ada multikolinieritas tinggi. Hapus variabel yang bermasalah.
# 3. Pemeriksaan dan Pembersihan Data
# Pertama, periksa apakah variabel ada dan bersihkan NA
if (!all(c("Y..IPM.", "X2..Indeks.Pengetahuan.", "X4..Pengeluaran.per.kapita.", "X3..UMR.", "X1..PDRB.") %in% names(jabar_merged))) {
stop("Salah satu variabel tidak ada di data. Periksa nama kolom!")
}
# Hapus baris dengan NA di variabel yang relevan
jabar_merged_clean <- jabar_merged %>%
dplyr::select(Y..IPM., X2..Indeks.Pengetahuan., X4..Pengeluaran.per.kapita., X3..UMR., X1..PDRB.) %>%
na.omit()
# Standardisasi data untuk mengurangi error numerik
jabar_merged_scaled <- jabar_merged_clean %>%
mutate(
Y..IPM._scaled = scale(Y..IPM.),
X2..Indeks.Pengetahuan._scaled = scale(X2..Indeks.Pengetahuan.),
X4..Pengeluaran.per.kapita._scaled = scale(X4..Pengeluaran.per.kapita.),
X3..UMR._scaled = scale(X3..UMR.),
X1..PDRB._scaled = scale(X1..PDRB.)
)
Hasil & Pembahasan
# OLS dan LMtest
model_ols <- lm(Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. + X3..UMR. + X4..Pengeluaran.per.kapita.,
data = jabar_merged)
summary(model_ols)
##
## Call:
## lm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9655 -3.0425 -0.2963 2.8723 8.7721
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.412e+01 9.123e+00 9.221 5.17e-09 ***
## X1..PDRB. 1.270e-05 8.428e-06 1.507 0.146
## X2..Indeks.Pengetahuan. -2.146e-01 1.992e-01 -1.077 0.293
## X3..UMR. -1.011e-06 1.111e-06 -0.910 0.373
## X4..Pengeluaran.per.kapita. 5.620e-04 6.657e-04 0.844 0.408
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.309 on 22 degrees of freedom
## Multiple R-squared: 0.1632, Adjusted R-squared: 0.01107
## F-statistic: 1.073 on 4 and 22 DF, p-value: 0.3937
lm_tests <- lm.LMtests(model_ols, listw = lwW, test=c("LMerr","LMlag","RLMerr","RLMlag","SARMA"))
## Please update scripts to use lm.RStests in place of lm.LMtests
print(lm_tests)
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged)
## test weights: listw
##
## RSerr = 0.42197, df = 1, p-value = 0.516
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged)
## test weights: listw
##
## RSlag = 1.9706, df = 1, p-value = 0.1604
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged)
## test weights: listw
##
## adjRSerr = 7.7336, df = 1, p-value = 0.00542
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged)
## test weights: listw
##
## adjRSlag = 9.2822, df = 1, p-value = 0.002314
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged)
## test weights: listw
##
## SARMA = 9.7042, df = 2, p-value = 0.007812
# SAR
slag_model <- lagsarlm(
Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. + X3..UMR. + X4..Pengeluaran.per.kapita.,
data = jabar_merged,
listw = lwB,
method = "eigen",
zero.policy = TRUE
)
summary(slag_model)
##
## Call:lagsarlm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged,
## listw = lwB, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.39191 -2.65992 0.41443 1.89333 6.61670
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.0270e+01 7.0072e+00 11.4553 < 2.2e-16
## X1..PDRB. 1.8285e-05 6.5632e-06 2.7861 0.005335
## X2..Indeks.Pengetahuan. 4.4210e-04 1.6402e-01 0.0027 0.997849
## X3..UMR. -1.0228e-06 8.4081e-07 -1.2164 0.223826
## X4..Pengeluaran.per.kapita. 4.3850e-05 5.2674e-04 0.0832 0.933655
##
## Rho: -0.015977, LR test value: 9.5024, p-value: 0.0020521
## Asymptotic standard error: 0.0046911
## z-value: -3.4058, p-value: 0.0006598
## Wald statistic: 11.599, p-value: 0.0006598
##
## Log likelihood: -70.23619 for lag model
## ML residual variance (sigma squared): 10.631, (sigma: 3.2606)
## Number of observations: 27
## Number of parameters estimated: 7
## AIC: 154.47, (AIC for lm: 161.97)
## LM test for residual autocorrelation
## test value: 0.35207, p-value: 0.55294
# SEM
serr_model <- errorsarlm(Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. + X3..UMR. + X4..Pengeluaran.per.kapita.,
data = jabar_merged, listw = lwW, zero.policy = TRUE)
summary(serr_model)
##
## Call:errorsarlm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged,
## listw = lwW, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0967 -3.4076 -0.4608 2.6751 6.1171
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.0578e+01 7.1084e+00 11.3355 <2e-16
## X1..PDRB. 2.2834e-06 7.0327e-06 0.3247 0.7454
## X2..Indeks.Pengetahuan. -9.3649e-02 1.5563e-01 -0.6017 0.5473
## X3..UMR. -1.1969e-06 8.7896e-07 -1.3618 0.1733
## X4..Pengeluaran.per.kapita. 4.4927e-04 5.0238e-04 0.8943 0.3712
##
## Lambda: 0.57132, LR test value: 2.1199, p-value: 0.14539
## Asymptotic standard error: 0.16647
## z-value: 3.4319, p-value: 0.00059933
## Wald statistic: 11.778, p-value: 0.00059933
##
## Log likelihood: -73.9274 for error model
## ML residual variance (sigma squared): 12.7, (sigma: 3.5637)
## Number of observations: 27
## Number of parameters estimated: 7
## AIC: 161.85, (AIC for lm: 161.97)
# SDM
sdm_model <- lagsarlm(
Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
X3..UMR. + X4..Pengeluaran.per.kapita.,
data = jabar_merged,
listw = lwB,
Durbin = TRUE,
method = "Matrix",
zero.policy = TRUE
)
## Warning in lagsarlm(Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. + X3..UMR. + : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## reciprocal condition number = 2.04162e-17 - using numerical Hessian.
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
## Warning in log(1 - coef * eig): NaNs produced
summary(sdm_model)
##
## Call:lagsarlm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged,
## listw = lwB, Durbin = TRUE, method = "Matrix", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.49104 -3.16605 -0.20043 2.92779 7.12871
##
## Type: mixed
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.2714e+01 NaN NaN NaN
## X1..PDRB. 1.9394e-05 NaN NaN NaN
## X2..Indeks.Pengetahuan. -6.1817e-02 NaN NaN NaN
## X3..UMR. -2.4115e-06 NaN NaN NaN
## X4..Pengeluaran.per.kapita. 5.0134e-04 NaN NaN NaN
## lag.(Intercept) 4.0475e+01 NaN NaN NaN
## lag.X1..PDRB. 2.5153e-05 NaN NaN NaN
## lag.X2..Indeks.Pengetahuan. -3.8203e-01 NaN NaN NaN
## lag.X3..UMR. -1.9962e-06 NaN NaN NaN
## lag.X4..Pengeluaran.per.kapita. 1.6515e-03 NaN NaN NaN
##
## Rho: -0.43346, LR test value: 20.203, p-value: 6.9646e-06
## Approximate (numerical Hessian) standard error: NaN
## z-value: NaN, p-value: NA
## Wald statistic: NaN, p-value: NA
##
## Log likelihood: -53.15572 for mixed model
## ML residual variance (sigma squared): 13.256, (sigma: 3.6409)
## Number of observations: 27
## Number of parameters estimated: 12
## AIC: 130.31, (AIC for lm: 148.51)
# SLX
slx_model <- lmSLX(
Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
X3..UMR. + X4..Pengeluaran.per.kapita.,
data = jabar_merged,
listw = lwW,
zero.policy = TRUE
)
summary(slx_model)
##
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
## data = as.data.frame(x), weights = weights)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.216e+01 1.915e+01 4.291e+00 4.394e-04
## X1..PDRB. 1.182e-05 8.079e-06 1.463e+00 1.607e-01
## X2..Indeks.Pengetahuan. -2.801e-02 1.799e-01 -1.557e-01 8.780e-01
## X3..UMR. -3.021e-06 1.111e-06 -2.718e+00 1.411e-02
## X4..Pengeluaran.per.kapita. 8.266e-04 5.898e-04 1.402e+00 1.780e-01
## lag.X1..PDRB. 6.947e-05 2.450e-05 2.836e+00 1.096e-02
## lag.X2..Indeks.Pengetahuan. -4.610e-01 3.450e-01 -1.336e+00 1.982e-01
## lag.X3..UMR. -6.151e-06 3.236e-06 -1.901e+00 7.343e-02
## lag.X4..Pengeluaran.per.kapita. 3.069e-03 1.563e-03 1.964e+00 6.519e-02
# SDEM
sdem_model <- errorsarlm(
Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
X3..UMR. + X4..Pengeluaran.per.kapita.,
data = jabar_merged,
listw = lwW,
Durbin = TRUE,
method = "eigen",
zero.policy = TRUE
)
## Warning in errorsarlm(Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. + X3..UMR. + : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## reciprocal condition number = 8.63242e-17 - using numerical Hessian.
summary(sdem_model)
##
## Call:errorsarlm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged,
## listw = lwW, Durbin = TRUE, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.167580 -2.061075 -0.066034 1.612535 6.492518
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.9439e+01 1.7502e+01 4.5388 5.657e-06
## X1..PDRB. 1.3899e-05 6.8089e-06 2.0413 0.0412210
## X2..Indeks.Pengetahuan. -3.0124e-03 1.4933e-01 -0.0202 0.9839058
## X3..UMR. -3.2525e-06 9.4477e-07 -3.4427 0.0005760
## X4..Pengeluaran.per.kapita. 8.6555e-04 4.8272e-04 1.7931 0.0729611
## lag.X1..PDRB. 7.5727e-05 2.0871e-05 3.6284 0.0002852
## lag.X2..Indeks.Pengetahuan. -4.7916e-01 2.9977e-01 -1.5984 0.1099461
## lag.X3..UMR. -7.4437e-06 2.7471e-06 -2.7097 0.0067345
## lag.X4..Pengeluaran.per.kapita. 3.6047e-03 1.2815e-03 2.8128 0.0049112
##
## Lambda: 0.35465, LR test value: 1.1129, p-value: 0.29145
## Approximate (numerical Hessian) standard error: 0.30545
## z-value: 1.1611, p-value: 0.24562
## Wald statistic: 1.348, p-value: 0.24562
##
## Log likelihood: -66.54229 for error model
## ML residual variance (sigma squared): 7.827, (sigma: 2.7977)
## Number of observations: 27
## Number of parameters estimated: 11
## AIC: 155.08, (AIC for lm: 154.2)
# SAC
sac_model <- sacsarlm(
Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
X3..UMR. + X4..Pengeluaran.per.kapita.,
data = jabar_merged,
listw = lwB,
type = "sacmixed",
method = "Matrix",
zero.policy = TRUE
)
summary(sac_model)
##
## Call:sacsarlm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged,
## listw = lwB, type = "sacmixed", method = "Matrix", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.99010 -1.50631 -0.19694 1.50011 3.25037
##
## Type: sacmixed
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.5692e+01 NaN NaN NaN
## X1..PDRB. 8.8403e-06 NaN NaN NaN
## X2..Indeks.Pengetahuan. 8.9627e-02 NaN NaN NaN
## X3..UMR. -2.0998e-06 NaN NaN NaN
## X4..Pengeluaran.per.kapita. 4.3294e-04 NaN NaN NaN
## lag.(Intercept) -2.3636e+01 NaN NaN NaN
## lag.X1..PDRB. 1.4170e-05 NaN NaN NaN
## lag.X2..Indeks.Pengetahuan. 7.1868e-02 NaN NaN NaN
## lag.X3..UMR. -1.5803e-06 NaN NaN NaN
## lag.X4..Pengeluaran.per.kapita. 2.3137e-04 NaN NaN NaN
##
## Rho: 0.25
## Approximate (numerical Hessian) standard error: NaN
## z-value: NaN, p-value: NA
## Lambda: -0.5
## Approximate (numerical Hessian) standard error: NaN
## z-value: NaN, p-value: NA
##
## LR test value: 112.23, p-value: < 2.22e-16
##
## Log likelihood: -18.87118 for sacmixed model
## ML residual variance (sigma squared): 3.4192, (sigma: 1.8491)
## Number of observations: 27
## Number of parameters estimated: 13
## AIC: 63.742, (AIC for lm: 161.97)
# GNS
gns_model <- sacsarlm(
Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
X3..UMR. + X4..Pengeluaran.per.kapita.,
data = jabar_merged,
listw = lwB,
type = "sacmixed",
Durbin = TRUE,
method = "Matrix",
zero.policy = TRUE
)
summary(gns_model)
##
## Call:sacsarlm(formula = Y..IPM. ~ X1..PDRB. + X2..Indeks.Pengetahuan. +
## X3..UMR. + X4..Pengeluaran.per.kapita., data = jabar_merged,
## listw = lwB, Durbin = TRUE, type = "sacmixed", method = "Matrix",
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.99010 -1.50631 -0.19694 1.50011 3.25037
##
## Type: sacmixed
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.5692e+01 NaN NaN NaN
## X1..PDRB. 8.8403e-06 NaN NaN NaN
## X2..Indeks.Pengetahuan. 8.9627e-02 NaN NaN NaN
## X3..UMR. -2.0998e-06 NaN NaN NaN
## X4..Pengeluaran.per.kapita. 4.3294e-04 NaN NaN NaN
## lag.(Intercept) -2.3636e+01 NaN NaN NaN
## lag.X1..PDRB. 1.4170e-05 NaN NaN NaN
## lag.X2..Indeks.Pengetahuan. 7.1868e-02 NaN NaN NaN
## lag.X3..UMR. -1.5803e-06 NaN NaN NaN
## lag.X4..Pengeluaran.per.kapita. 2.3137e-04 NaN NaN NaN
##
## Rho: 0.25
## Approximate (numerical Hessian) standard error: NaN
## z-value: NaN, p-value: NA
## Lambda: -0.5
## Approximate (numerical Hessian) standard error: NaN
## z-value: NaN, p-value: NA
##
## LR test value: 112.23, p-value: < 2.22e-16
##
## Log likelihood: -18.87118 for sacmixed model
## ML residual variance (sigma squared): 3.4192, (sigma: 1.8491)
## Number of observations: 27
## Number of parameters estimated: 13
## AIC: 63.742, (AIC for lm: 161.97)
# Compare Model
model_comparison <- data.frame(
Model = c("OLS", "SAR", "SEM", "SDM", "SLX", "SDEM", "SAC", "GNS"),
AIC = c(
AIC(model_ols),
AIC(slag_model),
AIC(serr_model),
AIC(sdm_model),
AIC(slx_model),
AIC(sdem_model),
AIC(sac_model),
AIC(gns_model)
),
LogLik = c(
logLik(model_ols),
logLik(slag_model),
logLik(serr_model),
logLik(sdm_model),
logLik(slx_model),
logLik(sdem_model),
logLik(sac_model),
logLik(gns_model)
)
)
model_comparison <- model_comparison[order(model_comparison$AIC), ]
print(model_comparison)
## Model AIC LogLik
## 7 SAC 63.74235 -18.87118
## 8 GNS 63.74235 -18.87118
## 4 SDM 130.31144 -53.15572
## 5 SLX 154.19750 -67.09875
## 2 SAR 154.47237 -70.23619
## 6 SDEM 155.08459 -66.54229
## 3 SEM 161.85480 -73.92740
## 1 OLS 161.97473 -74.98737
# Compare OLS with spatial models
spatialreg::LR.Sarlm(slag_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 9.5024, df = 1, p-value = 0.002052
## sample estimates:
## Log likelihood of slag_model Log likelihood of model_ols
## -70.23619 -74.98737
spatialreg::LR.Sarlm(serr_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 2.1199, df = 1, p-value = 0.1454
## sample estimates:
## Log likelihood of serr_model Log likelihood of model_ols
## -73.92740 -74.98737
spatialreg::LR.Sarlm(sdm_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 43.663, df = 6, p-value = 8.62e-08
## sample estimates:
## Log likelihood of sdm_model Log likelihood of model_ols
## -53.15572 -74.98737
spatialreg::LR.Sarlm(sdem_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 16.89, df = 5, p-value = 0.004713
## sample estimates:
## Log likelihood of sdem_model Log likelihood of model_ols
## -66.54229 -74.98737
spatialreg::LR.Sarlm(gns_model, model_ols)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 112.23, df = 7, p-value < 2.2e-16
## sample estimates:
## Log likelihood of gns_model Log likelihood of model_ols
## -18.87118 -74.98737
# Nested spatial versions
spatialreg::LR.Sarlm(sdm_model, slag_model)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 34.161, df = 5, p-value = 2.212e-06
## sample estimates:
## Log likelihood of sdm_model Log likelihood of slag_model
## -53.15572 -70.23619
spatialreg::LR.Sarlm(sdem_model, serr_model)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 14.77, df = 4, p-value = 0.005202
## sample estimates:
## Log likelihood of sdem_model Log likelihood of serr_model
## -66.54229 -73.92740
spatialreg::LR.Sarlm(gns_model, sdm_model)
##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 68.569, df = 1, p-value < 2.2e-16
## sample estimates:
## Log likelihood of gns_model Log likelihood of sdm_model
## -18.87118 -53.15572