Input Data

Data yang digunakan adalah data kemiskinan di Jawa Barat.

setwd("D:\\Kuliah S2 IPB\\Bahan Kuliah\\Semester 1 SSD 2020\\01 Metode Penelitian Kuantitatif\\uts\\spasial\\SPASIAL DENGAN R\\REGRESI SPASIAL DENGAN r\\Spatial Regression with R") 
# untuk melakukan pengaturan working directory

datajabar <- read_excel("Jabar Data.xlsx", sheet="data") 
# memasukkan file excel ke dalam objek datajabar

Untuk melihat dimensi dari data dapat menggunakan fungsi dim. Terlihat bahwa dimensi dari data tersebut adalah 26 observasi dan 32 variabel/peubah.

dim(datajabar)

## [1] 26 32

Untuk menempilkan peta Jabar.

#Peta Dasar Batas-batas Wilayah dari Jabar
jabar<-readOGR(dsn="petaJabar2", layer="Jabar2")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Kuliah S2 IPB\Bahan Kuliah\Semester 1 SSD 2020\01 Metode Penelitian Kuantitatif\uts\spasial\SPASIAL DENGAN R\REGRESI SPASIAL DENGAN r\Spatial Regression with R\petaJabar2", layer: "Jabar2"
## with 26 features
## It has 7 fields

plot(jabar)

Eksplore Data

Membuat peta sebaran persentase penduduk miskin di Jawa Barat Tahun 2016.

colfunc<-colorRampPalette(c("green", "yellow","red"))
color<-colfunc(16)
jabar$miskin<-datajabar$p.miskin16
spplot(jabar, "miskin", col.regions=color,
       main="Peta Sebaran Persentase Penduduk Miskin 2016")

Menuju warna hijau -> kemiskinan semakin rendah
Menuju warna merah -> kemiskinan semakin tinggi

Hasil plot data ini dapat menjadi identifikasi awal adanya ketergantungan spasial.

Moran’s Indeks

Moran’s indeks merupakan cara eskplorasi untuk melihat adanya autokorelasi spasial. Nilainya antara -1 sampai dengan 1. Jika positif berarti memiliki autokorelasi yang positif. Jika indeks nya makin mendekati 1, maka makin berkorelasi (autokorelasi) sedangkan jika mendekati nol maka datanya diduga tidak berkelompok (menyebar). Perlu dilakukan uji Moran (Moran test) untuk melihat ada tidaknya autokorelasi spasial.

Ho: Tidak terdapat autokorelasi spasial
H1: Terdapat autokorelasi spasial

w<-poly2nb(jabar)
ww<-nb2listw(w) ##ww adalah bobot
moran(datajabar$p.miskin16, ww, n=length(ww$neighbours), S0=Szero(ww))

## $I
## [1] 0.3932657
## 
## $K
## [1] 2.403804

moran.test(datajabar$p.miskin16, ww, alternative="greater")

## 
##  Moran I test under randomisation
## 
## data:  datajabar$p.miskin16  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.0168, p-value = 0.001277
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##         0.3932657        -0.0400000         0.0206265

Nilai indeks moran-nya adalah 0.3932657. Didapatkan p-value= 0.001277 kurang dari alpha=0.05 sehingga Tolak H0 berarti terbukti adanya autokorelasi spasial.

Dapat juga dibuat Moran Plot. Dalam Moran plot, data akan tersebar dalam 4 jenis kuadran, yaitu:

• Kuadran I atau HH (High-High)

• Kuadran II atau LH (Low-High)

• Kuadran III atau LL (Low-Low)

• Kuadran IV atau HL (High-Low)

moran.plot(datajabar$p.miskin16, ww, labels=datajabar$KABKOT)

#labels untuk memberi nama titik dalam sebarannya

Dari Moran plot, terdapat Cirebon yang masuk ke dalam kuadran IV (High-Low) berarti Cirebon memiliki persentase penduduk miskin yang tinggi sedangkan wilayah sekitarnya memiliki persentase penduduk miskin yang rendah.

Regresi Klasik

Membuat model regresi klasik:

reg.klasik<-lm(p.miskin16~EYS2016,datajabar)
summary(reg.klasik)

## 
## Call:
## lm(formula = p.miskin16 ~ EYS2016, data = datajabar)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3627 -1.9604 -0.3899  1.5238  8.1125 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  39.5732     9.0973   4.350 0.000217 ***
## EYS2016      -2.3948     0.7209  -3.322 0.002854 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.788 on 24 degrees of freedom
## Multiple R-squared:  0.315,  Adjusted R-squared:  0.2865 
## F-statistic: 11.04 on 1 and 24 DF,  p-value: 0.002854

Uji Asumsi Klasik dari model

err.regklasik<-residuals(reg.klasik)

#heterogenitas
bptest(reg.klasik) #digunakan modelnya dalam perhitungan

## 
##  studentized Breusch-Pagan test
## 
## data:  reg.klasik
## BP = 1.0644, df = 1, p-value = 0.3022

#dependensi
RunsTest(err.regklasik)

## 
##  Runs Test for Randomness
## 
## data:  err.regklasik
## runs = 13, m = 13, n = 13, p-value = 0.8358
## alternative hypothesis: true number of runs is not equal the expected number
## sample estimates:
##  median(x) 
## -0.3898635

#Normalitas
ad.test(err.regklasik) 

## 
##  Anderson-Darling normality test
## 
## data:  err.regklasik
## A = 0.39412, p-value = 0.3495

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

Dalam uji klasik, seluruh asumsi klasik dalam model regresi terpenuhi.

Spatial Autocorrelation: Uji Moran Index

Selanjutnya, akan dilakukan uji terhadap error dari regresi klasik. Ingin menguji apakah terdapat autokorelasi.

#Moran Test pada error
moran.test(err.regklasik, ww,randomisation=T, alternative="greater")

## 
##  Moran I test under randomisation
## 
## data:  err.regklasik  
## weights: ww    
## 
## Moran I statistic standard deviate = 3.5604, p-value = 0.0001851
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.44799554       -0.04000000        0.01878561

Hasil dari uji tersebut adalah tolak H0 sehingga terdapat autokorelasi spasial. Dari hasil ini dapat diidentifikasi bahwa ada kemungkinan model ini mendapat ketergantungan spasial. Sehingga dapat dilanjutkan dengan uji efek spasial menggunakan uji lm.

Model SAR

Menggunakan fungsi lagsarlm.

#Model SAR
sar<-lagsarlm(p.miskin16~EYS2016,data=datajabar,nb2listw(w))
summary(sar)

## 
## Call:
## lagsarlm(formula = p.miskin16 ~ EYS2016, data = datajabar, listw = nb2listw(w))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.17670 -0.93185 -0.11318  0.91353  7.69131 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 28.22898    7.66250  3.6840 0.0002296
## EYS2016     -1.94997    0.57325 -3.4016 0.0006700
## 
## Rho: 0.59078, LR test value: 7.9343, p-value: 0.0048507
## Asymptotic standard error: 0.1559
##     z-value: 3.7894, p-value: 0.00015101
## Wald statistic: 14.36, p-value: 0.00015101
## 
## Log likelihood: -58.54132 for lag model
## ML residual variance (sigma squared): 4.7513, (sigma: 2.1798)
## Number of observations: 26 
## Number of parameters estimated: 4 
## AIC: 125.08, (AIC for lm: 131.02)
## LM test for residual autocorrelation
## test value: 0.036687, p-value: 0.8481

#asumsi SAR
err.sar<-residuals(sar)
ad.test(err.sar)

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

bptest.sarlm(sar)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.97078, df = 1, p-value = 0.3245

moran.test(err.sar, ww,randomisation=T, alternative="greater")

## 
##  Moran I test under randomisation
## 
## data:  err.sar  
## weights: ww    
## 
## Moran I statistic standard deviate = 0.20765, p-value = 0.4178
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.01331736       -0.04000000        0.01651197

Hasil pengujian asumsi didapatkan nilai p-value > alpha, sehingga semua asumsi tersebut terpenuhi.

Model SEM

Perintah menggunakan errorsarlm.

#Model SEM
sem<-errorsarlm(p.miskin16~EYS2016,data=datajabar,nb2listw(w))
summary(sem)

## 
## Call:errorsarlm(formula = p.miskin16 ~ EYS2016, data = datajabar, 
##     listw = nb2listw(w))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.24699 -1.36122 -0.13809  1.15579  7.03019 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 36.88515    7.80246  4.7274 2.274e-06
## EYS2016     -2.17498    0.60042 -3.6224 0.0002919
## 
## Lambda: 0.61793, LR test value: 8.7676, p-value: 0.0030663
## Asymptotic standard error: 0.1576
##     z-value: 3.9208, p-value: 8.8267e-05
## Wald statistic: 15.372, p-value: 8.8267e-05
## 
## Log likelihood: -58.12466 for error model
## ML residual variance (sigma squared): 4.5459, (sigma: 2.1321)
## Number of observations: 26 
## Number of parameters estimated: 4 
## AIC: 124.25, (AIC for lm: 131.02)

# Uji Asumsi
err.sem<-residuals(sem)
ad.test(err.sem)

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

bptest.sarlm(sem)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.38918, df = 1, p-value = 0.5327

moran.test(err.sem,ww,randomisation=T, alternative="greater")

## 
##  Moran I test under randomisation
## 
## data:  err.sem  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.096995, p-value = 0.5386
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.05288671       -0.04000000        0.01765155

Hasil pengujian asumsi didapatkan nilai p-value > alpha, sehingga semua asumsi tersebut terpenuhi.

Model SARMA/ GSM

Perintahnya adalah sacsarlm.

#Model GSM/SARMA
gsm<-sacsarlm(p.miskin16~EYS2016,data=datajabar,nb2listw(w))
summary(gsm)

## 
## Call:
## sacsarlm(formula = p.miskin16 ~ EYS2016, data = datajabar, listw = nb2listw(w))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0817 -1.2815 -0.1627  1.1317  6.7128 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 38.02312    8.19667  4.6388 3.504e-06
## EYS2016     -2.15311    0.59378 -3.6261 0.0002877
## 
## Rho: -0.14689
## Asymptotic standard error: 0.38983
##     z-value: -0.3768, p-value: 0.70633
## Lambda: 0.69352
## Asymptotic standard error: 0.23014
##     z-value: 3.0135, p-value: 0.0025829
## 
## LR test value: 8.8931, p-value: 0.011719
## 
## Log likelihood: -58.06189 for sac model
## ML residual variance (sigma squared): 4.3254, (sigma: 2.0797)
## Number of observations: 26 
## Number of parameters estimated: 5 
## AIC: 126.12, (AIC for lm: 131.02)

#uji asumsi
err.gsm<-residuals(gsm)
ad.test(err.gsm)

## 
##  Anderson-Darling normality test
## 
## data:  err.gsm
## A = 0.55307, p-value = 0.1388

bptest.sarlm(gsm)

## 
##  studentized Breusch-Pagan test
## 
## data:  
## BP = 0.31097, df = 1, p-value = 0.5771

moran.test(err.gsm,ww,randomisation=T, alternative="greater")

## 
##  Moran I test under randomisation
## 
## data:  err.gsm  
## weights: ww    
## 
## Moran I statistic standard deviate = -0.0061095, p-value = 0.5024
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.04081914       -0.04000000        0.01797655

Hasil pengujian asumsi didapatkan nilai p-value > alpha, sehingga semua asumsi tersebut terpenuhi.

Goodness of Fit

Selanjutnya akan dibandingkan semua model tersebut.

Goodness of Fit

Dari hasil perbandingan semua model tersebut ada kemungkinan Regresi klasik, SAR dan SEM adalah model terbaik.

Terakhir lihat AIC-nya yang terkecil adalah model SEM. Sehingga dapat disimpulkan untuk data kemiskinan 2016 ini model yang terbaik adalah model SEM.

Input Data

Data yang digunakan adalah data GWR Jawa Timur (Data PDRB). Data ini terdapat 38 observasi dan 6 peubah/variabel

datajatim <- read_xlsx("data gwr jatim.xlsx")
dim(datajatim)

## [1] 38  6

#Input Peta
jatim<-readOGR(dsn="Poligon Jatim", layer="poly3500")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Kuliah S2 IPB\Bahan Kuliah\Semester 1 SSD 2020\01 Metode Penelitian Kuantitatif\uts\spasial\SPASIAL DENGAN R\REGRESI SPASIAL DENGAN r\Spatial Regression with R\Poligon Jatim", layer: "poly3500"
## with 38 features
## It has 8 fields

plot(jatim)

Eksplore Data

Eksplore data menggunakan fungsi spplot, dengan warnanya menggunakan palet dari brewer.pal.

jatim@data<-cbind(pdrb=datajatim[,2],jatim@data)
library(RColorBrewer)
my.palette.1<-brewer.pal(n=9,name="YlOrRd")
spplot(jatim[,1], cuts=8, col.regions=my.palette.1, cex=c(0.3,1,0.3), main="PDRB jatim")

Jika dilihat pada peta bahwa semakin ke kanan semakin memerah, sehingga data PDRB untuk wilayah yang semakin ke kanan, PDRB nya semakin besar, semakin ke kiri maka PDRB nya semakin rendah.

Regresi Klasik

#Regresi Klasik
reglin<-lm(Y~x1+x2,data=datajatim)
summary (reglin)

## 
## Call:
## lm(formula = Y ~ x1 + x2, data = datajatim)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -64.616 -13.475  -0.679  10.545 164.572 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -72.45320   18.54899  -3.906 0.000409 ***
## x1            0.10788    0.01952   5.528 3.24e-06 ***
## x2           12.10245    3.51225   3.446 0.001497 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.82 on 35 degrees of freedom
## Multiple R-squared:  0.5671, Adjusted R-squared:  0.5424 
## F-statistic: 22.93 on 2 and 35 DF,  p-value: 4.332e-07

err.ols<-residuals(reglin)

# Uji Asumsi
vif(reglin)

##       x1       x2 
## 1.006723 1.006723

ad.test(err.ols)

## 
##  Anderson-Darling normality test
## 
## data:  err.ols
## A = 1.9608, p-value = 4.258e-05

bptest(reglin)

## 
##  studentized Breusch-Pagan test
## 
## data:  reglin
## BP = 13.836, df = 2, p-value = 0.0009896

Uji asumsi semua nilai vif < 5 sehingga tidak ada korineliaritas.

Uji efek spasial

###UJI EFEK SPASIAL###
#Uji BP
bptest(reglin)

## 
##  studentized Breusch-Pagan test
## 
## data:  reglin
## BP = 13.836, df = 2, p-value = 0.0009896

Hasil uji spasial p_value=0.0009896 < alpha=0.05, dapat dikatakan bahwa terdapat keragaman galat. Karena efek spasial pada data ini terpenuhi maka data ini dapat digunakan untuk model GWR.

Pembobot GWR yang digunakan dalam praktikum kali ini adalah Pembobot Kernel Gaussian.

Pembobot Kernel Gaussian

#coordinat
cord<-cbind(datajatim$Long,datajatim$Lat)
DM <-gw.dist(dp.locat=cord) #matriks jarak

#konvert data spasial
datajatim.sp<-datajatim
coordinates(datajatim.sp)<- ~Long+Lat #set bahwa long lat sebagai koordinat
class(datajatim.sp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

#View(datajatim.sp)

GWR - Bandwith

Untuk mencari bandwith menggunakan fungsi bw.gwr. Bandwith digunakan Untuk menghitung pembobot.

#mencari bandwitch 
bw.fg<-bw.gwr(Y~x1+x2,data=datajatim.sp, kernel = "gaussian",dMat=DM, adaptive = FALSE)

## Fixed bandwidth: 2.568552 CV score: 77431.75 
## Fixed bandwidth: 1.58777 CV score: 76850.32 
## Fixed bandwidth: 0.9816131 CV score: 74446.36 
## Fixed bandwidth: 0.6069877 CV score: 67449.85 
## Fixed bandwidth: 0.3754565 CV score: 57557.27 
## Fixed bandwidth: 0.2323623 CV score: 58990.41 
## Fixed bandwidth: 0.4638935 CV score: 61547.68 
## Fixed bandwidth: 0.3207994 CV score: 55973.74 
## Fixed bandwidth: 0.2870194 CV score: 56026.58 
## Fixed bandwidth: 0.3416765 CV score: 56403.23 
## Fixed bandwidth: 0.3078966 CV score: 55867.61 
## Fixed bandwidth: 0.2999222 CV score: 55876.05 
## Fixed bandwidth: 0.312825 CV score: 55891.64 
## Fixed bandwidth: 0.3048507 CV score: 55863.7 
## Fixed bandwidth: 0.3029682 CV score: 55865.64 
## Fixed bandwidth: 0.3060141 CV score: 55864.17 
## Fixed bandwidth: 0.3041316 CV score: 55864.04 
## Fixed bandwidth: 0.305295 CV score: 55863.73 
## Fixed bandwidth: 0.304576 CV score: 55863.77 
## Fixed bandwidth: 0.3050204 CV score: 55863.69 
## Fixed bandwidth: 0.3051253 CV score: 55863.7 
## Fixed bandwidth: 0.3049556 CV score: 55863.69 
## Fixed bandwidth: 0.3050605 CV score: 55863.69 
## Fixed bandwidth: 0.3049956 CV score: 55863.69

GWR - Model

Untuk model yang digunakan menggunakan fungsi gwr.basic.

###MODEL GWR###
gwr1<-gwr.basic(Y~x1+x2, data=datajatim.sp, bw=bw.fg,kernel = "gaussian", dMat=DM)
gwr1

##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2021-03-22 15:27:50 
##    Call:
##    gwr.basic(formula = Y ~ x1 + x2, data = datajatim.sp, bw = bw.fg, 
##     kernel = "gaussian", dMat = DM)
## 
##    Dependent (y) variable:  Y
##    Independent variables:  x1 x2
##    Number of data points: 38
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##     Min      1Q  Median      3Q     Max 
## -64.616 -13.475  -0.679  10.545 164.572 
## 
##    Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
##    (Intercept) -72.45320   18.54899  -3.906 0.000409 ***
##    x1            0.10788    0.01952   5.528 3.24e-06 ***
##    x2           12.10245    3.51225   3.446 0.001497 ** 
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 36.82 on 35 degrees of freedom
##    Multiple R-squared: 0.5671
##    Adjusted R-squared: 0.5424 
##    F-statistic: 22.93 on 2 and 35 DF,  p-value: 4.332e-07 
##    ***Extra Diagnostic information
##    Residual sum of squares: 47462.13
##    Sigma(hat): 36.30967
##    AIC:  386.7832
##    AICc:  387.9953
##    BIC:  369.8839
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: gaussian 
##    Fixed bandwidth: 0.3049956 
##    Regression points: the same locations as observations are used.
##    Distance metric: A distance matrix is specified for this model calibration.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                     Min.     1st Qu.      Median     3rd Qu.    Max.
##    Intercept -185.220437  -80.448573  -34.683457  -13.604522  4.2308
##    x1           0.010088    0.033628    0.062080    0.122049  0.2896
##    x2          -3.672824    3.685552    8.853189   13.036107 27.4911
##    ************************Diagnostic information*************************
##    Number of data points: 38 
##    Effective number of parameters (2trace(S) - trace(S'S)): 24.9044 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 13.0956 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 430.7659 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 345.5029 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 361.6579 
##    Residual sum of squares: 11516.94 
##    R-square value:  0.894956 
##    Adjusted R-square value:  0.6786741 
## 
##    ***********************************************************************
##    Program stops at: 2021-03-22 15:27:50

Biasanya ditampilkan dalam bentuk boxplot supaya terlihat sebarannya/perbedaannya.

hasil=as.data.frame(gwr1$SDF)
View(hasil) #melihat output model GWR
par(mfrow=c(1,3))
boxplot(hasil[,1],xlab="Intercept") 
boxplot(hasil[,2],xlab="X1")
boxplot(hasil[,3],xlab="X2")

par(mfrow=c(1,1))

GWR- Uji t

Di dalam package ada fungsi khusus uji t menggunakan fungsi gwr.t.adjust

#hasil uji t
ujit <- gwr.t.adjust(gwr1)
a <- ujit$results$t

a

##    Intercept_t      x1_t         x2_t
## 1   0.07500666 0.1013885  0.130699595
## 2  -0.14929048 0.2157656  0.812943235
## 3  -0.24382033 0.2121404  1.100402805
## 4  -0.77809332 0.5073455  1.817783449
## 5  -1.03150976 1.1678547  1.901986651
## 6  -1.28452274 1.1962923  2.077186689
## 7  -1.30932769 1.8582087  1.659964084
## 8  -0.94327543 1.9873844  1.284157922
## 9  -0.12244755 1.1741318  0.071774007
## 10 -0.04982797 1.0021626 -0.222302493
## 11 -0.15878387 1.0351021  0.006655903
## 12 -0.21904872 1.0101595  0.038605042
## 13 -1.48059172 2.2312068  1.622192589
## 14 -2.90120720 5.1851252  2.495743406
## 15 -3.57748166 7.1541738  2.833026493
## 16 -3.57084305 6.6247468  3.180363650
## 17 -3.10016529 4.8917619  2.976466370
## 18 -1.05069180 1.2991928  1.681710514
## 19 -0.43015507 0.6148005  0.919913579
## 20 -0.09618664 0.3231631  0.333303056
## 21 -0.23681017 0.6419804  0.325916469
## 22 -0.87810784 1.7396100  0.928864528
## 23 -1.33830841 3.7635588  0.802950342
## 24 -4.26880431 6.3071611  2.974985041
## 25 -4.09884082 6.1570722  2.909876338
## 26 -4.36032747 4.4683034  0.711817195
## 27 -0.14320952 0.2399921  0.109873140
## 28 -0.11969202 0.2650994 -0.015039046
## 29 -0.36947655 0.4607246  0.276731865
## 30 -1.11517121 0.9381891  1.981325867
## 31 -1.06567257 1.0890099  1.983133030
## 32 -1.69965880 2.5145669  1.951437978
## 33 -1.95150401 3.1030348  1.925933418
## 34 -3.19647847 5.8304793  2.534254030
## 35 -3.81335593 6.6067751  3.293334465
## 36 -0.28338527 0.4823571  0.660346626
## 37 -3.85246507 7.0155083  2.869742511
## 38 -2.13110280 3.3716669  2.424206100

Hasil dari uji-t ini setiap baris mewakili satu wilayah/daerah.

TERIMA KASIH