UAS REGSPAS

Import Library

library(spdep)        # Analisis data spasial: spatial autocorrelation, spatial regression

## 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')`

## Loading required package: 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 TRUE

library(sp)           # Struktur data untuk manipulasi data spasial (poligon, titik, dll)

## Warning: package 'sp' was built under R version 4.3.3

library(sf)           # Mengelola data spasial dengan standar Simple Features (sf)
library(readxl)       # Membaca file Excel (format .xls dan .xlsx)

## Warning: package 'readxl' was built under R version 4.3.2

library(openxlsx)     # Membuat, membaca, dan menulis file Excel (lebih cepat dari readxl)

## Warning: package 'openxlsx' was built under R version 4.3.3

library(corrplot)     # Membuat visualisasi korelasi dalam bentuk plot

## Warning: package 'corrplot' was built under R version 4.3.2

## corrplot 0.92 loaded

library(DescTools)    # Kumpulan fungsi statistik deskriptif dan inferensial

## Warning: package 'DescTools' was built under R version 4.3.3

library(nortest)      # Melakukan uji normalitas (Anderson-Darling, Lilliefors, dll)
library(car)          # Fungsi regresi linear dan diagnostik model (vif, durbinWatsonTest, dll)

## Warning: package 'car' was built under R version 4.3.2

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.3.2

## 
## Attaching package: 'car'

## The following object is masked from 'package:DescTools':
## 
##     Recode

library(spatialreg)   # Estimasi model regresi spasial (SAR, SEM, SLM)

## 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.ZeroPolicyOption

library(lmtest)       # Fungsi untuk menguji model regresi (heteroskedastisitas, serial correlation)

## 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.2

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(dplyr)        # Manipulasi data seperti filter, select, mutate, dan summarise

## Warning: package 'dplyr' was built under R version 4.3.2

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(RColorBrewer) # Skema palet warna yang dapat digunakan dalam visualisasi
library(ggplot2)      # Membuat grafik statistik dengan pendekatan Grammar of Graphics

## Warning: package 'ggplot2' was built under R version 4.3.3

Tahap pertama dalam pemodelan regresi spasial adalah mngimport library yang diperlukan dalam analisis

Import Data

Data yang digunakan terdiri dari dua jenis data:

Data Polygon (Peta Pulau Jawa, dengan extension .shp)
Dataframe (Peubah Y dan X yang dilibatkan)

Berikut ini adalah peubah-peubah yang digunakan:

Peubah Y
Peubah X1
Peubah X2

Data Frame

data <- read_excel("C:/Users/Admin/Downloads/data-40.xlsx")
head(data)

## # A tibble: 6 × 3
##       Y     X1    X2
##   <dbl>  <dbl> <dbl>
## 1 5.68  0.684  0.171
## 2 5.95  0.690  0.440
## 3 1.51  0.195  0.308
## 4 0.674 0.204  0.101
## 5 3.62  0.374  0.280
## 6 0.514 0.0962 0.345

dim(data)

## [1] 119   3

Data terdiri dari 119 kabupaten/kota di Pulau Jawa. Peubah yang dilibatkan diantaranya adalah satu peubah respon (Y) dan dua peubah penjelas (X1 dan X2)

Data Peta

peta <- st_read(dsn ="C:/Users/Admin/Downloads/Adriano Excel Putra/Jawamap", layer = "jawa")

## Reading layer `jawa' from data source 
##   `C:\Users\Admin\Downloads\Adriano Excel Putra\Jawamap' using driver `ESRI Shapefile'
## Simple feature collection with 119 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.0998 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## Geodetic CRS:  WGS 84

peta

## Simple feature collection with 119 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.0998 ymin: -8.78036 xmax: 116.2702 ymax: -5.048857
## Geodetic CRS:  WGS 84
## First 10 features:
##    PROVNO KABKOTNO    PROVINSI           KABKOT ID2013
## 1      31       01 DKI JAKARTA KEPULAUAN SERIBU   3101
## 2      31       71 DKI JAKARTA  JAKARTA SELATAN   3171
## 3      31       72 DKI JAKARTA    JAKARTA TIMUR   3172
## 4      31       73 DKI JAKARTA    JAKARTA PUSAT   3173
## 5      31       74 DKI JAKARTA    JAKARTA BARAT   3174
## 6      31       75 DKI JAKARTA    JAKARTA UTARA   3175
## 7      32       01  JAWA BARAT            BOGOR   3201
## 8      32       02  JAWA BARAT         SUKABUMI   3202
## 9      32       03  JAWA BARAT          CIANJUR   3203
## 10     32       04  JAWA BARAT          BANDUNG   3204
##                          geometry
## 1  MULTIPOLYGON (((106.3842 -5...
## 2  MULTIPOLYGON (((106.8491 -6...
## 3  MULTIPOLYGON (((106.9719 -6...
## 4  MULTIPOLYGON (((106.8788 -6...
## 5  MULTIPOLYGON (((106.711 -6....
## 6  MULTIPOLYGON (((106.969 -6....
## 7  MULTIPOLYGON (((106.994 -6....
## 8  MULTIPOLYGON (((106.9652 -6...
## 9  MULTIPOLYGON (((107.2843 -6...
## 10 MULTIPOLYGON (((107.75 -6.8...

plot(peta, main = "Peta Pulau Jawa")

Berikut adalah visulaisisa Peta Jawa yang dibagi menjadi 5 Bagian

Eksplorasi Data

Pengecekan Missing Value

data[!complete.cases(data),]

## # A tibble: 0 × 3
## # ℹ 3 variables: Y <dbl>, X1 <dbl>, X2 <dbl>

Hasil pengecekan missing value menunjukkan bahwa tidak ada missing value pada data

Ringkasan Statistik Deskrpitif

summary(data[,c(1,2,3)])

##        Y                 X1                X2          
##  Min.   : 0.5135   Min.   :0.01925   Min.   :0.004079  
##  1st Qu.: 3.4436   1st Qu.:0.26193   1st Qu.:0.241331  
##  Median : 5.2935   Median :0.51078   Median :0.451687  
##  Mean   : 5.4175   Mean   :0.51291   Mean   :0.480641  
##  3rd Qu.: 7.3440   3rd Qu.:0.78740   3rd Qu.:0.720443  
##  Max.   :11.0567   Max.   :0.99968   Max.   :0.997271

# Menunjukkan kab/kota dengan Y tertinggi
kab_max <- data %>%
  arrange(desc(Y)) %>%
  slice(1)

print(paste("Kabupaten/kota dengan Y tertinggi adalah", kab_max$Y))

## [1] "Kabupaten/kota dengan Y tertinggi adalah 11.0566924649"

Boxplot

par(mfrow=c(1,3))
boxplot(data$Y,main = "Sebaran Y", col="dodgerblue")
boxplot(data$X1,main = "Sebaran X1", col="dodgerblue")
boxplot(data$X2,main = "Sebaran X2", col="dodgerblue")

Berdasarkan Boxplot dari peubah respon dan peubah penjelas diatas diketahui bahwa tidak terdapat pencilan pada data.

Scatter plot

pairs(data[,1:3], pch = 19, lower.panel = NULL)

Berdasarkan scatter plot di atas, tampak adanya pola hubungan linear positif antara peubah X1 maupun X2 terhadap peubah Y, dengan X1 tampak memiliki hubungan yang lebih kuat. Selain itu, tampak juga bahwa antar peubah penjelas tidak memiliki hubungan karena tidak memiliki pola tertentu.

Corplot

corrplot(cor(data[,1:3]), method = "number")

Peubah X1 dan X2 berkorelasi positif dengan Y dengan peubah X1 memiliki korelasi yang lebih tinggi daripada X2 . Hal ini menginidkasikan bahwa Korelasi X1 dan Peubah Y terindikasi kuat, sedangkan X2 lemah. Peubah X1 dan X2 juga tidak memiliki korelasi sehingga dapat mengindikasikan bahwa tidak ada multikoliniearitas.

Uji Korelasi Antar Peubah

Ho: Tidak terdapat korelasi antara peubah Y dan X1

H1: Ada korelasi antara peubah Y dan X1

Keputusan: Tolak Ho jika p-value < alpha=0.05

cor.test(data$Y, data$X1)

## 
##  Pearson's product-moment correlation
## 
## data:  data$Y and data$X1
## t = 16.598, df = 117, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7746143 0.8844228
## sample estimates:
##       cor 
## 0.8377981

cor.test(data$Y, data$X2)

## 
##  Pearson's product-moment correlation
## 
## data:  data$Y and data$X2
## t = 2.8594, df = 117, p-value = 0.005027
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.07922337 0.41641402
## sample estimates:
##       cor 
## 0.2555745

Dari Kedua uji , dapat dikatakan p-valuenya lebih kecil dari 0.05. Hal ini mengindiaksikan Tolak H-0 , artinya terdapat korelasi antar peubah X1, X2 terhadap Y secara parsial.

Sebaran Spasial Data

Sebaran Spasial Y

k=102
colfunc <- colorRampPalette(c("green", "yellow","red"))

peta$Y <- data$Y

# Plot menggunakan ggplot2
ggplot(peta) +
  geom_sf(aes(fill = Y)) +
  scale_fill_gradientn(colors = colfunc(k)) +
  labs(title = "Peta Sebaran Spasial Peubah Respon Y") +
  theme_minimal()

Dari Plot Sebaran Spasial Data, ada kecenderungan bahwa ada pola pengerombolan pada Peubah Y. Hal ini dapat diliahat dari warna gradasi yang cenderung mengelompok ( Hijau ke kkuningan) disekitar Jawa Tengah ke Jawa Timar. Hal ini bisa mengnindikasikan adanya autokorelasi positif pada data.

Sebaran Spasial X1

peta$X1<- data$X1
# Plot menggunakan ggplot2
ggplot(peta) +
  geom_sf(aes(fill = X1)) +
  scale_fill_gradientn(colors = colfunc(k)) +
  labs(title = "Sebaran Spasial Peubah Penjelas X1") +
  theme_minimal()

Sebaran Spasial X2

peta$X2<- data$X2

# Plot menggunakan ggplot2
ggplot(peta) +
  geom_sf(aes(fill = X2)) +
  scale_fill_gradientn(colors = colfunc(k)) +
  labs(title = "Sebaran Spasial Peubah Penjelas X2") +
  theme_minimal()

Hal yang menarik dapat dilihat pada peubah X1 dimana ada beberapa daerah yang cenderung kontras perbedaanya seperti di Jawa Tengah yang dimana ada beberapa daerah yang terindikasi sebagai Hot-Spot dimana daerahnya memiliki nilai yang lebih tinggi dibandingkan dengan derah yang lain.

Pembentukan Matrix Pembobot Spasial

Dari Visualisasi, Dapat Dilihat Bahwa Terdapat Daerah yang tidak memilli Tetangga, Hal ini mengindikasikan bahwa Matrix Pembobot dengan pendekatan ketetanggan kurang cocok dalam analisis data ini. Namun untuk perbandingan, Kedua Matrix tetangga tetap dimasukkan ( Rook dan Queen Conquinitiy).

Queen Conquinity

# Memvalidasi geometri dalam objek sf
peta <- st_make_valid(peta)
# Konversi objek sf ke objek sp
sp.peta <- as(peta, "Spatial")

# Mengakses slot polygons dari objek sp
sp.polygons <- SpatialPolygons(sp.peta@polygons)
queen <- poly2nb(sp.peta, queen = T)

## Warning in poly2nb(sp.peta, queen = T): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

## Warning in poly2nb(sp.peta, queen = T): neighbour object has 3 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

W.queen <- nb2listw(queen, style='W',zero.policy=TRUE)
qc = moran.test(data$Y, W.queen,randomisation=T,
alternative="greater", zero.policy = TRUE)
qc

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.queen  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.2013, p-value = 0.1148
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.073056069      -0.008547009       0.004614014

Indeks Moran tidak Singifikan pada Taraf nyata 5%. Hal ini mengindikasikan bahwa tidak terdapat autkorelasi spasial dalam penggunaan matrix pembobot Queen

Rook Conquinity

rook <- poly2nb(sp.peta, queen = F)

## Warning in poly2nb(sp.peta, queen = F): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

## Warning in poly2nb(sp.peta, queen = F): neighbour object has 3 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

W.rook <- nb2listw(rook, style='W',zero.policy=TRUE)
rc = moran.test(data$Y, W.rook,randomisation=T,
alternative="greater", zero.policy = TRUE)
rc

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.rook  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.0827, p-value = 0.1395
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.065154200      -0.008547009       0.004633723

Indeks Moran tidak Singifikan pada Taraf nyata 5%. Hal ini mengindikasikan bahwa tidak terdapat autkorelasi spasial dalam penggunaan matrix pembobot Rook.

Matrix Pembobot Jarak

Sebelum Menentukan matrix yang digunakan terlebih dahulu ditentukan jarak yang akan digunakan dalam pembentukan Matrix ini.Pada Kasus ini jarak matrix akan dihitung sebagai centroid dan digunakan metode jarak Euclidean. dalam pembentukan jaraknya.

# Menghitung centroid dari setiap unit spasial
centroids <- st_centroid(peta)

## Warning: st_centroid assumes attributes are constant over geometries

# Ekstrak koordinat centroid
longlat <- st_coordinates(centroids)

# Tambahkan koordinat long dan lat ke objek peta
peta$long <- longlat[, 1]  # Kolom X (long)
peta$lat  <- longlat[, 2]  # Kolom Y (lat

coords <- peta[c("long","lat")]
#class(coords)
koord <- as.data.frame(coords)
jarak<-dist(longlat, method = "euclidean")
m.jarak<-as.matrix(jarak)

KNN

Matriks bobot k-nearest neighbor menentukan tetangga berdasarkan sejumlah tetangga terdekat (k) yang dipilih untuk setiap titik atau wilayah. Pada kali ini saya akan mencoba 2 Tetangga yaitu dengan 3 dan 5 Tetangga.

KNN 5 Tetangga

W.knn<-knn2nb(knearneigh(longlat,k=5,longlat=TRUE))
W.knn.s <- nb2listw(W.knn,style='W')

MI.knn5 <- moran(data$Y,W.knn.s,n=length(W.knn.s$neighbours),S0=Szero(W.knn.s))
mt1 = moran.test(data$Y,W.knn.s,randomisation = T,alternative = "greater")
mt1

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = 1.6616, p-value = 0.0483
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.080037530      -0.008474576       0.002837664

KNN 3 Tetangga

W.knn<-knn2nb(knearneigh(longlat,k=3,longlat=TRUE))
W.knn.s <- nb2listw(W.knn,style='W')

MI.knn <- moran(data$Y,W.knn.s,n=length(W.knn.s$neighbours),S0=Szero(W.knn.s))
mt1 = moran.test(data$Y,W.knn.s,randomisation = T,alternative = "greater")
mt1

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.knn.s    
## 
## Moran I statistic standard deviate = 0.68051, p-value = 0.2481
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.038205920      -0.008474576       0.004705478

Indeks moran yang didapat memiliki variasi dari tetangga yang digunakan. Indeks Moran dengan 3 Tetangga tidak signifikan pada taraf nyata 5%. Namun hal ini berbeda dengan Indeks Moran dengan jumlah ketetanggan 5, dimana indeks morannya singinifkan pada taraf nyata 5%. Namun Niali indeks morannya tergolong kecil dan mendekati 0. Sehingga menunjukkan adanya korelasi acak antara peubah di daerahnya.

RDW

Matriks bobot radial distance mendefinisikan tetangga berdasarkan jarak geografis. Setiap titik atau wilayah dianggap sebagai tetangga jika jaraknya dari titik/wilayah pusat berada dalam jarak radius tertentu. Pada dasarnya Matrix RDW ( Radial Distance) ini aka mengelompokan ketetanggan berdasarkan suatu jarak yang telah ditentukan sebelumnya. Karena Pada kali ini jarak yang ditentukan adalah 100 km.

# Set zero.policy, false

W.rdw <-dnearneigh(longlat,0,65,longlat=TRUE)

## Warning in dnearneigh(longlat, 0, 65, longlat = TRUE): neighbour object has 2
## sub-graphs

W.rdw.s <- nb2listw(W.rdw,style='W',zero.policy=TRUE)

MI.rdw <- moran(data$Y,W.rdw.s,n=length(W.rdw.s$neighbours),S0=Szero(W.rdw.s))
mt2 = moran.test(data$Y,W.rdw.s,randomisation = T,alternative = "greater")
mt2

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.rdw.s  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 1.8134, p-value = 0.03489
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.073472037      -0.008547009       0.002045766

Didapatkan dengan RDW, Indeks Morannya signifikan pada taraf nyata 5%. Sehingga dapat dikatakan ada autokorelasi spasial pada data.

IDW ( Inverse Distance Weight)

Matriks pembobot Inverse Distance Weight (IDW) menentukan bobot berdasarkan kebalikan dari jarak antar titik atau wilayah. Konsepnya seperti suara, semakin dekat semakin besar, Semakin dekat dua titik, semakin besar bobot yang diberikan. Pada IDW sendiri ada ada indikator yang disebut alpha yang mengatur seberapa cepat suatu bobot akan berubah nilainya menjadi lebih kecil atau lebih besar, pada kali ini akan diatur 2 bobot yaitu alpha 1 dan alpha 2. ( Misa indikator disini kita sebut alpha).

IDW Alpha 1

alpha1=1
W.idw <-1/(m.jarak^alpha1)
diag(W.idw)<-0
rtot<-rowSums(W.idw,na.rm=TRUE)
W.idw.sd<-W.idw/rtot #row-normalized
W.idw.s = mat2listw(W.idw.sd,style='W')

MI.idw <- moran(data$Y,W.idw.s,n=length(W.idw.s$neighbours),S0=Szero(W.idw.s))
mt3 = moran.test(data$Y,W.idw.s,randomisation = T,alternative = "greater")
mt3

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.idw.s    
## 
## Moran I statistic standard deviate = 1.9816, p-value = 0.02376
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0325361575     -0.0084745763      0.0004283172

DIdapatkan p-value yang lebih kecil dari 0,05, hal ini mengindikasikan tolak H-0 yang artinya pada taraf 5%, terdapat autokorelasi spasial pada peubah y dengan metode IDW Alpha 1.

IDW Alpha 2

alpha2=2
W.idw2<-1/(m.jarak^alpha2)
diag(W.idw2)<-0
rtot<-rowSums(W.idw2,na.rm=TRUE)
W.idw2.sd<-W.idw2/rtot #row-normalized
W.idw2.s = mat2listw(W.idw2.sd,style='W')

MI.idw2 <- moran(data$Y,W.idw2.s,n=length(W.idw2.s$neighbours),S0=Szero(W.idw2.s))
mt4 = moran.test(data$Y,W.idw2.s,randomisation = T,alternative = "greater")
mt4

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.idw2.s    
## 
## Moran I statistic standard deviate = 1.999, p-value = 0.02281
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.090355206      -0.008474576       0.002444346

DIdapatkan p-value yang lebih kecil dari 0,05, hal ini mengindikasikan tolak H-0 yang artinya pada taraf 5%, terdapat autokorelasi spasial pada peubah y dengan metode IDW Alpha 1. Disini hal yang menarik adalah indeks morannya yang semakin besar di penambahan alpha.

EDW ( Exponential Distance Weight)

Matriks exponential distance weight memberikan bobot yang menurun secara eksponensial seiring dengan peningkatan jarak antar unit spasial. Ini mirip dengan inverse distance weight, tetapi penurunan bobot terhadap jarak lebih tajam karena mengikuti sebaran eksponesnial yang seetiap berubah satu titik semakin tajam juga perubahannya.

EDW Alpha 1

alpha=1
W.edw<-exp((-alpha)*m.jarak)
diag(W.edw)<-0
rtot<-rowSums(W.edw,na.rm=TRUE)
W.edw.sd<-W.edw/rtot #row-normalized
W.edw.s = mat2listw(W.edw.sd,style='W')

MI.edw <- moran(data$Y,W.edw.s,n=length(W.edw.s$neighbours),S0=Szero(W.edw.s))
mt5 = moran.test(data$Y,W.edw.s,randomisation = T,alternative = "greater")
mt5

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.edw.s    
## 
## Moran I statistic standard deviate = 1.5475, p-value = 0.06087
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0133757825     -0.0084745763      0.0001993578

Didapatkan Indeks Morannya tidak singikan pada taraf 5%.

EDW Alpha 2

alpha2=2
W.edw2<-exp((-alpha2)*m.jarak)
diag(W.edw2)<-0
rtot<-rowSums(W.edw2,na.rm=TRUE)
W.edw2.sd<-W.edw2/rtot #row-normalized
W.edw2.s = mat2listw(W.edw2.sd,style='W')
MI.edw2 <- moran(data$Y,W.edw2.s,n=length(W.edw2.s$neighbours),S0=Szero(W.edw2.s))
mt6 = moran.test(data$Y,W.edw2.s,randomisation = T,alternative = "greater")
mt6

## 
##  Moran I test under randomisation
## 
## data:  data$Y  
## weights: W.edw2.s    
## 
## Moran I statistic standard deviate = 1.7002, p-value = 0.04455
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0290258598     -0.0084745763      0.0004864974

Didapatkan Indeks Morannya signifikan pada taraf 5%.

Evaluasi Matrix Pembobot

Di tahap ini yang dimasukkan hanya KNN dengan 5 tetangganya saja karrena dengan 3 Tetangga tidak signifikan.

MatriksBobot <- c("K-Nearest Neighbor (k=5)",
                  "Radial Distance Weight (dmax=60)",
                  "Invers Distance Weight (alpha=1)",
                  "Invers Distance Weight (alpha=2)",
                  "Exponential Distance Weight (alpha=1)",
                  "Exponential Distance Weight (alpha=2)",
                  "Rook Contiguity",
                  "Queen Contiguity")
IndeksMoran <- c(MI.knn5$I,MI.rdw$I,MI.idw$I,MI.idw2$I,MI.edw$I,MI.edw2$I,rc$estimate[1],qc$estimate[1])
pv = c(mt1$p.value,mt2$p.value,mt3$p.value,mt4$p.value,mt5$p.value,mt6$p.value,rc$p.value,qc$p.value)
Matriks = cbind.data.frame(MatriksBobot, IndeksMoran,"p-value"=pv)
colnames(Matriks) <-c("Matriks Bobot", "Indeks Moran", "p-value")
Matriks

##                           Matriks Bobot Indeks Moran    p-value
## 1              K-Nearest Neighbor (k=5)   0.08003753 0.24809124
## 2      Radial Distance Weight (dmax=60)   0.07409468 0.03488730
## 3      Invers Distance Weight (alpha=1)   0.03253616 0.02376233
## 4      Invers Distance Weight (alpha=2)   0.09035521 0.02280574
## 5 Exponential Distance Weight (alpha=1)   0.01337578 0.06086652
## 6 Exponential Distance Weight (alpha=2)   0.02902586 0.04454811
## 7                       Rook Contiguity   0.06515420 0.13947001
## 8                      Queen Contiguity   0.07305607 0.11480906

Didapatkan dari 8 Model Matrix Pembobot Spasial, Terdapat variasi dari p-value. Dimana ada beberapa matrix pembobot yang menunjukkan tidak ada autokorelasi spasial, dan ada beberapa pula yang menunjukkan adanya korelasi spasial. Untuk itu akan difilter yang memiliki p-value < 0.05

Matriks[Matriks$`p-value`< 0.05,]

##                           Matriks Bobot Indeks Moran    p-value
## 2      Radial Distance Weight (dmax=60)   0.07409468 0.03488730
## 3      Invers Distance Weight (alpha=1)   0.03253616 0.02376233
## 4      Invers Distance Weight (alpha=2)   0.09035521 0.02280574
## 6 Exponential Distance Weight (alpha=2)   0.02902586 0.04454811

Didapat 5 Matrix Pembobot yang signifikan. Untuk itu akan dpilih berdasarkan Matrix pembobt yang memiliki Nilai Indeks Moran tertinggi. Disini ada 2 calon matrix pembobot yang nilainya paling tinggi dibandingkan matrix pembobt lainnya yaitu KNn dan IDW. Untuk Itu selanjutnya Pemodelan akan memakai Matrix IDW dengan Alpha (2).

Pemodelan Regresi

Regresi OLS

Pemodelan yang akan dilakukan pertama kali adalah pemodelan OLS, dimana Kita akan memodelkan peubah yang ada dengan pendekataan MKT.

ols <- lm(Y~X1+X2, data=data)
summary(ols)

## 
## Call:
## lm(formula = Y ~ X1 + X2, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.78980 -0.79573  0.00662  0.78422  2.43897 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.04566    0.27560  -0.166    0.869    
## X1           7.66908    0.34146  22.460  < 2e-16 ***
## X2           3.18246    0.35171   9.049 4.01e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.097 on 116 degrees of freedom
## Multiple R-squared:  0.8253, Adjusted R-squared:  0.8222 
## F-statistic: 273.9 on 2 and 116 DF,  p-value: < 2.2e-16

AIC(ols)

## [1] 364.6135

Dari Pemodelan dapat diambil beberapa kesimpulan:

Model memiliki adjusted R-Squared yang tinggi yaitu 0.8222 , Hal ini menandkan bahwa Model Regresi Linier berganda mampu dengan baik menjelaskan keragaman Y oleh variabel X1 dan X2 nya.
Semua peubah X1 dan X2 signifikan dan memiliki tanda koefisien yang selaras dengan korelasi yang sudah dibentuk
AIC dari Model adalah 364.6135

Uji Diagnostik OLS

Uji Normalitas Sisaan

Pengujian Normalitas akan menggunakan Uji Anderson Darling Test, karena Pada pengujian ini, sifat lokasi akan diperhitungkan sebagai kriteria kenormalan data

H0 : sisaan model menyebar normal H1 : sisaan model tidak menyebar normal

Tolak H0 jika p-value < 0.05

err.ols <- residuals(ols)
ad.test(err.ols)

## 
##  Anderson-Darling normality test
## 
## data:  err.ols
## A = 0.32445, p-value = 0.5202

Karena p-value > 5%, maka gagal tolak H0. Artinya sisaan menyebar normal.

Histogram

hist(err.ols)

Histogram dari rseidual model menunjukan bahwa data menyebar normal

QQ-Plot

qqnorm(err.ols,datax=T)
qqline(rnorm(length(err.ols),mean(err.ols),sd(err.ols)),datax=T, col="blue")

Dapat dilihat bahwa data cenderung menyebar normal karena ssebagian besar data sudah berada digaris Quantil. Terlihat pula ada satu pencilan yang keluar dari garisnya, Namun tidak menyebabkan data menjadi tidak normal.

Uji Kehomogenan Ragam

lmtest::bptest(ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 16.66, df = 2, p-value = 0.0002412

Karena p-value < 0.05, maka gagal tolak H0. Artinya ragam sisaan tidak homogen.

Uji Multikolinearitas

car::vif(ols)

##       X1       X2 
## 1.012596 1.012596

Didapatkan Bahwa Nilai VIF antar Peubah memiliki niali 1 yang dibawah treshold 5 dan 10. Hal ini menunjukkan bahwa tidak terjadi multikolinearitas antar peubah.

Pengujian Efek Spasial

Uji Indeks Moran

Pengujian autokorelasi spasial menggunakan Indeks Moran yang memerlukan matriks pembobot spasial. Dalam kasus ini, digunakan matriks bobot terbaik yang telah diperoleh sebelumnya, yaitu k-nearest neighbor dengan IDW Alpha 2

H0 : Tidak ada autokorelasi spasial

H1 : Ada autokorelasi spasial

ww = W.idw2.s
lm.morantest(ols, listw=ww, alternative="two.sided", zero.policy = TRUE)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = Y ~ X1 + X2, data = data)
## weights: ww
## 
## Moran I statistic standard deviate = 4.1328, p-value = 3.584e-05
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.194533435     -0.009058664      0.002426796

Didapatkan Bahwa odel memiliki autokorelasi spasial dibuktian pada indeks morannya yang nyata pada taraf 5%.

err.ols <- residuals(ols)
merr <- moran.test(err.ols, ww,randomisation=T, 
           alternative="two.sided")
merr

## 
##  Moran I test under randomisation
## 
## data:  err.ols  
## weights: ww    
## 
## Moran I statistic standard deviate = 4.1105, p-value = 3.948e-05
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.194533435      -0.008474576       0.002439118

Dapat dilihat bahwa pada model ( dan sisaanya) terdapat autkorelasi spasial. Hal ini dibuktikan pada Indeks moran yang signifikan dengan nilai yang lebih kecil dari taraf 0.05. Hal ini mengindikasikan bahwa terdapat autokorelasi spasial dan itu berpengaruh nyata terhadap model

Pengujian Hetogenitas Spasial

lmtest::bptest(ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 16.66, df = 2, p-value = 0.0002412

Dapat dikatakan bahwa terdapat heterogenitas spasial disini. Namun dalam penanganannya akan dicobakan terlebih dahulu dengan pemodelan spasial dan dilihat apakah dapat tertangani.

Pemodelan Regresi Spasial

Uji Lagrange Multiplier

model <- lm.LMtests(ols,listw=ww,zero.policy = TRUE, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))

## Please update scripts to use lm.RStests in place of lm.LMtests

summary(model)

##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = Y ~ X1 + X2, data = data)
## test weights: listw
##  
##          statistic parameter   p.value    
## RSerr      14.6817         1 0.0001273 ***
## adjRSerr    5.2525         1 0.0219149 *  
## RSlag      12.2624         1 0.0004622 ***
## adjRSlag    2.8332         1 0.0923358 .  
## SARMA      17.5149         2 0.0001573 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan dari output bahwa semua Pengujian bersifat nyata, Dimana pada Adjs Robust Lag nyata pada taraf 10 %, sedangkan pengujian lainnya nyata pada taraf 5%. Untuk itu akan dilakukan ketiga model Spasial untuk melihat apakah ketiga model tersebut dapat mengatasi autkorelasi spasial dan dependensi spasial yang ada.

SAR ( Spatial Autoregressive )

sar <- lagsarlm(Y ~ X1 + X2, data = data, listw = ww)
summary(sar, Nagelkerke = T)

## 
## Call:lagsarlm(formula = Y ~ X1 + X2, data = data, listw = ww)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.865092 -0.758024  0.017451  0.848347  2.396457 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -1.36972    0.46220 -2.9635  0.003042
## X1           7.58734    0.32251 23.5262 < 2.2e-16
## X2           3.15235    0.33088  9.5272 < 2.2e-16
## 
## Rho: 0.25983, LR test value: 10.56, p-value: 0.0011558
## Asymptotic standard error: 0.075756
##     z-value: 3.4299, p-value: 0.00060391
## Wald statistic: 11.764, p-value: 0.00060391
## 
## Log likelihood: -173.0269 for lag model
## ML residual variance (sigma squared): 1.0614, (sigma: 1.0302)
## Nagelkerke pseudo-R-squared: 0.84009 
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 356.05, (AIC for lm: 364.61)
## LM test for residual autocorrelation
## test value: 3.2231, p-value: 0.072607

Output di atas memperlihatkan bahwa koefisien Rho pada model SAR signifikan, dengan nilai AIC sebesar 356.05. Selain itu, terlihat pula hasil uji autokorelasi pada sisaan model memperlihatkan nilai p-value sebesar 0.07, artinya tidak terdapat lagi autokorelasi pada sisaan.

Uji asumsi Model SAR

Normalitas

err.sar<-residuals(sar)
ad.test(err.sar)

## 
##  Anderson-Darling normality test
## 
## data:  err.sar
## A = 0.51059, p-value = 0.1929

Kehomogenan Ragam

bptest.Sarlm(sar)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 6.6394, df = 2, p-value = 0.03616

Autokorelasi Spasial

moran.test(err.sar, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.sar  
## weights: ww    
## 
## Moran I statistic standard deviate = 1.7711, p-value = 0.07654
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.078940701      -0.008474576       0.002435932

Dari ketiga asumsi dapat dikatakan bahwa Normalitas dan Autokorelasi spasial sudah terselesaikan, namun belum untuk heterogenitas,.

SEM ( Spatial Error Model)

sem<-errorsarlm(Y ~ X1 + X2,data=data, listw= ww)
summary(sem, Nagelkerke = T)

## 
## Call:errorsarlm(formula = Y ~ X1 + X2, data = data, listw = ww)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.769945 -0.813055  0.045476  0.735133  2.404195 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.11059    0.30356  0.3643   0.7156
## X1           7.50443    0.30676 24.4635   <2e-16
## X2           3.15265    0.32168  9.8007   <2e-16
## 
## Lambda: 0.51875, LR test value: 13.281, p-value: 0.00026806
## Asymptotic standard error: 0.11862
##     z-value: 4.3734, p-value: 1.2233e-05
## Wald statistic: 19.127, p-value: 1.2233e-05
## 
## Log likelihood: -171.6661 for error model
## ML residual variance (sigma squared): 1.0008, (sigma: 1.0004)
## Nagelkerke pseudo-R-squared: 0.84371 
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 353.33, (AIC for lm: 364.61)

Uji Asumsi Model SEM

Normalitas

err.sem<-residuals(sem)
ad.test(err.sem)

## 
##  Anderson-Darling normality test
## 
## data:  err.sem
## A = 0.42164, p-value = 0.3177

Kehomogenan Ragam

bptest.Sarlm(sem)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 11.528, df = 2, p-value = 0.003139

moran.test(err.sem, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.46766, p-value = 0.64
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.014611413      -0.008474576       0.002436894

Pada Model SEM , terdapat 2 asumsi yang terlanggar yaitu kehomogenan ragam dan masih adanya autokorelasi spasial

Model SARMA

SARMA <- sacsarlm(Y ~ X1+X2,data=data,ww, zero.policy = TRUE)
summary(SARMA)

## 
## Call:
## sacsarlm(formula = Y ~ X1 + X2, data = data, listw = ww, zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.800523 -0.766824  0.054183  0.737547  2.407304 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.64223    0.72636 -0.8842   0.3766
## X1           7.60663    0.31765 23.9464   <2e-16
## X2           3.17940    0.32507  9.7808   <2e-16
## 
## Rho: 0.12466
## Asymptotic standard error: 0.11644
##     z-value: 1.0705, p-value: 0.28437
## Lambda: 0.38561
## Asymptotic standard error: 0.16813
##     z-value: 2.2934, p-value: 0.021823
## 
## LR test value: 14.146, p-value: 0.0008476
## 
## Log likelihood: -171.2337 for sac model
## ML residual variance (sigma squared): 1.0135, (sigma: 1.0067)
## Number of observations: 119 
## Number of parameters estimated: 6 
## AIC: 354.47, (AIC for lm: 364.61)

err.SARMA<-residuals(SARMA)
ad.test(err.SARMA)

## 
##  Anderson-Darling normality test
## 
## data:  err.SARMA
## A = 0.41745, p-value = 0.3251

bptest.Sarlm(SARMA)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 9.5116, df = 2, p-value = 0.008602

moran.test(err.SARMA, ww, alternative="two.sided")

## 
##  Moran I test under randomisation
## 
## data:  err.SARMA  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.29006, p-value = 0.7718
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.005842557      -0.008474576       0.002436318

Asumsi yang sama terjadi pada Model SARMA , yang sejalan dengan Model SAR dimana masih terdepat heterogenitas spasial yang terjadi.

Model SDM

sdm <- lagsarlm(Y ~ X1 + X2, data=data, ww, zero.policy = TRUE, type="mixed");
summary(sdm)

## 
## Call:lagsarlm(formula = Y ~ X1 + X2, data = data, listw = ww, type = "mixed", 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.811812 -0.779087  0.025882  0.758124  2.401374 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.61366    0.63752 -0.9626  0.33576
## X1           7.55991    0.31174 24.2507  < 2e-16
## X2           3.16786    0.32198  9.8387  < 2e-16
## lag.X1      -2.85557    1.28718 -2.2185  0.02652
## lag.X2      -1.13669    0.94816 -1.1988  0.23059
## 
## Rho: 0.49368, LR test value: 11.698, p-value: 0.00062577
## Asymptotic standard error: 0.12119
##     z-value: 4.0735, p-value: 4.6318e-05
## Wald statistic: 16.593, p-value: 4.6318e-05
## 
## Log likelihood: -171.0136 for mixed model
## ML residual variance (sigma squared): 0.99474, (sigma: 0.99737)
## Number of observations: 119 
## Number of parameters estimated: 7 
## AIC: 356.03, (AIC for lm: 365.72)
## LM test for residual autocorrelation
## test value: 0.17282, p-value: 0.67762

bptest.Sarlm(sdm)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 14.009, df = 4, p-value = 0.007268

Goodnes of Fit

Dari Kelima Model terdapat beberapa model yang sudah memenuhi asumsi autokorelasinya dimana ada beberapa calon kandidat model yang sudah tidak terdapat lagi autokorelasi spasial pada modelnya.

df2 <- data.frame("Model" = c("OLS (Regresi KlasiK)","SAR","SEM","SARMA", "SDM"),
           "AIC" = c(AIC(ols),AIC(sar), AIC(sem),AIC(SARMA), AIC(sdm)),
           "Normalitas" = c(ad.test(err.ols)$p.value, ad.test(err.sar)$p.value, ad.test(err.sem)$p.value, ad.test(err.SARMA)$p.value, ad.test(err.SARMA)$p.value),
           "Heterogenitas" = c(lmtest::bptest(ols)$p.value, bptest.Sarlm(sar)$p.value, bptest.Sarlm(sem)$p.value, bptest.Sarlm(SARMA)$p.value, bptest.Sarlm(sdm)$p.value),
           "Autokorelasi Spasial" = c(moran.test(err.ols, ww, alternative="two.sided")$p.value, moran.test(err.sar, ww, alternative="two.sided")$p.value, moran.test(err.sem, ww, alternative="two.sided")$p.value, moran.test(err.SARMA, ww, alternative="two.sided")$p.value, moran.test(err.SARMA, ww, alternative="two.sided")$p.value))

df2

##                  Model      AIC Normalitas Heterogenitas Autokorelasi.Spasial
## 1 OLS (Regresi KlasiK) 364.6135  0.5202404  0.0002411682         3.947689e-05
## 2                  SAR 356.0538  0.1928857  0.0361627962         7.653613e-02
## 3                  SEM 353.3322  0.3176869  0.0031389239         6.400278e-01
## 4                SARMA 354.4673  0.3250679  0.0086017569         7.717697e-01
## 5                  SDM 356.0272  0.3250679  0.0072675707         7.717697e-01

Berdasarkan output diatas, Dari Kriteria AIC sendiri Model yang menjadi model Terbaik adalah MODEL SEM. Namun bila diperhatikan, Kelima model ini masih tidak memenuhi asumsi Kehomogenan Ragam. Hal ini mengindikasikan bahwa masih terdapat heterogenitas spasial pada data yang ada. Namun pada Penerapan Model Dependensi, bila asumsi ini dianggap setara untuk semua model, Pilihan Model terbaik akan tetap pada SEM karena Normalitas dan Autkorelasi yang sudah terpenuh, dibuktikan sudah tidak adalagi autkorelasi spasial pada sisaaan. Perlu diperhatikan bahwa hetoregenitas perlu diatas dengan beberapa cara dan tidak boleh diabaikan, namun untuk saat ini, Model pilihannya adalah SEM.

Kesimpulan dan Interpretasi

Model SEM sendiri adalah ( Spatial Error Model) yang artinya memodelkan galat spasial. Hal ini berimplikasi bahwa tidak ada efek Spill over atauapun Umpan Balik karena yang dimodelkan pada konteks ini adalah Errornya.

sem<-errorsarlm(Y ~ X1 + X2,data=data, listw= ww)
summary(sem, Nagelkerke = T)

## 
## Call:errorsarlm(formula = Y ~ X1 + X2, data = data, listw = ww)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.769945 -0.813055  0.045476  0.735133  2.404195 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.11059    0.30356  0.3643   0.7156
## X1           7.50443    0.30676 24.4635   <2e-16
## X2           3.15265    0.32168  9.8007   <2e-16
## 
## Lambda: 0.51875, LR test value: 13.281, p-value: 0.00026806
## Asymptotic standard error: 0.11862
##     z-value: 4.3734, p-value: 1.2233e-05
## Wald statistic: 19.127, p-value: 1.2233e-05
## 
## Log likelihood: -171.6661 for error model
## ML residual variance (sigma squared): 1.0008, (sigma: 1.0004)
## Nagelkerke pseudo-R-squared: 0.84371 
## Number of observations: 119 
## Number of parameters estimated: 5 
## AIC: 353.33, (AIC for lm: 364.61)

Interptasi Model :

Model Tersebut memiliki Pseudo R-Squared : 74% yang lebih kecil dari OLS namun telah mengatasi efek dependensi spasial yang ditimbulkan pada saat pemodelan OLS
AIC yang didapat merupakan AIC yang paling kecil sehingga menjadi indikator model terbaik dari Pemillihan Model
Koefisien Spasialnya bernilai 0.51 yang nyata pada taraf 5%.
$ Y= 0.11059+7.50443X1+3.15265X2+0.515585W$
Model tersebut memiliki signifikansi disemua paramternya dan parameternya searah dengan korelasi plotnyaa. Dimana dapat diinterpratsikan setiap kenaikan X1 dan peubah lainnya tetap, akan meningkatkan Y sebesar 7.50443 dan setiap kenaikan X2 dan peubah lainnya tetap akan meningkatan Y sebsesar 3.16265. denga 0.515585 sebagai efek kesalah spasialnya

Penangan Kehomogenan Ragam ( Tambahan)

Terdapat Permasalahan pada Asumsi Kehomogenan Ragam pada model OLSnya kemudian setelah dilakukan pemodelan Spasial, Masalah ini masih belom terselesaikan, untuk itu akan dicoba penanganan dengan transformasi, ataupun pemodelan dengan GWR ( Geogrpahic Weighted Regrression)

Transformasi Logaritam

# Lakukan Transformasi Log pada Model OLS pada peubah y
data$Y <- log(data$Y)

# Pemodelan OLS dengan data yang sudah di trasnformasi
ols2 <- lm(Y~X1+X2, data=data)

# Lakukan pengujian yang sama seperti yang sudah saya lakukan 
summary(ols2)

## 
## Call:
## lm(formula = Y ~ X1 + X2, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.39506 -0.15932  0.01382  0.20188  0.65834 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.30652    0.08395   3.651 0.000393 ***
## X1           1.69149    0.10400  16.264  < 2e-16 ***
## X2           0.75294    0.10713   7.029 1.54e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.334 on 116 degrees of freedom
## Multiple R-squared:  0.7157, Adjusted R-squared:  0.7108 
## F-statistic:   146 on 2 and 116 DF,  p-value: < 2.2e-16

AIC(ols2)

## [1] 81.68016

# Uji asumsi kehomogenan ragam
lmtest::bptest(ols2)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols2
## BP = 19.516, df = 2, p-value = 5.782e-05

Dapat dilihat bahwa pemodelan Regresi dengan Transformasi Log belum mengatasi keheterogenan dn tidak mengatasi masalahnya, sehingga pemodelan mungkin bisa dilanjutkan dengan GWR