library(rgdal)
## Loading required package: sp
## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/Ny 07/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/Ny 07/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
library(maptools)
## Checking rgeos availability: TRUE
library(raster)
library(spdep)
## Loading required package: spData
## 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.0.5
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(spatialreg)
## Loading required package: Matrix
## Registered S3 methods overwritten by 'spatialreg':
## method from
## residuals.stsls spdep
## deviance.stsls spdep
## coef.stsls spdep
## print.stsls spdep
## summary.stsls spdep
## print.summary.stsls spdep
## residuals.gmsar spdep
## deviance.gmsar spdep
## coef.gmsar spdep
## fitted.gmsar spdep
## print.gmsar spdep
## summary.gmsar spdep
## print.summary.gmsar spdep
## print.lagmess spdep
## summary.lagmess spdep
## print.summary.lagmess spdep
## residuals.lagmess spdep
## deviance.lagmess spdep
## coef.lagmess spdep
## fitted.lagmess spdep
## logLik.lagmess spdep
## fitted.SFResult spdep
## print.SFResult spdep
## fitted.ME_res spdep
## print.ME_res spdep
## print.lagImpact spdep
## plot.lagImpact spdep
## summary.lagImpact spdep
## HPDinterval.lagImpact spdep
## print.summary.lagImpact spdep
## print.sarlm spdep
## summary.sarlm spdep
## residuals.sarlm spdep
## deviance.sarlm spdep
## coef.sarlm spdep
## vcov.sarlm spdep
## fitted.sarlm spdep
## logLik.sarlm spdep
## anova.sarlm spdep
## predict.sarlm spdep
## print.summary.sarlm spdep
## print.sarlm.pred spdep
## as.data.frame.sarlm.pred spdep
## residuals.spautolm spdep
## deviance.spautolm spdep
## coef.spautolm spdep
## fitted.spautolm spdep
## print.spautolm spdep
## summary.spautolm spdep
## logLik.spautolm spdep
## print.summary.spautolm spdep
## print.WXImpact spdep
## summary.WXImpact spdep
## print.summary.WXImpact spdep
## predict.SLX spdep
##
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
##
## anova.sarlm, as.spam.listw, as_dgRMatrix_listw, as_dsCMatrix_I,
## as_dsCMatrix_IrW, as_dsTMatrix_listw, bptest.sarlm, can.be.simmed,
## cheb_setup, coef.gmsar, coef.sarlm, coef.spautolm, coef.stsls,
## create_WX, deviance.gmsar, deviance.sarlm, deviance.spautolm,
## deviance.stsls, do_ldet, eigen_pre_setup, eigen_setup, eigenw,
## errorsarlm, fitted.gmsar, fitted.ME_res, fitted.sarlm,
## fitted.SFResult, fitted.spautolm, get.ClusterOption,
## get.coresOption, get.mcOption, get.VerboseOption,
## get.ZeroPolicyOption, GMargminImage, GMerrorsar, griffith_sone,
## gstsls, Hausman.test, HPDinterval.lagImpact, impacts, intImpacts,
## Jacobian_W, jacobianSetup, l_max, lagmess, lagsarlm, lextrB,
## lextrS, lextrW, lmSLX, logLik.sarlm, logLik.spautolm, LR.sarlm,
## LR1.sarlm, LR1.spautolm, LU_prepermutate_setup, LU_setup,
## Matrix_J_setup, Matrix_setup, mcdet_setup, MCMCsamp, ME, mom_calc,
## mom_calc_int2, moments_setup, powerWeights, predict.sarlm,
## predict.SLX, print.gmsar, print.ME_res, print.sarlm,
## print.sarlm.pred, print.SFResult, print.spautolm, print.stsls,
## print.summary.gmsar, print.summary.sarlm, print.summary.spautolm,
## print.summary.stsls, residuals.gmsar, residuals.sarlm,
## residuals.spautolm, residuals.stsls, sacsarlm, SE_classic_setup,
## SE_interp_setup, SE_whichMin_setup, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption, similar.listw, spam_setup, spam_update_setup,
## SpatialFiltering, spautolm, spBreg_err, spBreg_lag, spBreg_sac,
## stsls, subgraph_eigenw, summary.gmsar, summary.sarlm,
## summary.spautolm, summary.stsls, trW, vcov.sarlm, Wald1.sarlm
library(readxl)
library(sp)
memanggil data shp
ind<-readOGR(dsn="D:/SHP Indonesia", layer="prov")
## OGR data source with driver: ESRI Shapefile
## Source: "D:\SHP Indonesia", layer: "prov"
## with 34 features
## It has 4 fields
ind
## class : SpatialPolygonsDataFrame
## features : 34
## extent : 95.0097, 141.0194, -11.00762, 6.076941 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 4
## names : NAME_0, NAME_1, KODE, ipm
## min values : Indonesia, Aceh, 11, NA
## max values : Indonesia, Yogyakarta, 92, NA
##input data
covid<-read.csv("D:\\total.csv",header=T, sep=";", stringsAsFactors = F)
summary(covid)
## Provinsi Total
## Length:34 Min. : 1941
## Class :character 1st Qu.: 3806
## Mode :character Median : 9207
## Mean : 21839
## 3rd Qu.: 18135
## Max. :183248
membuat peta
ind$ipm<-covid$Total
ind$NAME_1<-covid$Provinsi
plot(ind,col=c(1:5), main="Peta Indonesia")
legend("bottomleft", legend = c("c 1","c 2","c 3","c 4","c 5"),
fill=c(1:5))
text(ind,"NAME_1",cex=0.4,col="gray40")
#indikasi kelasterisasinya
spplot(ind, "ipm", col.regions = rainbow(99, start=.1))
MORAN InDEX
#jarak_EUC <- dist(coordinates(jawa), method="euclidean")
#w<-as.matrix(1/jarak_EUC)
#w<-poly2nb(ind, queen=F)
w<-as.matrix(1/dist(coordinates(ind)))
ww<-mat2listw(w)
moran(covid$Total,ww,n=length(ww$neighbours),S0=Szero(ww))
## $I
## [1] 0.09778682
##
## $K
## [1] 13.14362
moran.test(covid$Total, ww,randomisation=T,alternative="two.sided")
##
## Moran I test under randomisation
##
## data: covid$Total
## weights: ww
##
## Moran I statistic standard deviate = 3.6247, p-value = 0.0002893
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.097786815 -0.030303030 0.001248797
a<-moran.plot(covid$Total, ww, labels=F)
text(a,covid$Provinsi,cex=0.5, col="red")
LISA
s<-localmoran(covid$Total, ww)
moran.map <- cbind(ind, s)
quadrant <- vector(mode="numeric",length=nrow(s))
quadrant
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
m.qualification <- covid$Total - mean(covid$Total)
# centers the local Moran's around the mean
m.local <- s[,4] - mean(s[,4])
# significance threshold
signif <- 0.3
# builds a data quadrant
quadrant[m.qualification >0 & m.local>0 & (s[, 5] <= 0.3)] <- 4
quadrant[m.qualification <0 & m.local<0 & (s[, 5] <= 0.3)] <- 1
quadrant[m.qualification <0 & m.local>0 & (s[, 5] <= 0.3)] <- 2
quadrant[m.qualification >0 & m.local<0 & (s[, 5] <= 0.3)] <- 3
quadrant[s[,5]>signif] <- 0
# plot in r
brks <- c(0,1,2,3,4)
colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")
plot(ind,border="lightgray",main="Local Moran's",col=colors[findInterval(quadrant,brks,all.inside=FALSE)])
box()
legend("bottomleft", legend = c("insignificant","low-low","low-high","high-low","high-high"),
fill=colors,bty="n")
text(ind,"NAME_1",cex=0.45, col="gray40")
