5. Model Epidemiologi
Deskripsi
Regresi binomial negatif digunakan untuk mengatasi overdisversi pada regresi Poisson. Pada regresi binomial negatif dapat ditambahkan efek dependensi spasial dan dapat dimodelkan menggunakan GLMM (Generalized Linear Mixed Models). Model ini juga dikenal sebagai conditional autoregressive regression (CAR). Pemodelan CAR umumnya menggunakan metode pendugaan Bayes dengan pendugaan sebaran posterior yang efisien menggunakan Integrated Nested Laplace Aproximation (INLA).
GLMM adalah pengembangan GLM yang terdapat penambahan komponen acak pada bagian prediktor model. GLMM dapat digunakan untuk menangani situasi yang terdapat korelasi antar observasi atau terdapat beberapa jenis kelompok dalam data, maupun juga data longitudinal (McCulloch, 2003). Persamaan umum dari GLMM bisa dibentuk dengan memodifikasi model GLM dengan menambahkan efek acak \(\theta_i\) (Tufvesson et al., 2019) sebagai berikut:
\[ g(\mu_i) = X_i^T \beta + \theta_i \tag{1}\]
dengan \(\theta_i\) merupakan efek acak yang independen (unstructured random effect) dapat berupa suatu noise yang memiliki sebaran normal.
Pengembangan lain dari GLMM yaitu mempertimbangkan unsur efek acak spasial terstruktur (structured random effect) pada prediktor melalui prior CAR. Model ini sering digunakan pada beberapa ilmu seperti demografi, ekonomi, epidemiologi, dan geografi (Oliveira et al., 2020). Misalkan peubah respon \(Y = (y_1, \ldots, y_n)\) merupakan vektor univariat yang menyatakan jumlah kasus pada lokasi ke-\(i\) dan vektor predictor pada lokasi ke-\(i\) adalah \(X_i^T = (1, x_{1i}, x_{2i}, \ldots, x_{pi})\) dengan \(i = 1, 2, 3, \ldots, n\) dengan jumlah peubah bebas \(p\). Adapun unsur efek acak spasial terstruktur direpresentasikan oleh \(\phi = (\phi_1, \phi_2, \ldots, \phi_n)\), dengan formulasi model CAR sebagai berikut:
\[ y_i \mid \mu_i \sim f(y_i \mid \mu_i, v^2) \tag{2}\]
\[ g(\mu_i) = X_i^T \beta + \phi_i \tag{3}\]
Seperti halnya pada GLM, peubah respon \(y_i\) berasal dari anggota keluarga sebaran eksponensial \(f(y_i \mid \mu_i, v^2)\) dengan \(v^2\) merupakan parameter skala yang digunakan untuk sebaran normal. Untuk sebaran binomial negatif dengan dua parameter, dengan asumsi tambahan \(v^2 = \gamma\), maka dari Persamaan 2 dan Persamaan 3 dapat dibentuk model CAR dengan fungsi sebaran binomial negatif yaitu sebagai berikut (Sparks, 2020):
\[ y_i \mid \mu_i \sim \text{NB}(y_i \mid \mu_i, \gamma) \tag{4}\]
\[ \log(\mu_i) = X_i^T \beta + \phi_i \tag{5}\]
Efek acak \(\phi_i\) dimodelkan dengan kelas sebaran prior CAR, yang merupakan jenis sebaran Gaussian Markov Random Field (GMRF). Sebaran tersebut dapat ditulis dalam bentuk \(\phi \sim N(0, \tau^2 Q^{-1})\), dengan \(Q\) merupakan matriks presisi yang mungkin berupa matriks tunggal (noninvertible matrix). Korelasi spasial antara efek acak ditentukan oleh matriks ketetanggaan biner \(W\) berukuran \(n \times n\), dengan elemen \(w_{ji}\) sama dengan satu jika area \((j,i)\) didefinisikan sebagai tetangga, dan nol jika tidak. Kemudian apabila dua area didefinisikan sebagai tetangga, maka efek randomnya akan berkorelasi, sedangkan efek random di area yang tidak bertetangga akan bebas bersyarat dengan mempertimbangkan elemen-elemen lain dari \(\phi\). Pendekatan ini paling umum yaitu, \(w_{ji}=1\) jika dan hanya jika mereka memiliki batas yang sama, yaitu dilambangkan dengan \(j \sim i\) sebagai notasi selanjutnya. Namun, setiap area harus memiliki setidaknya satu elemen positif \(w_{ji}\), sehingga jumlah suatu baris dari matriks \(W\) tidak boleh nol.
Umumnya prior CAR dispesifikasikan oleh satu set data berjumlah \(n\) memiliki sebaran bersyarat penuh univariat \(f(\phi_i \mid \phi_{-i})\) di mana \(\phi_{-i} = (\phi_1, \ldots, \phi_{i-1}, \phi_{i+1}, \ldots, \phi_i)\). Berbagai sebaran prior CAR telah banyak diperkenalkan oleh berbagai ahli terutama dalam konteks pemetaan penyakit, yaitu seperti Intrinsic CAR (ICAR), Besag York Mollié (BYM) (Besag et al., 1991), (Stern & Cressie, 1999), dan (Leroux et al., 2000). Namun pada penelitian ini hanya akan membahas prior ICAR dan BYM.
1. ICAR
Prior CAR yang paling sederhana saat ini adalah Intrinsic autoregressive (IAR), yang diusulkan oleh Besag et al. (1991) dan memiliki sebaran bersyarat penuh yaitu,
\[ \phi_i = \psi_i \tag{6}\]
\[ \psi_i \mid \psi_{-i}, W, \tau^2 \sim N\left(\frac{\sum_{j \sim i} w_{ji} \psi_i}{\sum_{j \sim i} w_{ji}}, \frac{\tau^2}{\sum_{j \sim i} w_{ji}}\right) \tag{7}\]
\[ \tau^2 \sim \text{Inverse-Gamma}(a, b) \tag{8}\]
Nilai harapan bersyarat dari \(\psi_i\) sama dengan rata-rata dari efek acak di area tetangga, sedangkan ragam bersyarat berbanding terbalik dengan jumlah tetangga. Kemudian struktur ragam dari \(\psi_i\) ini akan memiliki korelasi spasial yang kuat, apabila semakin banyak tetangga yang dimiliki oleh suatu area. Parameter ragam \(\tau^2\) digunakan untuk mengontrol besarnya variasi antara efek acak, yang pada umumnya mengikuti proses sebaran invers gamma.
2. BYM
Prior konvolusi atau Besag-York-Mollie (BYM) pertama kali dijelaskan oleh Besag et al. (1991). Memiliki efek acak yang dari dua komponen, yaitu mengkombinasikan efek acak spasial terstruktur (structured random effect) pada Persamaan Persamaan 6 dengan efek acak yang independen (unstructured random effect) seperti pada Persamaan 7. Model ini memiliki bentuk sebagai berikut,
\[ \phi_i = \psi_i + \theta_i \tag{9}\]
\[ \psi_i \mid \psi_{-i}, W, \tau^2 \sim N\left(\frac{\sum_{j \sim i} w_{ji} \psi_i}{\sum_{j \sim i} w_{ji}}, \frac{\tau^2}{\sum_{j \sim i} w_{ji}}\right) \tag{10}\]
\[ \theta_i \sim N(0, \sigma^2) \tag{11}\]
\[ \tau^2, \sigma^2 \sim \text{Inverse-Gamma}(a, b) \tag{12}\]
Efek acak \(\theta = (\theta_1, \theta_2, \ldots, \theta_n)\) bersifat independen dengan mean nol dan ragam konstan, sedangkan efek acak autokorelasi spasial dimodelkan melalui \(\psi\). Model ini merupakan model CAR yang paling sering digunakan dalam praktiknya.
Data
Penelitian ini menggunakan data Profil Kesehatan 2021 dari setiap provinsi di Pulau Jawa, yaitu Provinsi Banten, DKI Jakarta, Jawa Barat, DI Yogyakarta, Jawa Tengah, dan Jawa Timur. Satuan pengamatan pada penelitian ini adalah kabupaten/ kota di seluruh Pulau Jawa, yakni 119 kabupaten/kota. Peubah yang digunakan pada penelitian ini dapat dilihat pada Tabel 1. Penelitian selengkapnya terdapat (Fitri et al., 2024).
| Peubah | Keterangan | Faktor | Literatur |
|---|---|---|---|
| \(Y\) | Banyaknya stunting per 100.000 balita | ||
| \(X_1\) | Persentase balita yang mendapatkan ASI eksklusif | Faktor asupan gizi | Hasiru et al, (2022) |
| \(X_2\) | Persentase balita yang mendapatkan imunisasi dasar lengkap | Faktor pelayanan kesehatan | Manaf et al, (2022) |
| \(X_3\) | Persentase keluarga dengan akses sanitasi layak | Faktor pola asuh | Fadliana dan Drajat (2021) |
| \(X_4\) | Persentase penduduk miskin | Faktor ekonomi | Bele et al, (2022) |
Tahapan Analisis Data
Metode yang digunakan dalam analisis data stunting Negative Binomial Conditional Autoregressive (CAR) dengan pendekatan INLA. Pengolahan data menggunakan software R Studio. Tahapan analisis data yang dilakukan adalah sebagai berikut:
Melakukan eksplorasi pada data penelitian
Medeteksi multikolinieritas antar peubah penjelas menggunakan nilai VIF (Variance Inflation Factor), jika didapatkan nilai
VIF > 10maka dapat dikatakan bahwa terjadi multikolinearitas pada data (Kutner MH, 2005), sehingga perlu dilakukan penanganan lebih lanjut sebelum dilakukannya pemodelan.Memeriksa overdispersi, dengan melihat nilai devians dan chi-kuadrat pearson dibagi dengan derajat bebasnya. Jika nilai tersebut bernilai lebih dari satu maka terdapat overdispersi pada data, sehingga dapat menggunakan sebaran yang lebih fleksibel seperti sebaran binomial negatif.
-
Menghitung Matriks Pembobotan Spasial (\(W\)) dengan menggunakan queen contiguity, eksponensial jarak (exponential weight), kebalikan jarak (inverse distance weight), dan pembobot tetangga terdekat (k-Nearest Neighbor Weight) berikut rincian setiap pembobotan yang dilakukan (Djuraidah et al, 2022):
-
Untuk queen contiguity, pembobotan dilakukan dengan menggunakan rumus berikut:
\[ w_{ij} = \begin{cases} 1 & \text{jika } i \text{ dan } j \text{ bersinggungan} \\ 0 & \text{jika } i \text{ dan } j \text{ tidak bersinggungan} \end{cases} \]
-
Untuk bobot eksponensial (exponential weight), pembobotan dilakukan dengan menggunakan rumus berikut:
\[ w_{ij} = \exp(-d_{ij}) \]
dengan \(d_{ij}\) adalah jarak antara luas ke-\(i\) dan luas ke-\(j\).
-
Untuk bobot jarak terbalik (inverse distance weight), pembobotan dilakukan dengan menggunakan rumus berikut:
\[ w_{ij} = d_{ij}^{-1} \]
-
Untuk pembobot KNN:
- Hitung jarak pusat antara unit ke-\(i\) terhadap seluruh unit lainnya \(j \neq i\).
- Beri peringkat sebagai berikut \(d_{ij}(1) \leq d_{ij}(2) \leq \dots \leq d_{ij}(n-1)\).
- Kemudian untuk setiap \(k=1,\dots,n-1\), atur \(N_k (i) = \{j(1),j(2),\dots,j(k)\}\) yang berisi \(k\) unit terdekat terhadap \(i\).
- Untuk setiap \(k\), matriks pembobot \(W\) memiliki elemen \(w_{ij}\) bernilai 1 jika daerah \(i\) berdekatan dengan daerah \(j\), sedangkan elemen diagonal utama akan selalu bernilai nol.
-
-
Melakukan uji autokorelasi spasial dengan menghitung nilai indeks Moran pada peubah respon (Anselin, 1988), dengan langkah-langkah berikut:
Hipotesis
\[ H_0 : \text{Tidak terdapat autokorelasi spasial pada data} \]
\[ H_1 : \text{Terdapat autokorelasi spasial pada data} \]
Statistik Uji \[ z = \frac{(I - E(I))}{\sqrt{\text{var}(I)}} \]
dengan \[ I = \frac{\sum_{i=1}^{n}\sum_{j=1}^{n}W_{ij}(x_i-\bar{x})(x_j-\bar{x})}{S_0 \sum_{i=1}^{n}(x_i-\bar{x})^2} \]
Kriteria pengujiannya \(H_0\) ditolak jika nilai \(|Z| > Z_{(\alpha/2)}\)
-
Melakukan pendugaan model Negative Binomial CAR dengan pendekatan INLA,
6.1 Membuat model yang terdapat empat model yang akan dilihat (Sparks, 2018):
Model GLM: \[ \eta_i = X_i^T \beta \]
Model GLMM: \[ \eta_i = X_i^T \beta + \theta_i \]
Model ICAR: \[ \eta_i = X_i^T \beta + \psi_i \]
Model BYM: \[ \eta_i = X_i^T \beta + \psi_i + \theta_i \]
6.2 Menentukan sebaran prior untuk hiperparameter pada model dengan ketentuan nilai prior:
Model GLM, prior ~ Normal Structured Random Effect, \(\psi_i|\psi_{(-i)} \sim N\left(\frac{\sum_{j\sim i} w_{ji} \psi_i}{\sum_{j\sim i} w_{ji}}, \frac{\tau^2}{\sum_{j\sim i} w_{ji}}\right)\) Unstructured Random Effect, \(\theta_i \sim N(0, \sigma^2)\) Ragam \(\psi_i\), \(\tau^2 \sim \text{Inverse-Gamma}(1, 0.0005)\) Ragam \(\theta_i\), \(\sigma^2 \sim \text{Inverse-Gamma}(1, 0.0005)\)
-
Memilih model terbaik dengan mempertimbangkan beberapa kriteria sebagai berikut:
Defiance Information Criterion (DIC) (Wang et al, 2018) \[ \text{DIC} = \overline{D} + 2PD \]
dengan
\[ PD = E(D(\theta)) - D(E(\theta)) = \overline{D} - D(\theta_{\overline{}}) \]
dan
\[ D(\theta) = -2\log(p(y^*|\theta)) \]
Semakin kecil nilai AIC, maka model yang dibangun semakin baik.
Mean Absolute Deviance (MAD)
\[ \text{MAD} = \frac{\sum_{i=1}^{n}|y_i-\hat{y}_i|}{N} \]
dengan \(y_i\) adalah respon pada lokasi ke \(i\), \(\hat{y}_i\) adalah hasil prediksi pada lokasi ke \(i\) dan \(N\) adalah jumlah lokasi pengamatan. Semakin kecil nilai MAD, maka semakin kecil kesalahan hasil pendugaan.
-
Menghitung nilai resiko relatif dari model terbaik berdasarkan kriteria dari kebaikan model pada langkah ke-7. Risiko relatif ditentukan dengan cara mengeksponensialkan komponen acak terstruktur spasial \(\psi_i\) dari model terbaik menggunakan rumus sebagai berikut (Blangiardo dan Cameletti, 2015) dengan formula:
\[ \text{RR}_i = e^{\psi_i} \]
dengan \(\psi_i\) adalah efek acak terstruktur spasial.
Menginterpretasikan hasil penelitian dan menarik kesimpulan
Package
Data
jawa <- readOGR(dsn = "data/shp", layer = "jawa2")
#> Warning in readOGR(dsn = "data/shp", layer = "jawa2"): OGR support
#> is provided by the sf and terra packages among others
#> Warning in ogrInfo(dsn = dsn, layer = layer, encoding = encoding,
#> use_iconv = use_iconv, : OGR support is provided by the sf and
#> terra packages among others
#> Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is
#> provided by the sf and terra packages among others
#> Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI,
#> dumpSRS = dumpSRS, : OGR support is provided by the sf and terra
#> packages among others
#> Warning in ogrListLayers(dsn): OGR support is provided by the sf
#> and terra packages among others
#> Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is
#> provided by the sf and terra packages among others
#> OGR data source with driver: ESRI Shapefile
#> Source: "C:\Users\anugraha\Documents\Materi_Orasi\book\data\shp", layer: "jawa2"
#> with 119 features
#> It has 1 fields
jawa_data <- read_excel("data/data_stunting_epidemiologi.xlsx")
jawa_data <- jawa_data %>%
mutate(logStunting = log(Stunting), E_d = mean(Stunting))
jawa@data <- left_join(jawa@data, jawa_data)
#> Joining with `by = join_by(KabKota)`
jawa@data$ID <- 1:nrow(jawa@data)
summary(jawa_data)
#> Provinsi KabKota Stunting
#> Length:119 Length:119 Min. : 1
#> Class :character Class :character 1st Qu.: 22
#> Mode :character Mode :character Median : 40
#> Mean : 57
#> 3rd Qu.: 67
#> Max. :368
#> Imunisasi IMD ASI BBLR
#> Min. : 0.0 Min. : 30.6 Min. : 41.1 Min. : 0.0
#> 1st Qu.: 80.8 1st Qu.: 75.4 1st Qu.: 61.9 1st Qu.: 2.8
#> Median : 93.9 Median : 82.7 Median : 72.6 Median : 4.0
#> Mean : 88.3 Mean : 80.7 Mean : 70.4 Mean : 4.7
#> 3rd Qu.: 99.0 3rd Qu.: 89.6 3rd Qu.: 79.3 3rd Qu.: 5.3
#> Max. :138.5 Max. :110.3 Max. :100.0 Max. :66.4
#> Sanitasi TPM Miskin RLS
#> Min. : 4.8 Min. : 9.2 Min. : 2.57 Min. : 3.48
#> 1st Qu.: 90.2 1st Qu.: 51.4 1st Qu.: 7.47 1st Qu.: 7.18
#> Median : 97.7 Median : 62.8 Median :10.21 Median : 8.10
#> Mean : 92.9 Mean : 61.7 Mean :10.34 Mean : 8.47
#> 3rd Qu.:100.0 3rd Qu.: 74.5 3rd Qu.:12.69 3rd Qu.:10.01
#> Max. :232.2 Max. :100.0 Max. :23.76 Max. :11.89
#> kepadatan kb_aktif posyandu Pengeluaran
#> Min. : 280 Min. : 40.8 Min. : 0.1 Min. : 7829
#> 1st Qu.: 782 1st Qu.: 68.8 1st Qu.: 67.4 1st Qu.: 9830
#> Median : 1124 Median : 72.9 Median : 81.1 Median :10942
#> Mean : 3242 Mean : 73.6 Mean : 76.7 Mean :11749
#> 3rd Qu.: 3312 3rd Qu.: 78.3 3rd Qu.: 90.8 3rd Qu.:12762
#> Max. :20360 Max. :145.5 Max. :118.5 Max. :23888
#> Pengangguran Gini Kode Lat
#> Min. : 2.04 Min. :0.264 Min. : 1.0 Min. :-8.22
#> 1st Qu.: 4.79 1st Qu.:0.320 1st Qu.: 30.5 1st Qu.:-7.66
#> Median : 6.55 Median :0.340 Median : 60.0 Median :-7.19
#> Mean : 6.86 Mean :0.348 Mean : 60.0 Mean :-7.18
#> 3rd Qu.: 9.10 3rd Qu.:0.370 3rd Qu.: 89.5 3rd Qu.:-6.83
#> Max. :13.07 Max. :0.489 Max. :119.0 Max. :-5.75
#> Long kode_prov logStunting E_d
#> Min. :106 Min. :31.0 Min. :0.00 Min. :56.6
#> 1st Qu.:108 1st Qu.:32.0 1st Qu.:3.11 1st Qu.:56.6
#> Median :110 Median :33.0 Median :3.69 Median :56.6
#> Mean :110 Mean :33.6 Mean :3.62 Mean :56.6
#> 3rd Qu.:112 3rd Qu.:35.0 3rd Qu.:4.20 3rd Qu.:56.6
#> Max. :114 Max. :36.0 Max. :5.91 Max. :56.6palette(terrain.colors(6))
plot(jawa, col=jawa$kode_prov-30)
text(jawa,'Kode',cex=0.5)
Exploratory Data Analysis
Stunting
my.palette.1<-brewer.pal(n=4,name="OrRd")
spplot(jawa,'Stunting', cuts=3, col.regions=my.palette.1, cex=c(0.3,1,0.3), main="Banyaknya Stunting per 100.000 di Pulau Jawa", at = c(min(jawa$Stunting) - 1, 50, 100, 150, max(jawa$Stunting) + 1))
Insidensi stunting tinggi banyak terjadi di Jawa Barat, dengan kasus tertinggi di Kabupaten Bogor.
ASI
my.palette.1<-brewer.pal(n=4,name="OrRd")
spplot(jawa,'ASI', cuts=3, col.regions=my.palette.1, cex=c(0.3,1,0.3), main="Persentase tingkat balita yang mendapatkan ASI eksklusif di Pulau Jawa")
Imunisasi
Sanitasi
Miskin
my.palette.1<-brewer.pal(n=4,name="OrRd")
spplot(jawa,'Miskin', cuts=3, col.regions=my.palette.1, cex=c(0.3,1,0.3), main="Persentase Penduduk Miskin di Pulau Jawa")
Correlation plot
eda_data <- jawa_data %>%
dplyr::select(Stunting, ASI,Imunisasi, Sanitasi, Miskin)
PerformanceAnalytics::chart.Correlation(eda_data, histogram = TRUE, pch = 19)
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
Correlation Matrix
round(cor(eda_data), 4)
#> Stunting ASI Imunisasi Sanitasi Miskin
#> Stunting 1.0000 -0.0661 -0.0606 -0.1665 0.1187
#> ASI -0.0661 1.0000 0.1349 -0.0172 0.0292
#> Imunisasi -0.0606 0.1349 1.0000 -0.0314 -0.0364
#> Sanitasi -0.1665 -0.0172 -0.0314 1.0000 -0.0176
#> Miskin 0.1187 0.0292 -0.0364 -0.0176 1.0000Correlation plot
eda_data <- jawa_data %>%
dplyr::select(ASI, Imunisasi, Sanitasi, Miskin)
PerformanceAnalytics::chart.Correlation(eda_data, histogram = TRUE, pch = 19)
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
#> Warning in par(usr): argument 1 does not name a graphical parameter
Correlation Matrix
Multicollinearity Test
Dari nilai VIF di atas yang kurang dari 10 menunjukkan bahwa tidak terdapat multikolinieritas antar peubah penjelas.
Overdispersion Test
#Cek Overdispersi
poisson<-glm(jawa_data$Stunting ~ jawa_data$ASI + jawa_data$Imunisasi + jawa_data$Sanitasi + jawa_data$Miskin, data=jawa_data, family = poisson())
model1 <- glm(jawa_data$Stunting~1, data=jawa_data, family = poisson())
dispersiontest(poisson)
#>
#> Overdispersion test
#>
#> data: poisson
#> z = 3, p-value = 0.002
#> alternative hypothesis: true dispersion is greater than 1
#> sample estimates:
#> dispersion
#> 55.2
summary(poisson)
#>
#> Call:
#> glm(formula = jawa_data$Stunting ~ jawa_data$ASI + jawa_data$Imunisasi +
#> jawa_data$Sanitasi + jawa_data$Miskin, family = poisson(),
#> data = jawa_data)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 4.988258 0.096277 51.81 < 2e-16 ***
#> jawa_data$ASI -0.003785 0.000967 -3.91 9.1e-05 ***
#> jawa_data$Imunisasi -0.002410 0.000581 -4.15 3.4e-05 ***
#> jawa_data$Sanitasi -0.008468 0.000601 -14.09 < 2e-16 ***
#> jawa_data$Miskin 0.028068 0.002810 9.99 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for poisson family taken to be 1)
#>
#> Null deviance: 5011.7 on 118 degrees of freedom
#> Residual deviance: 4685.6 on 114 degrees of freedom
#> AIC: 5346
#>
#> Number of Fisher Scoring iterations: 5
anova(model1, poisson, test="Chisq")
#> Analysis of Deviance Table
#>
#> Model 1: jawa_data$Stunting ~ 1
#> Model 2: jawa_data$Stunting ~ jawa_data$ASI + jawa_data$Imunisasi + jawa_data$Sanitasi +
#> jawa_data$Miskin
#> Resid. Df Resid. Dev Df Deviance Pr(>Chi)
#> 1 118 5012
#> 2 114 4686 4 326 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Diperoleh p-value < 0.05, yang artinya bahwa nilai dispersi lebih dari 1. Sehingga dapat disimpulkan bahwa terdapat dispersi pada data.
#Cek Sebaran
poisson <- fitdist(jawa_data$Stunting, "pois", method = c("mle", "mme", "qme", "mge", "mse"),
start=NULL, fix.arg=NULL, T, keepdata = TRUE)
binom.negatif <- fitdist(jawa_data$Stunting, "nbinom", method = c("mle", "mme", "qme", "mge", "mse"),
start=NULL, fix.arg=NULL, T, keepdata = TRUE)
qqcomp(list(poisson, binom.negatif), legendtext = c("Pois", "Neg Binom"), fitpch="o", fitcol = c("red","green"))
descdist(jawa_data$Stunting, discrete = TRUE, boot=1000)
#> summary statistics
#> ------
#> min: 1 max: 368
#> median: 40
#> mean: 56.6
#> estimated sd: 56.6
#> estimated skewness: 2.66
#> estimated kurtosis: 12.6
Spatial Weighted Matrix
Queen continguity
wm_q <- poly2nb(jawa, queen = TRUE)
#> Warning in poly2nb(jawa, queen = TRUE): some observations have no neighbours;
#> if this seems unexpected, try increasing the snap argument.
#> Warning in poly2nb(jawa, queen = TRUE): neighbour object has 3 sub-graphs;
#> if this sub-graph count seems unexpected, try increasing the snap argument.
summary(wm_q)
#> Neighbour list object:
#> Number of regions: 119
#> Number of nonzero links: 520
#> Percentage nonzero weights: 3.67
#> Average number of links: 4.37
#> 1 region with no links:
#> 0
#> 3 disjoint connected subgraphs
#> Link number distribution:
#>
#> 0 1 2 3 4 5 6 7 8 9 11
#> 1 15 10 13 22 20 20 11 4 2 1
#> 15 least connected regions:
#> 3 45 46 47 49 51 53 54 55 56 58 59 61 62 109 with 1 link
#> 1 most connected region:
#> 12 with 11 links- Terdapat 119 Kabupaten/Kota di Pulau Jawa.
- Terdapat satu daerah yang tidak bersinggungan dengan daerah lain, yaitu Kepulauan Seribu.
- Daerah yang paling banyak bersinggungan dengan daerah lain, yaitu Bogor yang bersinggungan dengan 11 daerah.
plot(jawa, borders = 'lightgrey', main = "Queen Contiguity Based Neighbours Maps")
#> Warning in title(...): "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
#> Warning in polypath(x = mcrds[, 1], y = mcrds[, 2], border =
#> border, col = col, : "borders" is not a graphical parameter
plot(wm_q, coordinates(jawa), pch = 19, cex = 0.6, add = TRUE, col = "red")
Exponential Weight
coords <- coordinates(jawa)
dist <- nbdists(wm_q, coords, longlat = TRUE)
eds <- lapply(dist, function(x) exp(-1*x))
head(eds)
#> [[1]]
#> numeric(0)
#>
#> [[2]]
#> [1] 2.42e-14 2.36e-22 5.84e-16 2.15e-09 2.36e-11 9.44e-31 6.80e-23
#>
#> [[3]]
#> [1] 2.42e-14 1.40e-17 1.88e-11 7.20e-07 6.35e-15 1.88e-25
#>
#> [[4]]
#> [1] 1.94e-16
#>
#> [[5]]
#> [1] 2.43e-24 2.23e-19 1.95e-15 5.68e-15 8.69e-13 4.01e-13
#>
#> [[6]]
#> [1] 8.34e-14 1.60e-05 1.45e-10 1.55e-09Inverse Distance Weight
ids <- lapply(dist, function(x) 1/(x))
head(ids)
#> [[1]]
#> numeric(0)
#>
#> [[2]]
#> [1] 0.0319 0.0201 0.0285 0.0501 0.0409 0.0145 0.0196
#>
#> [[3]]
#> [1] 0.0319 0.0258 0.0405 0.0707 0.0306 0.0176
#>
#> [[4]]
#> [1] 0.0276
#>
#> [[5]]
#> [1] 0.0184 0.0233 0.0295 0.0305 0.0360 0.0350
#>
#> [[6]]
#> [1] 0.0332 0.0905 0.0441 0.0493Spatial Autocorrelation Test
rswm_q <- nb2listw(wm_q, style = "W", zero.policy = TRUE)
rswm_eds <- nb2listw(wm_q, glist = eds, style = "B", zero.policy = TRUE)
rswm_ids <- nb2listw(wm_q, glist = ids, style = "B", zero.policy = TRUE)
moran_q <- moran.test(jawa$logStunting, listw = rswm_q, zero.policy = TRUE)
moran_q
#>
#> Moran I test under randomisation
#>
#> data: jawa$logStunting
#> weights: rswm_q
#> n reduced by no-neighbour observations
#>
#> Moran I statistic standard deviate = 1, p-value = 0.2
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.05915 -0.00855 0.00454
moran_eds <- moran.test(jawa$logStunting, listw = rswm_eds, zero.policy = TRUE)
moran_eds
#>
#> Moran I test under randomisation
#>
#> data: jawa$logStunting
#> weights: rswm_eds
#> n reduced by no-neighbour observations
#>
#> Moran I statistic standard deviate = -0.6, p-value = 0.7
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> -0.55751 -0.00855 0.79378
moran_ids <- moran.test(jawa$logStunting, listw = rswm_ids, zero.policy = TRUE)
moran_ids
#>
#> Moran I test under randomisation
#>
#> data: jawa$logStunting
#> weights: rswm_ids
#> n reduced by no-neighbour observations
#>
#> Moran I statistic standard deviate = -0.3, p-value = 0.6
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> -0.06538 -0.00855 0.03140| Weighted Matrix | Moran’s Index | E(I) | Var(I) | p-value |
|---|---|---|---|---|
| Queen Continuity | 0.0591 | -0.0085 | 0.0045 | 0.1576 |
| Exponential Distance | -0.5575 | -0.0085 | 0.7938 | 0.7311 |
| Inverse Distance | -0.0654 | -0.0085 | 0.0314 | 0.6258 |
Matriks Bobot
# Distance Matrix
longlat<-cbind(jawa_data$Long ,jawa_data$Lat)
gdist<-pointDistance(longlat,lonlat=TRUE)
m.gdist<-as.matrix(gdist)
djarak<-dist(longlat)
m.djarak<-as.matrix(djarak)#K-Nearest Neighbour Weight dengan k=4
koord <- coordinates(jawa)
W1<-knn2nb(knearneigh(longlat,k=4,longlat=TRUE))
WW1<- nb2listw(W1,style='W')
WW1
#> Characteristics of weights list object:
#> Neighbour list object:
#> Number of regions: 119
#> Number of nonzero links: 476
#> Percentage nonzero weights: 3.36
#> Average number of links: 4
#> Non-symmetric neighbours list
#>
#> Weights style: W
#> Weights constants summary:
#> n nn S0 S1 S2
#> W 119 14161 119 52 496MI1 <- moran.test(jawa_data$Stunting,WW1) Uji Moran
moran.test(jawa_data$Stunting, WW1,randomisation=T,
alternative="two.sided")
#>
#> Moran I test under randomisation
#>
#> data: jawa_data$Stunting
#> weights: WW1
#>
#> Moran I statistic standard deviate = 2, p-value = 0.05
#> alternative hypothesis: two.sided
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.10526 -0.00847 0.00323Berdasarkan uji Indeks Moran di atas, diperoleh p-value < 0.05, artinya bahwa \(H_0\) ditolak. Sehingga dapat disimpulkan bahwa terdapat autokorelasi spasial pada data.
Model Estimation
GLM
glm_mod <- inla(form_fit, data = jawa@data, family = "nbinomial", E = E_d,
control.predictor = list(compute = TRUE, link = 1),
control.compute = list(dic = TRUE))
summary(glm_mod)
#> Time used:
#> Pre = 0.35, Running = 0.246, Post = 24.9, Total = 25.5
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> (Intercept) 1.109 0.692 -0.248 1.109 2.468 1.109 0
#> ASI -0.005 0.006 -0.018 -0.005 0.007 -0.005 0
#> Imunisasi -0.003 0.004 -0.011 -0.003 0.005 -0.003 0
#> Sanitasi -0.009 0.004 -0.017 -0.009 -0.002 -0.009 0
#> Miskin 0.037 0.020 -0.003 0.037 0.077 0.037 0
#>
#> Model hyperparameters:
#> mean sd
#> size for the nbinomial observations (1/overdispersion) 1.51 0.184
#> 0.025quant
#> size for the nbinomial observations (1/overdispersion) 1.18
#> 0.5quant
#> size for the nbinomial observations (1/overdispersion) 1.50
#> 0.975quant
#> size for the nbinomial observations (1/overdispersion) 1.90
#> mode
#> size for the nbinomial observations (1/overdispersion) 1.49
#>
#> Deviance Information Criterion (DIC) ...............: 1196.11
#> Deviance Information Criterion (DIC, saturated) ....: -6412.71
#> Effective number of parameters .....................: 6.00
#>
#> Marginal log-Likelihood: -633.55
#> is computed
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')glm1<-glm.nb(form_fit, data = jawa@data)
summary(glm1)
#>
#> Call:
#> glm.nb(formula = form_fit, data = jawa@data, init.theta = 1.502309197,
#> link = log)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 5.08798 0.65683 7.75 9.5e-15 ***
#> ASI -0.00472 0.00601 -0.79 0.432
#> Imunisasi -0.00293 0.00384 -0.76 0.446
#> Sanitasi -0.00939 0.00382 -2.46 0.014 *
#> Miskin 0.03701 0.01812 2.04 0.041 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for Negative Binomial(1.5) family taken to be 1)
#>
#> Null deviance: 140.70 on 118 degrees of freedom
#> Residual deviance: 130.81 on 114 degrees of freedom
#> AIC: 1196
#>
#> Number of Fisher Scoring iterations: 1
#>
#>
#> Theta: 1.502
#> Std. Err.: 0.186
#>
#> 2 x log-likelihood: -1184.083Residuals autocorrelation
moran.test(residuals(glm_mod)[[1]], listw = WW1, na.action = na.omit, zero.policy = T)
#>
#> Moran I test under randomisation
#>
#> data: residuals(glm_mod)[[1]]
#> weights: WW1
#>
#> Moran I statistic standard deviate = 0.6, p-value = 0.3
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.02844 -0.00847 0.00348Diperoleh p-value < 0.05, yang artinya bahwa terdapat autokorelasi spasial pada galat model GLM.
GLMM
form_fit2 <- Stunting ~ Imunisasi + ASI + Sanitasi + Miskin +
f(ID, model = "iid", param = c(1, .5))
glmm_mod <- inla(form_fit2, data = jawa@data, family = "nbinomial", E = E_d,
control.predictor = list(compute = TRUE, link = 1),
control.compute = list(dic = TRUE))
summary(glmm_mod)
#> Time used:
#> Pre = 0.222, Running = 0.291, Post = 31.4, Total = 31.9
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> (Intercept) 0.540 0.708 -0.853 0.541 1.930 0.541 0
#> Imunisasi -0.005 0.004 -0.013 -0.005 0.003 -0.005 0
#> ASI -0.004 0.006 -0.017 -0.004 0.009 -0.004 0
#> Sanitasi -0.008 0.004 -0.016 -0.008 0.000 -0.008 0
#> Miskin 0.056 0.020 0.017 0.056 0.095 0.056 0
#>
#> Random effects:
#> Name Model
#> ID IID model
#>
#> Model hyperparameters:
#> mean
#> size for the nbinomial observations (1/overdispersion) 1973.16
#> Precision for ID 1.35
#> sd
#> size for the nbinomial observations (1/overdispersion) 2.04e+04
#> Precision for ID 2.05e-01
#> 0.025quant
#> size for the nbinomial observations (1/overdispersion) 5.98
#> Precision for ID 0.99
#> 0.5quant
#> size for the nbinomial observations (1/overdispersion) 144.59
#> Precision for ID 1.34
#> 0.975quant
#> size for the nbinomial observations (1/overdispersion) 12474.58
#> Precision for ID 1.79
#> mode
#> size for the nbinomial observations (1/overdispersion) 11.33
#> Precision for ID 1.32
#>
#> Deviance Information Criterion (DIC) ...............: 1003.91
#> Deviance Information Criterion (DIC, saturated) ....: 253.96
#> Effective number of parameters .....................: 129.50
#>
#> Marginal log-Likelihood: -625.83
#> is computed
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Komponen acak yang digunakan adalah Independent Gaussian Random Effect. Komponennya adalah sebagai berikut:
rand_glmm <- cbind(jawa@data$Provinsi,jawa@data$KabKota, glmm_mod$summary.random$ID[,c('ID', 'mean')])
names(rand_glmm) <- c("Provinsi","KabKota", "ID", "iid")
rand_glmm
#> Provinsi KabKota ID iid
#> 1 DKI Jakarta Kepulauan Seribu 1 -1.92061
#> 2 Jawa Barat Bandung 2 1.58686
#> 3 Jawa Barat Bandung Barat 3 0.99858
#> 4 Jawa Timur Bangkalan 4 -1.31739
#> 5 Jawa Tengah Banjarnegara 5 0.71354
#> 6 DI Yogyakarta Bantul 6 -0.18243
#> 7 Jawa Tengah Banyumas 7 0.92673
#> 8 Jawa Timur Banyuwangi 8 0.05498
#> 9 Jawa Tengah Batang 9 0.18985
#> 10 Jawa Barat Bekasi 10 0.19760
#> 11 Jawa Timur Blitar 11 0.40268
#> 12 Jawa Tengah Blora 12 -0.03287
#> 13 Jawa Barat Bogor 13 2.24594
#> 14 Jawa Timur Bojonegoro 14 -0.56576
#> 15 Jawa Timur Bondowoso 15 -0.12135
#> 16 Jawa Tengah Boyolali 16 0.22575
#> 17 Jawa Tengah Brebes 17 0.74695
#> 18 Jawa Barat Ciamis 18 -0.31841
#> 19 Jawa Barat Cianjur 19 0.86343
#> 20 Jawa Tengah Cilacap 20 0.09518
#> 21 Jawa Barat Cirebon 21 1.00910
#> 22 Jawa Tengah Demak 22 -0.28924
#> 23 Jawa Barat Garut 23 2.07900
#> 24 Jawa Timur Gresik 24 0.74706
#> 25 Jawa Tengah Grobogan 25 0.33579
#> 26 DI Yogyakarta Gunungkidul 26 0.00408
#> 27 Jawa Barat Indramayu 27 0.06413
#> 28 Jawa Timur Jember 28 0.84319
#> 29 Jawa Tengah Jepara 29 0.80847
#> 30 Jawa Timur Jombang 30 0.50972
#> 31 Jawa Tengah Karanganyar 31 -0.64370
#> 32 Jawa Barat Karawang 32 -0.20285
#> 33 Jawa Tengah Kebumen 33 0.39124
#> 34 Jawa Timur Kediri 34 0.82172
#> 35 Jawa Tengah Kendal 35 0.37053
#> 36 Jawa Tengah Klaten 36 0.77489
#> 37 DKI Jakarta Jakarta Barat 37 0.85209
#> 38 DKI Jakarta Jakarta Pusat 38 -0.39096
#> 39 DKI Jakarta Jakarta Selatan 39 -1.17365
#> 40 DKI Jakarta Jakarta Timur 40 -0.17625
#> 41 DKI Jakarta Jakarta Utara 41 -0.40299
#> 42 Jawa Barat Kota Bandung 42 0.86226
#> 43 Jawa Barat Kota Banjar 43 -1.54385
#> 44 Jawa Timur Kota Batu 44 -0.31482
#> 45 Jawa Barat Kota Bekasi 45 0.47286
#> 46 Jawa Timur Kota Blitar 46 -1.48461
#> 47 Jawa Barat Kota Bogor 47 -0.05834
#> 48 Banten Kota Cilegon 48 -0.60054
#> 49 Jawa Barat Kota Cimahi 49 0.27902
#> 50 Jawa Barat Kota Cirebon 50 -0.39427
#> 51 Jawa Barat Kota Depok 51 -0.00110
#> 52 Jawa Timur Kota Kediri 52 -0.94249
#> 53 Jawa Timur Kota Madiun 53 -1.26011
#> 54 Jawa Tengah Kota Magelang 54 -1.58740
#> 55 Jawa Timur Kota Malang 55 0.30368
#> 56 Jawa Timur Kota Mojokerto 56 -1.89639
#> 57 Jawa Timur Kota Pasuruan 57 -0.24426
#> 58 Jawa Tengah Kota Pekalongan 58 -0.69005
#> 59 Jawa Timur Kota Probolinggo 59 -0.53377
#> 60 Jawa Tengah Kota Salatiga 60 -1.05816
#> 61 Jawa Tengah Kota Semarang 61 -0.42559
#> 62 Banten Kota Serang 62 0.34668
#> 63 Jawa Barat Kota Sukabumi 63 -0.71138
#> 64 Jawa Timur Kota Surabaya 64 -0.16827
#> 65 Jawa Tengah Kota Surakarta 65 -0.89246
#> 66 Banten Kota Tangerang 66 0.38752
#> 67 Banten Kota Tangerang Selatan 67 -1.05295
#> 68 Jawa Barat Kota Tasikmalaya 68 0.47174
#> 69 Jawa Tengah Kota Tegal 69 -0.97056
#> 70 DI Yogyakarta Kota Yogyakarta 70 -0.64684
#> 71 Jawa Tengah Kudus 71 0.25512
#> 72 DI Yogyakarta Kulon Progo 72 -0.73537
#> 73 Jawa Barat Kuningan 73 0.26467
#> 74 Jawa Timur Lamongan 74 -0.12476
#> 75 Banten Lebak 75 -0.23336
#> 76 Jawa Timur Lumajang 76 0.23570
#> 77 Jawa Timur Madiun 77 0.23037
#> 78 Jawa Tengah Magelang 78 1.05132
#> 79 Jawa Timur Magetan 79 -0.02828
#> 80 Jawa Barat Majalengka 80 -0.35300
#> 81 Jawa Timur Malang 81 1.20800
#> 82 Jawa Timur Mojokerto 82 -0.64610
#> 83 Jawa Timur Nganjuk 83 0.11795
#> 84 Jawa Timur Ngawi 84 0.03229
#> 85 Jawa Timur Pacitan 85 -0.24726
#> 86 Jawa Timur Pamekasan 86 -0.38758
#> 87 Banten Pandeglang 87 -0.71738
#> 88 Jawa Barat Pangandaran 88 -1.54383
#> 89 Jawa Timur Pasuruan 89 -0.65686
#> 90 Jawa Tengah Pati 90 0.36788
#> 91 Jawa Tengah Pekalongan 91 -0.78304
#> 92 Jawa Tengah Pemalang 92 0.40435
#> 93 Jawa Timur Ponorogo 93 0.37056
#> 94 Jawa Timur Probolinggo 94 0.54299
#> 95 Jawa Tengah Purbalingga 95 0.56835
#> 96 Jawa Barat Purwakarta 96 -0.30752
#> 97 Jawa Tengah Purworejo 97 0.02227
#> 98 Jawa Tengah Rembang 98 0.02017
#> 99 Jawa Timur Sampang 99 -0.96164
#> 100 Jawa Tengah Semarang 100 0.02064
#> 101 Banten Serang 101 0.50251
#> 102 Jawa Timur Sidoarjo 102 0.63501
#> 103 Jawa Timur Situbondo 103 -0.51098
#> 104 DI Yogyakarta Sleman 104 0.25168
#> 105 Jawa Tengah Sragen 105 -0.05524
#> 106 Jawa Barat Subang 106 -0.51864
#> 107 Jawa Barat Sukabumi 107 1.37433
#> 108 Jawa Tengah Sukoharjo 108 0.23425
#> 109 Jawa Barat Sumedang 109 0.49675
#> 110 Jawa Timur Sumenep 110 -0.77524
#> 111 Banten Tangerang 111 0.32874
#> 112 Jawa Barat Tasikmalaya 112 0.93012
#> 113 Jawa Tengah Tegal 113 1.70899
#> 114 Jawa Tengah Temanggung 114 0.51006
#> 115 Jawa Timur Trenggalek 115 -0.18662
#> 116 Jawa Timur Tuban 116 0.14917
#> 117 Jawa Timur Tulungagung 117 -0.19299
#> 118 Jawa Tengah Wonogiri 118 0.14849
#> 119 Jawa Tengah Wonosobo 119 0.41827
view(rand_glmm)Residuals autocorrelation
moran.test(residuals(glmm_mod)[[1]], listw = WW1, na.action = na.omit, zero.policy = T)
#>
#> Moran I test under randomisation
#>
#> data: residuals(glmm_mod)[[1]]
#> weights: WW1
#>
#> Moran I statistic standard deviate = 0.5, p-value = 0.3
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.02245 -0.00847 0.00356ICAR
wm_q.mat <- as(nb2mat(wm_q, style = "B", zero.policy = T), "Matrix")
# Fit model
form_fit3 <- Stunting ~ Imunisasi + ASI + Sanitasi + Miskin +
f(ID, model = "besag", graph = wm_q.mat)
icar_mod <- inla(form_fit3,
data = jawa@data, E = E_d, family ="nbinomial",
control.predictor = list(compute = TRUE, link = 1),
control.compute = list(dic = TRUE))
summary(icar_mod)
#> Time used:
#> Pre = 0.239, Running = 0.268, Post = 30, Total = 30.5
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> (Intercept) 0.988 0.698 -0.393 0.991 2.351 0.991 0
#> Imunisasi -0.002 0.004 -0.010 -0.002 0.005 -0.002 0
#> ASI -0.005 0.006 -0.018 -0.005 0.007 -0.005 0
#> Sanitasi -0.009 0.004 -0.017 -0.009 -0.002 -0.009 0
#> Miskin 0.042 0.021 0.001 0.042 0.084 0.042 0
#>
#> Random effects:
#> Name Model
#> ID Besags ICAR model
#>
#> Model hyperparameters:
#> mean
#> size for the nbinomial observations (1/overdispersion) 1.46e+00
#> Precision for ID 5.74e+05
#> sd
#> size for the nbinomial observations (1/overdispersion) 2.08e-01
#> Precision for ID 4.86e+06
#> 0.025quant
#> size for the nbinomial observations (1/overdispersion) 1.08
#> Precision for ID 1992.98
#> 0.5quant
#> size for the nbinomial observations (1/overdispersion) 1.45
#> Precision for ID 60176.24
#> 0.975quant
#> size for the nbinomial observations (1/overdispersion) 1.90e+00
#> Precision for ID 3.84e+06
#> mode
#> size for the nbinomial observations (1/overdispersion) 1.43
#> Precision for ID 3335.22
#>
#> Deviance Information Criterion (DIC) ...............: 1189.92
#> Deviance Information Criterion (DIC, saturated) ....: -3122.74
#> Effective number of parameters .....................: 6.73
#>
#> Marginal log-Likelihood: -697.07
#> is computed
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Komponen acak spasial yang digunakan pada model ICAR \(u\) adalah sebagai berikut:
rand_icar <- cbind(jawa@data$Provinsi,jawa@data$KabKota, icar_mod$summary.random$ID[,c('ID', 'mean')])
names(rand_icar) <- c("Provinsi","KabKota", "ID", "u")
rand_icar
#> Provinsi KabKota ID u
#> 1 DKI Jakarta Kepulauan Seribu 1 -3.118377
#> 2 Jawa Barat Bandung 2 0.046111
#> 3 Jawa Barat Bandung Barat 3 0.043445
#> 4 Jawa Timur Bangkalan 4 -0.002839
#> 5 Jawa Tengah Banjarnegara 5 0.005548
#> 6 DI Yogyakarta Bantul 6 -0.012046
#> 7 Jawa Tengah Banyumas 7 0.009261
#> 8 Jawa Timur Banyuwangi 8 -0.013732
#> 9 Jawa Tengah Batang 9 -0.001868
#> 10 Jawa Barat Bekasi 10 0.024956
#> 11 Jawa Timur Blitar 11 -0.016111
#> 12 Jawa Tengah Blora 12 -0.014657
#> 13 Jawa Barat Bogor 13 0.042305
#> 14 Jawa Timur Bojonegoro 14 -0.018918
#> 15 Jawa Timur Bondowoso 15 -0.013710
#> 16 Jawa Tengah Boyolali 16 -0.012239
#> 17 Jawa Tengah Brebes 17 0.013023
#> 18 Jawa Barat Ciamis 18 0.012135
#> 19 Jawa Barat Cianjur 19 0.044356
#> 20 Jawa Tengah Cilacap 20 0.009606
#> 21 Jawa Barat Cirebon 21 0.017525
#> 22 Jawa Tengah Demak 22 -0.010668
#> 23 Jawa Barat Garut 23 0.040641
#> 24 Jawa Timur Gresik 24 -0.013499
#> 25 Jawa Tengah Grobogan 25 -0.012133
#> 26 DI Yogyakarta Gunungkidul 26 -0.012004
#> 27 Jawa Barat Indramayu 27 0.020742
#> 28 Jawa Timur Jember 28 -0.014058
#> 29 Jawa Tengah Jepara 29 -0.013726
#> 30 Jawa Timur Jombang 30 -0.014916
#> 31 Jawa Tengah Karanganyar 31 -0.016995
#> 32 Jawa Barat Karawang 32 0.027443
#> 33 Jawa Tengah Kebumen 33 0.005510
#> 34 Jawa Timur Kediri 34 -0.014172
#> 35 Jawa Tengah Kendal 35 -0.008357
#> 36 Jawa Tengah Klaten 36 -0.011183
#> 37 DKI Jakarta Jakarta Barat 37 0.021133
#> 38 DKI Jakarta Jakarta Pusat 38 0.016619
#> 39 DKI Jakarta Jakarta Selatan 39 0.016794
#> 40 DKI Jakarta Jakarta Timur 40 0.018262
#> 41 DKI Jakarta Jakarta Utara 41 0.017579
#> 42 Jawa Barat Kota Bandung 42 0.037940
#> 43 Jawa Barat Kota Banjar 43 0.010670
#> 44 Jawa Timur Kota Batu 44 -0.016533
#> 45 Jawa Barat Kota Bekasi 45 0.025222
#> 46 Jawa Timur Kota Blitar 46 -0.023358
#> 47 Jawa Barat Kota Bogor 47 0.020434
#> 48 Banten Kota Cilegon 48 0.018558
#> 49 Jawa Barat Kota Cimahi 49 0.035170
#> 50 Jawa Barat Kota Cirebon 50 0.017615
#> 51 Jawa Barat Kota Depok 51 0.021677
#> 52 Jawa Timur Kota Kediri 52 -0.019472
#> 53 Jawa Timur Kota Madiun 53 -0.019470
#> 54 Jawa Tengah Kota Magelang 54 -0.021931
#> 55 Jawa Timur Kota Malang 55 -0.016615
#> 56 Jawa Timur Kota Mojokerto 56 -0.026302
#> 57 Jawa Timur Kota Pasuruan 57 -0.017077
#> 58 Jawa Tengah Kota Pekalongan 58 -0.005129
#> 59 Jawa Timur Kota Probolinggo 59 -0.016880
#> 60 Jawa Tengah Kota Salatiga 60 -0.019160
#> 61 Jawa Tengah Kota Semarang 61 -0.011137
#> 62 Banten Kota Serang 62 0.017764
#> 63 Jawa Barat Kota Sukabumi 63 0.021266
#> 64 Jawa Timur Kota Surabaya 64 -0.014904
#> 65 Jawa Tengah Kota Surakarta 65 -0.015362
#> 66 Banten Kota Tangerang 66 0.019273
#> 67 Banten Kota Tangerang Selatan 67 0.019199
#> 68 Jawa Barat Kota Tasikmalaya 68 0.015573
#> 69 Jawa Tengah Kota Tegal 69 0.014959
#> 70 DI Yogyakarta Kota Yogyakarta 70 -0.014405
#> 71 Jawa Tengah Kudus 71 -0.011546
#> 72 DI Yogyakarta Kulon Progo 72 -0.010793
#> 73 Jawa Barat Kuningan 73 0.012958
#> 74 Jawa Timur Lamongan 74 -0.016837
#> 75 Banten Lebak 75 0.023481
#> 76 Jawa Timur Lumajang 76 -0.014091
#> 77 Jawa Timur Madiun 77 -0.018714
#> 78 Jawa Tengah Magelang 78 -0.008455
#> 79 Jawa Timur Magetan 79 -0.018646
#> 80 Jawa Barat Majalengka 80 0.017791
#> 81 Jawa Timur Malang 81 -0.012918
#> 82 Jawa Timur Mojokerto 82 -0.020649
#> 83 Jawa Timur Nganjuk 83 -0.017408
#> 84 Jawa Timur Ngawi 84 -0.017222
#> 85 Jawa Timur Pacitan 85 -0.016466
#> 86 Jawa Timur Pamekasan 86 0.001917
#> 87 Banten Pandeglang 87 0.018359
#> 88 Jawa Barat Pangandaran 88 0.019974
#> 89 Jawa Timur Pasuruan 89 -0.018175
#> 90 Jawa Tengah Pati 90 -0.011457
#> 91 Jawa Tengah Pekalongan 91 0.002149
#> 92 Jawa Tengah Pemalang 92 0.009390
#> 93 Jawa Timur Ponorogo 93 -0.018003
#> 94 Jawa Timur Probolinggo 94 -0.013404
#> 95 Jawa Tengah Purbalingga 95 0.008448
#> 96 Jawa Barat Purwakarta 96 0.034539
#> 97 Jawa Tengah Purworejo 97 -0.007723
#> 98 Jawa Tengah Rembang 98 -0.013231
#> 99 Jawa Timur Sampang 99 -0.000746
#> 100 Jawa Tengah Semarang 100 -0.009623
#> 101 Banten Serang 101 0.017486
#> 102 Jawa Timur Sidoarjo 102 -0.013735
#> 103 Jawa Timur Situbondo 103 -0.014709
#> 104 DI Yogyakarta Sleman 104 -0.010300
#> 105 Jawa Tengah Sragen 105 -0.013590
#> 106 Jawa Barat Subang 106 0.033298
#> 107 Jawa Barat Sukabumi 107 0.037783
#> 108 Jawa Tengah Sukoharjo 108 -0.013020
#> 109 Jawa Barat Sumedang 109 0.033975
#> 110 Jawa Timur Sumenep 110 0.001665
#> 111 Banten Tangerang 111 0.023613
#> 112 Jawa Barat Tasikmalaya 112 0.025878
#> 113 Jawa Tengah Tegal 113 0.015364
#> 114 Jawa Tengah Temanggung 114 -0.008955
#> 115 Jawa Timur Trenggalek 115 -0.016945
#> 116 Jawa Timur Tuban 116 -0.014537
#> 117 Jawa Timur Tulungagung 117 -0.017670
#> 118 Jawa Tengah Wonogiri 118 -0.016048
#> 119 Jawa Tengah Wonosobo 119 -0.003386
view(rand_icar)Residuals autocorrelation
moran.test(residuals(icar_mod)[[1]], listw = WW1, na.action = na.omit, zero.policy = T)
#>
#> Moran I test under randomisation
#>
#> data: residuals(icar_mod)[[1]]
#> weights: WW1
#>
#> Moran I statistic standard deviate = 1, p-value = 0.1
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.05999 -0.00847 0.00348BYM
form_fit4 <- Stunting ~ ASI + Imunisasi + Sanitasi + Miskin +
f(ID, model = "bym", graph = wm_q.mat)
bym_mod <- inla(form_fit4,
data = jawa@data, E = E_d, family ="nbinomial",
control.predictor = list(compute = TRUE, link = 1),
control.compute = list(dic = TRUE))
summary(bym_mod)
#> Time used:
#> Pre = 0.878, Running = 0.523, Post = 34.4, Total = 35.8
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> (Intercept) 0.382 0.692 -0.979 0.382 1.742 0.382 0
#> ASI -0.004 0.006 -0.017 -0.004 0.008 -0.004 0
#> Imunisasi -0.003 0.004 -0.011 -0.003 0.005 -0.003 0
#> Sanitasi -0.009 0.004 -0.017 -0.009 -0.001 -0.009 0
#> Miskin 0.062 0.019 0.024 0.062 0.100 0.062 0
#>
#> Random effects:
#> Name Model
#> ID BYM model
#>
#> Model hyperparameters:
#> mean
#> size for the nbinomial observations (1/overdispersion) 1651.85
#> Precision for ID (iid component) 1.47
#> Precision for ID (spatial component) 2115.72
#> sd
#> size for the nbinomial observations (1/overdispersion) 1.60e+04
#> Precision for ID (iid component) 2.24e-01
#> Precision for ID (spatial component) 2.39e+03
#> 0.025quant
#> size for the nbinomial observations (1/overdispersion) 5.99
#> Precision for ID (iid component) 1.06
#> Precision for ID (spatial component) 126.01
#> 0.5quant
#> size for the nbinomial observations (1/overdispersion) 134.71
#> Precision for ID (iid component) 1.45
#> Precision for ID (spatial component) 1358.25
#> 0.975quant
#> size for the nbinomial observations (1/overdispersion) 10621.81
#> Precision for ID (iid component) 1.95
#> Precision for ID (spatial component) 8446.44
#> mode
#> size for the nbinomial observations (1/overdispersion) 11.63
#> Precision for ID (iid component) 1.43
#> Precision for ID (spatial component) 335.78
#>
#> Deviance Information Criterion (DIC) ...............: 986.94
#> Deviance Information Criterion (DIC, saturated) ....: 247.20
#> Effective number of parameters .....................: 124.01
#>
#> Marginal log-Likelihood: -584.28
#> is computed
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')Komponen acak spasial yang digunakan pada model ICAR \(u\) adalah sebagai berikut:
rand_bym <- cbind(jawa@data$Provinsi,jawa@data$KabKota, bym_mod$summary.random$ID[1:119, c('ID', 'mean')], bym_mod$summary.random$ID[120:238, c('mean')])
names(rand_bym) <- c("Provinsi","KabKota", "ID", "u", "v")
rand_bym
#> Provinsi KabKota ID u v
#> 1 DKI Jakarta Kepulauan Seribu 1 -2.86740 -2.867204
#> 2 Jawa Barat Bandung 2 1.56561 0.012690
#> 3 Jawa Barat Bandung Barat 3 0.98639 0.012841
#> 4 Jawa Timur Bangkalan 4 -1.38053 -0.001172
#> 5 Jawa Tengah Banjarnegara 5 0.65107 0.001118
#> 6 DI Yogyakarta Bantul 6 -0.23511 -0.005563
#> 7 Jawa Tengah Banyumas 7 0.91692 0.003352
#> 8 Jawa Timur Banyuwangi 8 0.04782 -0.008689
#> 9 Jawa Tengah Batang 9 0.15991 -0.000519
#> 10 Jawa Barat Bekasi 10 0.24293 0.011898
#> 11 Jawa Timur Blitar 11 0.36698 -0.009416
#> 12 Jawa Tengah Blora 12 -0.09052 -0.006879
#> 13 Jawa Barat Bogor 13 2.27172 0.012393
#> 14 Jawa Timur Bojonegoro 14 -0.55080 -0.007962
#> 15 Jawa Timur Bondowoso 15 -0.18632 -0.008803
#> 16 Jawa Tengah Boyolali 16 0.18830 -0.005782
#> 17 Jawa Tengah Brebes 17 0.71340 0.005249
#> 18 Jawa Barat Ciamis 18 -0.41788 0.006199
#> 19 Jawa Barat Cianjur 19 0.90190 0.012726
#> 20 Jawa Tengah Cilacap 20 0.11894 0.004508
#> 21 Jawa Barat Cirebon 21 0.99172 0.007392
#> 22 Jawa Tengah Demak 22 -0.32442 -0.005222
#> 23 Jawa Barat Garut 23 2.06530 0.012058
#> 24 Jawa Timur Gresik 24 0.71635 -0.008563
#> 25 Jawa Tengah Grobogan 25 0.31125 -0.005868
#> 26 DI Yogyakarta Gunungkidul 26 -0.07076 -0.005872
#> 27 Jawa Barat Indramayu 27 0.00804 0.009212
#> 28 Jawa Timur Jember 28 0.91922 -0.008250
#> 29 Jawa Tengah Jepara 29 0.84859 -0.004963
#> 30 Jawa Timur Jombang 30 0.49922 -0.008454
#> 31 Jawa Tengah Karanganyar 31 -0.67046 -0.006941
#> 32 Jawa Barat Karawang 32 -0.22692 0.011771
#> 33 Jawa Tengah Kebumen 33 0.32833 0.001165
#> 34 Jawa Timur Kediri 34 0.80931 -0.008744
#> 35 Jawa Tengah Kendal 35 0.43495 -0.002672
#> 36 Jawa Tengah Klaten 36 0.75319 -0.005169
#> 37 DKI Jakarta Jakarta Barat 37 0.86261 0.011753
#> 38 DKI Jakarta Jakarta Pusat 38 -0.39493 0.011228
#> 39 DKI Jakarta Jakarta Selatan 39 -1.22001 0.011119
#> 40 DKI Jakarta Jakarta Timur 40 -0.19032 0.011492
#> 41 DKI Jakarta Jakarta Utara 41 -0.41376 0.011448
#> 42 Jawa Barat Kota Bandung 42 0.88428 0.013372
#> 43 Jawa Barat Kota Banjar 43 -1.62291 0.003777
#> 44 Jawa Timur Kota Batu 44 -0.29799 -0.009368
#> 45 Jawa Barat Kota Bekasi 45 0.47862 0.012047
#> 46 Jawa Timur Kota Blitar 46 -1.55538 -0.012319
#> 47 Jawa Barat Kota Bogor 47 -0.07520 0.012049
#> 48 Banten Kota Cilegon 48 -0.61314 0.010296
#> 49 Jawa Barat Kota Cimahi 49 0.29941 0.013090
#> 50 Jawa Barat Kota Cirebon 50 -0.41743 0.006498
#> 51 Jawa Barat Kota Depok 51 0.02535 0.011659
#> 52 Jawa Timur Kota Kediri 52 -0.97194 -0.010542
#> 53 Jawa Timur Kota Madiun 53 -1.28079 -0.008994
#> 54 Jawa Tengah Kota Magelang 54 -1.60634 -0.007198
#> 55 Jawa Timur Kota Malang 55 0.35170 -0.007977
#> 56 Jawa Timur Kota Mojokerto 56 -1.86718 -0.013005
#> 57 Jawa Timur Kota Pasuruan 57 -0.20937 -0.009737
#> 58 Jawa Tengah Kota Pekalongan 58 -0.73355 -0.000612
#> 59 Jawa Timur Kota Probolinggo 59 -0.53714 -0.009802
#> 60 Jawa Tengah Kota Salatiga 60 -1.09778 -0.006821
#> 61 Jawa Tengah Kota Semarang 61 -0.42236 -0.004486
#> 62 Banten Kota Serang 62 0.41396 0.012207
#> 63 Jawa Barat Kota Sukabumi 63 -0.70396 0.011129
#> 64 Jawa Timur Kota Surabaya 64 -0.17561 -0.008823
#> 65 Jawa Tengah Kota Surakarta 65 -0.97069 -0.006925
#> 66 Banten Kota Tangerang 66 0.40536 0.011655
#> 67 Banten Kota Tangerang Selatan 67 -0.95515 0.011336
#> 68 Jawa Barat Kota Tasikmalaya 68 0.46733 0.007771
#> 69 Jawa Tengah Kota Tegal 69 -0.96972 0.003999
#> 70 DI Yogyakarta Kota Yogyakarta 70 -0.67398 -0.006015
#> 71 Jawa Tengah Kudus 71 0.19470 -0.005389
#> 72 DI Yogyakarta Kulon Progo 72 -0.82383 -0.004696
#> 73 Jawa Barat Kuningan 73 0.27451 0.006328
#> 74 Jawa Timur Lamongan 74 -0.17800 -0.008480
#> 75 Banten Lebak 75 -0.22161 0.011789
#> 76 Jawa Timur Lumajang 76 0.21520 -0.008438
#> 77 Jawa Timur Madiun 77 0.18357 -0.007930
#> 78 Jawa Tengah Magelang 78 1.02592 -0.004177
#> 79 Jawa Timur Magetan 79 -0.04509 -0.007676
#> 80 Jawa Barat Majalengka 80 -0.40680 0.007882
#> 81 Jawa Timur Malang 81 1.20668 -0.008676
#> 82 Jawa Timur Mojokerto 82 -0.70086 -0.009528
#> 83 Jawa Timur Nganjuk 83 0.08571 -0.008240
#> 84 Jawa Timur Ngawi 84 -0.03795 -0.007123
#> 85 Jawa Timur Pacitan 85 -0.28264 -0.007908
#> 86 Jawa Timur Pamekasan 86 -0.41430 0.000705
#> 87 Banten Pandeglang 87 -0.60675 0.011042
#> 88 Jawa Barat Pangandaran 88 -1.66738 0.005354
#> 89 Jawa Timur Pasuruan 89 -0.68818 -0.009381
#> 90 Jawa Tengah Pati 90 0.31988 -0.005860
#> 91 Jawa Tengah Pekalongan 91 -0.78127 0.000718
#> 92 Jawa Tengah Pemalang 92 0.35073 0.002887
#> 93 Jawa Timur Ponorogo 93 0.34109 -0.007872
#> 94 Jawa Timur Probolinggo 94 0.47410 -0.008806
#> 95 Jawa Tengah Purbalingga 95 0.53490 0.002248
#> 96 Jawa Barat Purwakarta 96 -0.35043 0.012045
#> 97 Jawa Tengah Purworejo 97 0.03056 -0.002269
#> 98 Jawa Tengah Rembang 98 -0.03864 -0.006766
#> 99 Jawa Timur Sampang 99 -1.02696 -0.000329
#> 100 Jawa Tengah Semarang 100 0.03919 -0.004753
#> 101 Banten Serang 101 0.47487 0.011599
#> 102 Jawa Timur Sidoarjo 102 0.63767 -0.008770
#> 103 Jawa Timur Situbondo 103 -0.51520 -0.009091
#> 104 DI Yogyakarta Sleman 104 0.23969 -0.005183
#> 105 Jawa Tengah Sragen 105 -0.08010 -0.006463
#> 106 Jawa Barat Subang 106 -0.55312 0.011306
#> 107 Jawa Barat Sukabumi 107 1.39501 0.012642
#> 108 Jawa Tengah Sukoharjo 108 0.22632 -0.006225
#> 109 Jawa Barat Sumedang 109 0.44168 0.010418
#> 110 Jawa Timur Sumenep 110 -0.87545 0.000795
#> 111 Banten Tangerang 111 0.29245 0.011787
#> 112 Jawa Barat Tasikmalaya 112 1.07078 0.008605
#> 113 Jawa Tengah Tegal 113 1.73440 0.004663
#> 114 Jawa Tengah Temanggung 114 0.53561 -0.003005
#> 115 Jawa Timur Trenggalek 115 -0.24421 -0.008272
#> 116 Jawa Timur Tuban 116 0.09672 -0.007478
#> 117 Jawa Timur Tulungagung 117 -0.20143 -0.008593
#> 118 Jawa Tengah Wonogiri 118 0.13474 -0.007048
#> 119 Jawa Tengah Wonosobo 119 0.34513 -0.001394
view(rand_bym)Residuals autocorrelation
moran.test(residuals(bym_mod)[[1]], listw = WW1, na.action = na.omit, zero.policy = T)
#>
#> Moran I test under randomisation
#>
#> data: residuals(bym_mod)[[1]]
#> weights: WW1
#>
#> Moran I statistic standard deviate = 0.6, p-value = 0.3
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.02819 -0.00847 0.00356Model Selection
mad <- function(y, yhat){
return(mean(abs(y-yhat)))
}
mad_glm <- mad(jawa@data$Stunting, glm_mod$summary.fitted.values$mean*jawa@data$E_d)
mad_glmm <- mad(jawa@data$Stunting, glmm_mod$summary.fitted.values$mean*jawa@data$E_d)
mad_icar <- mad(jawa@data$Stunting, icar_mod$summary.fitted.values$mean*jawa@data$E_d)
mad_bym <- mad(jawa@data$Stunting, bym_mod$summary.fitted.values$mean*jawa@data$E_d)
mad_glm
#> [1] 35.7
mad_glmm
#> [1] 3.34
mad_icar
#> [1] 35.1
mad_bym
#> [1] 2.7| Model | DIC | MAD |
|---|---|---|
| GLM | 1196.11 | 35.7369 |
| GLMM | 1003.91 | 3.3447 |
| ICAR | 1189.92 | 35.114 |
| BYM | 986.94 | 2.7039 |
Diperoleh nilai DIC dan MAD model ICAR lebih kecil daripada model BYM.
mape <- function(y, yhat){
return(sum(abs(y-yhat)/y))
}
y <- jawa@data$Stunting/jawa@data$E_d
mape_glm <- mape(y, glm_mod$summary.fitted.values$mean)
mape_glmm <- mape(y, glmm_mod$summary.fitted.values$mean)
mape_icar <- mape(y, icar_mod$summary.fitted.values$mean)
mape_bym <- mape(y, bym_mod$summary.fitted.values$mean)
mape_glm
#> [1] 214
mape_glmm
#> [1] 18.3
mape_icar
#> [1] 192
mape_bym
#> [1] 14.7Relative Risk
jawa@data$RR <- bym_mod$summary.fitted.values$mean
bym_mod$summary.fitted.values
#> mean sd 0.025quant 0.5quant 0.975quant
#> fitted.Predictor.001 0.0638 0.0428 0.0161 0.0530 0.175
#> fitted.Predictor.002 3.7502 0.6819 2.2912 3.7950 5.098
#> fitted.Predictor.003 1.7946 0.3486 1.1119 1.7884 2.536
#> fitted.Predictor.004 0.4388 0.1437 0.2482 0.4120 0.797
#> fitted.Predictor.005 2.2652 0.4340 1.4343 2.2498 3.224
#> fitted.Predictor.006 0.6324 0.1582 0.3828 0.6125 1.004
#> fitted.Predictor.007 2.1613 0.4117 1.3506 2.1554 3.043
#> fitted.Predictor.008 0.5611 0.1405 0.3337 0.5450 0.886
#> fitted.Predictor.009 0.9470 0.2112 0.5891 0.9266 1.432
#> fitted.Predictor.010 0.7820 0.1798 0.4785 0.7645 1.193
#> fitted.Predictor.011 1.0599 0.2291 0.6591 1.0410 1.579
#> fitted.Predictor.012 0.6594 0.1611 0.4001 0.6404 1.035
#> fitted.Predictor.013 5.9045 1.0799 3.4899 6.0608 7.846
#> fitted.Predictor.014 0.4677 0.1301 0.2724 0.4488 0.778
#> fitted.Predictor.015 0.7924 0.1885 0.4898 0.7698 1.236
#> fitted.Predictor.016 0.7844 0.1806 0.4806 0.7665 1.198
#> fitted.Predictor.017 2.2573 0.4316 1.4250 2.2442 3.204
#> fitted.Predictor.018 0.4962 0.1339 0.2919 0.4776 0.814
#> fitted.Predictor.019 1.3235 0.2685 0.8146 1.3133 1.902
#> fitted.Predictor.020 0.9149 0.2063 0.5680 0.8942 1.391
#> fitted.Predictor.021 2.7166 0.5056 1.7006 2.7170 3.791
#> fitted.Predictor.022 0.6016 0.1535 0.3626 0.5815 0.964
#> fitted.Predictor.023 4.7193 0.8624 2.8132 4.8209 6.315
#> fitted.Predictor.024 1.4418 0.2905 0.8968 1.4286 2.077
#> fitted.Predictor.025 1.2258 0.2596 0.7708 1.2043 1.818
#> fitted.Predictor.026 0.9465 0.2156 0.5916 0.9227 1.451
#> fitted.Predictor.027 0.7799 0.1827 0.4797 0.7598 1.204
#> fitted.Predictor.028 2.4345 0.4592 1.5243 2.4302 3.418
#> fitted.Predictor.029 1.3955 0.2812 0.8624 1.3847 2.003
#> fitted.Predictor.030 1.0434 0.2234 0.6449 1.0273 1.542
#> fitted.Predictor.031 0.3624 0.1083 0.2029 0.3458 0.621
#> fitted.Predictor.032 0.5781 0.1480 0.3469 0.5589 0.927
#> fitted.Predictor.033 1.4989 0.3075 0.9488 1.4762 2.199
#> fitted.Predictor.034 1.8638 0.3621 1.1654 1.8534 2.648
#> fitted.Predictor.035 1.2776 0.2667 0.7997 1.2586 1.878
#> fitted.Predictor.036 1.6113 0.3194 1.0048 1.5990 2.307
#> fitted.Predictor.037 1.1033 0.2307 0.6743 1.0920 1.604
#> fitted.Predictor.038 0.3498 0.1028 0.1950 0.3349 0.593
#> fitted.Predictor.039 0.1340 0.0565 0.0582 0.1232 0.273
#> fitted.Predictor.040 0.4311 0.1180 0.2482 0.4154 0.708
#> fitted.Predictor.041 0.3875 0.1109 0.2198 0.3716 0.650
#> fitted.Predictor.042 1.2737 0.2603 0.7829 1.2630 1.836
#> fitted.Predictor.043 0.2363 0.0910 0.1166 0.2189 0.460
#> fitted.Predictor.044 0.3259 0.0967 0.1794 0.3120 0.554
#> fitted.Predictor.045 0.9066 0.1996 0.5564 0.8907 1.354
#> fitted.Predictor.046 0.1264 0.0562 0.0530 0.1152 0.265
#> fitted.Predictor.047 0.4974 0.1298 0.2920 0.4812 0.800
#> fitted.Predictor.048 0.2672 0.0866 0.1409 0.2535 0.475
#> fitted.Predictor.049 0.6014 0.1459 0.3586 0.5864 0.934
#> fitted.Predictor.050 0.5165 0.1386 0.3051 0.4972 0.846
#> fitted.Predictor.051 0.7431 0.1756 0.4552 0.7237 1.150
#> fitted.Predictor.052 0.2254 0.0795 0.1135 0.2118 0.418
#> fitted.Predictor.053 0.1356 0.0572 0.0590 0.1246 0.276
#> fitted.Predictor.054 0.1348 0.0598 0.0568 0.1229 0.283
#> fitted.Predictor.055 0.6853 0.1608 0.4131 0.6700 1.049
#> fitted.Predictor.056 0.0834 0.0440 0.0289 0.0739 0.194
#> fitted.Predictor.057 0.5036 0.1324 0.2966 0.4865 0.814
#> fitted.Predictor.058 0.2909 0.0931 0.1558 0.2760 0.514
#> fitted.Predictor.059 0.3371 0.1012 0.1865 0.3219 0.578
#> fitted.Predictor.060 0.1676 0.0651 0.0783 0.1557 0.326
#> fitted.Predictor.061 0.3135 0.0951 0.1713 0.2994 0.539
#> fitted.Predictor.062 0.3140 0.0932 0.1673 0.3021 0.530
#> fitted.Predictor.063 0.2900 0.0926 0.1555 0.2753 0.512
#> fitted.Predictor.064 0.4110 0.1132 0.2353 0.3960 0.676
#> fitted.Predictor.065 0.2064 0.0750 0.1014 0.1934 0.388
#> fitted.Predictor.066 0.8416 0.1888 0.5152 0.8253 1.267
#> fitted.Predictor.067 0.2210 0.0802 0.1101 0.2069 0.416
#> fitted.Predictor.068 1.1174 0.2370 0.6950 1.1001 1.649
#> fitted.Predictor.069 0.2450 0.0845 0.1258 0.2307 0.449
#> fitted.Predictor.070 0.2878 0.0916 0.1542 0.2734 0.507
#> fitted.Predictor.071 0.7131 0.1680 0.4326 0.6960 1.097
#> fitted.Predictor.072 0.4646 0.1339 0.2698 0.4436 0.788
#> fitted.Predictor.073 0.9564 0.2110 0.5934 0.9375 1.437
#> fitted.Predictor.074 0.6656 0.1638 0.4049 0.6455 1.050
#> fitted.Predictor.075 0.9060 0.2124 0.5663 0.8799 1.410
#> fitted.Predictor.076 0.7646 0.1766 0.4671 0.7473 1.168
#> fitted.Predictor.077 0.8023 0.1841 0.4922 0.7841 1.223
#> fitted.Predictor.078 1.9078 0.3675 1.1812 1.9039 2.685
#> fitted.Predictor.079 0.6025 0.1496 0.3617 0.5849 0.950
#> fitted.Predictor.080 0.6630 0.1684 0.4035 0.6402 1.065
#> fitted.Predictor.081 2.2282 0.4220 1.3691 2.2331 3.100
#> fitted.Predictor.082 0.3451 0.1048 0.1914 0.3288 0.596
#> fitted.Predictor.083 0.7734 0.1800 0.4746 0.7545 1.188
#> fitted.Predictor.084 0.9440 0.2146 0.5895 0.9206 1.445
#> fitted.Predictor.085 0.6353 0.1598 0.3846 0.6148 1.012
#> fitted.Predictor.086 0.5892 0.1529 0.3541 0.5684 0.952
#> fitted.Predictor.087 0.6588 0.1742 0.3991 0.6330 1.081
#> fitted.Predictor.088 0.1725 0.0726 0.0783 0.1582 0.351
#> fitted.Predictor.089 0.3632 0.1088 0.2033 0.3464 0.623
#> fitted.Predictor.090 0.8284 0.1871 0.5075 0.8116 1.252
#> fitted.Predictor.091 0.3873 0.1159 0.2187 0.3691 0.666
#> fitted.Predictor.092 1.4969 0.3067 0.9470 1.4746 2.194
#> fitted.Predictor.093 1.0033 0.2185 0.6224 0.9850 1.498
#> fitted.Predictor.094 2.3412 0.4522 1.4941 2.3186 3.362
#> fitted.Predictor.095 1.7903 0.3546 1.1322 1.7708 2.584
#> fitted.Predictor.096 0.4577 0.1249 0.2663 0.4405 0.753
#> fitted.Predictor.097 0.8135 0.1887 0.5017 0.7932 1.251
#> fitted.Predictor.098 0.9259 0.2110 0.5776 0.9030 1.418
#> fitted.Predictor.099 0.6681 0.1866 0.4051 0.6367 1.132
#> fitted.Predictor.100 0.6695 0.1614 0.4060 0.6515 1.043
#> fitted.Predictor.101 1.8720 0.3738 1.1842 1.8499 2.717
#> fitted.Predictor.102 0.9675 0.2087 0.5919 0.9536 1.429
#> fitted.Predictor.103 0.6136 0.1599 0.3706 0.5911 0.996
#> fitted.Predictor.104 0.7101 0.1663 0.4307 0.6935 1.089
#> fitted.Predictor.105 0.8211 0.1919 0.5083 0.7993 1.270
#> fitted.Predictor.106 0.3575 0.1064 0.1994 0.3415 0.611
#> fitted.Predictor.107 2.9093 0.5370 1.7866 2.9301 3.995
#> fitted.Predictor.108 0.7109 0.1665 0.4316 0.6942 1.091
#> fitted.Predictor.109 1.2507 0.2614 0.7817 1.2321 1.838
#> fitted.Predictor.110 0.6383 0.1748 0.3860 0.6104 1.068
#> fitted.Predictor.111 1.8238 0.3706 1.1610 1.7958 2.679
#> fitted.Predictor.112 2.7569 0.5147 1.7151 2.7607 3.842
#> fitted.Predictor.113 3.3201 0.6112 2.0061 3.3650 4.504
#> fitted.Predictor.114 1.1283 0.2381 0.6997 1.1123 1.658
#> fitted.Predictor.115 0.6147 0.1549 0.3709 0.5950 0.979
#> fitted.Predictor.116 1.1006 0.2410 0.6918 1.0769 1.659
#> fitted.Predictor.117 0.4668 0.1247 0.2723 0.4505 0.759
#> fitted.Predictor.118 0.9485 0.2121 0.5902 0.9277 1.437
#> fitted.Predictor.119 1.7261 0.3477 1.0973 1.7021 2.519
#> mode
#> fitted.Predictor.001 0.0369
#> fitted.Predictor.002 3.8879
#> fitted.Predictor.003 1.7923
#> fitted.Predictor.004 0.3783
#> fitted.Predictor.005 2.2476
#> fitted.Predictor.006 0.5853
#> fitted.Predictor.007 2.1638
#> fitted.Predictor.008 0.5214
#> fitted.Predictor.009 0.9019
#> fitted.Predictor.010 0.7414
#> fitted.Predictor.011 1.0193
#> fitted.Predictor.012 0.6145
#> fitted.Predictor.013 6.3341
#> fitted.Predictor.014 0.4218
#> fitted.Predictor.015 0.7411
#> fitted.Predictor.016 0.7429
#> fitted.Predictor.017 2.2446
#> fitted.Predictor.018 0.4510
#> fitted.Predictor.019 1.3061
#> fitted.Predictor.020 0.8689
#> fitted.Predictor.021 2.7402
#> fitted.Predictor.022 0.5538
#> fitted.Predictor.023 5.0068
#> fitted.Predictor.024 1.4189
#> fitted.Predictor.025 1.1813
#> fitted.Predictor.026 0.8942
#> fitted.Predictor.027 0.7337
#> fitted.Predictor.028 2.4442
#> fitted.Predictor.029 1.3777
#> fitted.Predictor.030 1.0090
#> fitted.Predictor.031 0.3209
#> fitted.Predictor.032 0.5320
#> fitted.Predictor.033 1.4556
#> fitted.Predictor.034 1.8526
#> fitted.Predictor.035 1.2398
#> fitted.Predictor.036 1.5925
#> fitted.Predictor.037 1.0808
#> fitted.Predictor.038 0.3117
#> fitted.Predictor.039 0.1058
#> fitted.Predictor.040 0.3916
#> fitted.Predictor.041 0.3474
#> fitted.Predictor.042 1.2545
#> fitted.Predictor.043 0.1942
#> fitted.Predictor.044 0.2900
#> fitted.Predictor.045 0.8710
#> fitted.Predictor.046 0.0977
#> fitted.Predictor.047 0.4572
#> fitted.Predictor.048 0.2317
#> fitted.Predictor.049 0.5646
#> fitted.Predictor.050 0.4700
#> fitted.Predictor.051 0.6980
#> fitted.Predictor.052 0.1903
#> fitted.Predictor.053 0.1071
#> fitted.Predictor.054 0.1044
#> fitted.Predictor.055 0.6487
#> fitted.Predictor.056 0.0589
#> fitted.Predictor.057 0.4615
#> fitted.Predictor.058 0.2529
#> fitted.Predictor.059 0.2985
#> fitted.Predictor.060 0.1367
#> fitted.Predictor.061 0.2771
#> fitted.Predictor.062 0.2823
#> fitted.Predictor.063 0.2524
#> fitted.Predictor.064 0.3729
#> fitted.Predictor.065 0.1726
#> fitted.Predictor.066 0.8043
#> fitted.Predictor.067 0.1850
#> fitted.Predictor.068 1.0812
#> fitted.Predictor.069 0.2083
#> fitted.Predictor.070 0.2507
#> fitted.Predictor.071 0.6727
#> fitted.Predictor.072 0.4146
#> fitted.Predictor.073 0.9147
#> fitted.Predictor.074 0.6183
#> fitted.Predictor.075 0.8488
#> fitted.Predictor.076 0.7242
#> fitted.Predictor.077 0.7605
#> fitted.Predictor.078 1.9122
#> fitted.Predictor.079 0.5601
#> fitted.Predictor.080 0.6103
#> fitted.Predictor.081 2.2566
#> fitted.Predictor.082 0.3043
#> fitted.Predictor.083 0.7297
#> fitted.Predictor.084 0.8927
#> fitted.Predictor.085 0.5871
#> fitted.Predictor.086 0.5400
#> fitted.Predictor.087 0.6005
#> fitted.Predictor.088 0.1370
#> fitted.Predictor.089 0.3214
#> fitted.Predictor.090 0.7898
#> fitted.Predictor.091 0.3427
#> fitted.Predictor.092 1.4546
#> fitted.Predictor.093 0.9636
#> fitted.Predictor.094 2.3090
#> fitted.Predictor.095 1.7572
#> fitted.Predictor.096 0.4152
#> fitted.Predictor.097 0.7673
#> fitted.Predictor.098 0.8752
#> fitted.Predictor.099 0.5991
#> fitted.Predictor.100 0.6267
#> fitted.Predictor.101 1.8354
#> fitted.Predictor.102 0.9369
#> fitted.Predictor.103 0.5612
#> fitted.Predictor.104 0.6709
#> fitted.Predictor.105 0.7718
#> fitted.Predictor.106 0.3172
#> fitted.Predictor.107 2.9830
#> fitted.Predictor.108 0.6714
#> fitted.Predictor.109 1.2128
#> fitted.Predictor.110 0.5756
#> fitted.Predictor.111 1.7742
#> fitted.Predictor.112 2.7891
#> fitted.Predictor.113 3.4567
#> fitted.Predictor.114 1.0950
#> fitted.Predictor.115 0.5680
#> fitted.Predictor.116 1.0501
#> fitted.Predictor.117 0.4263
#> fitted.Predictor.118 0.9027
#> fitted.Predictor.119 1.6831
summary(jawa@data$RR)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.06 0.43 0.71 1.00 1.18 5.90
view(jawa@data$RR)my.palette.1<-brewer.pal(n=2,name="OrRd")
#> Warning in brewer.pal(n = 2, name = "OrRd"): minimal value for n is 3, returning requested palette with 3 different levels
spplot(jawa,'RR', cuts=1, col.regions=my.palette.1, cex=c(0.3,1,0.3), main="Peta risiko relatif kasus stunting di Pulau Jawa",at = c(min(jawa$RR) - 1 ,1.0 , max(jawa$RR) + 1))
Dari grafik tersebut juga dapat dilihat bahwa risiko terjadinya stunting yang tinggi berada di Jawa Barat, dengan risiko paling tinggi di Kabupaten Bogor.
Relative Risk yang hanya menggunakan komponen acak adalah sebagai berikut.
Prediction
yhat_icar <- icar_mod$summary.fitted.values$mean*jawa@data$E_d
yhat_bym <- bym_mod$summary.fitted.values$mean*jawa@data$E_d
yhat_glmm <- glmm_mod$summary.fitted.values$mean*jawa@data$E_d
pred <- data.frame(jawa@data$Provinsi,jawa@data$KabKota, jawa@data$Stunting, yhat_icar, yhat_bym, yhat_glmm)
names(pred) <- c("Provinsi","KabKota", "Stunting", "yhat ICAR", "yhat BYM", "yhat GLMM")
pred
#> Provinsi KabKota Stunting yhat ICAR
#> 1 DKI Jakarta Kepulauan Seribu 3 4.89
#> 2 Jawa Barat Bandung 227 68.76
#> 3 Jawa Barat Bandung Barat 106 52.53
#> 4 Jawa Timur Bangkalan 20 109.84
#> 5 Jawa Tengah Banjarnegara 131 83.10
#> 6 DI Yogyakarta Bantul 34 55.66
#> 7 Jawa Tengah Banyumas 127 63.07
#> 8 Jawa Timur Banyuwangi 31 41.75
#> 9 Jawa Tengah Batang 53 64.11
#> 10 Jawa Barat Bekasi 44 50.13
#> 11 Jawa Timur Blitar 60 59.05
#> 12 Jawa Tengah Blora 36 52.30
#> 13 Jawa Barat Bogor 368 52.02
#> 14 Jawa Timur Bojonegoro 24 55.86
#> 15 Jawa Timur Bondowoso 43 66.45
#> 16 Jawa Tengah Boyolali 44 48.84
#> 17 Jawa Tengah Brebes 131 75.76
#> 18 Jawa Barat Ciamis 26 57.44
#> 19 Jawa Barat Cianjur 78 47.35
#> 20 Jawa Tengah Cilacap 51 60.41
#> 21 Jawa Barat Cirebon 160 76.96
#> 22 Jawa Tengah Demak 32 59.54
#> 23 Jawa Barat Garut 292 50.64
#> 24 Jawa Timur Gresik 84 50.94
#> 25 Jawa Tengah Grobogan 69 66.55
#> 26 DI Yogyakarta Gunungkidul 52 66.68
#> 27 Jawa Barat Indramayu 43 58.48
#> 28 Jawa Timur Jember 143 75.02
#> 29 Jawa Tengah Jepara 82 47.90
#> 30 Jawa Timur Jombang 60 47.81
#> 31 Jawa Tengah Karanganyar 18 52.03
#> 32 Jawa Barat Karawang 31 58.13
#> 33 Jawa Tengah Kebumen 85 72.31
#> 34 Jawa Timur Kediri 109 62.66
#> 35 Jawa Tengah Kendal 73 63.87
#> 36 Jawa Tengah Klaten 94 53.67
#> 37 DKI Jakarta Jakarta Barat 65 41.77
#> 38 DKI Jakarta Jakarta Pusat 18 45.03
#> 39 DKI Jakarta Jakarta Selatan 5 38.75
#> 40 DKI Jakarta Jakarta Timur 23 47.51
#> 41 DKI Jakarta Jakarta Utara 20 48.20
#> 42 Jawa Barat Kota Bandung 75 48.15
#> 43 Jawa Barat Kota Banjar 9 87.18
#> 44 Jawa Timur Kota Batu 17 36.96
#> 45 Jawa Barat Kota Bekasi 52 49.31
#> 46 Jawa Timur Kota Blitar 4 44.46
#> 47 Jawa Barat Kota Bogor 27 48.41
#> 48 Banten Kota Cilegon 13 43.44
#> 49 Jawa Barat Kota Cimahi 34 41.62
#> 50 Jawa Barat Kota Cirebon 27 67.97
#> 51 Jawa Barat Kota Depok 41 56.40
#> 52 Jawa Timur Kota Kediri 10 46.10
#> 53 Jawa Timur Kota Madiun 5 38.81
#> 54 Jawa Tengah Kota Magelang 4 50.94
#> 55 Jawa Timur Kota Malang 39 40.72
#> 56 Jawa Timur Kota Mojokerto 1 39.84
#> 57 Jawa Timur Kota Pasuruan 27 49.79
#> 58 Jawa Tengah Kota Pekalongan 14 49.23
#> 59 Jawa Timur Kota Probolinggo 17 45.31
#> 60 Jawa Tengah Kota Salatiga 7 40.57
#> 61 Jawa Tengah Kota Semarang 16 40.22
#> 62 Banten Kota Serang 18 18.99
#> 63 Jawa Barat Kota Sukabumi 14 47.67
#> 64 Jawa Timur Kota Surabaya 22 40.86
#> 65 Jawa Tengah Kota Surakarta 9 41.67
#> 66 Banten Kota Tangerang 48 48.76
#> 67 Banten Kota Tangerang Selatan 10 52.28
#> 68 Jawa Barat Kota Tasikmalaya 64 57.72
#> 69 Jawa Tengah Kota Tegal 11 50.75
#> 70 DI Yogyakarta Kota Yogyakarta 14 44.23
#> 71 Jawa Tengah Kudus 40 49.47
#> 72 DI Yogyakarta Kulon Progo 23 67.85
#> 73 Jawa Barat Kuningan 54 55.69
#> 74 Jawa Timur Lamongan 36 55.54
#> 75 Banten Lebak 49 91.56
#> 76 Jawa Timur Lumajang 43 46.50
#> 77 Jawa Timur Madiun 45 48.53
#> 78 Jawa Tengah Magelang 113 50.44
#> 79 Jawa Timur Magetan 33 46.34
#> 80 Jawa Barat Majalengka 35 84.04
#> 81 Jawa Timur Malang 132 50.41
#> 82 Jawa Timur Mojokerto 17 51.52
#> 83 Jawa Timur Nganjuk 43 51.54
#> 84 Jawa Timur Ngawi 52 68.63
#> 85 Jawa Timur Pacitan 34 56.64
#> 86 Jawa Timur Pamekasan 31 60.70
#> 87 Banten Pandeglang 34 97.11
#> 88 Jawa Barat Pangandaran 6 70.00
#> 89 Jawa Timur Pasuruan 18 54.72
#> 90 Jawa Tengah Pati 47 45.83
#> 91 Jawa Tengah Pekalongan 19 63.36
#> 92 Jawa Tengah Pemalang 85 73.86
#> 93 Jawa Timur Ponorogo 57 55.23
#> 94 Jawa Timur Probolinggo 134 97.22
#> 95 Jawa Tengah Purbalingga 102 72.93
#> 96 Jawa Barat Purwakarta 24 53.08
#> 97 Jawa Tengah Purworejo 45 56.72
#> 98 Jawa Tengah Rembang 51 65.25
#> 99 Jawa Timur Sampang 33 113.79
#> 100 Jawa Tengah Semarang 37 51.59
#> 101 Banten Serang 107 114.49
#> 102 Jawa Timur Sidoarjo 56 42.66
#> 103 Jawa Timur Situbondo 32 73.83
#> 104 DI Yogyakarta Sleman 40 43.34
#> 105 Jawa Tengah Sragen 45 62.70
#> 106 Jawa Barat Subang 18 48.76
#> 107 Jawa Barat Sukabumi 175 55.55
#> 108 Jawa Tengah Sukoharjo 40 44.38
#> 109 Jawa Barat Sumedang 71 60.82
#> 110 Jawa Timur Sumenep 32 101.37
#> 111 Banten Tangerang 103 131.01
#> 112 Jawa Barat Tasikmalaya 163 74.08
#> 113 Jawa Tengah Tegal 203 47.14
#> 114 Jawa Tengah Temanggung 65 49.74
#> 115 Jawa Timur Trenggalek 33 58.06
#> 116 Jawa Timur Tuban 61 67.14
#> 117 Jawa Timur Tulungagung 25 45.12
#> 118 Jawa Tengah Wonogiri 53 62.51
#> 119 Jawa Tengah Wonosobo 98 82.69
#> yhat BYM yhat GLMM
#> 1 3.61 7.95
#> 2 212.28 209.38
#> 3 101.58 100.10
#> 4 24.84 26.32
#> 5 128.22 127.52
#> 6 35.80 35.98
#> 7 122.34 120.77
#> 8 31.76 31.79
#> 9 53.61 53.54
#> 10 44.27 44.42
#> 11 60.00 59.61
#> 12 37.33 37.36
#> 13 334.23 326.08
#> 14 26.47 27.06
#> 15 44.85 45.00
#> 16 44.40 44.28
#> 17 127.77 127.35
#> 18 28.09 28.24
#> 19 74.92 74.23
#> 20 51.79 52.25
#> 21 153.77 151.56
#> 22 34.05 34.40
#> 23 267.13 262.30
#> 24 81.61 80.71
#> 25 69.39 68.70
#> 26 53.58 53.71
#> 27 44.15 44.05
#> 28 137.81 137.00
#> 29 78.99 78.90
#> 30 59.06 58.72
#> 31 20.51 20.97
#> 32 32.72 32.95
#> 33 84.84 84.18
#> 34 105.50 104.38
#> 35 72.32 73.08
#> 36 91.21 90.20
#> 37 62.45 61.79
#> 38 19.80 19.85
#> 39 7.59 8.14
#> 40 24.40 24.58
#> 41 21.93 22.03
#> 42 72.10 71.37
#> 43 13.37 14.18
#> 44 18.45 18.52
#> 45 51.32 51.06
#> 46 7.15 7.83
#> 47 28.15 28.11
#> 48 15.13 15.42
#> 49 34.04 34.05
#> 50 29.24 29.59
#> 51 42.06 42.34
#> 52 12.76 13.18
#> 53 7.68 8.02
#> 54 7.63 8.15
#> 55 38.79 38.93
#> 56 4.72 4.82
#> 57 28.51 28.99
#> 58 16.47 16.79
#> 59 19.08 19.27
#> 60 9.49 10.13
#> 61 17.75 18.07
#> 62 17.78 17.70
#> 63 16.42 16.86
#> 64 23.27 23.49
#> 65 11.69 11.56
#> 66 47.64 47.52
#> 67 12.51 13.59
#> 68 63.25 63.18
#> 69 13.87 14.40
#> 70 16.29 16.64
#> 71 40.36 40.14
#> 72 26.30 26.79
#> 73 54.14 54.11
#> 74 37.67 37.77
#> 75 51.28 51.93
#> 76 43.28 43.22
#> 77 45.41 45.28
#> 78 107.99 106.37
#> 79 34.11 34.27
#> 80 37.53 37.91
#> 81 126.13 123.09
#> 82 19.53 19.89
#> 83 43.78 43.83
#> 84 53.43 53.40
#> 85 35.96 36.28
#> 86 33.35 33.77
#> 87 37.29 38.85
#> 88 9.76 10.50
#> 89 20.56 20.82
#> 90 46.89 46.59
#> 91 21.92 22.62
#> 92 84.73 84.11
#> 93 56.79 56.55
#> 94 132.52 131.26
#> 95 101.34 99.58
#> 96 25.91 26.10
#> 97 46.05 46.30
#> 98 52.41 52.42
#> 99 37.82 39.14
#> 100 37.90 38.17
#> 101 105.97 105.30
#> 102 54.77 54.17
#> 103 34.73 35.48
#> 104 40.19 40.14
#> 105 46.48 46.71
#> 106 20.23 20.51
#> 107 164.68 161.84
#> 108 40.24 40.22
#> 109 70.79 69.64
#> 110 36.13 36.92
#> 111 103.23 102.89
#> 112 156.05 155.32
#> 113 187.94 184.20
#> 114 63.87 63.64
#> 115 34.80 34.98
#> 116 62.30 61.81
#> 117 26.42 26.69
#> 118 53.69 53.77
#> 119 97.71 96.88Conclusion
Model yand dipilih adalah ICAR yang memiliki nilai DIC dan MAD lebih rendah daripada BYM.
yhat_bym <- bym_mod$summary.fitted.values$mean*jawa@data$E_d
jawa@data$yhat_bym <- bym_mod$summary.fitted.values$mean*jawa@data$E_d
my.palette.1<-brewer.pal(n=4,name="OrRd")
spplot(jawa,'yhat_bym', cuts=3, col.regions=my.palette.1, cex=c(0.3,1,0.3), main="Hasil dugaan menggunakan CAR BYM",at = c(min(jawa$yhat_bym) - 1, 50, 100, 150 , max(jawa$yhat_bym) + 1))