email:

Mata Kuliah: Spatial Statistics

Dosen Pengampu: I Gede Nyoman Mindra Jaya, M.Si., Ph.D.

Library

suppressPackageStartupMessages({
  library(ggplot2)
  library(maps)
  library(raster)
  library(maptools)
  library(sp)
  library(spdep)
  library(gstat)
  library(dplyr)
  library(sf)
  library(spatialreg)
  library(lwgeom)
  library(spgwr)
  library(car)
  library(lmtest)
  library(GWmodel)
  })
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'maps' was built under R version 4.2.3
## Warning: package 'raster' was built under R version 4.2.3
## Warning: package 'sp' was built under R version 4.2.3
## Warning: package 'spdep' was built under R version 4.2.3
## Warning: package 'spData' was built under R version 4.2.3
## Warning: package 'sf' was built under R version 4.2.3
## Warning: package 'gstat' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'spatialreg' was built under R version 4.2.3
## Warning: package 'Matrix' was built under R version 4.2.3
## Warning: package 'lwgeom' was built under R version 4.2.3
## Warning: package 'spgwr' was built under R version 4.2.3
## Warning: package 'car' was built under R version 4.2.3
## Warning: package 'carData' was built under R version 4.2.2
## Warning: package 'lmtest' was built under R version 4.2.3
## Warning: package 'zoo' was built under R version 4.2.3
## Warning: package 'GWmodel' was built under R version 4.2.3
## Warning: package 'robustbase' was built under R version 4.2.3
## Warning: package 'Rcpp' was built under R version 4.2.3

Data

Input Data

Data yang digunakan adalah Produk Domestik Regional Bruto (PDRB) sebagai variabel y, Indeks Pembangunan Manusia (IPM) sebagai variabel x1, Jumlah Penduduk sebagai variabel x2, Upah Minimum Kabupaten/Kota (UMK) sebagai variabel x3, dan Investasi Penanaman Modal Asing sebagai variabel x4.

setwd("C:/Users/User/Documents/Stat/sem 5/Spatial Statistics")
data <- read.csv("data spasial.csv") ; data
##      Kabupaten.Kota         y    x1      x2      x3           x4
## 1           Bandung 153950.74 73.74 3721.11 3492466 2.968322e+12
## 2     Bandung Barat  56945.18 69.61 1859.64 3480795 5.146667e+12
## 3            Bekasi 393822.98 75.76 3237.42 5137575 4.046329e+13
## 4             Bogor 289106.15 71.78 5627.02 4520212 4.310496e+12
## 5            Ciamis  39841.04 72.05 1251.54 2021657 1.458834e+10
## 6           Cianjur  58391.22 66.55 2558.14 2893229 1.011157e+12
## 7           Cirebon  61309.76 70.95 2360.44 2430781 1.564737e+12
## 8             Garut  72229.23 68.11 2683.67 2117318 3.547528e+11
## 9         Indramayu 103559.94 69.25 1894.33 2541997 2.963209e+11
## 10         Karawang 290916.32 72.35 2526.00 5176179 4.024312e+13
## 11     Kota Bandung 351284.45 83.04 2506.60 4048463 4.176886e+12
## 12      Kota Banjar   5246.65 73.08  207.51 1998119 1.479998e+06
## 13      Kota Bekasi 118963.02 83.03 2627.21 5158248 3.433265e+12
## 14       Kota Bogor  57003.76 77.85 1070.72 4639429 2.266825e+11
## 15      Kota Cimahi  40499.39 79.46  590.78 3514093 1.288216e+12
## 16     Kota Cirebon  28772.59 76.46  341.98 2456517 9.092816e+10
## 17       Kota Depok  87568.62 82.38 2145.40 4694494 1.724696e+12
## 18    Kota Sukabumi  15349.71 76.32  360.64 2747775 6.759153e+09
## 19 Kota Tasikmalaya  27411.43 74.47  741.76 2533341 1.674914e+10
## 20         Kuningan  32544.89 70.72 1201.76 2101734 4.163267e+11
## 21       Majalengka  41713.45 69.13 1340.62 2180603 2.271929e+12
## 22      Pangandaran  14256.42 69.38  431.46 2018389 1.656562e+10
## 23       Purwakarta  82129.26 72.09 1037.07 4464675 7.598054e+12
## 24           Subang  49692.22 70.70 1649.82 3273811 9.539257e+11
## 25         Sukabumi  82317.31 68.49 2802.40 3351883 3.381206e+12
## 26         Sumedang  43327.15 73.18 1178.24 3471134 5.934779e+11
## 27      Tasikmalaya  46097.96 67.76 1907.05 2499954 3.020084e+10

Deskriptif Data

Deskriptif data dilakukan untuk melihat karakteristik dari data yang digunakan.

summary(data)
##  Kabupaten.Kota           y                x1              x2        
##  Length:27          Min.   :  5247   Min.   :66.55   Min.   : 207.5  
##  Class :character   1st Qu.: 40170   1st Qu.:69.50   1st Qu.:1053.9  
##  Mode  :character   Median : 57004   Median :72.09   Median :1859.6  
##                     Mean   : 97935   Mean   :73.25   Mean   :1846.7  
##                     3rd Qu.: 95564   3rd Qu.:76.04   3rd Qu.:2542.1  
##                     Max.   :393823   Max.   :83.04   Max.   :5627.0  
##        x3                x4           
##  Min.   :1998119   Min.   :1.480e+06  
##  1st Qu.:2443649   1st Qu.:1.588e+11  
##  Median :3273811   Median :1.011e+12  
##  Mean   :3294995   Mean   :4.541e+12  
##  3rd Qu.:4256569   3rd Qu.:3.407e+12  
##  Max.   :5176179   Max.   :4.046e+13

Peta

Peta Jawa Barat

#Plot peta Jawa Barat
Jabar<-readShapePoly("RBI_50K_2023_Jawa Barat.shp")
## Warning: shapelib support is provided by GDAL through the sf and terra packages
## among others
plot(Jabar)

Peta Persebaran PDRB

Menggabungkan Data Geospasial dengan Data PDRB

PDRB<-read.csv("pdrb harga berlaku.csv")
Jabar$id<-c(1:27)
Jabar_sf<-st_as_sf(Jabar)
Jabar_merged <- Jabar_sf %>%
  left_join(PDRB, by = "id")

Kontinu

ggplot() +
  geom_sf(data=Jabar_merged, aes(fill = pdrb),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(),  # Menghapus label x-axis
        axis.text.y = element_blank())+  # Menghapus label y-axis   
  labs(title = "Peta Persebaran PDRB Jawa Barat 2023",
       fill = "PDRB")

Diskrit

summary(PDRB$pdrb)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5247   40170   57004   97935   95564  393823
breaks <- c(-Inf, 40170, 57004, 95564,Inf)

# Memberi label pada setiap interval
labels <- c("Very Low", "Low", "High", "Very High")

# Membuat kolom baru dengan PDRB diskrit
Jabar_merged$PDRB_Discrete <- cut(Jabar_merged$pdrb, 
                                  breaks = breaks, labels = labels, 
                                  right = TRUE)

ggplot() +
  geom_sf(data=Jabar_merged, aes(fill = PDRB_Discrete),color=NA) +
  theme_bw() +
  scale_fill_manual(values = c("Very Low" = "yellow", 
                               "Low" = "orange",
                               "High" = "red",
                               "Very High" = "red3"))+
  labs(fill = "PDRB")+theme(legend.position = "right",
                          axis.text.x = element_blank(),  # Menghapus label x-axis
                          axis.text.y = element_blank())+  # Menghapus label y-axis
  labs(title = "Peta Persebaran PDRB Jawa Barat 2023",
       fill = "PDRB")

Regresi Klasik

Ordinary Least Square (OLS)

reglin <- lm(y ~ x1 + x2 + x3 + x4, data=data)
summary(reglin)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -64624 -18189  -6412  10225 146336 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.309e+05  1.685e+05  -3.744  0.00113 ** 
## x1           9.139e+03  2.614e+03   3.496  0.00204 ** 
## x2           5.376e+01  8.907e+00   6.036 4.48e-06 ***
## x3          -2.078e-02  1.491e-02  -1.394  0.17723    
## x4           6.305e-09  1.084e-09   5.816 7.51e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43430 on 22 degrees of freedom
## Multiple R-squared:  0.8574, Adjusted R-squared:  0.8314 
## F-statistic: 33.06 on 4 and 22 DF,  p-value: 5.182e-09

Dari hasil analisis OLS, didapatkan bahwa variabel x1, x2, dan x4 signifikan. Selain itu, nilai R-square sebesar 0.8314 menunjukkan bahwa sekitar 83.14% varians PDRB dapat dijelaskan oleh variabel prediktor. Kemudian dilakukan uji asumsi untuk residual dari model OLS.

Uji Asumsi

err.ols <- residuals(reglin)
vif(reglin)             # Multicollinearity test
##       x1       x2       x3       x4 
## 2.070497 1.617534 3.587403 1.787092
shapiro.test(err.ols)   # Saphiro test for normality of residuals
## 
##  Shapiro-Wilk normality test
## 
## data:  err.ols
## W = 0.85124, p-value = 0.001219
bptest(reglin)          # Breusch-Pagan test for heteroskedasticity
## 
##  studentized Breusch-Pagan test
## 
## data:  reglin
## BP = 9.9762, df = 4, p-value = 0.04083
durbinWatsonTest(reglin)# Non-autocorrelation
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.05697496      2.082851   0.948
##  Alternative hypothesis: rho != 0

Dari hasil uji asumsi diatas, didapatkan bahwa nilai VIF semua variabel x kurang dari 5 sehingga memenuhi uji asumsi multikolinearitas. Namun, dari hasil uji normalitas, didapatkan p-value sebesar 0.001219 yang artinya data tidak berdistribusi normal. Oleh karena itu, dilakukan transformasi data menggunakan log.

Transformasi Data

# Transformasi log pada variabel
y_trans = log(data$y)
x1_trans = log(data$x1)
x2_trans = log(data$x2)
x3_trans = log(data$x3)
x4_trans = log(data$x4)
data_trans = data.frame(y_trans, x1_trans, x2_trans, x3_trans, x4_trans)

# Model regresi dengan variabel yang ditransformasi
ols_trans <- lm(y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans,
                  data = data_trans)
summary(ols_trans)
## 
## Call:
## lm(formula = y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans, 
##     data = data_trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56766 -0.24565 -0.01751  0.16982  0.69605 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -15.08852    5.33075  -2.830  0.00974 ** 
## x1_trans      2.55838    1.60921   1.590  0.12614    
## x2_trans      0.72878    0.14423   5.053 4.63e-05 ***
## x3_trans      0.51928    0.39318   1.321  0.20017    
## x4_trans      0.07747    0.03720   2.082  0.04914 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3709 on 22 degrees of freedom
## Multiple R-squared:  0.8822, Adjusted R-squared:  0.8608 
## F-statistic:  41.2 on 4 and 22 DF,  p-value: 6.48e-10

Dari hasil OLS pada data yang sudah ditransformasi, didapatkan bahwa variabel x2 dan x4 signifikan. Selain itu, nilai R-square sebesar 0.8608 menunjukkan bahwa sekitar 86.08% variasi PDRB dapat dijelaskan oleh variabel prediktor. Kemudian, dilakukan kembali uji asumsi pada residual model.

err.trans <- residuals(ols_trans)

# Uji Asumsi
vif(ols_trans)              # Multicollinearity test
## x1_trans x2_trans x3_trans x4_trans 
## 1.923126 2.557573 3.101231 3.074185
shapiro.test(err.trans)     # Saphiro test for normality of residuals
## 
##  Shapiro-Wilk normality test
## 
## data:  err.trans
## W = 0.9553, p-value = 0.2874
bptest(ols_trans)           # Breusch-Pagan test for heteroskedasticity
## 
##  studentized Breusch-Pagan test
## 
## data:  ols_trans
## BP = 11.148, df = 4, p-value = 0.02495
durbinWatsonTest(ols_trans) # Non-autocorrelation
##  lag Autocorrelation D-W Statistic p-value
##    1       0.2493034      1.501267   0.112
##  Alternative hypothesis: rho != 0

Dari hasil uji asumsi diatas, didapatkan bahwa seluruh variabel x memiliki nilai VIF kurang dari 5. Selain itu, didapatkan p-value untuk uji normalitas dan autokorelasi lebih besar dari alpha (0.05) sehingga data sudah berdistribusi normal dan data tidak berkorelasi. Namun, nilai p-value dari uji heterokedastisitas sebesar 0.02495 yang dibawah 0.05 menandakan bahwa terdapat variasi dari kabupaten/kota di Jawa Barat secara spasial, sehingga dilakukan analisis yang mempertimbangkan efek spasial, yaitu analisis GWR.

Autokorelasi Spasial

Untuk melihat autokorelasi spasial, kita dapat menggunakan nilai Moran’s I dengan memilih matriks bobot yang memiliki p-value signifikan atau nilai Moran’s I terbesar.

ID<-c(1:27)
Jabar$ID <-c(1:27)
coords<-coordinates(Jabar) #membuat koordinat centroid

Rook

WR <- poly2nb(Jabar, row.names=ID,
             queen=FALSE)
WBR <- nb2mat(WR, style='B', 
              zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"
WBR[c(1:5),c(1:5)]
##   [,1] [,2] [,3] [,4] [,5]
## 1    0    1    0    0    0
## 2    1    0    0    0    0
## 3    0    0    0    1    0
## 4    0    0    1    0    0
## 5    0    0    0    0    0
WLR<-nb2listw(WR) ;WLR #List neighbours
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 106 
## Percentage nonzero weights: 14.54047 
## Average number of links: 3.925926 
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 27 729 27 16.73552 121.5371

Queen

WQ <- poly2nb(Jabar, row.names=ID, 
             queen=TRUE) #Mendapatkan W
WBQ <- nb2mat(WQ, style='B', 
              zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"
WBQ[1:5,1:5]
##   [,1] [,2] [,3] [,4] [,5]
## 1    0    1    0    0    0
## 2    1    0    0    0    0
## 3    0    0    0    1    0
## 4    0    0    1    0    0
## 5    0    0    0    0    0
WLQ<-nb2listw(WQ);WLQ
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 106 
## Percentage nonzero weights: 14.54047 
## Average number of links: 3.925926 
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 27 729 27 16.73552 121.5371

K-Nearest Neighbours (KNN)

#k=5
WK <- knn2nb(knearneigh(coords, k = 5), row.names = ID); WK
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 135 
## Percentage nonzero weights: 18.51852 
## Average number of links: 5 
## Non-symmetric neighbours list
WBK <- nb2mat(WK, style='B', 
              zero.policy = TRUE) #menyajikan dalam bentuk matrix biner "B"
WBK[1:5,1:5]
##   [,1] [,2] [,3] [,4] [,5]
## 1    0    1    0    0    0
## 2    1    0    0    0    0
## 3    0    0    0    1    0
## 4    0    0    1    0    0
## 5    0    0    0    0    0
WLK<-nb2listw(WK) ; WLK
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 135 
## Percentage nonzero weights: 18.51852 
## Average number of links: 5 
## Non-symmetric neighbours list
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0    S1     S2
## W 27 729 27 10.28 109.36

Inverse Distance

#Matriks jarak
Dist.mat<-as.matrix(dist(coords, method="euclidean"))
Dist.mat[1:5,1:5]
##           0         1         2         3         4
## 0 0.0000000 0.2836346 1.0038865 0.9971936 0.8422494
## 1 0.2836346 0.0000000 0.7320059 0.7258399 1.0896619
## 2 1.0038865 0.7320059 0.0000000 0.4913495 1.6893934
## 3 0.9971936 0.7258399 0.4913495 0.0000000 1.8135196
## 4 0.8422494 1.0896619 1.6893934 1.8135196 0.0000000
#Inverse Distance
Dist.mat.inv<-1/Dist.mat
diag(Dist.mat.inv)<-0
Dist.mat.invs<-mat2listw(Dist.mat.inv, style="W")
summary(Dist.mat.invs)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 27 
## Number of nonzero links: 702 
## Percentage nonzero weights: 96.2963 
## Average number of links: 26 
## Link number distribution:
## 
## 26 
## 27 
## 27 least connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 with 26 links
## 27 most connected regions:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 with 26 links
## 
## Weights style: W 
## Weights constants summary:
##    n  nn S0       S1       S2
## W 27 729 27 6.455601 109.5919

Moran’s I

#Menghitung Indeks Moran tiap matriks bobot
moran.test(y_trans,WLR) #Rook
## 
##  Moran I test under randomisation
## 
## data:  y_trans  
## weights: WLR    
## 
## Moran I statistic standard deviate = 2.5015, p-value = 0.006183
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.30678806       -0.03846154        0.01904868
moran.test(y_trans,WLQ) #Queen
## 
##  Moran I test under randomisation
## 
## data:  y_trans  
## weights: WLQ    
## 
## Moran I statistic standard deviate = 2.5015, p-value = 0.006183
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.30678806       -0.03846154        0.01904868
moran.test(y_trans,WLK) #KNN (k=5)
## 
##  Moran I test under randomisation
## 
## data:  y_trans  
## weights: WLK    
## 
## Moran I statistic standard deviate = 3.635, p-value = 0.000139
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.34212052       -0.03846154        0.01096204
moran.test(y_trans,Dist.mat.invs) #Inverse Distance
## 
##  Moran I test under randomisation
## 
## data:  y_trans  
## weights: Dist.mat.invs    
## 
## Moran I statistic standard deviate = 2.0468, p-value = 0.02034
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.117539052      -0.038461538       0.005809265

Dari hasil diatas, didapatkan nilai Moran’s I terbesar oleh matriks bobot KNN dengan k=5 sehingga dipilih matriks bobot KNN dengan k=5. Selain itu, hasil uji Moran’s I dari matriks bobot KNN (k=5) yang signifikan menunjukkan bahwa terdapat autokorelasi spasial pada kabupaten/kota di Jawa Barat.

#Membuat plot Moran's I
moran.plot(y_trans,WLK)

Local Moran’s I

locm<-localmoran(y_trans,WLK)
locm
##              Ii          E.Ii       Var.Ii       Z.Ii Pr(z != E(Ii))
## 1   0.257923195 -3.380441e-02 1.481533e-01  0.7579173    0.448500497
## 2  -0.041384415 -2.579599e-04 1.169804e-03 -1.2024436    0.229191703
## 3   1.709890904 -1.388825e-01 5.424790e-01  2.5101068    0.012069465
## 4   0.451648936 -9.643094e-02 3.952308e-01  0.8718032    0.383315746
## 5   0.520819587 -7.719897e-03 3.474712e-02  2.8354240    0.004576490
## 6   0.007928303 -1.214473e-04 5.508179e-04  0.3429876    0.731607796
## 7   0.002176063 -1.478588e-06 6.706867e-06  0.8408270    0.400444882
## 8   0.037524444 -1.006919e-03 4.562786e-03  0.5704260    0.568388789
## 9  -0.188262908 -1.084927e-02 4.867837e-02 -0.8041166    0.421329665
## 10  0.955351320 -9.721176e-02 3.980868e-01  1.6682428    0.095267542
## 11  0.126890853 -1.222845e-01 4.868535e-01  0.3571133    0.721007020
## 12  1.886687793 -2.454090e-01 8.399919e-01  2.3263195    0.020001515
## 13  0.720613696 -1.743281e-02 7.769672e-02  2.6477830    0.008102152
## 14 -0.024026297 -2.513635e-04 1.139898e-03 -0.7041842    0.481318046
## 15 -0.221655344 -7.151813e-03 3.220861e-02 -1.1952216    0.232000466
## 16  0.153103030 -2.349753e-02 1.040803e-01  0.5474037    0.584101391
## 17  0.190401649 -4.962294e-03 2.239727e-02  1.3054098    0.191753358
## 18 -0.473400116 -7.817271e-02 3.268720e-01 -0.6912867    0.489385408
## 19  0.766804879 -2.657911e-02 1.173584e-01  2.3159352    0.020561809
## 20  0.545612455 -1.651848e-02 7.369013e-02  2.0707746    0.038379867
## 21  0.102167391 -6.182833e-03 2.787193e-02  0.6490024    0.516336833
## 22  1.429068401 -8.669896e-02 3.591707e-01  2.5291938    0.011432488
## 23  0.154844951 -3.312304e-03 1.497484e-02  1.2924328    0.196207317
## 24 -0.071829853 -1.887721e-03 8.546540e-03 -0.7565606    0.449313170
## 25  0.022476488 -3.365434e-03 1.521423e-02  0.2095076    0.834052005
## 26 -0.091269461 -5.041170e-03 2.275147e-02 -0.5716695    0.567545882
## 27  0.307147978 -3.426962e-03 1.549143e-02  2.4952880    0.012585495
## attr(,"call")
## localmoran(x = y_trans, listw = WLK)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##         mean    median     pysal
## 1  High-High High-High High-High
## 2   Low-High  Low-High  Low-High
## 3  High-High High-High High-High
## 4  High-High High-High High-High
## 5    Low-Low   Low-Low   Low-Low
## 6    Low-Low  High-Low   Low-Low
## 7    Low-Low  High-Low   Low-Low
## 8  High-High  High-Low High-High
## 9   High-Low  High-Low  High-Low
## 10 High-High High-High High-High
## 11 High-High  High-Low High-High
## 12   Low-Low   Low-Low   Low-Low
## 13 High-High High-High High-High
## 14  Low-High  Low-High  Low-High
## 15  Low-High  Low-High  Low-High
## 16   Low-Low   Low-Low   Low-Low
## 17 High-High High-High High-High
## 18  Low-High  Low-High  Low-High
## 19   Low-Low   Low-Low   Low-Low
## 20   Low-Low   Low-Low   Low-Low
## 21   Low-Low   Low-Low   Low-Low
## 22   Low-Low   Low-Low   Low-Low
## 23 High-High High-High High-High
## 24  Low-High  Low-High  Low-High
## 25 High-High  High-Low High-High
## 26  Low-High  Low-High  Low-High
## 27   Low-Low   Low-Low   Low-Low

Pada hasil uji Moran’s I secara lokal, didapatkan bahwa lokasi 3 (Kabupaten Bekasi), 5 (Kabupaten Ciamis), 12 (Kota Banjar), 13 (Kota Bekasi), 19 (Kota Tasikmalaya), 20 (Kabupaten Kuningan), 22 (Kabupaten Pangandaran), dan 27 (Kabupaten Tasikmalaya) signifikan secara lokal yang dapat dilihat dari nilai p-value nya yang lebih kecil dari 0.05. Lokasi 3 (Kabupaten Bekasi) dan 13 (Kota Bekasi) signifikan secara positif (High-High), sedangkan lokasi lainnya signifikan secara negatif (Low-Low).

GWR

Pemodelan GWR dapat dilakukan dengan menggunakan berbagai kernel, kemudian dipilih kernel yang menghasilkan model GWR terbaik. Pada analisis ini, digunakan kernel Fixed Gaussian, Adaptive Gaussian, dan Bi-square untuk membuat pemodelan GWR.

#Konvert data spasial
data.sp<-data_trans
coordinates(data.sp) <- coords

Fixed Gaussian

#Bandwith
bw.fg <- gwr.sel(y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans, 
                data=data.sp)
## Bandwidth: 0.8934582 CV score: 4.210253 
## Bandwidth: 1.444202 CV score: 4.494197 
## Bandwidth: 0.5530796 CV score: 4.231259 
## Bandwidth: 0.7709021 CV score: 4.152226 
## Bandwidth: 0.7358279 CV score: 4.142287 
## Bandwidth: 0.6660243 CV score: 4.140491 
## Bandwidth: 0.6956871 CV score: 4.137611 
## Bandwidth: 0.6967208 CV score: 4.137619 
## Bandwidth: 0.6950549 CV score: 4.137609 
## Bandwidth: 0.6949504 CV score: 4.137609 
## Bandwidth: 0.6949911 CV score: 4.137609 
## Bandwidth: 0.6949097 CV score: 4.137609 
## Bandwidth: 0.6949504 CV score: 4.137609
bw.fg
## [1] 0.6949504
#Model GWR
gwr.fg <- gwr(y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans, 
              data=data.sp, bandwidth = bw.fg,
              hatmatrix=TRUE, se.fit=TRUE)
gwr.fg
## Call:
## gwr(formula = y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans, 
##     data = data.sp, bandwidth = bw.fg, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Fixed bandwidth: 0.6949504 
## Summary of GWR coefficient estimates at data points:
##                    Min.    1st Qu.     Median    3rd Qu.       Max.   Global
## X.Intercept. -24.290112 -23.545905 -19.974348 -15.866137 -11.497890 -15.0885
## x1_trans       1.043983   2.823517   3.838368   4.952251   5.378982   2.5584
## x2_trans       0.636045   0.729941   0.790811   0.819540   0.868241   0.7288
## x3_trans       0.266458   0.330826   0.415180   0.468025   0.645028   0.5193
## x4_trans       0.054796   0.058582   0.076390   0.118178   0.156336   0.0775
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 10.62476 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 16.37524 
## Sigma (residual: 2traceS - traceS'S): 0.3414295 
## Effective number of parameters (model: traceS): 8.646768 
## Effective degrees of freedom (model: traceS): 18.35323 
## Sigma (model: traceS): 0.3225067 
## Sigma (ML): 0.2658968 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 36.94633 
## AIC (GWR p. 96, eq. 4.22): 13.73851 
## Residual sum of squares: 1.90893 
## Quasi-global R2: 0.9257207

Dari model GWR dengan Adaptive Gaussian, didapatkan nilai bandwith sebesar 0.6949504 dengan CV score 4.137049, dan didapatkan nilai AIC sebesar 13.73851 dengan R-square besar 0.9257207.

Adaptive Gaussian

#Bandwith
bw.ag <- gwr.sel(y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans,
                 data=data.sp, adapt=T)
## Adaptive q: 0.381966 CV score: 4.159289 
## Adaptive q: 0.618034 CV score: 4.313665 
## Adaptive q: 0.236068 CV score: 4.127492 
## Adaptive q: 0.2135536 CV score: 4.172644 
## Adaptive q: 0.3007631 CV score: 4.146316 
## Adaptive q: 0.2607793 CV score: 4.09379 
## Adaptive q: 0.2649003 CV score: 4.099304 
## Adaptive q: 0.2557008 CV score: 4.095478 
## Adaptive q: 0.2591551 CV score: 4.09195 
## Adaptive q: 0.2586316 CV score: 4.092452 
## Adaptive q: 0.2593859 CV score: 4.09201 
## Adaptive q: 0.2591958 CV score: 4.091911 
## Adaptive q: 0.2592498 CV score: 4.09186 
## Adaptive q: 0.2593018 CV score: 4.091905 
## Adaptive q: 0.2592498 CV score: 4.09186
bw.ag
## [1] 0.2592498
#Model GWR
gwr.ag <- gwr(y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans,
              data=data.sp, adapt=bw.ag,
              hatmatrix=TRUE, se.fit=TRUE)
gwr.ag
## Call:
## gwr(formula = y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans, 
##     data = data.sp, adapt = bw.ag, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Adaptive quantile: 0.2592498 (about 6 of 27 data points)
## Summary of GWR coefficient estimates at data points:
##                    Min.    1st Qu.     Median    3rd Qu.       Max.   Global
## X.Intercept. -25.270172 -24.448243 -22.902531 -15.824516 -11.338237 -15.0885
## x1_trans       0.737403   2.829374   5.010015   5.378364   5.780349   2.5584
## x2_trans       0.618861   0.731023   0.803543   0.844025   0.907612   0.7288
## x3_trans       0.209288   0.285916   0.380861   0.457725   0.724769   0.5193
## x4_trans       0.052851   0.058126   0.075181   0.133335   0.162513   0.0775
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 11.84455 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 15.15545 
## Sigma (residual: 2traceS - traceS'S): 0.333029 
## Effective number of parameters (model: traceS): 9.611612 
## Effective degrees of freedom (model: traceS): 17.38839 
## Sigma (model: traceS): 0.3109115 
## Sigma (ML): 0.2495081 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 38.89405 
## AIC (GWR p. 96, eq. 4.22): 11.26803 
## Residual sum of squares: 1.680865 
## Quasi-global R2: 0.934595

Dari model GWR dengan Adaptive Gaussian, didapatkan nilai bandwith sebesar 0.2592498 dengan CV score 4.09186, dan didapatkan nilai AIC sebesar 11.26803 dengan R-square besar 0.934595.

Bi-square

#Bandwith
bw.bi <- gwr.sel(y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans,
                 data=data.sp, gweight=gwr.bisquare)
## Bandwidth: 0.8934582 CV score: 8.120757 
## Bandwidth: 1.444202 CV score: 4.758151 
## Bandwidth: 1.784581 CV score: 4.433669 
## Bandwidth: 1.696832 CV score: 4.52975 
## Bandwidth: 1.994947 CV score: 4.311958 
## Bandwidth: 2.056773 CV score: 4.307671 
## Bandwidth: 2.044389 CV score: 4.308087 
## Bandwidth: 2.073824 CV score: 4.307374 
## Bandwidth: 2.173709 CV score: 4.310761 
## Bandwidth: 2.081143 CV score: 4.307337 
## Bandwidth: 2.082486 CV score: 4.307336 
## Bandwidth: 2.082428 CV score: 4.307336 
## Bandwidth: 2.082387 CV score: 4.307336 
## Bandwidth: 2.082428 CV score: 4.307336
bw.bi
## [1] 2.082428
#Model GWR
gwr.bi <- gwr(y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans,
              data=data.sp, gweight=gwr.bisquare,
              bandwidth = bw.bi, hatmatrix=TRUE)
gwr.bi
## Call:
## gwr(formula = y_trans ~ x1_trans + x2_trans + x3_trans + x4_trans, 
##     data = data.sp, bandwidth = bw.bi, gweight = gwr.bisquare, 
##     hatmatrix = TRUE)
## Kernel function: gwr.bisquare 
## Fixed bandwidth: 2.082428 
## Summary of GWR coefficient estimates at data points:
##                    Min.    1st Qu.     Median    3rd Qu.       Max.   Global
## X.Intercept. -25.368526 -23.105185 -17.528975 -15.280609 -13.112879 -15.0885
## x1_trans       2.079578   2.758271   3.331637   4.531751   5.750618   2.5584
## x2_trans       0.681993   0.743344   0.775315   0.799828   0.842552   0.7288
## x3_trans       0.268953   0.382490   0.414897   0.462624   0.531917   0.5193
## x4_trans       0.054774   0.059364   0.077081   0.112921   0.155629   0.0775
## Number of data points: 27 
## Effective number of parameters (residual: 2traceS - traceS'S): 8.43506 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 18.56494 
## Sigma (residual: 2traceS - traceS'S): 0.3533185 
## Effective number of parameters (model: traceS): 7.19463 
## Effective degrees of freedom (model: traceS): 19.80537 
## Sigma (model: traceS): 0.3420752 
## Sigma (ML): 0.2929754 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 35.18131 
## AIC (GWR p. 96, eq. 4.22): 17.52332 
## Residual sum of squares: 2.317534 
## Quasi-global R2: 0.9098212

Dari model GWR dengan Bi-square, didapatkan nilai bandwith sebesar 2.082428 dengan CV score 4.307336, dan didapatkan nilai AIC sebesar 17.52332 dengan R-square besar 0.9098212.

Model GWR

Untuk menentukan fungsi kernel yang memberikan model GWR terbaik dapat dilihat dari nilai Cross-Validation (CV), dimana nilai CV didapatkan dari penentuan nilai bandwith.

#Perbandingan 
kernel_data <- data.frame(
  Kernel = c("Fixed Gaussian", "Adaptive Gaussian", "Bi-square"),
  Bandwidth = c(bw.fg, bw.ag, bw.bi),
  CV_Score = c(4.137609, 4.09186, 4.307336)
)
kernel_data
##              Kernel Bandwidth CV_Score
## 1    Fixed Gaussian 0.6949504 4.137609
## 2 Adaptive Gaussian 0.2592498 4.091860
## 3         Bi-square 2.0824278 4.307336

Dari hasil perbandingan, didapatkan bahwa kernel Adaptive Gaussian memiliki nilai CV paling rendah, sehingga kernel terbaik untuk model GWR adalah Adaptive Gaussian.

# Memilih GWR dengan Kernel Adaptive Gaussian sebagai model GWR terbaik
model_gwr <- as.data.frame(gwr.ag$SDF)
head(model_gwr)
##       sum.w X.Intercept. x1_trans  x2_trans  x3_trans   x4_trans
## 1 11.553189    -24.18224 5.330128 0.9014073 0.2678061 0.06621268
## 2 11.628671    -21.14387 4.438276 0.8271203 0.2763216 0.10931242
## 3 10.972877    -11.35792 1.225085 0.6546603 0.5293555 0.16251322
## 4 10.681044    -12.80735 1.221347 0.6400496 0.6663721 0.14359040
## 5  9.946806    -24.45749 5.598711 0.8204582 0.2663583 0.05679727
## 6 11.075312    -19.44412 4.204185 0.8113370 0.2476483 0.10332172
##   X.Intercept._se x1_trans_se x2_trans_se x3_trans_se x4_trans_se       gwr.e
## 1        5.556757    1.579673   0.1448292   0.3759599  0.03644430 -0.14336170
## 2        5.059874    1.548637   0.1461335   0.4102873  0.04653511 -0.32582945
## 3        5.260473    1.562932   0.1508875   0.4568047  0.05627121  0.37720164
## 4        5.701016    1.635048   0.1580837   0.5386595  0.06341549  0.24642586
## 5        7.282240    2.135494   0.1642611   0.5015814  0.03595363  0.05419554
## 6        5.149545    1.512876   0.1422136   0.4255514  0.04733511 -0.13703976
##       pred   pred.se   localR2 X.Intercept._se_EDF x1_trans_se_EDF
## 1 12.08775 0.1315209 0.9268561            5.952051        1.692047
## 2 11.27567 0.1173041 0.9160320            5.419821        1.658803
## 3 12.50646 0.1399716 0.8959979            5.634690        1.674115
## 4 12.32812 0.1933138 0.9113134            6.106572        1.751361
## 5 10.53846 0.1572440 0.9512875            7.800280        2.287408
## 6 11.11196 0.1364504 0.9268461            5.515871        1.620498
##   x2_trans_se_EDF x3_trans_se_EDF x4_trans_se_EDF pred.se.1 coords.x1 coords.x2
## 1       0.1551320       0.4027048      0.03903685 0.1408770  107.6087 -7.097588
## 2       0.1565290       0.4394741      0.04984551 0.1256488  107.4150 -6.890438
## 3       0.1616213       0.4893007      0.06027420 0.1499289  107.1217 -6.219747
## 4       0.1693294       0.5769784      0.06792671 0.2070657  106.7682 -6.561007
## 5       0.1759462       0.5372627      0.03851129 0.1684300  108.4286 -7.290250
## 6       0.1523304       0.4558241      0.05070241 0.1461571  107.1592 -7.130889

Perbandingan OLS dengan GWR

data.frame("Model" = c("OLS","GWR"),
           "AIC" = c(AIC(ols_trans),gwr.ag$results$AICh),
           "R2" = c(0.8314, 0.934595))%>% arrange(AIC)
##   Model      AIC       R2
## 1   GWR 11.26803 0.934595
## 2   OLS 29.53951 0.831400

Dari hasil perbandingan model OLS dengan GWR, ditunjukkan bahwa nilai AIC model GWR lebih kecil dari model OLS dengan koefisien determinasi model GWR sebesar 93.45% yang menandakan bahwa model GWR lebih baik dibandingkan dengan model OLS.

Peta Sebaran Penduga Parameter

Peta sebaran penduga parameter untuk variabel x menunjukkan kemiripan antar kabupaten/kota yang berdekatan. Kemiripan ini ditunjukkan dengan pola kecenderungan warna-warna yang sama untuk mengelompok dengan wilayah tetangganya.

#beta1(x1) 
b.x1<- gwr.ag$SDF$x1_trans 
Jabar@data$b.x1<- b.x1  
spplot(Jabar, zcol="b.x1",        
       main="Peta Sebaran Penduga Parameter beta1 (x1)")

Berdasarkan peta sebaran penduga parameter di atas dapat diketahui bahwa seluruh nilai penduga parameter untuk X1 memiliki nilai positif. Sehingga dapat diartikan bahwa IPM berkontribusi positif terhadap PDRB di seluruh kabupaten/kota di Jawa Barat.

#beta2(x2)} 
b.x2<- gwr.ag$SDF$x2_trans  
Jabar@data$b.x2<- b.x2  
spplot(Jabar, zcol="b.x2",                
       main="Peta Sebaran Penduga Parameter beta2 (x2)")

Berdasarkan peta sebaran penduga parameter di atas dapat diketahui bahwa seluruh nilai penduga parameter untuk X2 memiliki nilai positif. Sehingga dapat diartikan bahwa jumlah penduduk berkontribusi positif terhadap PDRB di seluruh kabupaten/kota di Jawa Barat.

#beta2(x3)} 
b.x3<- gwr.ag$SDF$x3_trans  
Jabar@data$b.x3<- b.x3 
spplot(Jabar, zcol="b.x3",                
       main="Peta Sebaran Penduga Parameter beta3 (x3)")

Berdasarkan peta sebaran penduga parameter di atas dapat diketahui bahwa seluruh nilai penduga parameter untuk X3 memiliki nilai positif. Sehingga dapat diartikan bahwa Upah Minimum Kabupaten/Kota (UMK) berkontribusi positif terhadap PDRB di seluruh kabupaten/kota di Jawa Barat.

#beta4(x4) 
b.x4<- gwr.ag$SDF$x4_trans 
Jabar@data$b.x4<- b.x4 
spplot(Jabar, zcol="b.x4",        
       main="Peta Sebaran Penduga Parameter beta4 (x4)")

Berdasarkan peta sebaran penduga parameter di atas dapat diketahui bahwa seluruh nilai penduga parameter untuk X4 memiliki nilai positif. Sehingga dapat diartikan bahwa Investasi Penanaman Modal Asing berkontribusi positif terhadap PDRB di seluruh kabupaten/kota di Jawa Barat.

Uji F

Uji F dilakukan untuk melihat signifikansi estimasi parameter terhadap model.

# Uji signifikan secara global
LMZ.F3GWR.test(gwr.ag)
## 
## Leung et al. (2000) F(3) test
## 
##             F statistic Numerator d.f. Denominator d.f.   Pr(>)  
## (Intercept)     1.88224       11.79841           18.595 0.10756  
## x1_trans        3.14062       13.29632           18.595 0.01199 *
## x2_trans        0.89824       10.48150           18.595 0.55595  
## x3_trans        0.24043        9.65967           18.595 0.98600  
## x4_trans        1.67552        9.45965           18.595 0.16362  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari hasil diatas, dapat disimpulkan bahwa hanya variabel X1(IPM) yang berpengaruh secara signifikan terhadap model. Kemudian dilakukan uji T untuk melihat signifikansi variabel prediktor secara lokal.

Uji T

# Uji signifikan secara lokal
t1 = gwr.ag$SDF$x1_trans / gwr.ag$SDF$x1_trans_se
t2 = gwr.ag$SDF$x2_trans / gwr.ag$SDF$x2_trans_se
t3 = gwr.ag$SDF$x3_trans / gwr.ag$SDF$x3_trans_se
t4 = gwr.ag$SDF$x4_trans / gwr.ag$SDF$x4_trans_se

pvalue1 = 2*pt(q=t1, df=22, lower.tail=FALSE)
pvalue2 = 2*pt(q=t2, df=22, lower.tail=FALSE)
pvalue3 = 2*pt(q=t3, df=22, lower.tail=FALSE)
pvalue4 = 2*pt(q=t4, df=22, lower.tail=FALSE)
#Variabel x1: Indeks Pembangunan Manusia (IPM)
data.frame(data$Kabupaten.Kota, t=t1, pvalue=pvalue1) %>% 
  mutate(Keterangan=ifelse(pvalue1<0.05,"Reject", "Not Reject"))
##    data.Kabupaten.Kota         t      pvalue Keterangan
## 1              Bandung 3.3741969 0.002734504     Reject
## 2        Bandung Barat 2.8659247 0.008979802     Reject
## 3               Bekasi 0.7838373 0.441493579 Not Reject
## 4                Bogor 0.7469796 0.462985130 Not Reject
## 5               Ciamis 2.6217406 0.015575924     Reject
## 6              Cianjur 2.7789362 0.010946205     Reject
## 7              Cirebon 2.8301884 0.009743055     Reject
## 8                Garut 3.5376325 0.001849215     Reject
## 9            Indramayu 2.5403934 0.018642571     Reject
## 10            Karawang 1.2324967 0.230773179 Not Reject
## 11        Kota Bandung 3.2603907 0.003582998     Reject
## 12         Kota Banjar 2.5178176 0.019588855     Reject
## 13         Kota Bekasi 0.6044919 0.551700619 Not Reject
## 14          Kota Bogor 0.7655324 0.452089893 Not Reject
## 15         Kota Cimahi 3.1269047 0.004906266     Reject
## 16        Kota Cirebon 2.8318656 0.009705891     Reject
## 17          Kota Depok 0.4357332 0.667275819 Not Reject
## 18       Kota Sukabumi 1.9066723 0.069707553 Not Reject
## 19    Kota Tasikmalaya 2.9900408 0.006748525     Reject
## 20            Kuningan 2.7990392 0.010458355     Reject
## 21          Majalengka 3.1776378 0.004355417     Reject
## 22         Pangandaran 2.9207370 0.007918966     Reject
## 23          Purwakarta 2.1678785 0.041260264     Reject
## 24              Subang 2.3686939 0.027052876     Reject
## 25            Sukabumi 1.8889706 0.072149862 Not Reject
## 26            Sumedang 3.1949864 0.004181165     Reject
## 27         Tasikmalaya 3.2114314 0.004022237     Reject

Dilihat dari variabel IPM (X1), dapat disimpulkan bahwa IPM di Kabupaten Bekasi, Kabupaten Bogor, Kabupaten Karawang, Kota Bekasi, Kota Bogor, Kota Depok, Kota Sukabumi, dan Kabupaten Sukabumi tidak mempengaruhi PDRB di kabupaten/kota tersebut.

#Variabel x2: Jumlah Penduduk
data.frame(data$Kabupaten.Kota, t=t2, pvalue=pvalue2) %>% 
  mutate(Keterangan=ifelse(pvalue2<0.05,"Reject", "Not Reject"))
##    data.Kabupaten.Kota        t       pvalue Keterangan
## 1              Bandung 6.223933 2.899482e-06     Reject
## 2        Bandung Barat 5.660033 1.083307e-05     Reject
## 3               Bekasi 4.338730 2.637658e-04     Reject
## 4                Bogor 4.048801 5.354782e-04     Reject
## 5               Ciamis 4.994842 5.334624e-05     Reject
## 6              Cianjur 5.705058 9.738780e-06     Reject
## 7              Cirebon 4.905789 6.620115e-05     Reject
## 8                Garut 6.460037 1.687902e-06     Reject
## 9            Indramayu 5.832992 7.204209e-06     Reject
## 10            Karawang 4.883685 6.985023e-05     Reject
## 11        Kota Bandung 5.909066 6.026882e-06     Reject
## 12         Kota Banjar 4.813530 8.283229e-05     Reject
## 13         Kota Bekasi 4.071446 5.066980e-04     Reject
## 14          Kota Bogor 4.023149 5.700476e-04     Reject
## 15         Kota Cimahi 5.677315 1.039895e-05     Reject
## 16        Kota Cirebon 4.922171 6.362123e-05     Reject
## 17          Kota Depok 3.814907 9.462246e-04     Reject
## 18       Kota Sukabumi 4.886840 6.931731e-05     Reject
## 19    Kota Tasikmalaya 5.617811 1.197290e-05     Reject
## 20            Kuningan 4.899075 6.728869e-05     Reject
## 21          Majalengka 5.611825 1.214408e-05     Reject
## 22         Pangandaran 5.528980 1.478606e-05     Reject
## 23          Purwakarta 5.247677 2.898153e-05     Reject
## 24              Subang 5.781439 8.133045e-06     Reject
## 25            Sukabumi 5.297772 2.569579e-05     Reject
## 26            Sumedang 6.298714 2.441084e-06     Reject
## 27         Tasikmalaya 5.925240 5.802997e-06     Reject

Dilihat dari variabel Jumlah Penduduk(X2), dapat disimpulkan bahwa jumlah penduduk pada setiap kabupaten/kota di Jawa Barat mempengaruhi PDRB pada kabupaten/kota tersebut.

#Variabel x3: Upah Minimum Kabupaten/Kota (UMK)
data.frame(data$Kabupaten.Kota, t=t3, pvalue=pvalue3) %>% 
  mutate(Keterangan=ifelse(pvalue3<0.05,"Reject", "Not Reject"))
##    data.Kabupaten.Kota         t    pvalue Keterangan
## 1              Bandung 0.7123261 0.4837496 Not Reject
## 2        Bandung Barat 0.6734833 0.5076542 Not Reject
## 3               Bekasi 1.1588223 0.2589530 Not Reject
## 4                Bogor 1.2370933 0.2290953 Not Reject
## 5               Ciamis 0.5310371 0.6007138 Not Reject
## 6              Cianjur 0.5819467 0.5665189 Not Reject
## 7              Cirebon 1.0955974 0.2851043 Not Reject
## 8                Garut 0.5993386 0.5550697 Not Reject
## 9            Indramayu 1.4626840 0.1576942 Not Reject
## 10            Karawang 1.0987535 0.2837551 Not Reject
## 11        Kota Bandung 0.7882078 0.4389862 Not Reject
## 12         Kota Banjar 0.6210477 0.5409501 Not Reject
## 13         Kota Bekasi 1.2243421 0.2337729 Not Reject
## 14          Kota Bogor 1.2261398 0.2331091 Not Reject
## 15         Kota Cimahi 0.7175014 0.4806146 Not Reject
## 16        Kota Cirebon 1.1058814 0.2807250 Not Reject
## 17          Kota Depok 1.3095012 0.2038762 Not Reject
## 18       Kota Sukabumi 0.8444325 0.4075171 Not Reject
## 19    Kota Tasikmalaya 0.4477516 0.6587104 Not Reject
## 20            Kuningan 0.8314891 0.4146315 Not Reject
## 21          Majalengka 0.9366125 0.3591266 Not Reject
## 22         Pangandaran 0.8361543 0.4120583 Not Reject
## 23          Purwakarta 0.8052718 0.4292808 Not Reject
## 24              Subang 1.0427423 0.3083914 Not Reject
## 25            Sukabumi 0.9786833 0.3383768 Not Reject
## 26            Sumedang 1.0602487 0.3005334 Not Reject
## 27         Tasikmalaya 0.6128963 0.5462292 Not Reject

Dilihat dari variabel Upah Minimum Kabupaten/Kota (X3), dapat disimpulkan bahwa UMK pada setiap kabupaten/kota di Jawa Barat tidak mempengaruhi PDRB di kabupaten/kota tersebut.

#Variabel x4: Tingkat Pengangguran Terbuka (TPT)
data.frame(data$Kabupaten.Kota, t=t4, pvalue=pvalue4) %>% 
  mutate(Keterangan=ifelse(pvalue4<0.05,"Reject", "Not Reject"))
##    data.Kabupaten.Kota        t      pvalue Keterangan
## 1              Bandung 1.816819 0.082894749 Not Reject
## 2        Bandung Barat 2.349031 0.028213161     Reject
## 3               Bekasi 2.888035 0.008536495     Reject
## 4                Bogor 2.264280 0.033756714     Reject
## 5               Ciamis 1.579737 0.128438777 Not Reject
## 6              Cianjur 2.182771 0.040010052     Reject
## 7              Cirebon 1.618690 0.119762343 Not Reject
## 8                Garut 1.575015 0.129525192 Not Reject
## 9            Indramayu 1.963005 0.062411829 Not Reject
## 10            Karawang 3.135039 0.004813617     Reject
## 11        Kota Bandung 1.854087 0.077182648 Not Reject
## 12         Kota Banjar 1.617093 0.120108287 Not Reject
## 13         Kota Bekasi 2.649270 0.014650213     Reject
## 14          Kota Bogor 2.217515 0.037225908     Reject
## 15         Kota Cimahi 2.070292 0.050363570 Not Reject
## 16        Kota Cirebon 1.621780 0.119095540 Not Reject
## 17          Kota Depok 2.338842 0.028832268     Reject
## 18       Kota Sukabumi 2.165289 0.041481233     Reject
## 19    Kota Tasikmalaya 1.489033 0.150675070 Not Reject
## 20            Kuningan 1.639859 0.115256625 Not Reject
## 21          Majalengka 1.509150 0.145489423 Not Reject
## 22         Pangandaran 1.681404 0.106827745 Not Reject
## 23          Purwakarta 2.893859 0.008423257     Reject
## 24              Subang 2.533670 0.018919840     Reject
## 25            Sukabumi 2.405448 0.025001307     Reject
## 26            Sumedang 1.675819 0.107929575 Not Reject
## 27         Tasikmalaya 1.562014 0.132555934 Not Reject

Dilihat dari variabel Investasi Penanaman Modal Asing (X4), dapat disimpulkan bahwa Investasi Penanaman Modal Asing di Kabupaten Bandung Barat, Kabupaten Bekasi, Kabupaten Bogor, Kabupaten Karawang, Kota Bekasi, Kota Bogor, Kota Depok, Kota Sukabumi, Kabupaten Purwakarta, Kabupaten Subang, dan Kabupaten Sukabumi mempengaruhi PDRB di kabupaten/kota tersebut.