Spasial Rabies
Andini Dwi Kurnia Putri
2025-10-14
Statistik Deskriptif
file_path <- "D:/SEMESTER 5/Spasial/Data Rabies.csv"
rabies_data <- read.csv(file_path)
head(rabies_data)## No Provinsi GHPR22 VAR22 LYSSA22 GHPR23 VAR23 LYSSA23 GHPR24 VAR24
## 1 1 Aceh 1494 1130 0 1434 1036 0 1493 1046
## 2 2 Sumatera Utara 6883 5720 6 10360 7848 20 11377 8448
## 3 3 Sumatera Barat 4504 2883 8 6638 4636 7 7335 5206
## 4 4 Riau 2330 1627 0 4858 4065 2 5548 4823
## 5 5 Jambi 1031 784 0 1511 1278 1 1641 1165
## 6 6 Sumatera Selatan 2432 2197 1 4183 3820 4 4674 4489
## LYSSA24
## 1 0
## 2 14
## 3 3
## 4 5
## 5 0
## 6 5str(rabies_data)## 'data.frame': 38 obs. of 11 variables:
## $ No : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Provinsi: chr "Aceh" "Sumatera Utara" "Sumatera Barat" "Riau" ...
## $ GHPR22 : int 1494 6883 4504 2330 1031 2432 1380 2188 50 11 ...
## $ VAR22 : int 1130 5720 2883 1627 784 2197 1133 1987 37 NA ...
## $ LYSSA22 : int 0 6 8 0 0 1 1 0 0 0 ...
## $ GHPR23 : int 1434 10360 6638 4858 1511 4183 2147 3043 168 17 ...
## $ VAR23 : int 1036 7848 4636 4065 1278 3820 1945 2656 94 4 ...
## $ LYSSA23 : int 0 20 7 2 1 4 1 0 0 0 ...
## $ GHPR24 : int 1493 11377 7335 5548 1641 4674 2122 3152 286 116 ...
## $ VAR24 : int 1046 8448 5206 4823 1165 4489 1912 2702 193 70 ...
## $ LYSSA24 : int 0 14 3 5 0 5 1 0 0 0 ...summary(rabies_data$GHPR24)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 536 1573 4878 4604 57943summary(rabies_data$VAR24)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 276 1298 3664 4105 32565summary(rabies_data$LYSSA24)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 3.211 3.750 46.000Visualisasi Awal
Boxplot
par(mfrow = c(1, 3))
boxplot(rabies_data$GHPR24, main = "Boxplot GHPR24", col = "skyblue")
# Boxplot untuk VAR24
boxplot(rabies_data$VAR24, main = "Boxplot VAR24", col = "lightcoral")
# Boxplot untuk LYSSA24
boxplot(rabies_data$LYSSA24, main = "Boxplot LYSSA24", col = "yellowgreen")
# Mengatur ulang layout plot ke default
par(mfrow = c(1, 1))# Muat semua paket
library(geodata)## Loading required package: terra## Warning: package 'terra' was built under R version 4.3.3## terra 1.8.29library(sf)## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(ggplot2)
library(dplyr)## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(spdep)## Warning: package 'spdep' was built under R version 4.3.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.3.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(sp)## Warning: package 'sp' was built under R version 4.3.3library(terra)
library(GWmodel)## Warning: package 'GWmodel' was built under R version 4.3.3## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 4.3.3## Loading required package: Rcpp## Warning: package 'Rcpp' was built under R version 4.3.3## Welcome to GWmodel version 2.4-2.Peta Indo
peta_spat <- gadm(country = "IDN", level = 1, path = ".")
peta_indonesia <- st_as_sf(peta_spat)
print(names(peta_indonesia))## [1] "GID_1" "GID_0" "COUNTRY" "NAME_1" "VARNAME_1" "NL_NAME_1"
## [7] "TYPE_1" "ENGTYPE_1" "CC_1" "HASC_1" "ISO_1" "geometry"peta_prov <- "NAME_1"nama_kolom_peta <- "NAME_1"
nama_kolom_data <- "Provinsi"
peta_rabies_lengkap <- left_join(
peta_indonesia,
rabies_data,
by = setNames(nama_kolom_data, nama_kolom_peta)
)
head(peta_rabies_lengkap)## Simple feature collection with 6 features and 21 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95.00971 ymin: -8.849262 xmax: 123.5523 ymax: 6.076941
## Geodetic CRS: WGS 84
## GID_1 GID_0 COUNTRY NAME_1 VARNAME_1 NL_NAME_1 TYPE_1
## 1 IDN.1_1 IDN Indonesia Aceh <NA> <NA> Propinisi
## 2 IDN.2_1 IDN Indonesia Bali <NA> <NA> Propinisi
## 3 IDN.3_1 IDN Indonesia Bangka Belitung <NA> <NA> Propinisi
## 4 IDN.4_1 IDN Indonesia Banten <NA> <NA> Propinisi
## 5 IDN.5_1 IDN Indonesia Bengkulu <NA> <NA> Propinisi
## 6 IDN.6_1 IDN Indonesia Gorontalo <NA> <NA> Propinisi
## ENGTYPE_1 CC_1 HASC_1 ISO_1 No GHPR22 VAR22 LYSSA22 GHPR23 VAR23 LYSSA23
## 1 Province 11 ID.AC ID-AC 1 1494 1130 0 1434 1036 0
## 2 Province 51 ID.BA ID-BA 17 38009 20755 22 72522 47043 9
## 3 Province 19 ID.BB <NA> 9 50 37 0 168 94 0
## 4 Province 36 ID.BT ID-BT 16 254 233 0 932 869 0
## 5 Province 17 ID.BE ID-BE 7 1380 1133 1 2147 1945 1
## 6 Province 75 ID.GO ID-GO 29 812 712 2 1015 925 0
## GHPR24 VAR24 LYSSA24 geometry
## 1 1493 1046 0 MULTIPOLYGON (((98.14903 2....
## 2 57943 32565 7 MULTIPOLYGON (((115.5257 -8...
## 3 286 193 0 MULTIPOLYGON (((107.9614 -3...
## 4 790 666 0 MULTIPOLYGON (((106.3868 -6...
## 5 2122 1912 1 MULTIPOLYGON (((103.573 -4....
## 6 1075 994 2 MULTIPOLYGON (((123.5491 0....names(peta_rabies_lengkap)## [1] "GID_1" "GID_0" "COUNTRY" "NAME_1" "VARNAME_1" "NL_NAME_1"
## [7] "TYPE_1" "ENGTYPE_1" "CC_1" "HASC_1" "ISO_1" "No"
## [13] "GHPR22" "VAR22" "LYSSA22" "GHPR23" "VAR23" "LYSSA23"
## [19] "GHPR24" "VAR24" "LYSSA24" "geometry"Peta
GHPR
ggplot(data = peta_rabies_lengkap) +
geom_sf(aes(fill = GHPR24),
color = "white",
linewidth = 0.15) +
scale_fill_viridis_c(option = "magma",
name = "GHPR 2024\n(Gigitan Hewan Penular Rabies)") +
labs(
title = "Peta Sebaran Kasus GHPR (Gigitan Hewan Penular Rabies) di Indonesia, 2024",
subtitle = "Menunjukkan kepadatan kasus GHPR di setiap provinsi",
caption = "Sumber: Data Rabies dari Profil Kesehatan Indonesia 2024 "
) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
plot.subtitle = element_text(hjust = 0.5, size = 10),
legend.position = "bottom",
legend.title = element_text(face = "bold")
)
### VAR
ggplot(data = peta_rabies_lengkap) +
geom_sf(aes(fill = VAR24),
color = "white",
linewidth = 0.15) +
scale_fill_viridis_c(option = "magma",
name = "VAR 2024\n(Vaksin Anti Rabies)") +
labs(
title = "Peta Sebaran Penerima Vaksin Anti Rabies di Indonesia, 2024",
subtitle = "Menunjukkan Sebaran Vaksin di setiap provinsi",
caption = "Sumber: Data Rabies dari Profil Kesehatan Indonesia 2024 "
) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
plot.subtitle = element_text(hjust = 0.5, size = 10),
legend.position = "bottom",
legend.title = element_text(face = "bold")
)
### LYSSA
ggplot(data = peta_rabies_lengkap) +
geom_sf(aes(fill = LYSSA24),
color = "white",
linewidth = 0.15) +
scale_fill_viridis_c(option = "magma",
name = "LYSSA 2024\n(Positif Rabies dan Meninggal)") +
labs(
title = "Peta Penderita Rabies yang Meninggal di Indonesia, 2024",
subtitle = "Menunjukkan Penderita Rabies yang Meninggal di Setiap Provinsi",
caption = "Sumber: Data Rabies dari Profil Kesehatan Indonesia 2024 "
) +
theme_void() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
plot.subtitle = element_text(hjust = 0.5, size = 10),
legend.position = "bottom",
legend.title = element_text(face = "bold")
)
# Autokorelasi Spasial ## Moran’s I
peta_rabies_analisis <- peta_rabies_lengkap %>%
filter(!is.na(LYSSA24))
neighbors <- poly2nb(peta_rabies_analisis, queen = TRUE)## Warning in poly2nb(peta_rabies_analisis, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in poly2nb(peta_rabies_analisis, queen = TRUE): neighbour object has 10 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.listw <- nb2listw(neighbors, style = "W", zero.policy = TRUE)
moran_test <- moran.test(peta_rabies_analisis$LYSSA24, listw, zero.policy = TRUE)
print(moran_test)##
## Moran I test under randomisation
##
## data: peta_rabies_analisis$LYSSA24
## weights: listw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 0.6998, p-value = 0.242
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.035329903 -0.034482759 0.009952289moran_i_value <- moran_test$estimate[1]
p_value <- moran_test$p.value
cat("\n--- Hasil Moran's I Global (LYSSA24) ---\n")##
## --- Hasil Moran's I Global (LYSSA24) ---cat(paste("Moran's I Observed:", round(moran_i_value, 4), "\n"))## Moran's I Observed: 0.0353cat(paste("P-value:", round(p_value, 4), "\n"))## P-value: 0.242LISA LYSSA24
x <- scale(peta_rabies_analisis$LYSSA24)[,1]
lagx <- spdep::lag.listw(listw, x, zero.policy = TRUE)
lisa <- spdep::localmoran(x, listw, 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)"
)
peta_LISA <- dplyr::bind_cols(peta_rabies_analisis, lisa_df) %>%
dplyr::mutate(
quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant")
)
ggplot(peta_LISA) +
geom_sf(aes(fill = quad), color="white", linewidth=0.2) +
scale_fill_manual(
values=c(
"High-High"="#d73027", # Merah Gelap (Hotspot)
"Low-Low"="#4575b4", # Biru Gelap (Coldspot)
"High-Low (Outlier)"="#fdae61", # Orange
"Low-High (Outlier)"="#74add1", # Biru Muda
"Not significant"="grey85" # Abu-abu
),
drop = FALSE
) +
labs(
title = "Peta Klaster LISA Kasus Rabies Meninggal (LYSSA24), 2024",
fill = "Kategori Klaster"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 14)
)
## Sensitivitas Bobot Spasial
nb_rook <- spdep::poly2nb(peta_rabies_analisis, queen = FALSE) ## Warning in spdep::poly2nb(peta_rabies_analisis, queen = FALSE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.## Warning in spdep::poly2nb(peta_rabies_analisis, queen = FALSE): neighbour object has 10 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.lw_rook <- spdep::nb2listw(nb_rook, style="W", zero.policy = TRUE)
coords <- sf::st_coordinates(sf::st_centroid(peta_rabies_analisis))## Warning: st_centroid assumes attributes are constant over geometriescoords <- sf::st_coordinates(sf::st_centroid(peta_rabies_analisis))## Warning: st_centroid assumes attributes are constant over geometriesnb_knn4 <- spdep::knn2nb(spdep::knearneigh(coords, k = 4))
lw_knn4 <- spdep::nb2listw(nb_knn4, style="W")
lyssa24_data <- peta_rabies_analisis$LYSSA24
c(
# 1. Queen (listw)
mean_z_queen = mean(spdep::localG(lyssa24_data, listw, zero.policy = TRUE)),
# 2. Rook (lw_rook)
mean_z_rook = mean(spdep::localG(lyssa24_data, lw_rook, zero.policy = TRUE)),
# 3. kNN (lw_knn4)
# Menggunakan lw_knn4 (k=4) yang telah dibuat sebelumnya
mean_z_knn4 = mean(spdep::localG(lyssa24_data, lw_knn4, zero.policy = TRUE))
)## mean_z_queen mean_z_rook mean_z_knn4
## NaN NaN -0.257687Estimasi Model Spatial Econometrics
SEM
library(spatialreg)## Warning: package 'spatialreg' was built under R version 4.3.3## Loading required package: Matrix##
## 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.ZeroPolicyOptiondat_model <- peta_rabies_analisis
dat_model$X1 <- dat_model$GHPR24
dat_model$X2 <- dat_model$VAR24
lw_model <- listw
sem_rabies <- spatialreg::errorsarlm(
formula = LYSSA24 ~ X1 + X2,
data = dat_model,
listw = lw_model,
# type = "error" tidak diperlukan karena sudah menggunakan errorsarlm
method = "eigen",
zero.policy = TRUE
)
summary(sem_rabies)##
## Call:spatialreg::errorsarlm(formula = LYSSA24 ~ X1 + X2, data = dat_model,
## listw = lw_model, method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.62623 -1.21228 0.12367 0.74974 9.17157
##
## Type: error
## Regions with no neighbours included:
## 2 3 21 22
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.38477558 0.80004422 -0.4809 0.6306
## X1 -0.00163727 0.00019927 -8.2161 2.22e-16
## X2 0.00320350 0.00029927 10.7046 < 2.2e-16
##
## Lambda: 0.20637, LR test value: 1.4843, p-value: 0.22311
## Asymptotic standard error: 0.16571
## z-value: 1.2454, p-value: 0.213
## Wald statistic: 1.5509, p-value: 0.213
##
## Log likelihood: -89.30386 for error model
## ML residual variance (sigma squared): 10.989, (sigma: 3.315)
## Number of observations: 34
## Number of parameters estimated: 5
## AIC: 188.61, (AIC for lm: 188.09)SAR(SLM)
sar_rabies <- spatialreg::lagsarlm(
formula = LYSSA24 ~ X1 + X2,
data = dat_model,
listw = lw_model,
type = "lag", # Menentukan jenis model: Spatial Lag / SAR
method = "eigen",
zero.policy = TRUE
)
summary(sar_rabies)##
## Call:spatialreg::lagsarlm(formula = LYSSA24 ~ X1 + X2, data = dat_model,
## listw = lw_model, type = "lag", method = "eigen", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.828726 -1.015926 -0.086591 0.838446 9.551572
##
## Type: lag
## Regions with no neighbours included:
## 2 3 21 22
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.51745676 0.76156658 -0.6795 0.4968
## X1 -0.00160690 0.00020471 -7.8496 4.219e-15
## X2 0.00315736 0.00030673 10.2937 < 2.2e-16
##
## Rho: -0.036197, LR test value: 0.047154, p-value: 0.82809
## Asymptotic standard error: 0.13347
## z-value: -0.2712, p-value: 0.78624
## Wald statistic: 0.073548, p-value: 0.78624
##
## Log likelihood: -90.02242 for lag model
## ML residual variance (sigma squared): 11.67, (sigma: 3.4161)
## Number of observations: 34
## Number of parameters estimated: 5
## AIC: 190.04, (AIC for lm: 188.09)
## LM test for residual autocorrelation
## test value: 4.0642, p-value: 0.043801Implementasi Analisis Dampak
imp_rabies <- spatialreg::impacts(sar_rabies, listw = lw_model, R = 500)
summary(imp_rabies, zstats = TRUE)## Impact measures (lag, exact):
## Direct Indirect Total
## X1 -0.001607767 5.700268e-05 -0.001550764
## X2 0.003159065 -1.120033e-04 0.003047062
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## X1 -0.001623 0.0002095 9.369e-06 9.369e-06
## X2 0.003187 0.0003136 1.402e-05 1.402e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 -0.002028 -0.001763 -0.001625 -0.001492 -0.001203
## X2 0.002559 0.002986 0.003207 0.003384 0.003789
##
## ========================================================
## Indirect:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## X1 3.156e-05 0.0002321 1.038e-05 1.038e-05
## X2 -5.967e-05 0.0004540 2.030e-05 2.030e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 -0.0004607 -0.0001128 5.154e-05 0.0001906 0.0004656
## X2 -0.0008690 -0.0003722 -1.011e-04 0.0002225 0.0009045
##
## ========================================================
## Total:
##
## Iterations = 1:500
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 500
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## X1 -0.001592 0.0002983 1.334e-05 1.334e-05
## X2 0.003128 0.0005360 2.397e-05 2.397e-05
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## X1 -0.002213 -0.001763 -0.001571 -0.001385 -0.001096
## X2 0.002217 0.002751 0.003080 0.003458 0.004330
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## X1 0.0002095038 0.0002320797 0.0002983030
## X2 0.0003135843 0.0004540242 0.0005360489
##
## Simulated z-values:
## Direct Indirect Total
## X1 -7.748608 0.1359829 -5.336199
## X2 10.164540 -0.1314284 5.834855
##
## Simulated p-values:
## Direct Indirect Total
## X1 9.3259e-15 0.89183 9.4915e-08
## X2 < 2.22e-16 0.89544 5.3837e-09Pemilihan Model Terbaik
Memiliki AIC terendah
ols_rabies <- lm(
formula = LYSSA24 ~ X1 + X2,
data = dat_model
)
aic_ols <- AIC(ols_rabies)
aic_sem <- 188.61
aic_sar <- 190.04
aic_comparison <- data.frame(
Model = c("OLS (Benchmark)", "Spatial Error Model (SEM)", "Spatial Lag Model (SAR)"),
AIC_Value = c(aic_ols, aic_sem, aic_sar)
)
aic_comparison_sorted <- aic_comparison[order(aic_comparison$AIC_Value), ]
print(aic_comparison_sorted)## Model AIC_Value
## 1 OLS (Benchmark) 188.092
## 2 Spatial Error Model (SEM) 188.610
## 3 Spatial Lag Model (SAR) 190.040Uji Asumsi
Normalitas
library(car) ## Warning: package 'car' was built under R version 4.3.3## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(nortest)
qqPlot(ols_rabies$residuals, main="Plot Q-Q Normalitas Residual OLS")
## [1] 27 30shapiro.test(ols_rabies$residuals)##
## Shapiro-Wilk normality test
##
## data: ols_rabies$residuals
## W = 0.87428, p-value = 0.001027Heteroskedastisitas
plot(ols_rabies$fitted.values, ols_rabies$residuals,
xlab = "Fitted Values (Prediksi)",
ylab = "Residuals (Sisaan)",
main = "Plot Heteroskedastisitas (Residual vs. Fitted)")
abline(h = 0, col = "red") 
library(lmtest)## Warning: package 'lmtest' was built under R version 4.3.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.3##
## Attaching package: 'zoo'## The following object is masked from 'package:terra':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(ols_rabies)##
## studentized Breusch-Pagan test
##
## data: ols_rabies
## BP = 4.6526, df = 2, p-value = 0.09765Auotokorelasi
moran_test_res <- spdep::moran.test(ols_rabies$residuals, lw_model, zero.policy = TRUE)
print(moran_test_res)##
## Moran I test under randomisation
##
## data: ols_rabies$residuals
## weights: lw_model
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 1.4436, p-value = 0.07443
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.20502422 -0.03448276 0.02752764lm_test <- spdep::lm.RStests(
model = ols_rabies,
listw = lw_model,
test = "all", # Menggunakan "all" untuk menjalankan semua uji yang tersedia
zero.policy = TRUE
)
print(lm_test)##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = LYSSA24 ~ X1 + X2, data = dat_model)
## test weights: lw_model
##
## RSerr = 1.536, df = 1, p-value = 0.2152
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = LYSSA24 ~ X1 + X2, data = dat_model)
## test weights: lw_model
##
## RSlag = 0.029705, df = 1, p-value = 0.8632
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = LYSSA24 ~ X1 + X2, data = dat_model)
## test weights: lw_model
##
## adjRSerr = 4.192, df = 1, p-value = 0.04062
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = LYSSA24 ~ X1 + X2, data = dat_model)
## test weights: lw_model
##
## adjRSlag = 2.6857, df = 1, p-value = 0.1013
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = LYSSA24 ~ X1 + X2, data = dat_model)
## test weights: lw_model
##
## SARMA = 4.2217, df = 2, p-value = 0.1211Plot Residual Spasial (Surface Plot)
# Muat paket ggplot2 dan sf
library(ggplot2)
library(sf)
# 1. Tambahkan Residual OLS ke data frame spasial Anda
dat_plot_res <- dat_model %>%
mutate(residual_ols = ols_rabies$residuals)
# 2. Plot Peta Residual
ggplot(data = dat_plot_res) +
geom_sf(aes(fill = residual_ols),
color = "black",
linewidth = 0.15) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0,
name = "Residual OLS\n(Kasus Terlalu Rendah/Tinggi)") +
labs(title = "Peta Sebaran Residual Model OLS",
subtitle = "Merah: Model meremehkan kasus; Biru: Model melebih-lebihkan kasus") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))