library(sp)
library(rgdal) # membaca shapefile
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision (unknown))
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
## Path to GDAL shared files: C:/Users/htauf/AppData/Local/R/win-library/4.4/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
## Path to PROJ shared files: C:/Users/htauf/AppData/Local/R/win-library/4.4/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(spData)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(readxl) # membaca file excel
library(raster) # text pada peta
library(spdep) # pembobot data spasial
## Loading required package: sf
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(lmtest) # regresi klasik
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(DescTools) # Uji dependensi
library(nortest) # uji normalitas
library(spatialreg) # running model SAR
## 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.ZeroPolicyOption
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:DescTools':
## 
##     Recode
library(nortest)
library(GWmodel) # Model GWR
## Loading required package: robustbase
## Loading required package: Rcpp
## Welcome to GWmodel version 2.4-2.
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# untuk melakukan pengaturan working directory

datasulsel <- read_excel("DATA.xlsx")
datasulsel <- datasulsel %>%
  mutate(across(c(Y, X1, X2, X3, X4, X5, X6, X7), as.numeric))
data<-datasulsel[3:9]
data
## # A tibble: 24 × 7
##       X1    X2    X3    X4    X5    X6    X7
##    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1  72.6  13.6  8.54 11712  8.46 0.371   156
##  2  73.2  13    7.54  9682 10.5  0.389   180
##  3  74.0  13.4  8.26 11392  7.22 0.379   383
##  4  73.9  13.9  8.94 11636 12.7  0.365   130
##  5  74.0  13.7  8.41 10233  7.42 0.376   442
##  6  73.8  12.1  7     9781 13.1  0.34    521
##  7  74.3  13.0  8.93 13451  6.93 0.393    46
##  8  73.6  12.6  8.14 12513 12.7  0.342    45
##  9  73.3  13.4  8.73 10691 12.7  0.341   131
## 10  75.2  15.6 11.6  17889  5.07 0.387  8227
## # ℹ 14 more rows
summary(datasulsel)
##  Kabupaten / Kota         Y               X1              X2       
##  Length:24          Min.   :68.95   Min.   :72.57   Min.   :12.12  
##  Class :character   1st Qu.:71.42   1st Qu.:73.27   1st Qu.:12.95  
##  Mode  :character   Median :73.27   Median :73.69   Median :13.29  
##                     Mean   :73.84   Mean   :73.75   Mean   :13.39  
##                     3rd Qu.:74.46   3rd Qu.:74.32   3rd Qu.:13.64  
##                     Max.   :84.85   Max.   :75.15   Max.   :15.61  
##        X3               X4                 X5               X6        
##  Min.   : 7.000   Min.   :    9.83   Min.   : 5.070   Min.   :0.3360  
##  1st Qu.: 7.973   1st Qu.: 9756.25   1st Qu.: 7.370   1st Qu.:0.3450  
##  Median : 8.305   Median :11514.00   Median : 8.725   Median :0.3635  
##  Mean   : 8.519   Mean   :10339.58   Mean   : 9.332   Mean   :0.3618  
##  3rd Qu.: 8.633   3rd Qu.:12604.00   3rd Qu.:12.322   3rd Qu.:0.3767  
##  Max.   :11.560   Max.   :17889.00   Max.   :13.400   Max.   :0.3930  
##        X7        
##  Min.   :  45.0  
##  1st Qu.: 146.5  
##  Median : 216.5  
##  Mean   : 666.8  
##  3rd Qu.: 461.5  
##  Max.   :8227.0
dim(datasulsel)
## [1] 24  9
#menampilkan peta
#Peta Dasar Batas-batas Wilayah dari Sulsel
sulsel<-readOGR(dsn=".",layer = "sul-sel")
## Warning in readOGR(dsn = ".", layer = "sul-sel"): OGR support is provided by
## the sf and terra packages among others
## Warning in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv =
## use_iconv, : OGR support is provided by the sf and terra packages among others
## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : OGR support is provided by the sf and terra packages among others
## Warning in ogrListLayers(dsn): OGR support is provided by the sf and terra
## packages among others
## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\htauf\Documents\Documents\freelancer\C4\revisi", layer: "sul-sel"
## with 24 features
## It has 11 fields
## Integer64 fields read as strings:  Y
kode=sulsel$ID
koordinat=coordinates(sulsel)
plot(sulsel)

sulsel$Kabupaten
##  [1] "BARRU"        "BONE"         "BULUKUMBA"    "ENREKANG"     "GOWA"        
##  [6] "JENEPONTO"    "LUWU TIMUR"   "LUWU UTARA"   "LUWU"         "MAKASSAR"    
## [11] "MAROS"        "PALOPO"       "PANGKEP"      "PARE-PARE"    "PINRANG"     
## [16] "SELAYAR"      "SIDRAP"       "SINJAI"       "SOPPENG"      "TAKALAR"     
## [21] "TORAJA UTARA" "TORAJA"       "WAJO"         "BANTAENG"
## eksplor data
colfunc<-colorRampPalette(c("green", "yellow", "orange","red"))
color<-colfunc(16)
sulsel$pertumbuhanUMKM<-datasulsel$Y
spplot(sulsel, "pertumbuhanUMKM", col.regions=color,
       main="Peta Sebaran Pertumbuhan UMKM")

sulsel$pertumbuhanUMKM<-datasulsel$X1
spplot(sulsel, "X1", col.regions=color,
       main="Peta Sebaran X1")

sulsel$pertumbuhanUMKM<-datasulsel$X2
spplot(sulsel, "X2", col.regions=color,
       main="Peta Sebaran X2")

sulsel$pertumbuhanUMKM<-datasulsel$X3
spplot(sulsel, "X3", col.regions=color,
       main="Peta Sebaran X3")

sulsel$pertumbuhanUMKM<-datasulsel$X4
spplot(sulsel, "X4", col.regions=color,
       main="Peta Sebaran X4")

sulsel$pertumbuhanUMKM<-datasulsel$X5
spplot(sulsel, "X5", col.regions=color,
       main="Peta Sebaran X5")

sulsel$pertumbuhanUMKM<-datasulsel$X6
spplot(sulsel, "X6", col.regions=color,
       main="Peta Sebaran X6")

sulsel$pertumbuhanUMKM<-datasulsel$X7
spplot(sulsel, "X7", col.regions=color,
       main="Peta Sebaran X7")

## plot
ggplot(datasulsel, aes(x = `Kabupaten / Kota` , y = Y, group = 1)) +
  geom_line(color = "blue", size = 1) + 
  geom_point(color = "blue", size = 3) +
  geom_text(aes(label = Y), position = 'dodge') +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + # Rotasi teks kabupaten
  labs(title = "Plot",
       x = "Kabupaten/Kota", 
       y = "variabel Y")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`

## Moran's index
w = gw.dist(dp.locat = koordinat)
ww = gw.weight(vdist = w, bw = 2, kernel = 'gaussian', adaptive = T)
#dist = mat2listw(ww)

# Membuat matriks bobot spasial menggunakan Queen contiguity
nb <- poly2nb(sulsel, queen = TRUE) # Tetangga berbasis batas wilayah
dist <- nb2listw(nb, style = "W")   # Konversi ke matriks bobot berbobot W

#nb <- poly2nb(sulsel)
#dist <- nb2listw(nb, style = "S")
moran(datasulsel$Y, dist, n=length(dist$neighbours), S0=Szero(dist))
## $I
## [1] 0.007111038
## 
## $K
## [1] 5.286209
moran.test(datasulsel$Y, dist, randomisation = F)
## 
##  Moran I test under normality
## 
## data:  datasulsel$Y  
## weights: dist    
## 
## Moran I statistic standard deviate = 0.35248, p-value = 0.3622
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.007111038      -0.043478261       0.020598751
moran.plot(datasulsel$Y, dist, labels=datasulsel$`Kabupaten / Kota`)

## regresi klasik
reg.klasik<-lm(Y~X1 + X2 + X3 + X4 + X5 + X6 +X7,datasulsel)
summary(reg.klasik)
## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7, data = datasulsel)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06155 -0.56936 -0.03539  0.52839  1.57532 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.159e+02  2.701e+01   4.292  0.00056 ***
## X1          -6.733e-01  3.558e-01  -1.892  0.07669 .  
## X2          -5.711e-02  5.767e-01  -0.099  0.92235    
## X3           2.340e+00  4.174e-01   5.605 3.94e-05 ***
## X4           1.988e-04  5.766e-05   3.448  0.00331 ** 
## X5          -3.337e-01  8.835e-02  -3.777  0.00165 ** 
## X6          -2.978e+01  1.373e+01  -2.168  0.04558 *  
## X7           3.620e-04  1.599e-04   2.263  0.03788 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9277 on 16 degrees of freedom
## Multiple R-squared:  0.9541, Adjusted R-squared:  0.934 
## F-statistic: 47.49 on 7 and 16 DF,  p-value: 1.582e-09
## uji asumsi klasik dari model
err.regklasik<-residuals(reg.klasik)

#heterogenitas
bptest(reg.klasik) #digunakan modelnya dalam perhitungan
## 
##  studentized Breusch-Pagan test
## 
## data:  reg.klasik
## BP = 2.6302, df = 7, p-value = 0.917
#dependensi
RunsTest(err.regklasik)
## 
##  Runs Test for Randomness
## 
## data:  err.regklasik
## runs = 10, m = 12, n = 12, p-value = 0.3009
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
##   median(x) 
## -0.03539317
#Normalitas
ad.test(err.regklasik)
## 
##  Anderson-Darling normality test
## 
## data:  err.regklasik
## A = 0.31394, p-value = 0.523
qqnorm(err.regklasik,datax=T)
qqline(rnorm(length(err.regklasik),mean(err.regklasik),sd(err.regklasik)),datax=T, col="red")

## Spatial Autocorrelation: Uji Moran Index
#Moran Test pada error
moran.test(err.regklasik, dist,randomisation=T, alternative="greater")
## 
##  Moran I test under randomisation
## 
## data:  err.regklasik  
## weights: dist    
## 
## Moran I statistic standard deviate = 0.2645, p-value = 0.3957
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.00498414       -0.04347826        0.02118052
## uji efek spasial
#Uji LM
LM<-lm.LMtests(reg.klasik, dist,test=c("LMerr", "LMlag","RLMerr","RLMlag","SARMA"))
## Please update scripts to use lm.RStests in place of lm.LMtests
summary(LM)
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7, data =
## datasulsel)
## test weights: listw
##  
##          statistic parameter p.value
## RSerr    0.0010114         1  0.9746
## RSlag    0.5969856         1  0.4397
## adjRSerr 0.1402328         1  0.7080
## adjRSlag 0.7362071         1  0.3909
## SARMA    0.7372185         2  0.6917
#Model SEM
sem<-errorsarlm(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 ,data=datasulsel, dist)
summary(sem)
## 
## Call:
## errorsarlm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7, data = datasulsel, 
##     listw = dist)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.064248 -0.570503 -0.038926  0.526611  1.570754 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  1.1560e+02  2.2064e+01  5.2393 1.612e-07
## X1          -6.6790e-01  2.9038e-01 -2.3001  0.021440
## X2          -6.6468e-02  4.6961e-01 -0.1415  0.887444
## X3           2.3447e+00  3.4037e-01  6.8886 5.634e-12
## X4           1.9811e-04  4.6987e-05  4.2163 2.484e-05
## X5          -3.3343e-01  7.1964e-02 -4.6334 3.598e-06
## X6          -2.9747e+01  1.1219e+01 -2.6515  0.008013
## X7           3.6161e-04  1.3028e-04  2.7757  0.005508
## 
## Lambda: -0.018088, LR test value: 0.002161, p-value: 0.96292
## Asymptotic standard error: 0.26693
##     z-value: -0.067765, p-value: 0.94597
## Wald statistic: 0.0045921, p-value: 0.94597
## 
## Log likelihood: -27.38695 for error model
## ML residual variance (sigma squared): 0.57366, (sigma: 0.7574)
## Number of observations: 24 
## Number of parameters estimated: 10 
## AIC: 74.774, (AIC for lm: 72.776)
# Uji Asumsi
err.sem<-residuals(sem)
ad.test(err.sem)
## 
##  Anderson-Darling normality test
## 
## data:  err.sem
## A = 0.31268, p-value = 0.5255
moran.test(err.sem,dist,randomisation=T, alternative="greater")
## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: dist    
## 
## Moran I statistic standard deviate = 0.27912, p-value = 0.3901
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.002848914      -0.043478261       0.021188896
#Model GSM/SARMA
gsm<-sacsarlm(Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 ,data=datasulsel, dist)
summary(gsm)
## 
## Call:
## sacsarlm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7, data = datasulsel, 
##     listw = dist)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.078060 -0.543789 -0.074447  0.442193  1.525376 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  1.2076e+02  2.3345e+01  5.1727 2.307e-07
## X1          -5.8819e-01  2.8558e-01 -2.0596  0.039435
## X2          -1.5068e-01  4.5410e-01 -0.3318  0.740018
## X3           2.3966e+00  3.3200e-01  7.2186 5.254e-13
## X4           2.0866e-04  4.8517e-05  4.3007 1.702e-05
## X5          -3.3486e-01  6.9758e-02 -4.8004 1.584e-06
## X6          -3.1383e+01  1.1286e+01 -2.7806  0.005425
## X7           3.3220e-04  1.2992e-04  2.5569  0.010561
## 
## Rho: -0.13408
## Asymptotic standard error: 0.13782
##     z-value: -0.97286, p-value: 0.33062
## Lambda: -0.10607
## Asymptotic standard error: 0.30354
##     z-value: -0.34945, p-value: 0.72675
## 
## LR test value: 0.842, p-value: 0.65639
## 
## Log likelihood: -26.96704 for sac model
## ML residual variance (sigma squared): 0.54982, (sigma: 0.7415)
## Number of observations: 24 
## Number of parameters estimated: 11 
## AIC: 75.934, (AIC for lm: 72.776)
#uji asumsi
err.gsm<-residuals(gsm)
ad.test(err.gsm)
## 
##  Anderson-Darling normality test
## 
## data:  err.gsm
## A = 0.45677, p-value = 0.2432
moran.test(err.gsm,dist,randomisation=T, alternative="greater")
## 
##  Moran I test under randomisation
## 
## data:  err.gsm  
## weights: dist    
## 
## Moran I statistic standard deviate = 0.19239, p-value = 0.4237
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.01550888       -0.04347826        0.02113443
## moran test
moran.test(datasulsel$Y,nb2listw(nb, style = "W")) 
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$Y  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.37756, p-value = 0.3529
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.007111038      -0.043478261       0.017953794
moran.test(datasulsel$X1,nb2listw(nb, style = "W")) 
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X1  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 1.1758, p-value = 0.1198
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.12770161       -0.04347826        0.02119626
moran.test(datasulsel$X2,nb2listw(nb, style = "W")) 
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X2  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.50943, p-value = 0.3052
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.02687401       -0.04347826        0.01907127
moran.test(datasulsel$X3,nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X3  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.64342, p-value = 0.26
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.04453481       -0.04347826        0.01871139
moran.test(datasulsel$X4,nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X4  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 1.3837, p-value = 0.08323
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.14702262       -0.04347826        0.01895530
moran.test(datasulsel$X5,nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X5  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.9175, p-value = 0.1794
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.09174598       -0.04347826        0.02172193
moran.test(datasulsel$X6,nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X6  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = -0.098225, p-value = 0.5391
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.05796377       -0.04347826        0.02174814
moran.test(datasulsel$X7,nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X7  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.79781, p-value = 0.2125
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.004067919      -0.043478261       0.002440180
## Autokorelasi spasial pada residual
# Uji autokorelasi spasial dengan Moran's I pada residual model regresi klasik
moran.test(err.regklasik, dist, randomisation = TRUE, alternative = "greater")
## 
##  Moran I test under randomisation
## 
## data:  err.regklasik  
## weights: dist    
## 
## Moran I statistic standard deviate = 0.2645, p-value = 0.3957
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.00498414       -0.04347826        0.02118052
moran.plot(err.regklasik, dist, labels=datasulsel$`Kabupaten / Kota`)

# Uji autokorelasi spasial dengan Moran's I pada residual model SEM (Spatial Error Model)
moran.test(err.sem, dist, randomisation = TRUE, alternative = "greater")
## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: dist    
## 
## Moran I statistic standard deviate = 0.27912, p-value = 0.3901
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.002848914      -0.043478261       0.021188896
moran.plot(err.sem, dist, labels=datasulsel$`Kabupaten / Kota`)

# Uji autokorelasi spasial dengan Moran's I pada residual model GSM/SARMA (Generalized Spatial Model)
moran.test(err.gsm, dist, randomisation = TRUE, alternative = "greater")
## 
##  Moran I test under randomisation
## 
## data:  err.gsm  
## weights: dist    
## 
## Moran I statistic standard deviate = 0.19239, p-value = 0.4237
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.01550888       -0.04347826        0.02113443
moran.plot(err.gsm, dist, labels=datasulsel$`Kabupaten / Kota`)

# Anda bisa menguji juga setiap variabel independen, jika diperlukan
# Moran's Test dan Plot untuk variabel X1
moran.test(datasulsel$X1, nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X1  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 1.1758, p-value = 0.1198
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.12770161       -0.04347826        0.02119626
moran.plot(datasulsel$X1, dist, labels=datasulsel$`Kabupaten / Kota`)

# Moran's Test dan Plot untuk variabel X2
moran.test(datasulsel$X2, nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X2  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.50943, p-value = 0.3052
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.02687401       -0.04347826        0.01907127
moran.plot(datasulsel$X2, dist, labels=datasulsel$`Kabupaten / Kota`)

# Moran's Test dan Plot untuk variabel X3
moran.test(datasulsel$X3, nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X3  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.64342, p-value = 0.26
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.04453481       -0.04347826        0.01871139
moran.plot(datasulsel$X3, dist, labels=datasulsel$`Kabupaten / Kota`)

# Moran's Test dan Plot untuk variabel X4
moran.test(datasulsel$X4, nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X4  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 1.3837, p-value = 0.08323
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.14702262       -0.04347826        0.01895530
moran.plot(datasulsel$X4, dist, labels=datasulsel$`Kabupaten / Kota`)

# Moran's Test dan Plot untuk variabel X5
moran.test(datasulsel$X5, nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X5  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.9175, p-value = 0.1794
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.09174598       -0.04347826        0.02172193
moran.plot(datasulsel$X5, dist, labels=datasulsel$`Kabupaten / Kota`)

# Moran's Test dan Plot untuk variabel X6
moran.test(datasulsel$X6, nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X6  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = -0.098225, p-value = 0.5391
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.05796377       -0.04347826        0.02174814
moran.plot(datasulsel$X6, mat2listw(ww, style="B"), labels = datasulsel$`Kabupaten / Kota`)

# Moran's Test dan Plot untuk variabel X7
moran.test(datasulsel$X7, nb2listw(nb, style = "W"))
## 
##  Moran I test under randomisation
## 
## data:  datasulsel$X7  
## weights: nb2listw(nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.79781, p-value = 0.2125
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.004067919      -0.043478261       0.002440180
moran.plot(datasulsel$X7, dist, labels=datasulsel$`Kabupaten / Kota`)

## peta baru (model SARMA)

coef_gsm <- coef(gsm)
# Hitung prediksi Y secara manual dari koefisien SARMA
pred_gsm <- coef_gsm[1] + coef_gsm[2]*datasulsel$X1 + coef_gsm[3]*datasulsel$X2 + 
  coef_gsm[4]*datasulsel$X3 + coef_gsm[5]*datasulsel$X4 + 
  coef_gsm[6]*datasulsel$X5 + coef_gsm[7]*datasulsel$X6 + 
  coef_gsm[8]*datasulsel$X7

# Menyimpan prediksi ke shapefile Sulsel
sulsel$pred_SARMA <- pred_gsm

# Menggunakan 4 warna untuk gradasi
colfunc <- colorRampPalette(c("green", "yellow", "orange", "red"))
color<-colfunc(16)

# Visualisasi peta prediksi SARMA dengan 4 warna
# Menampilkan peta prediksi SARMA
spplot(sulsel, "pred_SARMA", col.regions = color,
       main = "Peta Prediksi Pertumbuhan UMKM dengan Model SARMA (4 Warna)",
       sp.layout = list("sp.text", coordinates(sulsel), sulsel$Kabupaten, cex = 0.5, col = "black"))