Indo_Kec<-readRDS('gadm36_IDN_3_sp.rds')
Bandung<-Indo_Kec[Indo_Kec$NAME_2 == "Kota Bandung",]
plot(Bandung)
Bandung$id<-c(1:30)
Bandung_sf<-st_as_sf(Bandung)
Bandung_merged <- Bandung_sf %>%
left_join(Diare, by = "id")
ggplot() +
geom_sf(data=Bandung_merged, aes(fill = Diare),color=NA) +
theme_bw() +
scale_fill_gradient(low = "yellow", high = "red") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
theme(legend.position = "right",
axis.text.x = element_blank(), # Remove x-axis labels
axis.text.y = element_blank())+ # Remove y-axis labels
labs(title = "",
fill = "Diare")
summary(Diare$Diare)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235.0 529.0 775.0 939.1 1119.2 2898.0
breaks <- c(-Inf, 550, 800, 1200,Inf)
# Define labels for each interval
labels <- c("Very Low", "Low", "High", "Very High")
# Create a new column with discretized Diare
Bandung_merged$Diare_Discrete <- cut(Bandung_merged$Diare, breaks = breaks, labels = labels, right = TRUE)
ggplot() +
geom_sf(data=Bandung_merged, aes(fill = Diare_Discrete),color=NA) +
theme_bw() +
scale_fill_manual(values = c("Very Low" = "yellow",
"Low" = "orange",
"High" = "red",
"Very High" = "red3"))+
labs(fill = "Diare")+theme(legend.position = "right",
axis.text.x = element_blank(), # Remove x-axis labels
axis.text.y = element_blank())+ # Remove y-axis labels
labs(title = "",
fill = "Diare")
Spatial
Atocorrelation
Matriks Bobot
Spasial
#Membuat matrix bobot
row.names(Bandung) <- as.character(1:30)
W <- poly2nb(Bandung, row.names(Bandung), queen=FALSE)
WB <- nb2mat(W, style='B', zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"
WL<-nb2listw(W)#List neighbours
CoordK<-coordinates(Bandung)
plot(Bandung, axes=T, col="gray90")
text(CoordK[,1], CoordK[,2], row.names(Bandung), col="black", cex=0.8, pos=1.5)
points(CoordK[,1], CoordK[,2], pch=19, cex=0.7,col="blue")
plot.nb(W, CoordK, add=TRUE, col="red", lwd=2)
Global Moran’s I
Global_Moran<-moran.test(Diare$Diare,WL)
Global_Moran
##
## Moran I test under randomisation
##
## data: Diare$Diare
## weights: WL
##
## Moran I statistic standard deviate = 0.67226, p-value = 0.2507
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.04045578 -0.03448276 0.01242630
Local Moran’s I
Local_Moran <- localmoran(Diare$Diare, WL)
quadr_data <- attr(Local_Moran, "quadr")
mean_values <- quadr_data$mean
Bandung_sf <- st_as_sf(Bandung)
mean_values_char <- as.character(attr(Local_Moran, "quadr")$mean)
Bandung_sf$mean_values <- mean_values_char
ggplot() +
geom_sf(data = Bandung_sf, aes(fill = mean_values)) + # Fill based on 'mean_values'
scale_fill_manual(values = c("Low-Low" = "blue", "High-Low" = "green",
"Low-High" = "yellow", "High-High" = "red")) +
labs(title = "Local Moran's I Mean Values in Bandung", fill = "Mean Type") +
theme_minimal()
Spatial Regression
LM Test
LM<-lm.LMtests(Diare.ols, WL, test="all")
## Please update scripts to use lm.RStests in place of lm.LMtests
print(LM)
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
##
## RSerr = 0.22911, df = 1, p-value = 0.6322
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
##
## RSlag = 0.14934, df = 1, p-value = 0.6992
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
##
## adjRSerr = 0.27935, df = 1, p-value = 0.5971
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
##
## adjRSlag = 0.19957, df = 1, p-value = 0.6551
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare)
## test weights: listw
##
## SARMA = 0.42869, df = 2, p-value = 0.8071
Spatial lag model
(SAR)
sar.Diare<-lagsarlm(Prevalensi~PHBS+Kepadatan, data=Diare, WL)
summary(sar.Diare)
##
## Call:lagsarlm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare,
## listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02698 -0.56529 -0.12428 0.30811 3.21625
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.8112361 1.0066961 2.7925 0.00523
## PHBS -0.0138353 0.0104439 -1.3247 0.18526
## Kepadatan -0.0033585 0.0022850 -1.4698 0.14163
##
## Rho: -0.097136, LR test value: 0.14535, p-value: 0.70302
## Asymptotic standard error: 0.26437
## z-value: -0.36743, p-value: 0.7133
## Wald statistic: 0.135, p-value: 0.7133
##
## Log likelihood: -36.8056 for lag model
## ML residual variance (sigma squared): 0.67957, (sigma: 0.82436)
## Number of observations: 30
## Number of parameters estimated: 5
## AIC: 83.611, (AIC for lm: 81.757)
## LM test for residual autocorrelation
## test value: 0.15074, p-value: 0.69783
sar2sls.Diare<-stsls(Prevalensi~PHBS+Kepadatan, data=Diare, WL)
summary(sar2sls.Diare)
##
## Call:stsls(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare,
## listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.95775 -0.86809 -0.14850 0.79392 3.03386
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Rho -2.4145328 3.1173479 -0.7745 0.4386
## (Intercept) 7.1811683 6.0915490 1.1789 0.2384
## PHBS -0.0337948 0.0321396 -1.0515 0.2930
## Kepadatan -0.0041641 0.0040115 -1.0380 0.2992
##
## Residual variance (sigma squared): 1.9444, (sigma: 1.3944)
Cek residual
Diare$Diare.ols.res<-resid(Diare.ols) #residuals ols
Diare$Diare.sar.res<-resid(sar.Diare) #residual sar
Bandung@data<-Diare
spplot(Bandung,"Diare.ols.res", at=seq(min(Bandung$Diare.ols.res,na.rm=TRUE),max(Bandung$Diare.ols.res,na.rm=TRUE),length=12),col.regions=rev(brewer.pal(11,"RdBu")))
spplot(Bandung,"Diare.sar.res", at=seq(min(Bandung$Diare.ols.res,na.rm=TRUE),max(Bandung$Diare.ols.res,na.rm=TRUE),length=12),col.regions=rev(brewer.pal(11,"RdBu")))
Menghitung
Impact
impacts(sar.Diare, listw=WL)
## Impact measures (lag, exact):
## Direct Indirect Total
## PHBS -0.013864508 0.0012541191 -0.012610389
## Kepadatan -0.003365544 0.0003044315 -0.003061113
Spatial Error
model
errorsalm.Diare<-errorsarlm(Prevalensi~PHBS+Kepadatan, data=Diare, WL)
summary(errorsalm.Diare)
##
## Call:errorsarlm(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare,
## listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.04284 -0.55891 -0.12354 0.30761 3.19555
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.7387708 0.9262401 2.9569 0.003108
## PHBS -0.0147006 0.0101926 -1.4423 0.149221
## Kepadatan -0.0032921 0.0022216 -1.4819 0.138376
##
## Lambda: -0.1277, LR test value: 0.23418, p-value: 0.62845
## Asymptotic standard error: 0.27122
## z-value: -0.47083, p-value: 0.63776
## Wald statistic: 0.22168, p-value: 0.63776
##
## Log likelihood: -36.76118 for error model
## ML residual variance (sigma squared): 0.67652, (sigma: 0.82251)
## Number of observations: 30
## Number of parameters estimated: 5
## AIC: 83.522, (AIC for lm: 81.757)
fgls.Diare<-GMerrorsar(Prevalensi~PHBS+Kepadatan, data=Diare, WL)
summary(fgls.Diare)
##
## Call:GMerrorsar(formula = Prevalensi ~ PHBS + Kepadatan, data = Diare,
## listw = WL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.015843 -0.583648 -0.099188 0.305259 3.221005
##
## Type: GM SAR estimator
## Coefficients: (GM standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6857318 0.9368412 2.8668 0.004147
## PHBS -0.0138851 0.0103359 -1.3434 0.179148
## Kepadatan -0.0033074 0.0022558 -1.4661 0.142613
##
## Lambda: -0.063722 (standard error): 0.48559 (z-value): -0.13123
## Residual variance (sigma squared): 0.67996, (sigma: 0.8246)
## GM argmin sigma squared: 0.69511
## Number of observations: 30
## Number of parameters estimated: 5
