---
title: "GIZI BURUK SAR BN"
output: html_notebook
---
## Data Imputation
## library
```{r}
lapply(c("googlesheets4", "tidyverse", "car", "classInt", "ggplot2", "lmtest",
"spdep", "sp", "spgwr", "glmnet", "reshape", "GWmodel", "sf", "tmap",
"fields", "dplyr", "ggplot2", "corrplot", "fitdistrplus"),
function(pkg) {
library(pkg, character.only = TRUE)
})
```
## input data
```{r}
gizbur_jateng <- read_sheet("https://docs.google.com/spreadsheets/d/1EN4HTa8sZDrxvK_YLlmuX-5NyUnl3ELHxGgpsERoyls/edit?usp=sharing")
view(gizbur_jateng)
```
## statistika deskriptif
```{r}
# Fungsi untuk menghitung statistik deskriptif termasuk standar deviasi dan CV, dibulatkan ke 2 digit desimal
descriptive_stats <- function(x) {
c(
Min = round(min(x, na.rm = TRUE), 3),
`1st Qu.` = round(quantile(x, 0.25, na.rm = TRUE), 3),
Median = round(median(x, na.rm = TRUE), 3),
Mean = round(mean(x, na.rm = TRUE), 3),
`3rd Qu.` = round(quantile(x, 0.75, na.rm = TRUE), 3),
Max = round(max(x, na.rm = TRUE), 3),
SD = round(sd(x, na.rm = TRUE), 3),
CV = round(sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE) * 100, 3) # Koefisien variasi dalam %
)
}
# Menghitung statistika deskriptif untuk setiap kolom
summary_stats <- apply(gizbur_jateng[, 2:13], 2, descriptive_stats)
# Menampilkan hasil
print(summary_stats)
```
## input peta pulau jawa
```{r}
sf_use_s2(FALSE)
petajawa <- st_read(dsn = "jawa.shp", layer = "jawa")
```
```{r}
library(sf)
petajateng <- subset(petajawa, PROVINSI == "JAWA TENGAH")
centroids <- st_centroid(petajateng)
longlat <- st_coordinates(centroids)
petajateng$long <- longlat[, 1]
petajateng$lat <- longlat[, 2]
```
```{r}
petajateng$y = gizbur_jateng$y
petajateng$x1 = gizbur_jateng$x1
petajateng$x2 = gizbur_jateng$x2
petajateng$x3 = gizbur_jateng$x3
petajateng$x4 = gizbur_jateng$x4
petajateng$x5 = gizbur_jateng$x5
petajateng$x6 = gizbur_jateng$x6
petajateng$x7 = gizbur_jateng$x7
petajateng$x8 = gizbur_jateng$x8
petajateng$x9 = gizbur_jateng$x9
petajateng$x10 = gizbur_jateng$x10
petajateng$x11 = gizbur_jateng$x11
view(petajateng)
```
```{r}
sp.label <- function(x, label) {list("sp.text", coordinates(x),
label, cex = 0.7, col = "black",
fontface = "bold")}
NUMB.sp.label <- function(x) {sp.label(x, as.vector(petajateng$KABKOT))}
make.NUMB.sp.label <- function(x) {do.call("list", NUMB.sp.label(x))}
colfunc <- colorRampPalette(c("green", "orange", "red"))
color <- colfunc(1000)
petajateng$y <- gizbur_jateng$y
tm_shape(petajateng) +
tm_polygons("y", palette = color, title = "Kasus Malnutrisi per 1000 Balita") +
tm_text("KABKOT", size = 0.4, col = "black", fontface = "bold") +
tm_layout(main.title = "Sebaran Kasus Malnutrisi di Jawa Tengah 2021",
bg.color = "white", legend.outside = TRUE)
```
## corrplot
```{r}
sp.label <- function(x, label) {list("sp.text", coordinates(x),
label, cex = 0.7, col = "black",
fontface = "bold")}
NUMB.sp.label <- function(x) {sp.label(x, as.vector(petajateng$KABKOT))}
make.NUMB.sp.label <- function(x) {do.call("list", NUMB.sp.label(x))}
colfunc <- colorRampPalette(c("red", "orange", "green"))
color <- colfunc(1000)
petajateng$x6 <- gizbur_jateng$x6
tm_shape(petajateng) +
tm_polygons("x6", palette = color, title = "per Kapita per Bulan Makanan (Rp. 100rb)") +
tm_text("KABKOT", size = 0.4, col = "black", fontface = "bold") +
tm_layout(main.title = "Sebaran Mean Expend/Kapita/Bulan Makanan",
bg.color = "white", legend.outside = TRUE)
```
```{r}
colfunc <- colorRampPalette(c("green", "orange", "red"))
color <- colfunc(1000)
petajateng$x8 <- gizbur_jateng$x8
tm_shape(petajateng) +
tm_polygons("x8", palette = color, title = "% penduduk di bawah garis miskin") +
tm_text("KABKOT", size = 0.4, col = "black", fontface = "bold") +
tm_layout(main.title = "Sebaran penduduk di bawah garis miskin Jawa Tengah 2021",
bg.color = "white", legend.outside = TRUE)
```
## corrplot
```{r}
corrplot(cor(gizbur_jateng[, 2:13]),
method = "color",
type = "upper",
tl.pos = 'tp',
tl.cex = 0.5, # Mengatur ukuran teks label
tl.col = "black", # Mengubah warna label teks menjadi hitam
tl.offset = 0.5) # Menambahkan sedikit offset pada label teks di atas
corrplot(cor(gizbur_jateng[, 2:13]),
method = "number",
diag = FALSE,
add = TRUE,
type = "lower",
tl.pos = 'n',
cl.pos = 'n',
number.cex = 0.5) # Mengatur ukuran angka korelasi
```
```{r}
colnames(gizbur_jateng)
```
## hubungan x dan y
```{r}
# Reshape the data to long format for ggplot
plot_trend <- gizbur_jateng %>%
dplyr::select(y, starts_with("x")) %>%
pivot_longer(cols = starts_with("x"), names_to = "variable", values_to = "value")
# Create a scatter plot with facets for each variable
plot_relationships <- ggplot(plot_trend, aes(x = value, y = y)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE, color = "blue") + # Add a linear regression line
facet_wrap(~ variable, scales = "free_x") + # Create facets for each variable
labs(title = "Relationship between y and Each x Variable",
x = "Predictor Variables (x)",
y = "Response Variable (y)") +
theme_minimal()
print(plot_relationships)
```
```{r}
corrplot(corr = cor(gizbur_jateng[, 2:13]), method = "color", diag = F,
outline = T, addCoef.col = "black", number.cex = 0.7,
tl.offset = 0.5, tl.cex = 0.8)
```
```{r}
boxplot(gizbur_jateng$y, ylab="Banyaknya Gizi Buruk")
```
## Grafik Cullen-Frey
```{r}
mean(gizbur_jateng$y)
var(gizbur_jateng$y)
```
```{r}
# Membuat boxplot dengan outlier yang lebih jelas
boxplot(gizbur_jateng$y,
ylab = "Banyaknya Gizi Buruk",
main = "Distribusi Gizi Buruk",
col = "lightblue", # Warna boxplot
border = "darkblue", # Warna garis boxplot
outpch = 19, # Bentuk titik outlier
outcol = "red", # Warna titik outlier
outcex = 1.5) # Ukuran titik outlier
```
```{r}
poisson <- fitdist(gizbur_jateng$y, "pois", method = c('mle', 'mme', 'qme', 'mge', 'mse'), start=NULL, fix.arg=NULL, T, keepdata=T)
binom.negatif <- fitdist(gizbur_jateng$y, "nbinom", method = c('mle', 'mme', 'qme', 'mge', 'mse'), start=NULL, fix.arg=NULL, T, keepdata=T)
descdist(gizbur_jateng$y, discrete = T, boot=1000)
```
## Q-Q Plot
```{r}
qqcomp(list(poisson, binom.negatif), legendtext = c("Pois", "Neg Binom"), fitpch = "o", fitcol = c("red", "green"))
```
## periksa multikol
```{r}
reglin <- lm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11, data=gizbur_jateng)
as.matrix(vif(reglin))
```
## uji overdispersi
```{r}
poisson <- glm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11,
data = gizbur_jateng,
family = poisson())
# nilai deviance
dv_df <- deviance(poisson)/df.residual(poisson)
dv_df
# nilai chi sq
pr_df <- sum(residuals(poisson, "pearson")^2)/df.residual(poisson)
pr_df
# overdispersi
overdis <- cbind.data.frame(dv_df, pr_df)
overdis
```
```{r}
# Membuat tabel hasil overdispersi
hasil_overdispersi <- data.frame(
Uji = c("Khi-Kuadrat Pearson", "Devians"),
`Statistik Uji` = c(sum(residuals(poisson, type = "pearson")^2), deviance(poisson)),
db = c(df.residual(poisson), df.residual(poisson)),
`Statistik uji/db` = c(pr_df, dv_df)
)
# Cetak tabel hasil
print(hasil_overdispersi)
```
```{r}
summary(poisson)
```
## negative binomial
```{r}
bn_model <- glm.nb(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11,
data = gizbur_jateng)
summary(bn_model)
```
```{r}
# Model dasar (hanya intercept) untuk Poisson
null_poisson <- glm(y ~ 1, family = poisson(link = "log"), data = gizbur_jateng)
# Model penuh (dengan 11 prediktor) untuk Poisson
full_poisson <- glm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 +
x10 + x11, family = poisson(link = "log"),
data = gizbur_jateng)
G_stat_poisson <- 2 * (logLik(poisson) - logLik(null_poisson))
cat("Nilai Statistik Uji-G untuk Poisson:", as.numeric(G_stat_poisson), "\n")
```
```{r}
# Model dasar (hanya intercept) untuk Binomial Negatif
null_negbin <- glm.nb(y ~ 1, data = gizbur_jateng)
# Model penuh (dengan 6 prediktor) untuk Binomial Negatif
full_negbin <- glm.nb(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 +
x10 + x11, data=gizbur_jateng)
# Menghitung statistik uji-G untuk Binomial Negatif
G_stat_negbin <- 2 * (logLik(full_negbin) - logLik(null_negbin))
cat("Nilai Statistik Uji-G untuk Binomial Negatif:", as.numeric(G_stat_negbin), "\n")
```
```{r}
chi_sq_critical <- qchisq(0.95, df = 11)
cat("Nilai kritis Chi-Square untuk df = 11 dan α = 0.05 adalah:", chi_sq_critical, "\n")
```
## efek spasial
# selang kelas mean
```{r}
classIntervals(gizbur_jateng$y, n=4, style = 'kmeans')
```
## Membuat peta angka Gizi Buruk
```{r}
cc <- c("green", "lightgreen", "orange", "red")
plot.gizbur = ggplot(data=petajateng) +
geom_sf(aes(fill = cut(y, breaks = c(6, 36.5, 73, 106, 133))),
lwd = 0.1, col = "black") +
scale_fill_manual("GIZI BURUK", values = cc, labels = c("6 - 36.5",
"36.5 - 73",
"73 - 106",
"106 - 133"),
guide = guide_legend(reverse = T)) +
xlab("") +
ylab("")
plot.gizbur
```
```{r}
# Membuat peta tematik dengan cluster
cc <- c("green", "lightgreen", "orange", "red")
plot.gizbur <- ggplot(data=petajateng) +
geom_sf(aes(fill = cut(y, breaks = c(6, 36.5, 73, 106, 133))),
lwd = 0.1, col = "black") +
scale_fill_manual("Malnutrisi", values = cc, labels = c("6 - 36.5", "36.5 - 73", "73 - 106", "106 - 133"),
guide = guide_legend(reverse = T)) +
xlab("") + ylab("")
# Menambahkan label nama kabupaten
plot.gizbur +
geom_sf_text(aes(label = KABKOT), size = 2, col = "black", fontface = "bold", check_overlap = TRUE)
```
## plot x2
```{r}
classIntervals(gizbur_jateng$x2, n=4, style = 'kmeans')
```
```{r}
cc <- c("red", "lightgreen", "orange", "green")
plot.asieks = ggplot(data=petajateng) +
geom_sf(aes(fill = cut(x2, breaks = c(0.1824, 0.22745, 0.2687, 0.32885, 0.3916))),
lwd = 0.1, col = "black") +
scale_fill_manual("ASI EKSKLUSIF", values = cc, labels = c("0.1924 - 0.22745",
"0.22745 - 0.2687",
"0.2687 - 0.32885",
"0.32885 - 0.3916"),
guide = guide_legend(reverse = F)) +
xlab("") +
ylab("")
plot.asieks
```
```{r}
# Konversi petajateng ke format Spatial untuk fungsi spdep
library(spdep)
petajateng_sp <- as_Spatial(petajateng)
# Buat matriks bobot Queen menggunakan poly2nb
w.queen <- poly2nb(petajateng_sp, queen = TRUE)
# Plot peta dan konektivitas Queen
plot(st_geometry(petajateng), main = "Queen Contiguity")
plot(w.queen, coords = longlat, add = TRUE, col = "red")
```
## Uji efek spasial
## Uji dependensi spasial
## menghitung jarak antar 35 kabkot
```{r}
# menghitung jarak gcd
djarak <- dist(longlat)
m.djarak <- as.matrix(djarak)
```
## Queen
```{r}
queen.nb <- poly2nb(petajateng, queen = T)
queen.w <- nb2listw(queen.nb)
queen.jateng <- queen.w
```
## ilustrasi knn
```{r}
# Membuat objek matriks jarak untuk centroid
coords <- as.matrix(longlat)
knn_2 <- knn2nb(knearneigh(coords, k = 2))
knn_4 <- knn2nb(knearneigh(coords, k = 4))
par(mfrow=c(1,2))
# Plot peta dan konektivitas KNN untuk k = 2
plot(st_geometry(petajateng), main = "K-Nearest Neighbors (k = 2)")
plot(knn_2, coords = coords, add = TRUE, col = "blue")
# Plot peta dan konektivitas KNN untuk k = 4
plot(st_geometry(petajateng), main = "K-Nearest Neighbors (k = 4)")
plot(knn_4, coords = coords, add = TRUE, col = "red")
```
```{r}
euclid <- dist(coords)
gcd <- fields::rdist.earth(coords, miles=F)
library(ecodist)
gcd <- lower(gcd)
euclid <- lower(euclid)
```
```{r}
summary(gcd)
```
```{r}
summary(euclid)
```
```{r}
# Menghitung KNN untuk k = 3 dengan jarak Euclidean
euclid_distance <- dist(coords) # Jarak Euclidean
knn_euclid <- knn2nb(knearneigh(coords, k = 3, longlat = FALSE)) # Gunakan longlat = FALSE untuk Euclidean
# Menghitung KNN untuk k = 3 dengan jarak GCD
gcd_distance <- rdist.earth(coords, miles = FALSE) # Jarak GCD (Great Circle Distance)
knn_gcd <- knn2nb(knearneigh(coords, k = 3, longlat = TRUE)) # Gunakan longlat = TRUE untuk GCD
# Plot ilustrasi KNN dengan jarak Euclidean
par(mfrow=c(1,2)) # Membuat dua plot dalam satu baris
# Plot KNN dengan jarak Euclidean
plot(st_geometry(petajateng), main = "KNN3 - Euclidean")
plot(knn_euclid, coords = coords, add = TRUE, col = "blue")
# Plot KNN dengan jarak GCD
plot(st_geometry(petajateng), main = "KNN3 - GCD")
plot(knn_gcd, coords = coords, add = TRUE, col = "red")
```
## matriks pembobot knn 1
```{r}
W.knn1 <- knn2nb(knearneigh(longlat, k=1, longlat = T))
W.knn1.s <- nb2listw(W.knn1, style = "W")
knn1.jateng <- W.knn1.s
```
## k = 2
```{r}
W.knn2 <- knn2nb(knearneigh(longlat, k=2, longlat = T))
W.knn2.s <- nb2listw(W.knn2, style = "W")
knn2.jateng <- W.knn2.s
```
## k = 3
```{r}
W.knn3 <- knn2nb(knearneigh(longlat, k=3, longlat = T))
W.knn3.s <- nb2listw(W.knn3, style = "W")
knn3.jateng <- W.knn3.s
```
## k = 4
```{r}
W.knn4 <- knn2nb(knearneigh(longlat, k=4, longlat = T))
W.knn4.s <- nb2listw(W.knn4, style = "W")
knn4.jateng <- W.knn4.s
```
## k = 5
```{r}
W.knn5 <- knn2nb(knearneigh(longlat, k=5, longlat = T))
W.knn5.s <- nb2listw(W.knn5, style = "W")
knn5.jateng <- W.knn5.s
```
## matriks pembobot idw alpha = 1
```{r}
alpha1 = 1
W.idw1 <- 1/(m.djarak^alpha1)
```
## dinormalisasi
```{r}
diag(W.idw1) <- 0
rtot1 <- rowSums(W.idw1, na.rm = T)
W.idw.sd1 <- W.idw1/rtot1
W.idw.s1 = mat2listw(W.idw.sd1, style = "W")
idw1.jateng <- W.idw.s1
```
## idw alpha = 2
```{r}
alpha2 = 2
W.idw2 <- 1/(m.djarak^alpha2)
diag(W.idw2) <- 0
rtot2 <- rowSums(W.idw2, na.rm = T)
W.idw.sd2 <- W.idw2/rtot2
rowSums(W.idw.sd2, na.rm = T)
W.idw.s2 = mat2listw(W.idw.sd2, style = "W")
idw2.jateng <- W.idw.s2
```
## matriks pembobot expo alpha = 1
```{r}
alpha3 <- 1
W.eksp1 <- exp((-alpha3)*m.djarak)
diag(W.eksp1) <- 0
rtot3 <- rowSums(W.eksp1, na.rm = T)
W.eksp.sd1 <- W.eksp1/rtot3
W.eksp.s1 = mat2listw(W.eksp.sd1, style = "W")
eks1.jateng <- W.eksp.s1
```
## matriks pembobot expo alpha = 2
```{r}
alpha4 <- 2
W.eksp2 <- exp((-alpha4)*m.djarak)
diag(W.eksp2) <- 0
rtot4 <- rowSums(W.eksp2, na.rm = T)
W.eksp.sd2 <- W.eksp2/rtot4
W.eksp.s2 = mat2listw(W.eksp.sd2, style = "W")
eks2.jateng <- W.eksp.s2
```
# indeks moran untuk X1
```{r}
x1.k1 <- moran.test(petajateng$x1, knn1.jateng, randomisation=T)
x1.k2 <- moran.test(petajateng$x1, knn2.jateng, randomisation=T)
x1.k3 <- moran.test(petajateng$x1, knn3.jateng, randomisation=T)
x1.k4 <- moran.test(petajateng$x1, knn4.jateng, randomisation=T)
x1.k5 <- moran.test(petajateng$x1, knn5.jateng, randomisation=T)
x1.eks1 <- moran.test(petajateng$x1, eks1.jateng, randomisation = T)
x1.eks2 <- moran.test(petajateng$x1, eks2.jateng, randomisation = T)
x1.idw1 <- moran.test(petajateng$x1, idw1.jateng, randomisation = T)
x1.idw2 <- moran.test(petajateng$x1, idw2.jateng, randomisation = T)
x1.queen <- moran.test(petajateng$x1, queen.jateng, randomisation = T)
tab <- matrix(c(x1.k1$statistic, x1.k2$statistic, x1.k3$statistic,
x1.k4$statistic, x1.k5$statistic, x1.eks1$statistic,
x1.eks2$statistic, x1.idw1$statistic, x1.idw2$statistic,
x1.queen$statistic,
x1.k1$p.value, x1.k2$p.value, x1.k3$p.value, x1.k4$p.value,
x1.k5$p.value, x1.eks1$p.value, x1.eks2$p.value,
x1.idw1$p.value, x1.idw2$p.value, x1.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X1\n")
tab
```
## indeks moran untuk X2
```{r}
x2.k1 <- moran.test(petajateng$x2, knn1.jateng, randomisation=T)
x2.k2 <- moran.test(petajateng$x2, knn2.jateng, randomisation=T)
x2.k3 <- moran.test(petajateng$x2, knn3.jateng, randomisation=T)
x2.k4 <- moran.test(petajateng$x2, knn4.jateng, randomisation=T)
x2.k5 <- moran.test(petajateng$x2, knn5.jateng, randomisation=T)
x2.eks1 <- moran.test(petajateng$x2, eks1.jateng, randomisation = T)
x2.eks2 <- moran.test(petajateng$x2, eks2.jateng, randomisation = T)
x2.idw1 <- moran.test(petajateng$x2, idw1.jateng, randomisation = T)
x2.idw2 <- moran.test(petajateng$x2, idw2.jateng, randomisation = T)
x2.queen <- moran.test(petajateng$x2, queen.jateng, randomisation = T)
tab <- matrix(c(x2.k1$statistic, x2.k2$statistic, x2.k3$statistic,
x2.k4$statistic, x2.k5$statistic, x2.eks1$statistic,
x2.eks2$statistic, x2.idw1$statistic, x2.idw2$statistic,
x2.queen$statistic,
x2.k1$p.value, x2.k2$p.value, x2.k3$p.value, x2.k4$p.value,
x2.k5$p.value, x2.eks1$p.value, x2.eks2$p.value,
x2.idw1$p.value, x2.idw2$p.value, x2.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X2\n")
tab
```
## indeks moran untuk X3
```{r}
x3.k1 <- moran.test(petajateng$x3, knn1.jateng, randomisation=T)
x3.k2 <- moran.test(petajateng$x3, knn2.jateng, randomisation=T)
x3.k3 <- moran.test(petajateng$x3, knn3.jateng, randomisation=T)
x3.k4 <- moran.test(petajateng$x3, knn4.jateng, randomisation=T)
x3.k5 <- moran.test(petajateng$x3, knn5.jateng, randomisation=T)
x3.eks1 <- moran.test(petajateng$x3, eks1.jateng, randomisation = T)
x3.eks2 <- moran.test(petajateng$x3, eks2.jateng, randomisation = T)
x3.idw1 <- moran.test(petajateng$x3, idw1.jateng, randomisation = T)
x3.idw2 <- moran.test(petajateng$x3, idw2.jateng, randomisation = T)
x3.queen <- moran.test(petajateng$x3, queen.jateng, randomisation = T)
tab <- matrix(c(x3.k1$statistic, x3.k2$statistic, x3.k3$statistic,
x3.k4$statistic, x3.k5$statistic, x3.eks1$statistic,
x3.eks2$statistic, x3.idw1$statistic, x3.idw2$statistic,
x3.queen$statistic,
x3.k1$p.value, x3.k2$p.value, x3.k3$p.value, x3.k4$p.value,
x3.k5$p.value, x3.eks1$p.value, x3.eks2$p.value,
x3.idw1$p.value, x3.idw2$p.value, x3.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X3\n")
tab
```
## Indeks Moran Untuk X4
```{r}
x4.k1 <- moran.test(petajateng$x4, knn1.jateng, randomisation=T)
x4.k2 <- moran.test(petajateng$x4, knn2.jateng, randomisation=T)
x4.k3 <- moran.test(petajateng$x4, knn3.jateng, randomisation=T)
x4.k4 <- moran.test(petajateng$x4, knn4.jateng, randomisation=T)
x4.k5 <- moran.test(petajateng$x4, knn5.jateng, randomisation=T)
x4.eks1 <- moran.test(petajateng$x4, eks1.jateng, randomisation = T)
x4.eks2 <- moran.test(petajateng$x4, eks2.jateng, randomisation = T)
x4.idw1 <- moran.test(petajateng$x4, idw1.jateng, randomisation = T)
x4.idw2 <- moran.test(petajateng$x4, idw2.jateng, randomisation = T)
x4.queen <- moran.test(petajateng$x4, queen.jateng, randomisation = T)
tab <- matrix(c(x4.k1$statistic, x4.k2$statistic, x4.k3$statistic,
x4.k4$statistic, x4.k5$statistic, x4.eks1$statistic,
x4.eks2$statistic, x4.idw1$statistic, x4.idw2$statistic,
x4.queen$statistic,
x4.k1$p.value, x4.k2$p.value, x4.k3$p.value, x4.k4$p.value,
x4.k5$p.value, x4.eks1$p.value, x4.eks2$p.value,
x4.idw1$p.value, x4.idw2$p.value, x4.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X4\n")
tab
```
## Indeks Moran Untuk X5
```{r}
x5.k1 <- moran.test(petajateng$x5, knn1.jateng, randomisation=T)
x5.k2 <- moran.test(petajateng$x5, knn2.jateng, randomisation=T)
x5.k3 <- moran.test(petajateng$x5, knn3.jateng, randomisation=T)
x5.k4 <- moran.test(petajateng$x5, knn4.jateng, randomisation=T)
x5.k5 <- moran.test(petajateng$x5, knn5.jateng, randomisation=T)
x5.eks1 <- moran.test(petajateng$x5, eks1.jateng, randomisation = T)
x5.eks2 <- moran.test(petajateng$x5, eks2.jateng, randomisation = T)
x5.idw1 <- moran.test(petajateng$x5, idw1.jateng, randomisation = T)
x5.idw2 <- moran.test(petajateng$x5, idw2.jateng, randomisation = T)
x5.queen <- moran.test(petajateng$x5, queen.jateng, randomisation = T)
tab <- matrix(c(x5.k1$statistic, x5.k2$statistic, x5.k3$statistic,
x5.k4$statistic, x5.k5$statistic, x5.eks1$statistic,
x5.eks2$statistic, x5.idw1$statistic, x5.idw2$statistic,
x5.queen$statistic,
x5.k1$p.value, x5.k2$p.value, x5.k3$p.value, x5.k4$p.value,
x5.k5$p.value, x5.eks1$p.value, x5.eks2$p.value,
x5.idw1$p.value, x5.idw2$p.value, x5.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X5\n")
tab
```
## Indeks Moran Untuk X6
```{r}
x6.k1 <- moran.test(petajateng$x6, knn1.jateng, randomisation=T)
x6.k2 <- moran.test(petajateng$x6, knn2.jateng, randomisation=T)
x6.k3 <- moran.test(petajateng$x6, knn3.jateng, randomisation=T)
x6.k4 <- moran.test(petajateng$x6, knn4.jateng, randomisation=T)
x6.k5 <- moran.test(petajateng$x6, knn5.jateng, randomisation=T)
x6.eks1 <- moran.test(petajateng$x6, eks1.jateng, randomisation = T)
x6.eks2 <- moran.test(petajateng$x6, eks2.jateng, randomisation = T)
x6.idw1 <- moran.test(petajateng$x6, idw1.jateng, randomisation = T)
x6.idw2 <- moran.test(petajateng$x6, idw2.jateng, randomisation = T)
x6.queen <- moran.test(petajateng$x6, queen.jateng, randomisation = T)
tab <- matrix(c(x6.k1$statistic, x6.k2$statistic, x6.k3$statistic,
x6.k4$statistic, x6.k5$statistic, x6.eks1$statistic,
x6.eks2$statistic, x6.idw1$statistic, x6.idw2$statistic,
x6.queen$statistic,
x6.k1$p.value, x6.k2$p.value, x6.k3$p.value, x6.k4$p.value,
x6.k5$p.value, x6.eks1$p.value, x6.eks2$p.value,
x6.idw1$p.value, x6.idw2$p.value, x6.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X6\n")
tab
```
## Indeks Moran Untuk X7
```{r}
x7.k1 <- moran.test(petajateng$x7, knn1.jateng, randomisation=T)
x7.k2 <- moran.test(petajateng$x7, knn2.jateng, randomisation=T)
x7.k3 <- moran.test(petajateng$x7, knn3.jateng, randomisation=T)
x7.k4 <- moran.test(petajateng$x7, knn4.jateng, randomisation=T)
x7.k5 <- moran.test(petajateng$x7, knn5.jateng, randomisation=T)
x7.eks1 <- moran.test(petajateng$x7, eks1.jateng, randomisation = T)
x7.eks2 <- moran.test(petajateng$x7, eks2.jateng, randomisation = T)
x7.idw1 <- moran.test(petajateng$x7, idw1.jateng, randomisation = T)
x7.idw2 <- moran.test(petajateng$x7, idw2.jateng, randomisation = T)
x7.queen <- moran.test(petajateng$x7, queen.jateng, randomisation = T)
tab <- matrix(c(x7.k1$statistic, x7.k2$statistic, x7.k3$statistic,
x7.k4$statistic, x7.k5$statistic, x7.eks1$statistic,
x7.eks2$statistic, x7.idw1$statistic, x7.idw2$statistic,
x7.queen$statistic,
x7.k1$p.value, x7.k2$p.value, x7.k3$p.value, x7.k4$p.value,
x7.k5$p.value, x7.eks1$p.value, x7.eks2$p.value,
x7.idw1$p.value, x7.idw2$p.value, x7.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X7\n")
tab
```
## Indeks Moran Untuk X8
```{r}
x8.k1 <- moran.test(petajateng$x8, knn1.jateng, randomisation=T)
x8.k2 <- moran.test(petajateng$x8, knn2.jateng, randomisation=T)
x8.k3 <- moran.test(petajateng$x8, knn3.jateng, randomisation=T)
x8.k4 <- moran.test(petajateng$x8, knn4.jateng, randomisation=T)
x8.k5 <- moran.test(petajateng$x8, knn5.jateng, randomisation=T)
x8.eks1 <- moran.test(petajateng$x8, eks1.jateng, randomisation = T)
x8.eks2 <- moran.test(petajateng$x8, eks2.jateng, randomisation = T)
x8.idw1 <- moran.test(petajateng$x8, idw1.jateng, randomisation = T)
x8.idw2 <- moran.test(petajateng$x8, idw2.jateng, randomisation = T)
x8.queen <- moran.test(petajateng$x8, queen.jateng, randomisation = T)
tab <- matrix(c(x8.k1$statistic, x8.k2$statistic, x8.k3$statistic,
x8.k4$statistic, x8.k5$statistic, x8.eks1$statistic,
x8.eks2$statistic, x8.idw1$statistic, x8.idw2$statistic,
x8.queen$statistic,
x8.k1$p.value, x8.k2$p.value, x8.k3$p.value, x8.k4$p.value,
x8.k5$p.value, x8.eks1$p.value, x8.eks2$p.value,
x8.idw1$p.value, x8.idw2$p.value, x8.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X8\n")
tab
```
## Indeks Moran Untuk X9
```{r}
x9.k1 <- moran.test(petajateng$x9, knn1.jateng, randomisation=T)
x9.k2 <- moran.test(petajateng$x9, knn2.jateng, randomisation=T)
x9.k3 <- moran.test(petajateng$x9, knn3.jateng, randomisation=T)
x9.k4 <- moran.test(petajateng$x9, knn4.jateng, randomisation=T)
x9.k5 <- moran.test(petajateng$x9, knn5.jateng, randomisation=T)
x9.eks1 <- moran.test(petajateng$x9, eks1.jateng, randomisation = T)
x9.eks2 <- moran.test(petajateng$x9, eks2.jateng, randomisation = T)
x9.idw1 <- moran.test(petajateng$x9, idw1.jateng, randomisation = T)
x9.idw2 <- moran.test(petajateng$x9, idw2.jateng, randomisation = T)
x9.queen <- moran.test(petajateng$x9, queen.jateng, randomisation = T)
tab <- matrix(c(x9.k1$statistic, x9.k2$statistic, x9.k3$statistic,
x9.k4$statistic, x9.k5$statistic, x9.eks1$statistic,
x9.eks2$statistic, x9.idw1$statistic, x9.idw2$statistic,
x9.queen$statistic,
x9.k1$p.value, x9.k2$p.value, x9.k3$p.value, x9.k4$p.value,
x9.k5$p.value, x9.eks1$p.value, x9.eks2$p.value,
x9.idw1$p.value, x9.idw2$p.value, x9.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X9\n")
tab
```
## Indeks Moran Untuk X10
```{r}
x10.k1 <- moran.test(petajateng$x10, knn1.jateng, randomisation=T)
x10.k2 <- moran.test(petajateng$x10, knn2.jateng, randomisation=T)
x10.k3 <- moran.test(petajateng$x10, knn3.jateng, randomisation=T)
x10.k4 <- moran.test(petajateng$x10, knn4.jateng, randomisation=T)
x10.k5 <- moran.test(petajateng$x10, knn5.jateng, randomisation=T)
x10.eks1 <- moran.test(petajateng$x10, eks1.jateng, randomisation = T)
x10.eks2 <- moran.test(petajateng$x10, eks2.jateng, randomisation = T)
x10.idw1 <- moran.test(petajateng$x10, idw1.jateng, randomisation = T)
x10.idw2 <- moran.test(petajateng$x10, idw2.jateng, randomisation = T)
x10.queen <- moran.test(petajateng$x10, queen.jateng, randomisation = T)
tab <- matrix(c(x10.k1$statistic, x10.k2$statistic, x10.k3$statistic,
x10.k4$statistic, x10.k5$statistic, x10.eks1$statistic,
x10.eks2$statistic, x10.idw1$statistic, x10.idw2$statistic,
x10.queen$statistic,
x10.k1$p.value, x10.k2$p.value, x10.k3$p.value, x10.k4$p.value,
x10.k5$p.value, x10.eks1$p.value, x10.eks2$p.value,
x10.idw1$p.value, x10.idw2$p.value, x10.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X10\n")
tab
```
## Indeks Moran Untuk X11
```{r}
x11.k1 <- moran.test(petajateng$x11, knn1.jateng, randomisation=T)
x11.k2 <- moran.test(petajateng$x11, knn2.jateng, randomisation=T)
x11.k3 <- moran.test(petajateng$x11, knn3.jateng, randomisation=T)
x11.k4 <- moran.test(petajateng$x11, knn4.jateng, randomisation=T)
x11.k5 <- moran.test(petajateng$x11, knn5.jateng, randomisation=T)
x11.eks1 <- moran.test(petajateng$x11, eks1.jateng, randomisation = T)
x11.eks2 <- moran.test(petajateng$x11, eks2.jateng, randomisation = T)
x11.idw1 <- moran.test(petajateng$x11, idw1.jateng, randomisation = T)
x11.idw2 <- moran.test(petajateng$x11, idw2.jateng, randomisation = T)
x11.queen <- moran.test(petajateng$x11, queen.jateng, randomisation = T)
tab <- matrix(c(x11.k1$statistic, x11.k2$statistic, x11.k3$statistic,
x11.k4$statistic, x11.k5$statistic, x11.eks1$statistic,
x11.eks2$statistic, x11.idw1$statistic, x11.idw2$statistic,
x11.queen$statistic,
x11.k1$p.value, x11.k2$p.value, x11.k3$p.value, x11.k4$p.value,
x11.k5$p.value, x11.eks1$p.value, x11.eks2$p.value,
x11.idw1$p.value, x11.idw2$p.value, x11.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk X11\n")
tab
```
## Indeks Moran untuk Y
```{r}
y.k1 <- moran.test(petajateng$y, knn1.jateng, randomisation=T)
y.k2 <- moran.test(petajateng$y, knn2.jateng, randomisation=T)
y.k3 <- moran.test(petajateng$y, knn3.jateng, randomisation=T)
y.k4 <- moran.test(petajateng$y, knn4.jateng, randomisation=T)
y.k5 <- moran.test(petajateng$y, knn5.jateng, randomisation=T)
y.eks1 <- moran.test(petajateng$y, eks1.jateng, randomisation = T)
y.eks2 <- moran.test(petajateng$y, eks2.jateng, randomisation = T)
y.idw1 <- moran.test(petajateng$y, idw1.jateng, randomisation = T)
y.idw2 <- moran.test(petajateng$y, idw2.jateng, randomisation = T)
y.queen <- moran.test(petajateng$y, queen.jateng, randomisation = T)
tab <- matrix(c(y.k1$statistic, y.k2$statistic, y.k3$statistic,
y.k4$statistic, y.k5$statistic, y.eks1$statistic,
y.eks2$statistic, y.idw1$statistic, y.idw2$statistic,
y.queen$statistic,
y.k1$p.value, y.k2$p.value, y.k3$p.value, y.k4$p.value,
y.k5$p.value, y.eks1$p.value, y.eks2$p.value,
y.idw1$p.value, y.idw2$p.value, y.queen$p.value),
ncol = 2, byrow = F)
colnames(tab) <- c("Indeks Moran", "P-value")
rownames(tab) <- c("1-NN", "2-NN", "3-NN", "4-NN", "5-NN", "Eks a=1",
"Eks a=2", "IDW a=1", "IDW a=2", "Queen")
tab <- as.table(tab)
cat("Indeks Moran Untuk Y\n")
tab
```
## Moran Plot
```{r}
moran.plot(gizbur_jateng$y, queen.w, labels = gizbur_jateng$KABKOT,
xlim = c(min(gizbur_jateng$y) - 12, max(gizbur_jateng$y) + 12),
ylim = c(min(lag.listw(queen.w, gizbur_jateng$y)) - 10,
max(lag.listw(queen.w, gizbur_jateng$y)) + 10))
```
## Geary C asumsi acak
```{r}
geary.test(petajateng$y, queen.w)
```
```{r}
geary.mc(petajateng$y, queen.w, nsim = 500)
```
## Getis-Ord G
```{r}
globalG.test(petajateng$y, listw = nb2listw(queen.nb, style = "B"))
```
```{r}
coba <- globalG.test(petajateng$y, listw = nb2listw(queen.nb, style = "B"))
str(coba)
```
```{r}
coba$estimate[1]
```
## Autokorelasi spasial lokal (LISA)
```{r}
localmoran(petajateng$y, listw = nb2listw(queen.nb))
```
## peta LISA (Z score)
```{r}
lisa <- localmoran(petajateng$y, listw = nb2listw(queen.nb))
lm.palette <- colorRampPalette(c("yellow", "orange", "red"),
space = "rgb")
petajateng$z.li <- lisa[, 4]
petajateng$pvalue <- lisa[, 5]
ggplot(data = petajateng) +
geom_sf(aes(fill = z.li), color = NA) + # Menggunakan z.li untuk warna
scale_fill_gradientn(colors = lm.palette(20), na.value = "grey50") +
labs(title = "Local Moran's I", fill = "Local I") +
theme_minimal()
```
## Getis-Ord Gi (lokal)
```{r}
localG(petajateng$y, queen.w)
```
## skor Getis-Ord
```{r}
localg <- localG(petajateng$y, queen.w)
petajateng$localg <- as.numeric(localg)
lm.palette <- colorRampPalette(c("yellow", "orange", "red"),
space = "rgb")
ggplot(data = petajateng) +
geom_sf(aes(fill = localg), color = NA) + # Menggunakan localg untuk warna
scale_fill_gradientn(colors = lm.palette(20), na.value = "grey50") +
labs(title = "Local Gi", fill = "Local Gi") +
theme_minimal()
```
```{r}
petajateng$cluster <- NA
petajateng$cluster[lisa[, 1] > 0 & gizbur_jateng$y > mean(gizbur_jateng$y)] <- "High-High"
petajateng$cluster[lisa[, 1] > 0 & gizbur_jateng$y <= mean(gizbur_jateng$y)] <- "Low-Low"
petajateng$cluster[lisa[, 1] < 0 & gizbur_jateng$y > mean(gizbur_jateng$y)] <- "High-Low"
petajateng$cluster[lisa[, 1] < 0 & gizbur_jateng$y <= mean(gizbur_jateng$y)] <- "Low-High"
# Mengubah cluster menjadi faktor untuk visualisasi
petajateng$cluster <- factor(petajateng$cluster, levels = c("High-High", "High-Low", "Low-High", "Low-Low"))
```
```{r}
ggplot(data = petajateng) +
geom_sf(aes(fill = cluster), color = "black", lwd = 0.1) +
scale_fill_manual(values = c("High-High" = "red",
"Low-Low" = "blue",
"High-Low" = "pink",
"Low-High" = "lightblue",
"Tidak Signifikan" = "gray80"),
na.value = "white") +
labs(fill = "Cluster Moran's I") +
theme_minimal() +
theme(legend.position = "right") +
ggtitle("Moran Plot Cluster Map of Malnutrition in Jawa Tengah") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_sf_text(aes(label = KABKOT), size = 1.5, color = "black")
```
## uji LM Eksponensial
```{r}
reglinlog <- lm(log(y)~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11, data = gizbur_jateng)
LM_eksp <- lm.RStests(reglinlog, eks2.jateng,
test = c("LMlag", "LMerr", "RLMlag", "RLMerr", "SARMA"),
zero.policy = T)
summary(LM_eksp)
```
## Uji LM KNN
```{r}
LM_knn <- lm.RStests(reglinlog, knn5.jateng,
test = c("LMlag", "LMerr", "RLMlag", "RLMerr", "SARMA"),
zero.policy = TRUE)
summary(LM_knn)
```
## Uji LM IDW
```{r}
LM_idw <- lm.RStests(reglinlog, idw2.jateng,
test = c("LMlag", "LMerr", "RLMlag", "RLMerr", "SARMA"),
zero.policy = TRUE)
summary(LM_idw)
```
## Uji LM Queen
```{r}
LM_queen <- lm.RStests(reglinlog, queen.jateng,
test = c("LMlag", "LMerr", "RLMlag", "RLMerr", "SARMA"),
zero.policy = TRUE)
summary(LM_queen)
```
```{r}
bptest(reglinlog)
```
## pendugaan model SAR BN
```{r}
petajateng$lagY <- lag.listw(x=knn5.jateng, var = petajateng$y)
sar_model <- glm.nb(petajateng$y ~ lagY + x1 + x2 + x3 + x4 + x5 + x6 + x7 +
x8 + x9 + x10 + x11,
data = petajateng)
summary(sar_model, Nagelkerke=T)
```
```{r}
class(W.eksp.s1)
class(W.idw.s1)
W.eksp.s1
```
# SLX-BN
```{r}
library(spatialreg)
var.x5 <- as.matrix(gizbur_jateng[,7])
var.x8 <- as.matrix(gizbur_jateng[,10])
var.x9 <- as.matrix(gizbur_jateng[,11])
var.x10 <- as.matrix(gizbur_jateng[,12])
# Konversi W.eksp.s1 dan W.idw.s1 ke matriks
eks1.mat <- as.matrix(as_dgRMatrix_listw(W.eksp.s1))
idw1.mat <- as.matrix(as_dgRMatrix_listw(W.idw.s1))
W.knn1.mat <- nb2mat(W.knn1)
W.knn4.mat <- nb2mat(W.knn4)
petajateng$lag.x5 <- eks1.mat%*%var.x5
petajateng$lag.x8 <- W.knn1.mat%*%var.x8
petajateng$lag.x9 <- W.knn4.mat%*%var.x9
petajateng$lag.x10 <- idw1.mat%*%var.x10
slx_model <- glm.nb(y ~x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 +
lag.x5 + lag.x8 + lag.x9 + lag.x10, data = petajateng)
summary(slx_model, Nagelkerke=T)
```
## SDM
```{r}
sdm_model <- glm.nb(y ~x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 +
lagY + lag.x5 + lag.x8 + lag.x9 + lag.x10,
data = petajateng)
summary(sdm_model, Nagelkerke=T)
```
## SEM
```{r}
residuals_bn <- residuals(bn_model, type="response")
# Menghitung lag dari residuals
lag_residuals <- lag.listw(x=queen.jateng, var = residuals_bn)
# Membangun model SEM dengan residuals
sem_bn_model <- glm.nb(petajateng$y ~ lag_residuals + x1 + x2 + x3 + x4 + x5 +
x6 + x7 + x8 + x9 + x10 + x11,
data = petajateng)
# Menampilkan ringkasan model
summary(sem_bn_model)
```
```{r}
# Menggunakan SEM pada residuals
sem_model1 <- errorsarlm(residuals_bn ~ x1 + x2 + x3 + x4 + x5 +
x6 + x7 + x8 + x9 + x10 + x11,
listw = queen.jateng,
data = petajateng)
# Memasukkan hasil residual SEM ke dalam model Negatif Binomial
lag_residuals1 <- residuals(sem_model1, type = "response")
sem_bn_model1 <- glm.nb(petajateng$y ~ lag_residuals1 + x1 + x2 + x3 + x4 + x5 +
x6 + x7 + x8 + x9 + x10 + x11,
data = petajateng)
# Ringkasan model
summary(sem_bn_model1)
```
## Sarma
```{r}
# Tahap 2: Model pertama untuk menghitung residuals awal
sarma_initial <- sar_model
# Mendapatkan residual dari model awal
residuals_initial <- residuals(sarma_initial, type = "response")
# Tahap 3: Membuat lag dari residuals untuk MA (Moving Average)
petajateng$lagResiduals <- lag.listw(queen.jateng, residuals_initial)
# Tahap 4: Model SARMA-BN, memasukkan lag variabel dependen (lagY) dan lag residual (lagResiduals)
sarma_bn_model <- glm.nb(y ~ lagY + lagResiduals + x1 + x2 + x3 + x4 + x5 +
x6 + x7 + x8 + x9 + x10 + x11,
data = petajateng)
# Menampilkan ringkasan model
summary(sarma_bn_model)
```
## goodness of fit
```{r}
AIC.BN <- AIC(bn_model)
AIC.SAR <- AIC(sar_model)
AIC.SLX <- AIC(slx_model)
AIC.SDM <- AIC(sdm_model)
AIC.SEM <- AIC(sem_bn_model)
AIC.SARMA <- AIC(sarma_bn_model)
nilai.aic <- matrix(c(AIC.BN, AIC.SAR, AIC.SLX, AIC.SDM,
AIC.SEM, AIC.SARMA),
ncol = 1,
byrow = F)
colnames(nilai.aic) <- c("AIC")
rownames(nilai.aic) <- c("Reg Bin.Negatif", "SAR-BN", "SLX-BN",
"SDM-BN", "SEM-BN", "SARMA-BN")
nilai.aic <- as.table(nilai.aic)
# Mengurutkan tabel berdasarkan nilai AIC
nilai.aic.sorted <- nilai.aic[order(nilai.aic[, "AIC"]), , drop = FALSE]
nilai.aic.sorted
```
## asumsi residual
## sar
```{r}
sisaansar <- resid(sar_model)
normsar <- shapiro.test(sisaansar)
hetesar <- bptest(sar_model)
autosar <- moran.test(sisaansar, listw = W.knn5.s)
asumsisar <- data.frame(matrix(c("Kenormalan",
"Kehomogenan Ragam",
"Autokorelasi",
round(normsar$p.value, 2),
round(hetesar$p.value, 2),
round(autosar$p.value, 2)),
byrow = F, ncol = 2,
dimnames = list(1:3, c("Uji", "p-value"))))
asumsisar$Kesimpulan <- ifelse(asumsisar$p.value > 0.05,
"Terpenuhi", "Tidak Terpnuhi")
asumsisar
```
## residual sdm
```{r}
sisaansdm <- resid(sdm_model)
normsdm <- shapiro.test(sisaansdm)
hetesdm <- bptest(sdm_model)
autosdm <- moran.test(sisaansdm, listw = W.knn5.s)
asumsisdm <- data.frame(matrix(c("Kenormalan",
"Kehomogenan Ragam",
"Autokorelasi",
round(normsdm$p.value, 2),
round(hetesdm$p.value, 2),
round(autosdm$p.value, 2)),
byrow = F, ncol = 2,
dimnames = list(1:3, c("Uji", "p-value"))))
asumsisdm$Kesimpulan <- ifelse(asumsisdm$p.value > 0.05,
"Terpenuhi", "Tidak Terpnuhi")
asumsisdm
```
## residual sem
```{r}
sisaansem <- resid(sem_bn_model)
normsem <- shapiro.test(sisaansem)
hetesem <- bptest(sem_bn_model)
autosem <- moran.test(sisaansem, listw = queen.jateng)
asumsisem <- data.frame(matrix(c("Kenormalan",
"Kehomogenan Ragam",
"Autokorelasi",
round(normsem$p.value, 2),
round(hetesem$p.value, 2),
round(autosem$p.value, 2)),
byrow = F, ncol = 2,
dimnames = list(1:3, c("Uji", "p-value"))))
asumsisem$Kesimpulan <- ifelse(asumsisem$p.value >= 0.05,
"Terpenuhi", "Tidak Terpnuhi")
asumsisem
```
```{r}
sisaanslx <- resid(slx_model)
normslx <- shapiro.test(sisaanslx)
heteslx <- bptest(slx_model)
autoslx <- moran.test(sisaanslx, listw = knn5.jateng)
asumsislx <- data.frame(matrix(c("Kenormalan",
"Kehomogenan Ragam",
"Autokorelasi",
round(normslx$p.value, 2),
round(heteslx$p.value, 2),
round(autoslx$p.value, 2)),
byrow = F, ncol = 2,
dimnames = list(1:3, c("Uji", "p-value"))))
asumsislx$Kesimpulan <- ifelse(asumsislx$p.value > 0.05,
"Terpenuhi", "Tidak Terpnuhi")
asumsislx
```
```{r}
sisaansarma <- resid(sarma_bn_model)
normsarma <- shapiro.test(sisaansarma)
hetesarma <- bptest(sarma_bn_model)
autosarma <- moran.test(sisaansarma, listw = queen.jateng)
asumsisarma <- data.frame(matrix(c("Kenormalan",
"Kehomogenan Ragam",
"Autokorelasi",
round(normsarma$p.value, 2),
round(hetesarma$p.value, 2),
round(autosarma$p.value, 2)),
byrow = F, ncol = 2,
dimnames = list(1:3, c("Uji", "p-value"))))
asumsisarma$Kesimpulan <- ifelse(asumsisarma$p.value > 0.05,
"Terpenuhi", "Tidak Terpnuhi")
asumsisarma
```
```{r}
sisaanbn <- resid(bn_model)
normbn <- shapiro.test(sisaanbn)
hetebn <- bptest(bn_model)
autobn <- moran.test(sisaanbn, listw = queen.jateng)
asumsibn <- data.frame(matrix(c("Kenormalan",
"Kehomogenan Ragam",
"Autokorelasi",
round(normbn$p.value, 2),
round(hetebn$p.value, 2),
round(autobn$p.value, 2)),
byrow = F, ncol = 2,
dimnames = list(1:3, c("Uji", "p-value"))))
asumsibn$Kesimpulan <- ifelse(asumsibn$p.value > 0.05,
"Terpenuhi", "Tidak Terpnuhi")
asumsibn
```
## efek marjinal
## interpretasi koefisien
## matrix idntitas
```{r}
I <- matrix(0, 35, 35)
diag(I) <- 1
It <- matrix(1, 35, 1)
It_transpose <- t(It)
```
```{r}
queen <- nb2mat(queen.nb)
W <- as.matrix(queen)
rho <- sarma_bn_model$coefficients[2]
rhow <- rho*W
op <- I - rhow
invers <- solve(op)
```
## x1
```{r}
beta1 <- sarma_bn_model$coefficients[4]
b1 <- beta1*invers
de_x1 <- mean(diag(b1))
te_x1 <- (1/35)*It_transpose%*%b1%*%It
ie_x1 <- te_x1-de_x1
```
## x2
```{r}
beta2 <- sarma_bn_model$coefficients[5]
b2 <- beta2*invers
de_x2 <- mean(diag(b2))
te_x2 <- (1/35)*It_transpose%*%b2%*%It
ie_x2 <- te_x2-de_x2
```
## x4
```{r}
beta4 <- sarma_bn_model$coefficients[7]
b4 <- beta4*invers
de_x4 <- mean(diag(b4))
te_x4 <- (1/35)*It_transpose%*%b4%*%It
ie_x4 <- te_x4-de_x4
```
## x5
```{r}
beta5 <- sarma_bn_model$coefficients[8]
b5 <- beta5*invers
de_x5 <- mean(diag(b5))
te_x5 <- (1/35)*It_transpose%*%b5%*%It
ie_x5 <- te_x5-de_x5
```
## x7
```{r}
beta7 <- sarma_bn_model$coefficients[10]
b7 <- beta7*invers
de_x7 <- mean(diag(b7))
te_x7 <- (1/35)*It_transpose%*%b7%*%It
ie_x7 <- te_x7-de_x7
```
## x8
```{r}
beta8 <- sarma_bn_model$coefficients[11]
b8 <- beta8*invers
de_x8 <- mean(diag(b8))
te_x8 <- (1/35)*It_transpose%*%b8%*%It
ie_x8 <- te_x8-de_x8
```
## x9
```{r}
beta9 <- sarma_bn_model$coefficients[12]
b9 <- beta9*invers
de_x9 <- mean(diag(b9))
te_x9 <- (1/35)*It_transpose%*%b9%*%It
ie_x9 <- te_x9-de_x9
```
## x11
```{r}
beta11 <- sarma_bn_model$coefficients[14]
b11 <- beta11*invers
de_x11 <- mean(diag(b11))
te_x11 <- (1/35)*It_transpose%*%b11%*%It
ie_x11 <- te_x11-de_x11
```
## penyajian DE, IE, dan TE
```{r}
marginal <- matrix(c(de_x1, de_x2, de_x4, de_x5, de_x7, de_x8, de_x9, de_x11,
ie_x1, ie_x2, ie_x4, ie_x5, ie_x7, ie_x8, ie_x9, ie_x11,
te_x1, te_x2, te_x4, te_x5, te_x7, te_x8, te_x9, te_x11),
ncol = 3, byrow = F)
colnames(marginal) <- c("Efek Langsung", "Efek Tak Langsung", "Efek Total")
rownames(marginal) <- c("x1", "x2", "x4", "x5", "x7", "x8", "x9", "x11")
marginal <- as.table(marginal)
marginal
```
## efek DE dan IE ke provinsi untuk x8
```{r}
direct_b8 <- as.matrix(diag(b8))
direct_b8_kabkot <- cbind(petajateng$KABKOT, direct_b8)
b8.2 <- b8
diag(b8.2) <- 0
indirect_b8 <- as.matrix(rowSums(b8.2))
indirect_b8_kabkot <- cbind(petajateng$KABKOT, indirect_b8)
dat.map <- st_read("jawa.shp")
dat.map <- dat.map %>%filter(PROVNO == "33")
dat.map$long <- longlat[, 1]
dat.map$lat <- longlat[, 2]
```
## DE untuk x8
```{r}
dat.map$de_x8 <- as.numeric(direct_b8[,1])
mycol <- c("green", "yellow", "red")
peta.de8 <- ggplot() +
geom_sf(data = dat.map, mapping = aes(geometry = geometry,
fill = de_x8)) +
scale_fill_gradientn(colors = mycol,
name = NULL) + geom_text(data = dat.map,
aes(long, lat, label = KABKOT,
fontface = "bold"),
color = "black", size = 1.5)+
labs(title = "Efek Langsung Variabel 8") +
theme(panel.background = element_rect(fill = "white"))
peta.de8
```
## IE X8
```{r}
dat.map$ie_x8 <- as.numeric(indirect_b8)
peta.ie8 <- ggplot() +
geom_sf(data = dat.map, mapping = aes(geometry = geometry,
fill = ie_x8)) +
scale_fill_gradientn(colors = mycol,
name = NULL) + geom_text(data = dat.map,
aes(long, lat, label = KABKOT,
fontface = "bold"),
color = "black", size = 1.5)+
labs(title = "Efek Tidak Langsung Variabel 8") +
theme(panel.background = element_rect(fill = "white"))
peta.ie8
```
## efek DE dan IE ke provinsi untuk x4
```{r}
direct_b4 <- as.matrix(diag(b4))
direct_b4_kabkot <- cbind(petajateng$KABKOT, direct_b4)
b4.2 <- b4
diag(b4.2) <- 0
indirect_b4 <- as.matrix(rowSums(b4.2))
indirect_b4_kabkot <- cbind(petajateng$KABKOT, indirect_b4)
dat.map <- st_read("jawa.shp")
dat.map <- dat.map %>%filter(PROVNO == "33")
dat.map$long <- longlat[, 1]
dat.map$lat <- longlat[, 2]
```
## DE untuk x4
```{r}
dat.map$de_x4 <- as.numeric(direct_b4[,1])
mycol <- c("green", "yellow", "red")
peta.de4 <- ggplot() +
geom_sf(data = dat.map, mapping = aes(geometry = geometry,
fill = de_x4)) +
scale_fill_gradientn(colors = mycol,
name = NULL) + geom_text(data = dat.map,
aes(long, lat, label = KABKOT,
fontface = "bold"),
color = "black", size = 1.5)+
labs(title = "Efek Langsung Variabel 4") +
theme(panel.background = element_rect(fill = "white"))
peta.de4
```
## IE x4
```{r}
dat.map$ie_x4 <- as.numeric(indirect_b4)
peta.ie4 <- ggplot() +
geom_sf(data = dat.map, mapping = aes(geometry = geometry,
fill = ie_x4)) +
scale_fill_gradientn(colors = mycol,
name = NULL) + geom_text(data = dat.map,
aes(long, lat, label = KABKOT,
fontface = "bold"),
color = "black", size = 1.5)+
labs(title = "Efek Tidak Langsung Variabel 4") +
theme(panel.background = element_rect(fill = "white"))
peta.ie4
```
## x5
```{r}
direct_b5 <- as.matrix(diag(b5))
direct_b5_kabkot <- cbind(petajateng$KABKOT, direct_b5)
b5.2 <- b5
diag(b5.2) <- 0
indirect_b5 <- as.matrix(rowSums(b5.2))
indirect_b5_kabkot <- cbind(petajateng$KABKOT, indirect_b5)
dat.map <- st_read("jawa.shp")
dat.map <- dat.map %>%filter(PROVNO == "33")
dat.map$long <- longlat[, 1]
dat.map$lat <- longlat[, 2]
```
## de x5
```{r}
dat.map$de_x5 <- as.numeric(direct_b5[,1])
mycol <- c("green", "yellow", "red")
peta.de5 <- ggplot() +
geom_sf(data = dat.map, mapping = aes(geometry = geometry,
fill = de_x5)) +
scale_fill_gradientn(colors = mycol,
name = NULL) + geom_text(data = dat.map,
aes(long, lat, label = KABKOT,
fontface = "bold"),
color = "black", size = 1.5)+
labs(title = "Efek Langsung Variabel 5") +
theme(panel.background = element_rect(fill = "white"))
peta.de5
```
## ie x5
```{r}
dat.map$ie_x5 <- as.numeric(indirect_b5)
peta.ie5 <- ggplot() +
geom_sf(data = dat.map, mapping = aes(geometry = geometry,
fill = ie_x5)) +
scale_fill_gradientn(colors = mycol,
name = NULL) + geom_text(data = dat.map,
aes(long, lat, label = KABKOT,
fontface = "bold"),
color = "black", size = 1.5)+
labs(title = "Efek Tidak Langsung Variabel 5") +
theme(panel.background = element_rect(fill = "white"))
peta.ie5
```
```{r}
# Fungsi untuk menghitung MAD, MAPE, dan RMSE
calculate_metrics <- function(actual, predicted) {
mad <- mean(abs(actual - predicted)) # MAD
mape <- mean(abs((actual - predicted) / actual)) * 100 # MAPE
rmse <- sqrt(mean((actual - predicted)^2)) # RMSE
return(c(MAD = mad, MAPE = mape, RMSE = rmse))
}
# Mendapatkan metrik untuk masing-masing model
# Model BN
predicted_bn <- predict(bn_model, type = "response")
actual_bn <- petajateng$y # Nilai aktual
metrics_bn <- calculate_metrics(actual_bn, predicted_bn)
# Model SAR
predicted_sar <- predict(sar_model, type = "response")
actual_sar <- petajateng$y # Nilai aktual
metrics_sar <- calculate_metrics(actual_sar, predicted_sar)
# Model SEM
predicted_sem <- predict(sem_bn_model, type = "response")
actual_sem <- petajateng$y # Nilai aktual
metrics_sem <- calculate_metrics(actual_sem, predicted_sem)
# Model SARMA
predicted_sarma <- predict(sarma_bn_model, type = "response")
actual_sarma <- petajateng$y # Nilai aktual
metrics_sarma <- calculate_metrics(actual_sarma, predicted_sarma)
# Membuat data frame untuk menampilkan hasil
results <- data.frame(
Model = c("BN", "SAR", "SEM", "SARMA"),
MAD = c(metrics_bn["MAD"],metrics_sar["MAD"], metrics_sem["MAD"], metrics_sarma["MAD"]),
MAPE = c(metrics_bn["MAPE"], metrics_sar["MAPE"], metrics_sem["MAPE"], metrics_sarma["MAPE"]),
RMSE = c(metrics_bn["RMSE"], metrics_sar["RMSE"], metrics_sem["RMSE"], metrics_sarma["RMSE"])
)
# Menampilkan hasil
results
```