Geographically Weighted Regression
EKSPLORASI DATA
Import Data
#load package
#Packages Yang Digunakan
library(readxl)
library(openxlsx)
library(spdep)
library(rgdal)
library(raster)
library(corrplot)
library(DescTools)
library(nortest)
library(car)
library(spatialreg)
library(sf)
library(nortest)
library(DescTools)
library(lmtest)
library(tidyverse)
library(spgwr)
#input file excel
datauas <-read.csv("D:/MATERI KULIAH S2 IPB/SEMESTER 2/SPASIAL/MATERI BU ANIK/UAS2021-22.csv")
head(datauas)## X ID_KAB KAB y x1 x2 x3
## 1 1 3101 DKI JAKARTA 22.90460 85.54 72.22 4.08
## 2 2 3171 DKI JAKARTA 26.15227 98.01 74.05 4.75
## 3 3 3172 DKI JAKARTA 25.00793 96.62 79.49 4.88
## 4 4 3173 DKI JAKARTA 30.17872 93.10 67.41 4.68
## 5 5 3174 DKI JAKARTA 46.84015 90.66 67.52 5.94
## 6 6 3175 DKI JAKARTA 24.21561 86.20 67.16 4.20dim(datauas)## [1] 119 7Data terdiri dari 119 baris/amatan dan 7 kolom/peubah. Terdapat satu peubah respon Y dan 3 peubah bebas X1 X2 dan X3.
#Input peta shp
peta<-readOGR(dsn="D:/MATERI KULIAH S2 IPB/SEMESTER 2/SPASIAL/MATERI BU ANIK/Jawamap", layer="jawa") ## OGR data source with driver: ESRI Shapefile
## Source: "D:\MATERI KULIAH S2 IPB\SEMESTER 2\SPASIAL\MATERI BU ANIK\Jawamap", layer: "jawa"
## with 119 features
## It has 5 fieldsCek Missing Value Data
datauas[!complete.cases(datauas),]## [1] X ID_KAB KAB y x1 x2 x3
## <0 rows> (or 0-length row.names)Tidak terdapat missing value/nilai yang hilang pada data.
Summary Data
summary(datauas)## X ID_KAB KAB y
## Min. : 1.0 Min. :3101 Length:119 Min. : 6.282
## 1st Qu.: 30.5 1st Qu.:3276 Class :character 1st Qu.:18.781
## Median : 60.0 Median :3327 Mode :character Median :25.999
## Mean : 60.0 Mean :3386 Mean :26.069
## 3rd Qu.: 89.5 3rd Qu.:3516 3rd Qu.:32.926
## Max. :119.0 Max. :3674 Max. :46.840
## x1 x2 x3
## Min. :85.05 Min. :65.30 Min. :4.010
## 1st Qu.:88.75 1st Qu.:69.27 1st Qu.:4.395
## Median :92.02 Median :72.61 Median :4.970
## Mean :92.29 Mean :72.79 Mean :4.966
## 3rd Qu.:95.52 3rd Qu.:75.97 3rd Qu.:5.570
## Max. :99.96 Max. :79.70 Max. :6.000Sebaran Data
Sebaran dari data akan dieksplorasi menggunakan boxplot dan diperoleh hasil sebagai berikut:
boxplot(datauas$y,main="Sebaran Peubah Y") 
boxplot(datauas$x1,main="Sebaran Peubah X1") 
boxplot(datauas$x2,main="Sebaran Peubah X2") 
boxplot(datauas$x3,main="Sebaran Peubah X3")
Berdasrakan Boxplot di atas, terlihat bahwa tidak terdapat nilai pencilan pada peubah Y, X1, X2, dan X3.
Relationship of the data
- Scatter Plot X1 vs Y
plot(datauas$x1, datauas$y,
xlab="X1",
ylab="Y",
main="Scatter Plot of X1 and Y",
pch=20, col="orange", cex=2)
reg.klasik = lm(y~x1, data=datauas)
lines.lm(reg.klasik, col=2, add=T)
Scatter plot X1 vs Y di atas menunjukkan bahwa terdapat hubungan linier negatif antara X1 dengan Y.
- Scatter Plot X2 vs Y
plot(datauas$x2, datauas$y,
xlab="X2",
ylab="Y",
main="Scatter Plot of X2 and Y",
pch=20, col="orange", cex=2)
reg.klasik = lm(y~x2, data=datauas)
lines.lm(reg.klasik, col=2, add=T)
Scatter plot X2 vs Y di atas menunjukkan bahwa terdapat hubungan linier negatif antara X2 dengan Y.
- Scatter Plot X3 vs Y
plot(datauas$x3, datauas$y,
xlab="X3",
ylab="Y",
main="Scatter Plot of X3 and Y",
pch=20, col="orange", cex=2)
reg.klasik = lm(y~x3, data=datauas)
lines.lm(reg.klasik, col=2, add=T)
Scatter plot X3 vs Y di atas menunjukkan bahwa terdapat hubungan linier positif antara X3 dengan Y.
Correlation of the data
1. Matriks Korelasi
corrplot(cor(datauas[,4:7]), method = "number")
Berdasarkan matriks korelasi di atas terlihat bahwa Y memiliki korelasi negatif dengan X1 dan X2, dan memiliki korelasi positif dengan X3.
2. Pengujian Korelasi
- X1 vs Y
H0: Tidak terdapat korelasi antara peubah X1 dan Y
H1: Terdapat korelasi antara peubah X1 dan Y
cor.test(datauas$y, datauas$x1)##
## Pearson's product-moment correlation
##
## data: datauas$y and datauas$x1
## t = -2.9295, df = 117, p-value = 0.004083
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.42157021 -0.08543512
## sample estimates:
## cor
## -0.2614107Kesimpulan :
Karena diperoleh p-value < 0.05 maka dapat disimpulkan bahwa H0 ditolak yang berarti terdapat korelasi antara peubah X1 dan Y.
- X2 vs Y
H0: Tidak terdapat korelasi antara peubah X2 dan Y
H1: Terdapat korelasi antara peubah X2 dan Y
cor.test(datauas$y, datauas$x2)##
## Pearson's product-moment correlation
##
## data: datauas$y and datauas$x2
## t = -4.9965, df = 117, p-value = 2.065e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5572821 -0.2588973
## sample estimates:
## cor
## -0.419351Kesimpulan :
Karena diperoleh p-value < 0.05 maka dapat disimpulkan bahwa H0 ditolak yang berarti terdapat korelasi antara peubah X2 dan Y.
- X3 vs Y
H0: Tidak terdapat korelasi antara peubah X3 dan Y
H1: Terdapat korelasi antara peubah X3 dan Y
cor.test(datauas$y, datauas$x3)##
## Pearson's product-moment correlation
##
## data: datauas$y and datauas$x3
## t = 17.882, df = 117, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7986491 0.8974246
## sample estimates:
## cor
## 0.8556431Kesimpulan :
Karena diperoleh p-value < 0.05 maka dapat disimpulkan bahwa H0 ditolak yang berarti terdapat korelasi antara peubah X3 dan Y.
Sebaran Spasial Data
1. Sebaran Spasial Y
k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$y<- datauas$y
spplot(peta, "y", col.regions=color, main="Sebaran Spasial Peubah Y")
Berdasarkan plot di atas, dapat dilihat adanya kecenderungan pola bergerombol pada peubah Y. Hal ini tampak dari gradasi warna yang cenderung mengumpul, seperti pada warna kuning dan oranye.
2. Sebaran Spasial X1
k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$x1<- datauas$x1
spplot(peta, "x1", col.regions=color, main="Sebaran Spasial Peubah X1")
3. Sebaran Spasial X2
k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$x2<- datauas$x2
spplot(peta, "x2", col.regions=color, main="Sebaran Spasial Peubah X2")
3. Sebaran Spasial X3
k=16
colfunc <- colorRampPalette(c("green", "yellow","red"))
color <- colfunc(k)
peta$x3<- datauas$x3
spplot(peta, "x3", col.regions=color, main="Sebaran Spasial Peubah X3")
Matriks Pembobot Spasial
Matriks bobot berdasarkan jarak
longlat <- coordinates(peta)
head(longlat)## [,1] [,2]
## 0 106.4933 -5.796823
## 1 106.8100 -6.272549
## 2 106.9003 -6.255373
## 3 106.8351 -6.181230
## 4 106.7483 -6.165469
## 5 106.8695 -6.129185plot(longlat)
peta$long <- longlat[,1] #longitudinal yg nilainya ratusan
peta$lat <- longlat[,2] #latitude yang nilainya satuan
coords <- peta[c("long","lat")]
class(coords)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"# menghitung matriks jarak
djarak<-dist(longlat)
m.djarak<-as.matrix(djarak) #matriks jarak1. K-Nearest Neighbor Weight (Bobot k-Tetangga Terdekat)
#k=5, 5 tetangga terdekat
W.knn<-knn2nb(knearneigh(longlat,k=3,longlat=TRUE)) #matriks bobot dengan knn k=5 #knearneigh(x, k=1, longlat = NULL, use_kd_tree=TRUE)
W.knn## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 357
## Percentage nonzero weights: 2.521008
## Average number of links: 3
## Non-symmetric neighbours listclass(W.knn)## [1] "nb"W.knn1 <- nb2listw(W.knn,style='W')
W.knn1## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 357
## Percentage nonzero weights: 2.521008
## Average number of links: 3
## Non-symmetric neighbours list
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 119 14161 119 68.77778 495.1111plot(peta, col='gray', border='blue', main ="KNN dengan k=5")
plot(W.knn1, longlat, col='red', lwd=2, add=TRUE)
2. Radial Distance Weight
#d=70 ditetapkan nilai ambang dmax adalah 70 km
W.dmax<-dnearneigh(longlat,0,70,longlat=TRUE) #dnearneigh(x, d1, d2, row.names = NULL, longlat = NULL, bounds=c("GE", "LE"), use_kd_tree=TRUE, symtest=FALSE)
class(W.dmax) #nb## [1] "nb"W.dmax.s <- nb2listw(W.dmax,style='W') #W is row standardised (sums over all links to n). Standardisasi Baris
W.dmax.s## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 1300
## Percentage nonzero weights: 9.180143
## Average number of links: 10.92437
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 119 14161 119 25.99101 482.61023. Power Distance Weight / Invers Distance Weight (Bobot Jarak Invers)
#Alpha = 1
alpha1=1
W.idw <-1/(m.djarak^alpha1)
class(W.idw)## [1] "matrix" "array"#dinormalisasi
diag(W.idw)<-0
rtot<-rowSums(W.idw,na.rm=TRUE)
W.idw.sd<-W.idw/rtot #row-normalized
rowSums(W.idw.sd,na.rm=TRUE)## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1class(W.idw.sd)## [1] "matrix" "array"W.idw.1 = mat2listw(W.idw.sd,style='W')
summary(W.idw.1)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 14042
## Percentage nonzero weights: 99.15966
## Average number of links: 118
## Link number distribution:
##
## 118
## 119
## 119 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 119 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 119 14161 119 8.084636 479.8627#Alpha = 2
alpha1=2
W.idw2 <-2/(m.djarak^alpha1)
class(W.idw2)## [1] "matrix" "array"#dinormalisasi
diag(W.idw2)<-0
rtot<-rowSums(W.idw2,na.rm=TRUE)
W.idw.sd2<-W.idw2/rtot #row-normalized
rowSums(W.idw.sd2,na.rm=TRUE)## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1class(W.idw.sd2)## [1] "matrix" "array"W.idw.22 = mat2listw(W.idw.sd2,style='W')
summary(W.idw.22)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 14042
## Percentage nonzero weights: 99.15966
## Average number of links: 118
## Link number distribution:
##
## 118
## 119
## 119 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 119 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 119 14161 119 37.02077 490.92324. Exponential Distance Weight
alpha=1
W.e<-exp((-alpha)*m.djarak) #dinormalisasi
diag(W.e)<-0
rtot<-rowSums(W.e,na.rm=TRUE)
W.e.sd<-W.e/rtot #row-normalized
rowSums(W.e.sd,na.rm=TRUE)## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1class(W.e.sd)## [1] "matrix" "array"W.ed1 = mat2listw(W.e.sd,style='W')
summary(W.ed1)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 14042
## Percentage nonzero weights: 99.15966
## Average number of links: 118
## Link number distribution:
##
## 118
## 119
## 119 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 119 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 119 14161 119 4.839428 478.0838alpha=2
W.e2<-exp((-alpha)*m.djarak) #dinormalisasi
diag(W.e2)<-0
rtot<-rowSums(W.e2,na.rm=TRUE)
W.e.sd2<-W.e2/rtot #row-normalized
rowSums(W.e.sd2,na.rm=TRUE)## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1class(W.e.sd2)## [1] "matrix" "array"W.ed2 = mat2listw(W.e.sd2,style='W')
summary(W.ed2)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 14042
## Percentage nonzero weights: 99.15966
## Average number of links: 118
## Link number distribution:
##
## 118
## 119
## 119 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
## 119 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 with 118 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 119 14161 119 8.900784 479.4497List bobot berdasarkan jarak:
W.knn1 (KNN)
W.dmax.s (Radial Distance Weight)
W.idw.1 (Power Distance Weight)
W.idw.22 (Power Distance Weight)
W.ed1 (Exponential Distance Weight)
W.ed2 (Exponential Distance Weight)
Matriks Bobot berdasarkan ketetanggaan (Spatial Contiguity Weight)
1. Rook Contiguity
library(readr)
peta1 <- read_sf("D:/MATERI KULIAH S2 IPB/SEMESTER 2/SPASIAL/MATERI BU ANIK/Jawamap/jawa.shp")
sf::sf_use_s2(FALSE)## Spherical geometry (s2) switched off#Rook
W.rook<-poly2nb(peta1,queen=FALSE) ## although coordinates are longitude/latitude, st_intersects assumes that they are planarW.rook #matriks bobot Rook## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 518
## Percentage nonzero weights: 3.657934
## Average number of links: 4.352941
## 1 region with no links:
## 1class(W.rook) #nb## [1] "nb"W.rook.s <- nb2listw(W.rook,style='W', zero.policy = TRUE) #W is row standardised (sums over all links to n). Standardisasi Baris
print(W.rook.s, zero.policy=TRUE)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 518
## Percentage nonzero weights: 3.657934
## Average number of links: 4.352941
## 1 region with no links:
## 1
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 118 13924 118 66.41255 515.06332. Queen Contiguity
# queen
W.queen<-poly2nb(peta,queen=TRUE) ## although coordinates are longitude/latitude, st_intersects assumes that they are planarW.queen #matriks bobot Rook## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 522
## Percentage nonzero weights: 3.68618
## Average number of links: 4.386555
## 1 region with no links:
## 0class(W.queen) #nb## [1] "nb"W.queen.s <- nb2listw(W.queen,style='W', zero.policy = TRUE) #W is row standardised (sums over all links to n). Standardisasi Baris
print(W.queen.s, zero.policy=TRUE)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 119
## Number of nonzero links: 522
## Percentage nonzero weights: 3.68618
## Average number of links: 4.386555
## 1 region with no links:
## 0
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 118 13924 118 66.14482 515.6022Pemilihan Matriks Bobot
1. KNN
MI.knn <- moran(datauas$y, W.knn1, n=length(W.knn1$neighbours), S0=Szero(W.knn1))2. Radial Distance Weight
MI.radial <- moran(datauas$y, W.dmax.s, n=length(W.dmax.s$neighbours), S0=Szero(W.dmax.s))3. Power Distance Weight
MI.power1 <- moran(datauas$y, W.idw.1, n=length(W.idw.1$neighbours), S0=Szero(W.idw.1))
MI.power2 <- moran(datauas$y, W.idw.22, n=length(W.idw.22$neighbours), S0=Szero(W.idw.22))4. Exponential Distance Weight
MI.exp1 <- moran(datauas$y, W.ed1, n=length(W.ed1$neighbours), S0=Szero(W.ed1))
MI.exp2 <- moran(datauas$y, W.ed2, n=length(W.ed2$neighbours), S0=Szero(W.ed2))5. Spatial Contiguity
MI.rook <- moran.test(datauas$y, W.rook.s, randomisation = TRUE, zero.policy = TRUE)
MI.rook$estimate## Moran I statistic Expectation Variance
## -0.008530337 -0.008547009 0.004627167MI.queen <- moran.test(datauas$y, W.queen.s, randomisation = TRUE, zero.policy = TRUE)
MI.queen$estimate## Moran I statistic Expectation Variance
## -0.010844530 -0.008547009 0.004607486Perbandingan Nilai Indeks Moran’s
moranindeks<-data.frame(
"Matriks Bobot"=c("KNN (k=5)", "Radial distance weight", "Power distance weight (alpha=1)", "Power distance weight (alpha=2)", "Exponential distance weight (alpha=1)", "Exponential distance weight (alpha=2)", "Rook Contiguity", "Queen Contiguity"),
"Nilai Indeks Moran"=c(MI.knn$I, MI.radial$I, MI.power1$I, MI.power2$I, MI.exp1$I, MI.exp2$I, -0.008530337, -0.010844530))
moranindeks## Matriks.Bobot Nilai.Indeks.Moran
## 1 KNN (k=5) -0.069697688
## 2 Radial distance weight 0.041973476
## 3 Power distance weight (alpha=1) 0.003805885
## 4 Power distance weight (alpha=2) 0.027038800
## 5 Exponential distance weight (alpha=1) 0.011250108
## 6 Exponential distance weight (alpha=2) 0.017438536
## 7 Rook Contiguity -0.008530337
## 8 Queen Contiguity -0.010844530Kesimpulan:
Matriks bobot yang menghasilkan nilai Indeks Moran terbesar adalah matriks bobot radial distance weight dengan d=70.
Selanjutnya akan dilakukan pengujian indeks moran dengan menggunakan matriks bobot radial distance weight:
H0: I=0 (tidak ada autokorelasi spasial antarlokasi) H1: I>0 (ada autokorelasi spasial antar lokasi)
Woptimum <- W.dmax.s
moran.test(datauas$y, Woptimum, randomisation = TRUE, zero.policy = TRUE)##
## Moran I test under randomisation
##
## data: datauas$y
## weights: Woptimum
##
## Moran I statistic standard deviate = 1.2237, p-value = 0.1105
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.041973476 -0.008474576 0.001699493Kesimpulan :
Berdasarkan hasil di atas, diperoleh p-value < 0.05 berarti H0
ditolak, dengan alternative hypothesis: greater berarti
pada taraf nyata 5%, terdapat autokorelasi spasial positif.
Moran Plot
moran.plot(datauas$y, Woptimum, labels=datauas$KABKOT)
Pengujian Efek Dependensi Spasial
Regresi Linier
reg.klasik = lm(y~x1+x2+x3, data = datauas)
summary(reg.klasik)##
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = datauas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.0899 -2.1014 -0.1042 2.1141 6.4331
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.30956 8.14433 5.195 8.96e-07 ***
## x1 -0.27537 0.06860 -4.014 0.000107 ***
## x2 -0.66634 0.07084 -9.406 6.26e-16 ***
## x3 11.61586 0.48504 23.948 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.204 on 115 degrees of freedom
## Multiple R-squared: 0.8687, Adjusted R-squared: 0.8653
## F-statistic: 253.6 on 3 and 115 DF, p-value: < 2.2e-16Uji Multikolinearitas
car::vif(reg.klasik)## x1 x2 x3
## 1.033498 1.034220 1.013487Kesimpulan:
Karena nilai VIF pada x1, x2, dan x3 kurang dari 5, maka dapat disimpulkan bahwa tidak terdapat multikolinearitas antarpeubah penjelas.
Uji Normalitas
H0: sisaan model menyebar normal
H1: sisaan model tidak menyebar normal
err.regklasik<-reg.klasik$residuals
ad.test(err.regklasik) #anderson darling##
## Anderson-Darling normality test
##
## data: err.regklasik
## A = 0.54235, p-value = 0.1605hist(err.regklasik)
qqnorm(err.regklasik,datax=T)
qqline(rnorm(length(err.regklasik),mean(err.regklasik),sd(err.regklasik)),datax=T, col="red")
Kesimpulan:
Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).
Uji Autokorelasi Spasial
H0 : Tidak terdapat Autokorelasi pada sisaan
H1 : Terdapat Autokorelasi pada sisaan
moran.test(reg.klasik$residuals, listw=Woptimum, alternative="two.sided", zero.policy = TRUE)##
## Moran I test under randomisation
##
## data: reg.klasik$residuals
## weights: Woptimum
##
## Moran I statistic standard deviate = 19.214, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.783196536 -0.008474576 0.001697693Kesimpulan:
Berdasarkan hasil uji di atas diperoleh p-value<0.05 yang berarti tolak H0 sehingga tidak terdapat autokorelasi pada sisaan model regresi klasik pada taraf nyata 5%.
Uji Kehomogenan Ragam
H0: ragam sisaan homogen
H1: ragam sisaan tidak homogen
bptest(reg.klasik) #inputnya adalah modelnya##
## studentized Breusch-Pagan test
##
## data: reg.klasik
## BP = 9.3018, df = 3, p-value = 0.02554Kesimpulan:
Karena diperoleh p-value >0.05 maka dapat disimpulkan bahwa ragam sisaan tidak homogen.
Karena ragam sisaan tidak homogen dan terdapat autokorelasi spasial pada sisaan maka akan dilakukan pemodelan spasial dengan menggunakan Regresi Terboboti Geografis (RTG) atau Geographically Weighted Regression (GWR)
Geographically Weighted Regression (GWR)
Suatu pemodelan dapat bersifat global maupun lokal. Regresi linier klasik merupakan salah satu model global. Dikatakan global karena terdapat satu model yang berlaku umum untuk semua pengamatan. Suatu model lokal bersifat lebih fleksibel, yang dalam konteks spasial, artinya setiap daerah/lokasi dapat memiliki model masing-masing.
Geographically Weighted Regression (GWR) atau RTG merupakan salah satu model yang bersifat lokal. Beberapa keuntungan dengen menggunakan model ini, diantaranya adalah kita dapat:
menduga galat baku lokal
menghitung ukuran leverage lokal
melakukan pengujian terhadap signifikansi keragaman spasial pada penduga parameter lokal
menguji apakah model lokal lebih baik daripada model global
Menentukan Bandwith (Lebar Jendela) Optimum
Bandwidth (lebar jendela) merupakan lingkaran dengan radius r dari titik pusat lokasi yang berfungsi sebagai dasar penentuan bobot setiap pengamatan terhadap model regresi setiap lokasi.
Metode pemilihan lebar jendela sangat penting digunakan untuk pendugaan fungsi (kernel) yang tepat. Nilai lebar jendela yang sangat kecil akan mengakibatkan ragam membesar.Sebaliknya, jika nilai lebar jendela besar akan menimbulkan bias yang semakin besar pula.
Untuk itu pemilihan lebar jendela harus dilakukan secara hati-hati agar terpilih lebar jendela yang optimum sehingga ketepatan model yang dibentuk akan tinggi. Salah satu cara untuk memilih lebar jendela yang optimum dapat dilakukan dengan menggunakan metode validasi silang (Cross Validation, CV).
Pada pemilihan bandwith optimum ini digunakan fungsi pembobot kernel bi-square.
b.gwr <- gwr.sel(y~x1+x2+x3,data=datauas,coords = longlat,gweight=gwr.bisquare)## Bandwidth: 3.373553 CV score: 348.9663
## Bandwidth: 5.453074 CV score: 464.0503
## Bandwidth: 2.088339 CV score: 317.7549
## Bandwidth: 1.415407 CV score: 325.2034
## Bandwidth: 2.058408 CV score: 317.5577
## Bandwidth: 1.953423 CV score: 317.0788
## Bandwidth: 1.747919 CV score: 318.5092
## Bandwidth: 1.853998 CV score: 317.3191
## Bandwidth: 1.939104 CV score: 317.0588
## Bandwidth: 1.930691 CV score: 317.0544
## Bandwidth: 1.92806 CV score: 317.0541
## Bandwidth: 1.927527 CV score: 317.0541
## Bandwidth: 1.927652 CV score: 317.0541
## Bandwidth: 1.927693 CV score: 317.0541
## Bandwidth: 1.927611 CV score: 317.0541
## Bandwidth: 1.927652 CV score: 317.0541b.gwr## [1] 1.927652Berdasarkan output di atas, diperoleh bandwith optimum sebesar 1.927652.
Pemodelan RTG dengan Bandwith optimum
rtg.model <- gwr(y~x1+x2+x3, data=datauas, coords = longlat, gweight=gwr.bisquare,bandwidth = b.gwr,hatmatrix=TRUE)
rtg.model## Call:
## gwr(formula = y ~ x1 + x2 + x3, data = datauas, coords = longlat,
## bandwidth = b.gwr, gweight = gwr.bisquare, hatmatrix = TRUE)
## Kernel function: gwr.bisquare
## Fixed bandwidth: 1.927652
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max. Global
## X.Intercept. 23.22921 33.18121 36.49068 41.45227 63.51380 42.3096
## x1 -0.39207 -0.29222 -0.27736 -0.25998 -0.21696 -0.2754
## x2 -0.83797 -0.73152 -0.68524 -0.54073 -0.46302 -0.6663
## x3 10.53279 12.32362 12.44657 12.72839 12.91543 11.6159
## Number of data points: 119
## Effective number of parameters (residual: 2traceS - traceS'S): 25.16185
## Effective degrees of freedom (residual: 2traceS - traceS'S): 93.83815
## Sigma (residual: 2traceS - traceS'S): 1.492291
## Effective number of parameters (model: traceS): 19.97816
## Effective degrees of freedom (model: traceS): 99.02184
## Sigma (model: traceS): 1.452706
## Sigma (ML): 1.325164
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 456.1736
## AIC (GWR p. 96, eq. 4.22): 424.6912
## Residual sum of squares: 208.9712
## Quasi-global R2: 0.9767543Berdasarkan hasil pemodelan RTG di atas, diperoleh penduga parameter bagi peubah penjelas X1 dan X2 yang memiliki nilai minimum dan maksimum negatif sehingga peubah penjelas X1 dan X2 berkontribusi negatif terhadap peubah respon Y di semua lokasi amatan. Sedangkan untuk peubah penjelas X3, nilai minimum dan maksimum penduga parameter memiliki tanda positif, sehingga peubah penjelas X3 berkontribusi positif terhadap peubah respon Y di semua lokasi amatan.
Peta Sebaran Penduga Parameter RTG di Setiap Lokasi
1. Penduga parameter beta1 (X1)
b.x1<- rtg.model$SDF$x1
peta@data$b1<- b.x1
spplot(peta,zcol="b1",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 negatif. Sehingga dapat diartikan bahwa peubah penjelas X1 berkontribusi negatif terhadap peubah respon Y di semua lokasi amatan.
Selain itu, peta sebaran penduga parameter untuk peubah penjelas X1 menunjukkan kemiripan antarkabupaten/kota yang berdekatan. Kemiripan ini ditunjukkan dengan pola kecenderungan warna-warna yang sama untuk mengelompok dengan wilayah tetangganya.
2. Penduga parameter beta2 (X2)
b.x2<- rtg.model$SDF$x2
peta@data$b2<- b.x2
spplot(peta,zcol="b2",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 negatif. Sehingga dapat diartikan bahwa peubah penjelas X2 berkontribusi negatif terhadap peubah respon Y di semua lokasi amatan.
Selain itu, peta sebaran penduga parameter untuk peubah penjelas X2 menunjukkan kemiripan antarkabupaten/kota yang berdekatan. Kemiripan ini ditunjukkan dengan pola kecenderungan warna-warna yang sama untuk mengelompok dengan wilayah tetangganya
3. Penduga parameter beta3 (X3)
b.x3<- rtg.model$SDF$x3
peta@data$b3<- b.x3
spplot(peta,zcol="b3",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 peubah penjelas X3 berkontribusi negatif terhadap peubah respon Y di semua lokasi amatan.
Selain itu, peta sebaran penduga parameter untuk peubah penjelas X3 menunjukkan kemiripan antarkabupaten/kota yang berdekatan. Kemiripan ini ditunjukkan dengan pola kecenderungan warna-warna yang sama untuk mengelompok dengan wilayah tetangganya.
Plot Sebaran Dugaan Peubah Respon Y (Model RTG)
gwr.prediksi <- rtg.model$SDF$pred
peta@data$pred <- gwr.prediksi
peta@data$y <- datauas$y
spplot(peta,c("y","pred"),names.attr=c("Peta Sebaran Peubah Y pada Data","Peta Sebaran Dugaan Peubah Respon Y Model RTG"),as.table = TRUE, main = "Nilai Peubah Y VS Dugaan Y RTG")
#Statistik Amatan Peubah Y dan dugaan model RTG
summary(peta@data[,c("y","pred")])## y pred
## Min. : 6.282 Min. : 6.628
## 1st Qu.:18.781 1st Qu.:19.280
## Median :25.999 Median :26.391
## Mean :26.069 Mean :26.079
## 3rd Qu.:32.926 3rd Qu.:33.034
## Max. :46.840 Max. :45.468Berdasarkan hasil peta sebaran nilai di atas terlihat bahwa nilai dugaan yang dihasilkan dari model GWR tidak jauh berbeda dengan nilai pada data aktual. Selain itu, apabila dilihat dari ringkasan statistik di atas , ringkasan statistik dari nilai dugaan model GWR tidak jauh berbeda dengan nilai data aktual.
Plot Sisaan Model RTG
rtg.sisaan <- datauas$y-rtg.model$SDF$pred
peta@data$rtg.sisaan <- rtg.sisaan
spplot(peta,zcol="rtg.sisaan",main="Peta Sebaran Sisaan Model RTG")
Berdasarkan plot di atas terlihat bahwa model RTG memiliki nilai sisaan yang cenderung kecil (di bawah 2).
Uji Diagnostik Sisaan
Uji Normalitas
H0: sisaan model menyebar normal
H1: sisaan model tidak menyebar normal
err.rtg<-rtg.model$SDF$gwr.e
ad.test(err.rtg) #anderson darling##
## Anderson-Darling normality test
##
## data: err.rtg
## A = 0.64295, p-value = 0.09137hist(err.rtg)
qqnorm(err.rtg,datax=T)
qqline(rnorm(length(err.rtg),mean(err.rtg),sd(err.rtg)),datax=T, col="red")
Kesimpulan:
Berdasarkan QQ-Plot terlihat bahwa banyak titik-titik di sekitar garis yang menunjukkan secara eksploratif bahwa normalitas terpenuhi. Selain itu, berdasarkan uji Anderson-Darling diperoleh p-value>0.05 yang berarti asumsi normalitas terpenuhi (sisaan model menyebar normal).
Uji Kehomogenan Ragam
H0: ragam sisaan homogen
H1: ragam sisaan tidak homogen
bptest(rtg.model$lm, weights = rtg.model$gweight) ##
## studentized Breusch-Pagan test
##
## data: rtg.model$lm
## BP = 9.3018, df = 3, p-value = 0.02554Kesimpulan:
Karena diperoleh p-value >0.05 maka dapat disimpulkan bahwa ragam sisaan homogen.
Uji Autokorelasi Spasial
H0 : Tidak terdapat Autokorelasi pada sisaan
H1 : Terdapat Autokorelasi pada sisaan
gwr.morantest(rtg.model, Woptimum, zero.policy = TRUE)##
## Leung et al. 2000 three moment approximation for Moran's I
##
## data: GWR residuals
## statistic = 40.694, df = 108.21, p-value = 4.095e-10
## sample estimates:
## I
## 0.1406797Kesimpulan:
Berdasarkan hasil uji di atas diperoleh p-value>0.05 yang berarti tidak tolak H0 sehingga tidak terdapat autokorelasi pada sisaan model SEM pada taraf nyata 5%.
Evaluasi Model
Hipotesis:
H0: model RTG dan model RLB sama baik
H1: model RTG lebih baik dari model RLB
BFC99.gwr.test(rtg.model)##
## Brunsdon, Fotheringham & Charlton (1999) ANOVA
##
## data: rtg.model
## F = 20.615, df1 = 87.280, df2 = 99.673, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## SS GWR improvement SS GWR residuals
## 971.5162 208.9712BFC02.gwr.test(rtg.model)##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA
##
## data: rtg.model
## F = 5.649, df1 = 115.000, df2 = 93.838, p-value = 2.688e-16
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 1180.4874 208.9712BFC02.gwr.test(rtg.model,approx = T)##
## Brunsdon, Fotheringham & Charlton (2002, pp. 91-2) ANOVA (approximate
## degrees of freedom - only tr(S))
##
## data: rtg.model
## F = 5.649, df1 = 115.000, df2 = 99.022, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## SS OLS residuals SS GWR residuals
## 1180.4874 208.9712anova(rtg.model)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## OLS Residuals 4.000 1180.49
## GWR Improvement 21.162 971.52 45.909
## GWR Residuals 93.838 208.97 2.227 20.615anova(rtg.model,approx=T)## Analysis of Variance Table
## approximate degrees of freedom (only tr(S))
## Df Sum Sq Mean Sq F value
## OLS Residuals 4.000 1180.49
## GWR Improvement 15.978 971.52 60.803
## GWR Residuals 99.022 208.97 2.110 28.812Berdasarkan hasil pengujian di atas, diperoleh p-value<0.05 sehingga H0 ditolak yang berarti model RTG lebih sesuai untuk memodelkan data dibandingkan dengan regresi linier klasik.
Perbandingan Nilai Dugaan RTG dan Regresi Klasik
reg.klasik.pred <- Predict(reg.klasik,data=datauas)
peta@data$rlbpred <- reg.klasik.pred
spplot(peta,c("pred","rlbpred","y"),names.attr=c("Hasil Dugaan Model RTG","Hasil Dugaan Model Regresi Klasik","Peta Sebaran Nilai Aktual"),as.table = TRUE, main = "Nilai Y vs Dugaan Y Model RTG VS RLB")
#Statistik Amatan Peubah Y, dugaan Y model RTG,dan dugaan Y model RLB
summary(peta@data[,c("y","pred","rlbpred")])## y pred rlbpred
## Min. : 6.282 Min. : 6.628 Min. : 9.539
## 1st Qu.:18.781 1st Qu.:19.280 1st Qu.:19.174
## Median :25.999 Median :26.391 Median :25.723
## Mean :26.069 Mean :26.079 Mean :26.069
## 3rd Qu.:32.926 3rd Qu.:33.034 3rd Qu.:33.066
## Max. :46.840 Max. :45.468 Max. :42.114peta@data$err.regklasik <- err.regklasik
#Statistik Sisaan model RTG dan model RLB
summary(peta@data[,c("err.regklasik","rtg.sisaan")])## err.regklasik rtg.sisaan
## Min. :-8.0899 Min. :-3.047945
## 1st Qu.:-2.1014 1st Qu.:-1.053783
## Median :-0.1042 Median : 0.002010
## Mean : 0.0000 Mean :-0.009175
## 3rd Qu.: 2.1141 3rd Qu.: 1.004264
## Max. : 6.4331 Max. : 2.551854Hasil ringkasan statistik di atas dapat menunjukkan bahwa model RTG lebih sesuai dibandingkan dengan regresi klasik karena nilai prediksi yang dihasilkan oleh RTG lebih mendekati nilai Y aktual, sehingga nilai sisaan pada RTG juga lebih kecil dibandingkan dengan sisaan pada regresi klasik.
Nilai AIC RTG dan Regresi Klasik
data.frame("MODEL" = c("RTG","Regresi Klasik"),
"AIC" = c(rtg.model[["results"]][["AICh"]],AIC(reg.klasik)))%>% arrange(AIC)## MODEL AIC
## 1 RTG 424.6912
## 2 Regresi Klasik 620.7599Berdasarkan nilai AIC dapat disimpulkan bawa model RTG lebih baik dalam memodelkan data daripada regresi klasik karena RTG memiliki nilai AIC yang lebih kecil.