Tugas Epidem
Kiki Amelia
2026-05-20
#================ 2.1 Input Data Stunting (sesuai file) ================#
library(readxl); library(dplyr); library(stringr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2); library(forcats); library(sf)## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(raster)## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectStunting <- read.csv(“Stunting4.csv”, header = TRUE, sep = “;”, dec = “.”) head(Stunting)
# Filter hanya wilayah Jawa Barat
# Membaca file shapefile
jabar_sf <- st_read("Jabar.dbf")## Reading layer `Jabar' from data source
## `/Users/kikiamelia/Desktop/Amelia/Jabar.dbf' using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
## CRS: NA# Cek data
print(jabar_sf)## Simple feature collection with 27 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 106.3705 ymin: -7.823398 xmax: 108.8338 ymax: -5.91377
## CRS: NA
## First 10 features:
## PROVNO PROVINSI KABKOT PROV12_2 KABKOT12_2 Kabupaten
## 1 32 JAWA BARAT Kab Bogor JAWA BARAT BOGOR Kab Bogor
## 2 32 JAWA BARAT Kab Sukabumi JAWA BARAT SUKABUMI Kab Sukabumi
## 3 32 JAWA BARAT Kab Cianjur JAWA BARAT CIANJUR Kab Cianjur
## 4 32 JAWA BARAT Kab Bandung JAWA BARAT BANDUNG Kab Bandung
## 5 32 JAWA BARAT Kab Garut JAWA BARAT GARUT Kab Garut
## 6 32 JAWA BARAT Kab Tasikmalaya JAWA BARAT TASIKMALAYA Kab Tasikmalaya
## 7 32 JAWA BARAT Kab Pangandaran JAWA BARAT CIAMIS Kab Ciamis
## 8 32 JAWA BARAT Kab Kuningan JAWA BARAT KUNINGAN Kab Kuningan
## 9 32 JAWA BARAT Kab Cirebon JAWA BARAT CIREBON Kab Cirebon
## 10 32 JAWA BARAT Kab Majalengka JAWA BARAT MAJALENGKA Kab Majalengka
## longitude latitude ID geometry
## 1 106.7684 -6.561213 1 MULTIPOLYGON (((106.9909 -6...
## 2 106.7101 -7.074610 2 MULTIPOLYGON (((106.488 -7....
## 3 107.1595 -7.128689 3 MULTIPOLYGON (((106.8641 -7...
## 4 107.6104 -7.099844 4 MULTIPOLYGON (((107.4771 -7...
## 5 107.7887 -7.359583 5 MULTIPOLYGON (((107.8113 -7...
## 6 108.1412 -7.496875 6 MULTIPOLYGON (((108.0603 -7...
## 7 108.4287 -7.289809 7 MULTIPOLYGON (((108.21 -7.2...
## 8 108.5594 -7.003270 8 MULTIPOLYGON (((108.6072 -6...
## 9 108.5513 -6.745514 9 MULTIPOLYGON (((108.685 -6....
## 10 108.2570 -6.815748 10 MULTIPOLYGON (((108.1809 -6...plot(st_geometry(jabar_sf)) # Plot dasar untuk melihat peta
B=ggplot(data = jabar_sf) +
geom_sf(fill = "white", color = "black") +
# Menambahkan label ID secara otomatis di tengah poligon
geom_sf_text(aes(label = ID),
size = 3, # Ukuran teks (setara txt.cex)
color = "black",
fontface = "bold",
check_overlap = TRUE) +
theme_minimal()
B
#================ 2. Preprocessing Data Stunting ================#
# Anggap nama dataset Anda adalah 'data_stunting'
# Pastikan nama Kab/Kota seragam (Uppercase) untuk Join
Stunting <- read.csv("DataStuntingJabar5.csv", header = TRUE, sep = ";", dec = ".")
head(Stunting)## No id kode_provinsi nama_provinsi kode_kabupaten_kota Kab.Kota
## 1 1 1 32 JAWA BARAT 3201 Kab Bogor
## 2 2 2 32 JAWA BARAT 3202 Kab Sukabumi
## 3 3 3 32 JAWA BARAT 3203 Kab Cianjur
## 4 4 4 32 JAWA BARAT 3204 Kab Bandung
## 5 5 5 32 JAWA BARAT 3205 Kab Garut
## 6 6 6 32 JAWA BARAT 3206 Kab Tasikmalaya
## persentase_balita_stunting satuan tahun Tingkat.Pendidikan.Ibu PDRB
## 1 4.06 PERSEN 2019 30.04 237113
## 2 8.29 PERSEN 2019 13.48 67406
## 3 6.61 PERSEN 2019 13.25 46807
## 4 7.32 PERSEN 2019 30.23 124001
## 5 4.80 PERSEN 2019 17.79 57579
## 6 15.06 PERSEN 2019 14.93 37219
## Tingkat.Kemiskinan JumlahBalita JumlahKasusStunting Ei SMR
## 1 6.66 442162 132825.46 587303.73 22.61615
## 2 6.22 184865 24919.80 46067.99 54.09353
## 3 9.15 179399 23770.37 42643.80 55.74167
## 4 5.94 240699 72763.31 175140.55 41.54566
## 5 8.98 194656 34629.30 67408.01 51.37268
## 6 9.12 108122 16142.61 17453.72 92.48812
## lnSMR AirMinumLayak SanitasiLayak PHBS persentase_pemberian_asi IT IS
## 1 3.118664 91.02 56.00 53.91 53.12 1 1
## 2 3.990715 79.57 69.00 51.10 57.15 1 2
## 3 4.020728 81.50 62.00 62.00 68.78 1 3
## 4 3.726793 97.84 80.00 57.11 63.84 1 4
## 5 3.939106 84.17 76.09 60.02 74.32 1 5
## 6 4.527080 76.18 81.00 54.01 66.77 1 6
## ST
## 1 6.305517
## 2 9.354103
## 3 8.715372
## 4 8.664596
## 5 10.102140
## 6 11.205440# Join Data Peta dengan Data Stunting
jabar_stunting_join <- jabar_sf %>%
left_join(Stunting, by = c("Kabupaten" = "Kab.Kota"))## ================= Stunting ================= ##
## Interval (Menyesuaikan dengan persentase_balita_stunting)
# pretty_breaks diatur untuk skala persentase (misal: 10%, 20%, 30%)
pretty_breaks <- c(1.28, 5.37, 9.4, 13.55, 17.64)
# compute labels
labels <- c()
brks <- c(1.28, 5.37, 9.4, 13.55, 17.64)
format(brks, scientific=F)## [1] " 1.28" " 5.37" " 9.40" "13.55" "17.64"# Definisikan variabel baru pada dataset hasil join
# Catatan: Pastikan kolom di data CSV Anda bernama 'persentase_balita_stunting'
labels <- labels[1:length(labels)-1]
jabar_stunting_join$brks <- cut(jabar_stunting_join$persentase_balita_stunting,
breaks = brks,
include.lowest = TRUE,
labels = labels)
brks_scale <- levels(jabar_stunting_join$brks)
labels_scale <- rev(brks_scale)
# Simpan hasil
#ggsave("Stunting_Jabar.png", height = 6, width = 8, dpi = 600)
# Plotting menggunakan ggplot
# Menggunakan geom_sf karena jabar_stunting_join adalah objek 'sf'
A <- ggplot(jabar_stunting_join) +
geom_sf(aes(fill = brks), color="black", size=0.1) +
# Menggunakan Tahun sebagai pengganti Triwulan untuk facet
geom_sf_text(aes(label = id),
size = 1.5, # Ukuran teks (setara txt.cex)
color = "black",
fontface = "bold",
check_overlap = TRUE) +
facet_wrap(.~tahun, ncol=3) +
# Skala warna disesuaikan dengan label kategori stunting
scale_fill_manual(
labels = c("<5.36", "[5.37 – 9.45)", "[9.46 – 13.54)", "[13.55 – 17.63)", ">=17.64"),
values = c("grey95", "#FFFDEB", "#F97316", "#DC2626", "#7F1D1D")
) +
theme_bw() +
ylab("Northing") + xlab("Easting") +
theme(axis.text.x = element_text(angle = 90)) +
theme(legend.position = "bottom") +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
labs(fill = "Stunting (%)") +
theme(legend.position="right", text = element_text(size=25)) +
theme(text = element_text(size=13))
A
# Modeling #################
library(mclust)## Package 'mclust' version 6.1.2
## Type 'citation("mclust")' for citing this R package in publications.##
## Attaching package: 'mclust'## The following object is masked from 'package:dplyr':
##
## countlibrary(INLA)## Loading required package: Matrix## This is INLA_25.06.07 built 2025-06-11 19:15:03 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - List available models/likelihoods/etc with inla.list.models()
## - Use inla.doc(<NAME>) to access documentation
## - Consider upgrading R-INLA to testing[26.05.10] or stable[25.10.19].library(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(lattice)
library(RColorBrewer)
library(glmnet)## Loaded glmnet 4.1-10library(caret)
library(readxl); library(dplyr); library(stringr)
library(ggplot2); library(forcats); library(sf)
library(raster)
library(sf)
library(spdep)
library(spatialreg)##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionlibrary(spdep)
library(sp)
library(sf)
library(ggplot2)
library(dplyr)###Autokorelasi spasial
WB_nb <- poly2nb(jabar_sf, queen = TRUE)
lw <- nb2listw(WB_nb, style = "W", zero.policy = TRUE)
W_matrix <- listw2mat(lw)
library(Matrix)
W1 <- as(as_dgRMatrix_listw(lw), "CsparseMatrix")
print(W1)## 27 x 27 sparse Matrix of class "dgCMatrix"## [[ suppressing 27 column names '1', '2', '3' ... ]]##
## 1 . 0.1250000 0.1250000 . . . .
## 2 0.3333333 . 0.3333333 . . . .
## 3 0.1666667 0.1666667 . 0.1666667 0.1666667 . .
## 4 . . 0.1428571 . 0.1428571 . .
## 5 . . 0.2500000 0.2500000 . 0.2500000 .
## 6 . . . . 0.1666667 . 0.1666667
## 7 . . . . . 0.1666667 .
## 8 . . . . . . 0.3333333
## 9 . . . . . . .
## 10 . . . . . 0.1666667 0.1666667
## 11 . . . 0.1666667 0.1666667 0.1666667 .
## 12 . . . . . . .
## 13 . . . 0.1666667 . . .
## 14 0.2000000 . 0.2000000 . . . .
## 15 0.2500000 . . . . . .
## 16 0.3333333 . . . . . .
## 17 . . 0.1666667 0.1666667 . . .
## 18 . . . . . 0.5000000 0.5000000
## 19 1.0000000 . . . . . .
## 20 . 1.0000000 . . . . .
## 21 . . . 0.3333333 . . .
## 22 . . . . . . .
## 23 0.3333333 . . . . . .
## 24 0.5000000 . . . . . .
## 25 . . . 0.3333333 . . .
## 26 . . . . . 0.5000000 0.5000000
## 27 . . . . . . 1.0000000
##
## 1 . . . . . . 0.1250000
## 2 . . . . . . .
## 3 . . . . . . 0.1666667
## 4 . . . 0.1428571 . 0.1428571 .
## 5 . . . 0.2500000 . . .
## 6 . . 0.1666667 0.1666667 . . .
## 7 0.1666667 . 0.1666667 . . . .
## 8 . 0.3333333 0.3333333 . . . .
## 9 0.2500000 . 0.2500000 . 0.2500000 . .
## 10 0.1666667 0.1666667 . 0.1666667 0.1666667 . .
## 11 . . 0.1666667 . 0.1666667 0.1666667 .
## 12 . 0.2500000 0.2500000 0.2500000 . 0.2500000 .
## 13 . . . 0.1666667 0.1666667 . 0.1666667
## 14 . . . . . 0.2000000 .
## 15 . . . . . 0.2500000 0.2500000
## 16 . . . . . . .
## 17 . . . . . 0.1666667 0.1666667
## 18 . . . . . . .
## 19 . . . . . . .
## 20 . . . . . . .
## 21 . . . . . . .
## 22 . 1.0000000 . . . . .
## 23 . . . . . . .
## 24 . . . . . . .
## 25 . . . . . . .
## 26 . . . . . . .
## 27 . . . . . . .
##
## 1 0.1250000 0.1250000 . . 0.125 . . .
## 2 . . . . . 0.3333333 . .
## 3 . . 0.1666667 . . . . .
## 4 . . 0.1428571 . . . 0.1428571 .
## 5 . . . . . . . .
## 6 . . . 0.1666667 . . . .
## 7 . . . 0.1666667 . . . .
## 8 . . . . . . . .
## 9 . . . . . . . 0.25
## 10 . . . . . . . .
## 11 . . . . . . . .
## 12 . . . . . . . .
## 13 0.1666667 . 0.1666667 . . . . .
## 14 0.2000000 . 0.2000000 . . . . .
## 15 . 0.2500000 . . . . . .
## 16 0.3333333 . . . . . . .
## 17 . . . . . . 0.1666667 .
## 18 . . . . . . . .
## 19 . . . . . . . .
## 20 . . . . . . . .
## 21 . . 0.3333333 . . . . .
## 22 . . . . . . . .
## 23 . 0.3333333 . . . . . .
## 24 . . . . . . . .
## 25 . . 0.3333333 . . . 0.3333333 .
## 26 . . . . . . . .
## 27 . . . . . . . .
##
## 1 0.1250000 0.1250000 . . .
## 2 . . . . .
## 3 . . . . .
## 4 . . 0.1428571 . .
## 5 . . . . .
## 6 . . . 0.1666667 .
## 7 . . . 0.1666667 0.1666667
## 8 . . . . .
## 9 . . . . .
## 10 . . . . .
## 11 . . . . .
## 12 . . . . .
## 13 . . . . .
## 14 . . . . .
## 15 . . . . .
## 16 0.3333333 . . . .
## 17 . . 0.1666667 . .
## 18 . . . . .
## 19 . . . . .
## 20 . . . . .
## 21 . . 0.3333333 . .
## 22 . . . . .
## 23 . 0.3333333 . . .
## 24 0.5000000 . . . .
## 25 . . . . .
## 26 . . . . .
## 27 . . . . .Matriks bobot spasial
row.names(jabar_sf) <- as.character(1:27)
Plot jaringan tetangga
CoordK <- coordinates(jabar_sf) plot(jabar_sf, axes=T, col=“gray90”) text(CoordK[,1], CoordK[,2], row.names(jabar_sf), col=“black”, cex=0.8, pos=1.5) points(CoordK[,1], CoordK[,2], pch=19, cex=0.7, col=“blue”) plot.nb(WB_nb, CoordK, add=TRUE, col=“red”, lwd=2)
##### Open Your Data
DataL<-read.csv("DataStuntingJabar5.csv", header = TRUE, sep = ";", dec = ".")
head(DataL)## No id kode_provinsi nama_provinsi kode_kabupaten_kota Kab.Kota
## 1 1 1 32 JAWA BARAT 3201 Kab Bogor
## 2 2 2 32 JAWA BARAT 3202 Kab Sukabumi
## 3 3 3 32 JAWA BARAT 3203 Kab Cianjur
## 4 4 4 32 JAWA BARAT 3204 Kab Bandung
## 5 5 5 32 JAWA BARAT 3205 Kab Garut
## 6 6 6 32 JAWA BARAT 3206 Kab Tasikmalaya
## persentase_balita_stunting satuan tahun Tingkat.Pendidikan.Ibu PDRB
## 1 4.06 PERSEN 2019 30.04 237113
## 2 8.29 PERSEN 2019 13.48 67406
## 3 6.61 PERSEN 2019 13.25 46807
## 4 7.32 PERSEN 2019 30.23 124001
## 5 4.80 PERSEN 2019 17.79 57579
## 6 15.06 PERSEN 2019 14.93 37219
## Tingkat.Kemiskinan JumlahBalita JumlahKasusStunting Ei SMR
## 1 6.66 442162 132825.46 587303.73 22.61615
## 2 6.22 184865 24919.80 46067.99 54.09353
## 3 9.15 179399 23770.37 42643.80 55.74167
## 4 5.94 240699 72763.31 175140.55 41.54566
## 5 8.98 194656 34629.30 67408.01 51.37268
## 6 9.12 108122 16142.61 17453.72 92.48812
## lnSMR AirMinumLayak SanitasiLayak PHBS persentase_pemberian_asi IT IS
## 1 3.118664 91.02 56.00 53.91 53.12 1 1
## 2 3.990715 79.57 69.00 51.10 57.15 1 2
## 3 4.020728 81.50 62.00 62.00 68.78 1 3
## 4 3.726793 97.84 80.00 57.11 63.84 1 4
## 5 3.939106 84.17 76.09 60.02 74.32 1 5
## 6 4.527080 76.18 81.00 54.01 66.77 1 6
## ST
## 1 6.305517
## 2 9.354103
## 3 8.715372
## 4 8.664596
## 5 10.102140
## 6 11.205440DataLY<-data.frame(Case=DataL$persentase_balita_stunting)
Case=DataL$persentase_balita_stunting# Global Moran's I
Moran19 <- moran.test(DataL$persentase_balita_stunting[1:27], lw)
Moran20 <- moran.test(DataL$persentase_balita_stunting[28:54], lw)
Moran21 <- moran.test(DataL$persentase_balita_stunting[55:81], lw)
Moran22 <- moran.test(DataL$persentase_balita_stunting[82:108], lw)
Moran23 <- moran.test(DataL$persentase_balita_stunting[109:135], lw)
Moran24 <- moran.test(DataL$persentase_balita_stunting[136:162], lw)
Moran=c(Moran19$estimate[1],Moran20$estimate[1],Moran21$estimate[1],Moran22$estimate[1],Moran23$estimate[1],Moran24$estimate[1])
pv.Moran=c(Moran19$p.value,Moran20$p.value,Moran21$p.value,Moran22$p.value,Moran23$p.value,Moran24$p.value)
Moran.Indeks=data.frame(Moran,pv.Moran)
Moran.Indeks## Moran pv.Moran
## 1 0.19436344 0.038336774
## 2 0.03784775 0.293879138
## 3 0.02207477 0.332934871
## 4 0.26544701 0.014847020
## 5 0.30035373 0.008385097
## 6 0.18586762 0.050010172####AUTOKOREALASI library(MoranST) DBD<-DataBdg$DBD MoranST<-MoranST.MC(DBD,WB,1000) LMoranST<-LocalMoranST.MC(DBD,W,1000)
library(spdep) library(sf) library(Matrix)
# Membuat neighbour
nb <- poly2nb(jabar_sf, queen = TRUE)
# Spatial weight
lw <- nb2listw(nb,
style = "W",
zero.policy = TRUE)
# Matriks bobot spasial
W <- listw2mat(lw)
# Jumlah tahun
T <- length(unique(DataL$tahun))
# Membuat bobot spatio-temporal
W_st <- kronecker(Diagonal(T), W)
# Mengubah ke listw
lw_st <- mat2listw(as.matrix(W_st))## Warning in mat2listw(as.matrix(W_st)): style is M (missing); style should be
## set to a valid value## Warning in mat2listw(as.matrix(W_st)): neighbour object has 6 sub-graphs# Moran Spatio-Temporal
moran_st <- moran.test(
DataL$persentase_balita_stunting,
lw_st,
zero.policy = TRUE
)
moran_st##
## Moran I test under randomisation
##
## data: DataL$persentase_balita_stunting
## weights: lw_st
##
## Moran I statistic standard deviate = 5.5602, p-value = 1.347e-08
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.33306021 -0.00621118 0.00372314data.frame(
Pengujian = "Moran Spatio-temporal",
Statistik_I = round(moran_st$estimate[1],4),
P_value = round(moran_st$p.value,6)
)## Pengujian Statistik_I P_value
## Moran I statistic Moran Spatio-temporal 0.3331 0DataLE<-data.frame(E=DataL$E)
DataLX<-data.frame(Tingkat.Pendidikan.Ibu=DataL$Tingkat.Pendidikan.Ibu, PDRB=DataL$PDRB, Tingkat.Kemiskinan=DataL$Tingkat.Kemiskinan, AirMinumLayak=DataL$AirMinumLayak,
SanitasiLayak=DataL$SanitasiLayak, PHBS=DataL$PHBS, persentase_pemberian_asi=DataL$persentase_pemberian_asi)
DataLlnSMR<-DataL$lnSMR
DataLXlnSMR<-data.frame(DataLlnSMR,DataLX)
#FullDataL<-data.frame(DataLY, DataLE, DataLXR)
#head(FullDataL)
FullDataL<-data.frame(DataLY, DataLE, DataLX)
head(FullDataL)## Case E Tingkat.Pendidikan.Ibu PDRB Tingkat.Kemiskinan
## 1 4.06 587303.73 30.04 237113 6.66
## 2 8.29 46067.99 13.48 67406 6.22
## 3 6.61 42643.80 13.25 46807 9.15
## 4 7.32 175140.55 30.23 124001 5.94
## 5 4.80 67408.01 17.79 57579 8.98
## 6 15.06 17453.72 14.93 37219 9.12
## AirMinumLayak SanitasiLayak PHBS persentase_pemberian_asi
## 1 91.02 56.00 53.91 53.12
## 2 79.57 69.00 51.10 57.15
## 3 81.50 62.00 62.00 68.78
## 4 97.84 80.00 57.11 63.84
## 5 84.17 76.09 60.02 74.32
## 6 76.18 81.00 54.01 66.77DATAL1<-FullDataL
DATAL1S<-cbind(DATAL1[,1:2], scale(DATAL1[,c(-1,-2)]))
DATAL1S$IS1<-DataL$IS
DATAL1S$IS2<-DataL$IS
DATAL1S$IS3<-DataL$IS
DATAL1S$IT1<-DataL$IT
DATAL1S$IT2<-DataL$IT
DATAL1S$IT3<-DataL$IT###===4. MENGHITUNG SMR===###
#Estimasi SMR#
Ei=DataL$Ei
SMR_19=DataL$persentase_balita_stunting[1:27]/Ei[1:27]
SMR_20=DataL$persentase_balita_stunting[28:54]/Ei[28:54]
SMR_21=DataL$persentase_balita_stunting[55:81]/Ei[55:81]
SMR_22=DataL$persentase_balita_stunting[82:108]/Ei[82:108]
SMR_23=DataL$persentase_balita_stunting[109:135]/Ei[109:135]
SMR_24=DataL$persentase_balita_stunting[136:162]/Ei[136:162]
boxplot(SMR_19,SMR_20,SMR_21,SMR_22, SMR_23,SMR_24, main="Boxplot SMR",names=c("2019","2020","2021","2022","2023","2024"))
##### Define your own prior
##### (1) Inverse Gamma Prior
prior.c5=c(1,1*10^(-5))
igprior5=list(theta=list(param=prior.c5))
prior.c1=c(1,0.01)
igprior1=list(theta=list(param=prior.c1))
##### (2) Penalize Complexity Prior
pcprior <- list(prec = list(prior="pc.prec", param = c(3,0.01)))
##### (3) Half Cauchy Prior
halfcauchy = "expression:
lambda = 25;
precision = exp(log_precision);
logdens = -1.5*log_precision-log(pi*lambda)-log(1+1/(precision*lambda^2));
log_jacobian = log_precision;
return(logdens+log_jacobian);"
hcprior = list(prec = list(prior = halfcauchy))
hcprior<-hcprior
##### (4) Uniform Prior
sdunif="expression:
logdens=-log_precision/2;
return(logdens)"
uprior=list(prec=list(prior=sdunif))
control <- list(
predictor = list(compute = TRUE, cdf=c(log(1))),
results = list(return.marginals.random = TRUE, return.marginals.predictor=TRUE),
compute = list(hyperpar=TRUE, return.marginals=TRUE, dic=TRUE, mlik = TRUE, cpo = TRUE,
po = TRUE, waic=TRUE, graph=TRUE, gdensity=TRUE, openmp.strategy="huge"))##DATA EXPLORATION
library(ggcorrplot)
DataLX1<-DataLX[-c(1:60),]
VIF<-diag(solve(cor(DataLX1)))
DataLX1<-DataLX[-c(1:60),]
VIF<-diag(solve(cor(DataLX1)))
DataLXlnSMR1<-DataLXlnSMR[-c(1:60),]
Spearman.Correlation <- round(cor(DataLXlnSMR1,method="spearman"), 1)
r<-ggcorrplot(Spearman.Correlation,type = "lower",
lab = TRUE)## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggcorrplot package.
## Please report the issue at <https://github.com/kassambara/ggcorrplot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.library(extrafont)## Registering fonts with R#### Variable Selection
library(My.stepwise)
### Add Covariate
###(1) Parameteric Poisson Regression
#GLM0<-glm(Case~MeanTemperature+Rainfall+SunshineDuration+RelativeHumidity+Evaporation+
# MeanTemperatureL1+RainfallL1+SunshineDurationL1+RelativeHumidityL1+EvaporationL1+
# MeanTemperatureL2+RainfallL2+SunshineDurationL2+RelativeHumidityL2+EvaporationL2+offset(log(E)), family="poisson", data=FullDataL)
#summary(GLM0)
GLM0<-glm(Case~Tingkat.Pendidikan.Ibu+PDRB+Tingkat.Kemiskinan+AirMinumLayak+SanitasiLayak+PHBS+persentase_pemberian_asi +offset(log(E)), family="gaussian", data=FullDataL)
summary(GLM0)##
## Call:
## glm(formula = Case ~ Tingkat.Pendidikan.Ibu + PDRB + Tingkat.Kemiskinan +
## AirMinumLayak + SanitasiLayak + PHBS + persentase_pemberian_asi +
## offset(log(E)), family = "gaussian", data = FullDataL)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.638e+00 6.268e+00 1.378 0.1702
## Tingkat.Pendidikan.Ibu 9.021e-02 3.603e-02 2.503 0.0133 *
## PDRB -2.337e-05 3.862e-06 -6.052 1.04e-08 ***
## Tingkat.Kemiskinan 2.644e-01 1.940e-01 1.363 0.1748
## AirMinumLayak -1.020e-01 7.280e-02 -1.401 0.1633
## SanitasiLayak -2.658e-02 1.929e-02 -1.377 0.1704
## PHBS -6.549e-02 3.174e-02 -2.063 0.0408 *
## persentase_pemberian_asi 2.115e-02 2.309e-02 0.916 0.3612
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 16.82728)
##
## Null deviance: 3722.1 on 161 degrees of freedom
## Residual deviance: 2591.4 on 154 degrees of freedom
## AIC: 926.86
##
## Number of Fisher Scoring iterations: 2GLM1<-glm(Case~AirMinumLayak+PHBS+offset(log(E)), family="gaussian", data=FullDataL)
summary(GLM0)##
## Call:
## glm(formula = Case ~ Tingkat.Pendidikan.Ibu + PDRB + Tingkat.Kemiskinan +
## AirMinumLayak + SanitasiLayak + PHBS + persentase_pemberian_asi +
## offset(log(E)), family = "gaussian", data = FullDataL)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.638e+00 6.268e+00 1.378 0.1702
## Tingkat.Pendidikan.Ibu 9.021e-02 3.603e-02 2.503 0.0133 *
## PDRB -2.337e-05 3.862e-06 -6.052 1.04e-08 ***
## Tingkat.Kemiskinan 2.644e-01 1.940e-01 1.363 0.1748
## AirMinumLayak -1.020e-01 7.280e-02 -1.401 0.1633
## SanitasiLayak -2.658e-02 1.929e-02 -1.377 0.1704
## PHBS -6.549e-02 3.174e-02 -2.063 0.0408 *
## persentase_pemberian_asi 2.115e-02 2.309e-02 0.916 0.3612
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 16.82728)
##
## Null deviance: 3722.1 on 161 degrees of freedom
## Residual deviance: 2591.4 on 154 degrees of freedom
## AIC: 926.86
##
## Number of Fisher Scoring iterations: 2#moran residual
residual=resid(GLM1)
Moran19=moran.test(residual[1:27],listw=lw,zero.policy=TRUE)
Moran20=moran.test(residual[28:54],listw=lw,zero.policy=TRUE)
Moran21=moran.test(residual[55:81],listw=lw,zero.policy=TRUE)
Moran22=moran.test(residual[82:108],listw=lw,zero.policy=TRUE)
Moran23=moran.test(residual[109:135],listw=lw,zero.policy=TRUE)
Moran24=moran.test(residual[136:162],listw=lw,zero.policy=TRUE)
Moran=c(Moran19$estimate[1],Moran20$estimate[1],Moran21$estimate[1],Moran22$estimate[1],Moran23$estimate[1],Moran24$estimate[1])
pv.Moran=c(Moran19$p.value,Moran20$p.value,Moran21$p.value,Moran22$p.value,Moran23$p.value,Moran24$p.value)
Moran.Indeks=data.frame(Moran,pv.Moran)
Moran.Indeks## Moran pv.Moran
## 1 0.2172541 0.028584469
## 2 0.2662122 0.015237051
## 3 0.1124755 0.140993977
## 4 0.2966179 0.008304696
## 5 0.3274352 0.003820227
## 6 0.2567484 0.014645800###(2) Bayesian Poisson Regression
Formula.Cov<-Case~Tingkat.Pendidikan.Ibu+PDRB+Tingkat.Kemiskinan+AirMinumLayak+SanitasiLayak+PHBS+persentase_pemberian_asi
Model.Cov<- inla(Formula.Cov, data = FullDataL, family = "Gaussian", E=E,
control.fixed = list(mean.intercept = 0, prec.intercept = 0.000001, mean = 0, prec =0.000001),
control.predictor=list(compute=TRUE, cdf=c(log(1))),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE),
control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = "eb", strategy = "simplified.laplace"))
summary(Model.Cov)## Time used:
## Pre = 1.72, Running = 8.92, Post = 0.935, Total = 11.6
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 22.762 5.797 11.400 22.762 34.124 22.762 0
## Tingkat.Pendidikan.Ibu 0.083 0.033 0.018 0.083 0.149 0.083 0
## PDRB 0.000 0.000 0.000 0.000 0.000 0.000 0
## Tingkat.Kemiskinan 0.261 0.179 -0.091 0.261 0.613 0.261 0
## AirMinumLayak -0.133 0.067 -0.265 -0.133 -0.001 -0.133 0
## SanitasiLayak -0.017 0.018 -0.052 -0.017 0.018 -0.017 0
## PHBS -0.091 0.029 -0.148 -0.091 -0.033 -0.091 0
## persentase_pemberian_asi 0.001 0.021 -0.041 0.001 0.043 0.001 0
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.069 0.001 0.067 0.069
## 0.975quant mode
## Precision for the Gaussian observations 0.072 0.069
##
## Deviance Information Criterion (DIC) ...............: 902.06
## Deviance Information Criterion (DIC, saturated) ....: 172.31
## Effective number of parameters .....................: 8.00
##
## Watanabe-Akaike information criterion (WAIC) ...: 902.96
## Effective number of parameters .................: 8.29
##
## Marginal log-Likelihood: -547.74
## CPO, PIT 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)')#### BAYESIAN RIDGE REGRESION
#### (2a) Bayesian Poisson ridge regression "Copy" approch
#FullDataLS<-cbind(FullDataL[,1:2], scale(FullDataL[,c(-1,-2)]))
FullDataLS <- data.frame(
FullDataL[,1:2],
scale(FullDataL[, -c(1,2)])
)
head(FullDataLS)## Case E Tingkat.Pendidikan.Ibu PDRB Tingkat.Kemiskinan
## 1 4.06 587303.73 -0.1707676 1.5406327 -0.5899169
## 2 8.29 46067.99 -1.2147292 -0.2261874 -0.7506988
## 3 6.61 42643.80 -1.2292287 -0.4406436 0.3199622
## 4 7.32 175140.55 -0.1587898 0.3630233 -0.8530145
## 5 4.80 67408.01 -0.9430218 -0.3284963 0.2578419
## 6 15.06 17453.72 -1.1233195 -0.5404643 0.3089998
## AirMinumLayak SanitasiLayak PHBS persentase_pemberian_asi
## 1 -0.4534219 -1.6711889 -0.8356904 -1.2566722
## 2 -2.3562983 -0.9313206 -1.0758064 -0.9897552
## 3 -2.0355514 -1.3297113 -0.1443956 -0.2194709
## 4 0.6799944 -0.3052782 -0.5622487 -0.5466596
## 5 -1.5918239 -0.5278079 -0.3135876 0.1474573
## 6 -2.9196826 -0.2483653 -0.8271454 -0.3525983m <- nrow(FullDataLS)
FullDataLS$beta1 <- rep(1,m)
FullDataLS$beta2 <- rep(2,m)
FullDataLS$beta3 <- rep(3,m)
FullDataLS$beta4 <- rep(4,m)
FullDataLS$beta5 <- rep(5,m)
FullDataLS$beta6 <- rep(6,m)
FullDataLS$beta7 <- rep(7,m)###Model dengan effect fixed saja
param.beta = list(prec = list(param = c(1.0e-6, 1.0e-6)))
Formula.Ridge1 = Case ~ f(beta1, Tingkat.Pendidikan.Ibu, model="iid", values = c(1:7), hyper = param.beta)+
f(beta2, PDRB, copy="beta1", fixed=T)+
f(beta3, Tingkat.Kemiskinan, copy="beta1", fixed=T)+
f(beta4, AirMinumLayak, copy="beta1", fixed=T)+
f(beta5, SanitasiLayak, copy="beta1", fixed=T)+
f(beta6, PHBS, copy="beta1", fixed=T)+
f(beta7, persentase_pemberian_asi, copy="beta1", fixed=T)
Model.Ridge1 <- inla(Formula.Ridge1 , family="Gaussian", data=FullDataLS, E=E,
control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001),
control.predictor=list(compute=TRUE, cdf=c(log(1))),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE),
control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = "eb", strategy = "simplified.laplace"))##
## *** inla.core.safe: The inla program failed, but will rerun in case better initial values may help. try=1/1
##
## *** inla.core.safe: rerun with improved initial valuesCoeff.Ridge1 <- rbind(Model.Ridge1$summary.fixed, Model.Ridge1$summary.random$beta1[,-1])
round(Coeff.Ridge1,4) ## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 6.9491 0.2987 6.3636 6.9491 7.5346 6.9491 0
## 1 0.7270 0.3860 -0.0296 0.7270 1.4836 0.7270 0
## 2 -0.9368 0.2812 -1.4880 -0.9368 -0.3857 -0.9368 0
## 3 0.4188 0.3807 -0.3274 0.4188 1.1651 0.4188 0
## 4 -0.4853 0.3369 -1.1457 -0.4853 0.1750 -0.4853 0
## 5 -0.2921 0.3022 -0.8845 -0.2921 0.3003 -0.2921 0
## 6 -0.9070 0.2820 -1.4598 -0.9070 -0.3543 -0.9070 0
## 7 -0.0007 0.3135 -0.6152 -0.0007 0.6137 -0.0007 0summary(Model.Ridge1)## Time used:
## Pre = 1.52, Running = 0.992, Post = 0.04, Total = 2.55
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 6.949 0.299 6.364 6.949 7.535 6.949 0
##
## Random effects:
## Name Model
## beta1 IID model
## beta2 Copy
## beta3 Copy
## beta4 Copy
## beta5 Copy
## beta6 Copy
## beta7 Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.069 0.008 0.055 0.069
## Precision for beta1 2.539 2.070 0.463 1.969
## Beta for beta2 1.096 0.285 0.541 1.094
## Beta for beta3 0.967 0.309 0.360 0.966
## Beta for beta4 0.977 0.306 0.378 0.975
## Beta for beta5 0.935 0.317 0.313 0.934
## Beta for beta6 1.090 0.287 0.531 1.087
## Beta for beta7 0.912 0.325 0.273 0.912
## 0.975quant mode
## Precision for the Gaussian observations 0.086 0.068
## Precision for beta1 8.016 1.163
## Beta for beta2 1.662 1.086
## Beta for beta3 1.577 0.963
## Beta for beta4 1.583 0.969
## Beta for beta5 1.562 0.931
## Beta for beta6 1.661 1.078
## Beta for beta7 1.551 0.912
##
## Deviance Information Criterion (DIC) ...............: 900.89
## Deviance Information Criterion (DIC, saturated) ....: 170.41
## Effective number of parameters .....................: 6.47
##
## Watanabe-Akaike information criterion (WAIC) ...: 901.27
## Effective number of parameters .................: 6.49
##
## Marginal log-Likelihood: -485.98
## CPO, PIT 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)')#### (2c) Bayesian Poisson ridge regression Matrix Approach
n <- nrow(DataLX)
X <- matrix(1,nrow = n, ncol= 1)
Case<-DataLY$Case
E<-DataLE$E
Z <- as.matrix(DataLX)
pX = ncol(X); pZ = ncol(Z)
idx.X = c(1:pX, rep(NA, pZ))
idx.Z = c(rep(NA,pX), 1:pZ)
hyper.fixed = list(prec = list(initial = log(1.0e-9), fixed=TRUE))
param.data = list(prec = list(param = c(1.0e-3, 1.0e-3)))
param.Z <- list(prec = list(param = c(1.0e-3, 1.0e-3)))
Formula.Ridge2 = Case ~ -1 + f(idx.X, model="iid", hyper = hyper.fixed) + f(idx.Z, model="iid", hyper = param.Z)
Model.Ridge2 = inla(Formula.Ridge2, family="Gaussian", E=E, data = list(Case=Case, E=E, idx.X=idx.X, idx.Z=idx.Z),
control.predictor = list(A=cbind(X, Z),compute=TRUE),
control.fixed = list(mean.intercept = 0, prec.intercept = 0.000001, mean = 0, prec =0.000001),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE),
control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = "eb", strategy = "simplified.laplace"))
Coef.Ridge2 <- Model.Ridge2$summary.random$idx.Z
Coef.Ridge2## ID mean sd 0.025quant 0.5quant 0.975quant
## 1 1 2.742659e-02 1.991895e-02 -0.0116138281 2.742659e-02 6.646701e-02
## 2 2 -1.370985e-05 3.314062e-06 -0.0000202053 -1.370985e-05 -7.214411e-06
## 3 3 1.245215e-02 4.324985e-02 -0.0723159936 1.245215e-02 9.722030e-02
## 4 4 -3.213307e-02 3.626058e-02 -0.1032024934 -3.213307e-02 3.893635e-02
## 5 5 -2.121233e-02 1.648999e-02 -0.0535321211 -2.121233e-02 1.110747e-02
## 6 6 -7.799938e-02 2.348611e-02 -0.1240313076 -7.799938e-02 -3.196746e-02
## 7 7 -1.736978e-03 1.910878e-02 -0.0391894988 -1.736978e-03 3.571554e-02
## mode kld
## 1 2.742659e-02 0.000000e+00
## 2 -1.370985e-05 8.952919e-16
## 3 1.245215e-02 4.247592e-19
## 4 -3.213307e-02 5.396315e-18
## 5 -2.121233e-02 0.000000e+00
## 6 -7.799938e-02 2.262270e-17
## 7 -1.736978e-03 0.000000e+00summary(Model.Ridge2)## Time used:
## Pre = 1.27, Running = 9.44, Post = 1.35, Total = 12.1
## Random effects:
## Name Model
## idx.X IID model
## idx.Z IID model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.068 0.001 0.066 0.068
## Precision for idx.Z 502.541 7.412 492.025 501.578
## 0.975quant mode
## Precision for the Gaussian observations 0.071 0.068
## Precision for idx.Z 520.058 496.227
##
## Deviance Information Criterion (DIC) ...............: 901.82
## Deviance Information Criterion (DIC, saturated) ....: 168.88
## Effective number of parameters .....................: 5.63
##
## Watanabe-Akaike information criterion (WAIC) ...: 901.71
## Effective number of parameters .................: 5.28
##
## Marginal log-Likelihood: -498.62
## CPO, PIT 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)')### Add spatial and temporal component
###(1) Parametric Approach--Temporal Component
FullDataLS$IS1<-DataL$IS
FullDataLS$IS2<-DataL$IS
FullDataLS$IS3<-DataL$IS
FullDataLS$IT1<-DataL$IT
FullDataLS$IT2<-DataL$IT
FullDataLS$IT3<-DataL$ITparam.beta = list(prec = list(param = c(1.0e-6, 1.0e-6)))###(1) Poisson Formula.Ridge1t = Case ~ f(beta1, Tingkat.Pendidikan.Ibu, model=“iid”, values = c(1:7), hyper = param.beta)+ f(beta2, PDRB, copy=“beta1”, fixed=T)+ f(beta3, Tingkat.Kemiskinan, copy=“beta1”, fixed=T)+ f(beta4, AirMinumLayak, copy=“beta1”, fixed=T)+ f(beta5, SanitasiLayak, copy=“beta1”, fixed=T)+ f(beta6, PHBS, copy=“beta1”, fixed=T)+ f(beta7, persentase_pemberian_asi, copy=“beta1”, fixed=T)+ f(IS1, IT1, model = “iid”, hyper = hcprior, constr = T)+ IT2
Model.Ridge1t <- inla(Formula.Ridge1t , family=“Gaussian”, data=FullDataLS, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.Ridge1t <- rbind(Model.Ridge1t\(summary.fixed, Model.Ridge1t\)summary.random$beta1[,-1]) round(Coeff.Ridge1t,4)
summary(Model.Ridge1t)
(3) Spatial component—Besag
(a) Poisson
Formula.Ridge1s = Case ~f(beta1, Tingkat.Pendidikan.Ibu, model=“iid”, values = c(1:7), hyper = param.beta)+ f(beta2, PDRB, copy=“beta1”, fixed=T)+ f(beta3, Tingkat.Kemiskinan, copy=“beta1”, fixed=T)+ f(beta4, AirMinumLayak, copy=“beta1”, fixed=T)+ f(beta5, SanitasiLayak, copy=“beta1”, fixed=T)+ f(beta6, PHBS, copy=“beta1”, fixed=T)+ f(beta7, persentase_pemberian_asi, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior)
Model.Ridge1s <- inla(Formula.Ridge1s , family=“Gaussian”, data=FullDataLS, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.Ridge1s <- rbind(Model.Ridge1s\(summary.fixed, Model.Ridge1s\)summary.random$beta1[,-1]) round(Coeff.Ridge1s,4)
summary(Model.Ridge1s )
(4) Spatial component—BYM
(a) Poisson
Formula.Ridge1s2 = Case ~ f(beta1, Tingkat.Pendidikan.Ibu, model=“iid”, values = c(1:7), hyper = param.beta)+ f(beta2, PDRB, copy=“beta1”, fixed=T)+ f(beta3, Tingkat.Kemiskinan, copy=“beta1”, fixed=T)+ f(beta4, AirMinumLayak, copy=“beta1”, fixed=T)+ f(beta5, SanitasiLayak, copy=“beta1”, fixed=T)+ f(beta6, PHBS, copy=“beta1”, fixed=T)+ f(beta7, persentase_pemberian_asi, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior)+ f(IS2,model = “iid”, constr = T, hyper = hcprior)
Model.Ridge1s2 <- inla(Formula.Ridge1s2 , family=“Gaussian”, data=FullDataLS, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.Ridge1s2 <- rbind(Model.Ridge1s2\(summary.fixed, Model.Ridge1s2\)summary.random$beta1[,-1]) round(Coeff.Ridge1s2,4)
summary(Model.Ridge1s2)
(5) Spatial + Temporal component—Besag Parametric
Formula.Ridge1nbst = Case ~ f(beta1, Tingkat.Pendidikan.Ibu, model=“iid”, values = c(1:7), hyper = param.beta)+ f(beta2, PDRB, copy=“beta1”, fixed=T)+ f(beta3, Tingkat.Kemiskinan, copy=“beta1”, fixed=T)+ f(beta4, AirMinumLayak, copy=“beta1”, fixed=T)+ f(beta5, SanitasiLayak, copy=“beta1”, fixed=T)+ f(beta6, PHBS, copy=“beta1”, fixed=T)+ f(beta7, persentase_pemberian_asi, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior)+ f(IS3, IT1, model = “iid”, hyper = hcprior, constr = T)+ IT2
Model.Ridge1nbst <- inla(Formula.Ridge1nbst , family=“Gaussian”, data=FullDataLS, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.Ridge1nbst <- rbind(Model.Ridge1nbst\(summary.fixed, Model.Ridge1nbst\)summary.random$beta1[,-1]) round(Coeff.Ridge1nbst,4)
summary(Model.Ridge1nbst)
ST_19=Model.Ridge1nbst\(summary.fitted.values\)mean[1:27] ST_20=Model.Ridge1nbst\(summary.fitted.values\)mean[28:54] ST_21=Model.Ridge1nbst\(summary.fitted.values\)mean[55:81] ST_22=Model.Ridge1nbst\(summary.fitted.values\)mean[82:108] ST_23=Model.Ridge1nbst\(summary.fitted.values\)mean[109:135] ST_24=Model.Ridge1nbst\(summary.fitted.values\)mean[136:162]
boxplot(ST_19,ST_20,ST_21,ST_22, ST_23, ST_24, main=“Boxplot ST”,names=c(“2019”,“2020”,“2021”,“2022”,“2023”,“2024”))
ST=c(ST_19,ST_20,ST_21,ST_22, ST_23, ST_24)
#================ 2. Preprocessing Data Stunting ================#
# Anggap nama dataset Anda adalah 'data_stunting'
# Pastikan nama Kab/Kota seragam (Uppercase) untuk Join
Stunting <- read.csv("DataStuntingJabar5.csv", header = TRUE, sep = ";", dec = ".")
head(Stunting)## No id kode_provinsi nama_provinsi kode_kabupaten_kota Kab.Kota
## 1 1 1 32 JAWA BARAT 3201 Kab Bogor
## 2 2 2 32 JAWA BARAT 3202 Kab Sukabumi
## 3 3 3 32 JAWA BARAT 3203 Kab Cianjur
## 4 4 4 32 JAWA BARAT 3204 Kab Bandung
## 5 5 5 32 JAWA BARAT 3205 Kab Garut
## 6 6 6 32 JAWA BARAT 3206 Kab Tasikmalaya
## persentase_balita_stunting satuan tahun Tingkat.Pendidikan.Ibu PDRB
## 1 4.06 PERSEN 2019 30.04 237113
## 2 8.29 PERSEN 2019 13.48 67406
## 3 6.61 PERSEN 2019 13.25 46807
## 4 7.32 PERSEN 2019 30.23 124001
## 5 4.80 PERSEN 2019 17.79 57579
## 6 15.06 PERSEN 2019 14.93 37219
## Tingkat.Kemiskinan JumlahBalita JumlahKasusStunting Ei SMR
## 1 6.66 442162 132825.46 587303.73 22.61615
## 2 6.22 184865 24919.80 46067.99 54.09353
## 3 9.15 179399 23770.37 42643.80 55.74167
## 4 5.94 240699 72763.31 175140.55 41.54566
## 5 8.98 194656 34629.30 67408.01 51.37268
## 6 9.12 108122 16142.61 17453.72 92.48812
## lnSMR AirMinumLayak SanitasiLayak PHBS persentase_pemberian_asi IT IS
## 1 3.118664 91.02 56.00 53.91 53.12 1 1
## 2 3.990715 79.57 69.00 51.10 57.15 1 2
## 3 4.020728 81.50 62.00 62.00 68.78 1 3
## 4 3.726793 97.84 80.00 57.11 63.84 1 4
## 5 3.939106 84.17 76.09 60.02 74.32 1 5
## 6 4.527080 76.18 81.00 54.01 66.77 1 6
## ST
## 1 6.305517
## 2 9.354103
## 3 8.715372
## 4 8.664596
## 5 10.102140
## 6 11.205440# Join Data Peta dengan Data Stunting
jabar_stunting_join <- jabar_sf %>%
left_join(Stunting, by = c("Kabupaten" = "Kab.Kota"))
## ================= TB ================= ##
## Interval (Menyesuaikan dengan persentase_balita_stunting)
# pretty_breaks diatur untuk skala persentase (misal: 10%, 20%, 30%)
pretty_breaks <- c(-1.25, 1.79, 4.82, 7.85, 10.88,13.90)
# compute labels
labels <- c()
brks <- c(-1.25, 1.79, 4.82, 7.85, 10.88,13.90)
format(brks, scientific=F)## [1] "-1.25" " 1.79" " 4.82" " 7.85" "10.88" "13.90"levels(jabar_stunting_join$brks)## NULL# Definisikan variabel baru pada dataset hasil join
# Catatan: Pastikan kolom di data CSV Anda bernama 'persentase_balita_stunting'
labels <- labels[1:length(labels)-1]
jabar_stunting_join$brks <- cut(jabar_stunting_join$ST,
breaks = brks,
include.lowest = TRUE,
labels = labels)
brks_scale <- levels(jabar_stunting_join$brks)
labels_scale <- rev(brks_scale)
# Simpan hasil
#ggsave("Stunting_Jabar.png", height = 6, width = 8, dpi = 600)
# Plotting menggunakan ggplot
# Menggunakan geom_sf karena jabar_stunting_join adalah objek 'sf'
A <- ggplot(jabar_stunting_join) +
geom_sf(aes(fill = brks), color="black", size=0.1) +
# Menggunakan Tahun sebagai pengganti Triwulan untuk facet
geom_sf_text(aes(label = id),
size = 1.5, # Ukuran teks (setara txt.cex)
color = "black",
fontface = "bold",
check_overlap = TRUE) +
facet_wrap(.~tahun, ncol=3) +
# Skala warna disesuaikan dengan label kategori stunting
scale_fill_manual(
labels = c("<1.78", "[1.79 – 4.81)", "[4.82 – 7.84)", "[7.85 – 10.87)", ">=10.88"),
values = c("grey95", "#FFFF00", "#F97316","#FF0000","#7F1D1D")
) +
theme_bw() +
ylab("Northing") + xlab("Easting") +
theme(axis.text.x = element_text(angle = 90)) +
theme(legend.position = "bottom",legend.title = element_text(size = 8),
legend.text = element_text(size = 8)) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
labs(fill = "ST 2019-2024 (%)") +
theme(legend.position="right", text = element_text(size=25)) +
theme(text = element_text(size=13))
A
A <- ggplot(jabar_stunting_join) +
geom_sf(aes(fill = brks), color="black", size=0.1) +
# Menggunakan Tahun sebagai pengganti Triwulan untuk facet
geom_sf_text(aes(label = id),
size = 1.5, # Ukuran teks (setara txt.cex)
color = "black",
fontface = "bold",
check_overlap = TRUE) +
facet_wrap(.~tahun, ncol=3) +
# Skala warna disesuaikan dengan label kategori stunting
scale_fill_manual(
labels = c("<1.78", "[1.79 – 4.81)", "[4.82 – 7.84)", "[7.85 – 10.87)", ">=10.88"),
values = c("grey95", "#FFFDEB", "#F97316", "#DC2626", "#7F1D1D")
) +
theme_bw() +
ylab("Northing") + xlab("Easting") +
theme(axis.text.x = element_text(angle = 90)) +
theme(legend.position = "bottom",legend.title = element_text(size = 8),
legend.text = element_text(size = 8)) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
labs(fill = "ST 2019-2024 (%)") +
theme(legend.position="right", text = element_text(size=25)) +
theme(text = element_text(size=13))
A
(7) Spatial + Temporal component—Interaction
Formula.Ridge1nbstInt = Case ~ f(beta1, Tingkat.Pendidikan.Ibu, model=“iid”, values = c(1:7), hyper = param.beta)+ f(beta2, PDRB, copy=“beta1”, fixed=T)+ f(beta3, Tingkat.Kemiskinan, copy=“beta1”, fixed=T)+ f(beta4, AirMinumLayak, copy=“beta1”, fixed=T)+ f(beta5, SanitasiLayak, copy=“beta1”, fixed=T)+ f(beta6, PHBS, copy=“beta1”, fixed=T)+ f(beta7, persentase_pemberian_asi, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, scale.model = TRUE, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
Model.Ridge1nbstInt <- inla(Formula.Ridge1nbstInt , family=“Gaussian”, data=FullDataLS, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.Ridge1nbstInt <- rbind(Model.Ridge1nbstInt\(summary.fixed, Model.Ridge1nbstInt\)summary.random$beta1[,-1]) round(Coeff.Ridge1nbstInt,4)
summary(Model.Ridge1nbstInt )
(9) Spatial + Temporal component—Interaction Poisson MEB
Prepare Unic Data
FullDataLS1<-FullDataLS FullDataLS1[c(1:60),c(8:17)]<-0
set.seed(1) MeanTemperature1<-FullDataLS1\(MeanTemperature+runif(360,min=0.000000001, max=0.000000002) set.seed(2) Rainfall1<-FullDataLS1\)Rainfall+runif(360,min=0.000000001, max=0.000000002) set.seed(3) SunshineDuration1<-FullDataLS1\(SunshineDuration+runif(360,min=0.000000001, max=0.000000002) set.seed(4) Evaporation1<-FullDataLS1\)Evaporation+runif(360,min=0.000000001, max=0.000000002) set.seed(5) RelativeHumidity1<-FullDataLS1\(RelativeHumidity+runif(360,min=0.000000001, max=0.000000002) set.seed(7) MeanTemperatureL11<-FullDataLS1\)MeanTemperatureL1+runif(360,min=0.000000001, max=0.000000002) set.seed(8) RainfallL11<-FullDataLS1\(RainfallL1+runif(360,min=0.000000001, max=0.000000002) set.seed(9) SunshineDurationL11<-FullDataLS1\)SunshineDurationL1+runif(360,min=0.000000001, max=0.000000002) set.seed(10) EvaporationL11<-FullDataLS1\(EvaporationL1+runif(360,min=0.000000001, max=0.000000002) set.seed(11) RelativeHumidityL11<-FullDataLS1\)RelativeHumidityL1+runif(360,min=0.000000001, max=0.000000002) set.seed(12) MeanTemperatureL21<-FullDataLS1\(MeanTemperatureL2+runif(360,min=0.000000001, max=0.000000002) set.seed(13) RainfallL21<-FullDataLS1\)RainfallL2+runif(360,min=0.000000001, max=0.000000002) set.seed(14) SunshineDurationL21<-FullDataLS1\(SunshineDurationL2+runif(360,min=0.000000001, max=0.000000002) set.seed(15) EvaporationL21<-FullDataLS1\)EvaporationL2+runif(360,min=0.000000001, max=0.000000002) set.seed(16) RelativeHumidityL21<-FullDataLS1$RelativeHumidityL2+runif(360,min=0.000000001, max=0.000000002)
FullDataLS1<-data.frame(FullDataLS, MeanTemperature1,Rainfall1,SunshineDuration1,Evaporation1,RelativeHumidity1, MeanTemperatureL11,RainfallL11,SunshineDurationL11,EvaporationL11,RelativeHumidityL11, MeanTemperatureL21,RainfallL21,SunshineDurationL21,EvaporationL21,RelativeHumidityL21)
set.seed(17) s<-runif(360, min=0.1, max=1)
write.table(s, row.names=FALSE, “s_final.txt”) write.table(FullDataLS1, row.names=FALSE, “FullDataLS1.txt”)
prior.beta = c(0, 0.0001) prior.prec.u = c(10, 9) prior.prec.x = c(10, 9) prior.prec.y = c(10, 9) prec.u = 1 prec.x=1
prior.beta = c(0, 0.0001) prior.prec.u = c(1, 0.0005) prior.prec.x = c(1, 0.0005) prior.prec.y = c(1, 0.0005) prec.u = 1 prec.x=1
Formula.MEB1= Case ~ f(MeanTemperature,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(Rainfall,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(SunshineDuration,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(Evaporation,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RelativeHumidity,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(MeanTemperatureL1,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RainfallL1,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(SunshineDurationL1,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(EvaporationL1,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RelativeHumidityL1,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(MeanTemperatureL2,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RainfallL2,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(SunshineDurationL2,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(EvaporationL2,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RelativeHumidityL2,model=“meb”, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
Model.MEB1 <- inla(Formula.MEB1, data =FullDataLS1, family = “poisson”, E=E, control.fixed = list(mean.intercept = 0, prec.intercept = 0.00001, mean = 0, prec =0.00001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
summary(Model.MEB1)
Add Tiny values
Formula.MEB2= Case ~ f(MeanTemperature1,model=“meb”,scale=s, values=MeanTemperature1, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(Rainfall1,model=“meb”,scale=s, values= Rainfall1,constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(SunshineDuration1,model=“meb”,scale=s, values=SunshineDuration1, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(Evaporation1,model=“meb”,scale=s, values=Evaporation1, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RelativeHumidity1,model=“meb”,scale=s,values=RelativeHumidity1, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(MeanTemperatureL11,model=“meb”,scale=s, values= MeanTemperatureL11, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RainfallL11,model=“meb”,scale=s, values= RainfallL11, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(SunshineDurationL11,model=“meb”,scale=s, values= SunshineDurationL11, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(EvaporationL11,model=“meb”,scale=s,values=EvaporationL11, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RelativeHumidityL11,model=“meb”,scale=s,values=RelativeHumidityL11,constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(MeanTemperatureL21,model=“meb”,scale=s, values=MeanTemperatureL21 , constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RainfallL21,model=“meb”,scale=s, values= RainfallL21, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(SunshineDurationL21,model=“meb”,scale=s, values= SunshineDurationL21, constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(EvaporationL21,model=“meb”,scale=s,values=EvaporationL21,constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(RelativeHumidityL21,model=“meb”,scale=s, values= RelativeHumidityL21,constr=F, hyper =list(beta = list(param = prior.beta, fixed = FALSE), prec.u =list(param = prior.prec.u, initial = log(prec.u), fixed = FALSE)))+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
Model.MEB2 <- inla(Formula.MEB2, data =FullDataLS1, family = “poisson”, E=E, control.fixed = list(mean.intercept = 0, prec.intercept = 0.00001, mean = 0, prec =0.00001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
summary(Model.MEB2)
(10) Spatial + Temporal - Interaction MEC
prior.beta = c(0, 0.0001) prior.prec.u = c(10, 9) prior.prec.x = c(10, 9) prior.prec.y = c(10, 9) prec.u = 1 prec.x=1
Formula.MEC1 = Case ~ f(MeanTemperature, model=“mec”, constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(Rainfall, model=“mec”, constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDuration, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(Evaporation, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidity, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(MeanTemperatureL1, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RainfallL1, model=“mec”, constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDurationL1, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(EvaporationL1, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidityL1, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(MeanTemperatureL2, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RainfallL2, model=“mec”, constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDurationL2, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(EvaporationL2, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidityL2, model=“mec”,
constr=F,
hyper = list(beta = list(prior = “gaussian”, param = prior.beta, fixed =
FALSE ), prec.u = list( prior = “loggamma”, param = prior.prec.u,
initial = log(prec.u), fixed = FALSE ), prec.x = list( prior =
“loggamma”, param = prior.prec.x, initial = log(prec.x), fixed = FALSE
), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(IS1,model = “besag”, graph = W1, constr
= T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”,
cyclic = T, hyper = igprior5))
Model.MEC1 <- inla(Formula.MEC1 , family=“poisson”, data=FullDataLS1, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.00001, prec.intercept=0.00001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
summary(Model.MEC1)
Add Tiny values
Formula.MEC2 = Case ~ f(MeanTemperature1, model=“mec”,scale=s,
constr=F, values=MeanTemperature1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(Rainfall1, model=“mec”,scale=s,
constr=F, values=Rainfall1, hyper = list(beta = list(prior = “gaussian”,
param = prior.beta, fixed = FALSE ), prec.u = list( prior = “loggamma”,
param = prior.prec.u, initial = log(prec.u), fixed = FALSE ), prec.x =
list( prior = “loggamma”, param = prior.prec.x, initial = log(prec.x),
fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDuration1, model=“mec”,scale=s,
constr=F, values=SunshineDuration1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(Evaporation1, model=“mec”,scale=s,
constr=F, values=Evaporation1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidity1, model=“mec”,scale=s,
constr=F, values=RelativeHumidity1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(MeanTemperatureL11, model=“mec”,scale=s,
constr=F, values=MeanTemperatureL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RainfallL11, model=“mec”,scale=s,
constr=F, values=RainfallL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDurationL11,
model=“mec”,scale=s, constr=F, values=SunshineDurationL11, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(EvaporationL11, model=“mec”,scale=s,
constr=F, values=EvaporationL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidityL11,
model=“mec”,scale=s, constr=F, values=RelativeHumidityL11, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(MeanTemperatureL21, model=“mec”,scale=s,
constr=F, values=MeanTemperatureL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RainfallL21, model=“mec”,scale=s,
constr=F, values=RainfallL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDurationL21,
model=“mec”,scale=s, constr=F, values=SunshineDurationL21, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(EvaporationL21, model=“mec”,scale=s,
constr=F, values=EvaporationL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidityL21,
model=“mec”,scale=s, constr=F, values=RelativeHumidityL21, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(IS1,model = “besag”, graph = W1, constr
= T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”,
cyclic = T, hyper = igprior5))
Model.MEC2 <- inla(Formula.MEC2 , family=“poisson”, data=FullDataLS1, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.00001, prec.intercept=0.00001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
summary(Model.MEC2)
###CLUSTER setwd(“/Users/mindra/MINDRA/PENELITIAN/@DISEASE CLUSTERING/DATA”) source(“AHC_clusteringST_function.R”) ## Spatio-temporal AHC algorithm source(“AHC_function.R”)
Stage 1: applied the AHC algorithm to mortality data (in the same regions and time period)
SMR.prior1<-data.frame(logSMR=DataL$lnSMR) cluster.conf<-
clusteringST.function(data=SMR.prior1, W)
cluster.prior<-cluster.conf
Cluster Configuration Selection
model.selectionR <- function(cluster.prior, data, W, Qs, Qt, Carto=NULL, max.cluster=NULL, final.cluster=NULL, plot.dic=TRUE, strategy=“simplified.laplace”, lincomb=FALSE) {
#### Fit a sequence of models with upto “max.cluster” and choose the best #### by minimising the Deviance Information Criterion.
time <- system.time({
## Define our uniform prior distributions ##
n <- length(unique(data$IS1)) # n <- 30
t <- length(unique(data$IT1)) # t <- 12
t.from <- min(unique(data$IT1)) # t.from Jan-Dec
t.to <- max(unique(data$IT1)) # t.to <- Dec
Rs <- diag(apply(W,2,sum))-W ## Spatial neighborhood matrix
D1 <- diff(diag(t),differences=1) ## RW1 structure matrix
Rt <- t(D1)%*%D1
Qs=Rs
Rt <- t(D1)%*%D1
R.area.Leroux <- diag(dim(Qs)[1])-Qs
if(is.null(final.cluster)){
## Store the DIC values
n <- dim(Qs)[1]
t <- dim(Qt)[1]
if(is.null(max.cluster)) {max.cluster <- n*t-1}
dic.list <- rep(NA, max.cluster)
waic.list <- rep(NA, max.cluster)
## Fit a model with a single cluster
i <- 1
param.beta = list(prec = list(param = c(1.0e-6, 1.0e-6)))
formula = Case ~ f(beta1, MeanTemperature, model="iid", values = c(1:15), hyper = param.beta)+
f(beta2, Rainfall, copy="beta1", fixed=T)+
f(beta3, SunshineDuration, copy="beta1", fixed=T) +
f(beta4, Evaporation, copy="beta1", fixed=T)+
f(beta5, RelativeHumidity, copy="beta1", fixed=T)+
f(beta6, MeanTemperatureL1, copy="beta1", fixed=T)+
f(beta7, RainfallL1, copy="beta1", fixed=T)+
f(beta8, SunshineDurationL1, copy="beta1", fixed=T) +
f(beta9, EvaporationL1, copy="beta1", fixed=T)+
f(beta10, RelativeHumidityL1, copy="beta1", fixed=T)+
f(beta11, MeanTemperatureL2, copy="beta1", fixed=T)+
f(beta12, RainfallL2, copy="beta1", fixed=T)+
f(beta13, SunshineDurationL2, copy="beta1", fixed=T) +
f(beta14, EvaporationL2, copy="beta1", fixed=T)+
f(beta15, RelativeHumidityL2, copy="beta1", fixed=T)+
f(IS1,model = "besag", graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = "rw1", cyclic = T, hyper = igprior5))
model<- inla(formula, family="poisson", data=data, E=E,
control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001),
control.predictor=list(compute=TRUE, cdf=c(log(1))),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE),
control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = "eb", strategy = "simplified.laplace"))
dic.list[i] <- model$dic$dic
waic.list[i] <- model$waic$waic
if(plot.dic==TRUE) {
quartz()
y.lim <- dic.list[i]*c(1/1.01,1.005)
plot(1:max.cluster, dic.list, type="l", xlab="Number of spatial clusters", ylab="DIC",
main=paste("Model",i,"of",max.cluster), ylim=y.lim)
points(1,dic.list[i],pch=19)
mtext(paste("Minimun DIC value:",round(min(dic.list,na.rm=T),2)), side=3, line=-2)
}
## Fit separate models with between 2 and the max number of clusters.
for(i in 2:max.cluster){
j <- n*t-i+1
factor.clust <- as.numeric(as.factor(cluster.prior$cluster.store[j,]))
K <- length(unique(factor.clust))
## M: Cluster correspondence matrix ##
M <- matrix(0,n*t,K)
aux <- as.numeric(factor(factor.clust,labels=1:K))
for(ii in 1:(n*t)){
M[ii,aux[ii]] <- 1
}
ID.clust=as.numeric(factor(factor.clust,labels=1:K))
data$Group<-ID.clust
data$G0<-ID.clust
data$G1<-ID.clust
data$G2<-ID.clust
data$G3<-ID.clust
data$G4<-ID.clust
data$G5<-ID.clust
data$G6<-ID.clust
data$G7<-ID.clust
data$G8<-ID.clust
data$G9<-ID.clust
data$G10<-ID.clust
data$G11<-ID.clust
data$G12<-ID.clust
data$G13<-ID.clust
data$G14<-ID.clust
data$G15<-ID.clust
formula = Case ~ f(G1,MeanTemperature, model="iid", constr=T,hyper=pcprior)+
f(G2,Rainfall, model="iid", constr=T,hyper=pcprior)+
f(G3,SunshineDuration, model="iid", constr=T,hyper=pcprior)+
f(G4,Evaporation, model="iid", constr=T,hyper=pcprior)+
f(G5,RelativeHumidity, model="iid", constr=T,hyper=pcprior)+
f(G6,MeanTemperatureL1, model="iid", constr=T,hyper=pcprior)+
f(G7,RainfallL1, model="iid", constr=T,hyper=pcprior)+
f(G8,SunshineDurationL1, model="iid", constr=T,hyper=pcprior)+
f(G9,EvaporationL1, model="iid", constr=T,hyper=pcprior)+
f(G10,RelativeHumidityL1, model="iid", constr=T,hyper=pcprior)+
f(G11,MeanTemperatureL2, model="iid", constr=T,hyper=pcprior)+
f(G12,RainfallL2, model="iid", constr=T,hyper=pcprior)+
f(G13,SunshineDurationL2, model="iid", constr=T,hyper=pcprior)+
f(G14,EvaporationL2, model="iid", constr=T,hyper=pcprior)+
f(G15,RelativeHumidityL2, model="iid", constr=T,hyper=pcprior)+
f(beta1, MeanTemperature, model="iid", values = c(1:15), hyper = param.beta)+
f(beta2, Rainfall, copy="beta1", fixed=T)+
f(beta3, SunshineDuration, copy="beta1", fixed=T) +
f(beta4, Evaporation, copy="beta1", fixed=T)+
f(beta5, RelativeHumidity, copy="beta1", fixed=T)+
f(beta6, MeanTemperatureL1, copy="beta1", fixed=T)+
f(beta7, RainfallL1, copy="beta1", fixed=T)+
f(beta8, SunshineDurationL1, copy="beta1", fixed=T) +
f(beta9, EvaporationL1, copy="beta1", fixed=T)+
f(beta10, RelativeHumidityL1, copy="beta1", fixed=T)+
f(beta11, MeanTemperatureL2, copy="beta1", fixed=T)+
f(beta12, RainfallL2, copy="beta1", fixed=T)+
f(beta13, SunshineDurationL2, copy="beta1", fixed=T) +
f(beta14, EvaporationL2, copy="beta1", fixed=T)+
f(beta15, RelativeHumidityL2, copy="beta1", fixed=T)+
f(IS1,model = "besag", graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = "rw1", cyclic = T, hyper = igprior5))
model <- inla(formula, family="poisson", data=data, E=E,
control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001),
control.predictor=list(compute=TRUE, cdf=c(log(1))),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE),
control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = "eb", strategy = "simplified.laplace"))
dic.list[i] <- model$dic$dic
waic.list[i] <- model$waic$waic
if(plot.dic==TRUE) {
if(dic.list[i]<y.lim[1]) {
y.lim[1] <- y.lim[1]-abs(y.lim[1]-dic.list[i])
}
if(dic.list[i]>y.lim[2]) {
y.lim[2] <- y.lim[2]+abs(y.lim[2]-dic.list[i])
}
plot(1:max.cluster, dic.list, type="l", xlab="Number of spatial clusters", ylab="DIC",
main=paste("Model",i,"of",max.cluster), ylim=y.lim)
points(i,dic.list[i],pch=19)
mtext(paste("Minimun DIC value:",round(min(dic.list,na.rm=T),2)), side=3, line=-2)
}
}
}
cat("\n\n ************* RUNNING FINAL MODEL ************* \n\n")
###Fit the final model
if(is.null(final.cluster)){
best.model <- which(dic.list==min(dic.list))
}else{
best.model <- final.cluster
plot.dic <- FALSE
dic.list <- NULL
waic.list <- NULL
}
j <- n*t-best.model+1
factor.clust <- as.numeric(as.factor(cluster.prior$cluster.store[j,]))
K <- length(unique(factor.clust))
## M: Cluster correspondence matrix ##
M <- matrix(0,n*t,K)
aux <- as.numeric(factor(factor.clust,labels=1:K))
for(ii in 1:(n*t)){
M[ii,aux[ii]] <- 1
}
ID.clust=as.numeric(factor(factor.clust,labels=1:K))
data$Group<-ID.clust
if(lincomb) {
lc1 = inla.make.lincomb(Predictor = rep(1/(n*t),n*t))
names(lc1) <- "intercept"
M1 <- matrix(-1/(n*t),n,n*t)
for (i in 1:n) {
M1[i,seq(0,n*(t-1),by=n)+i] <- 1/(n*t)*(n-1)
}
lc2 = inla.make.lincombs(Predictor = M1)
names(lc2) <- paste("spatial.",as.character(seq(1:n)),sep="")
M2 <- matrix(-1/(n*t),t,n*t)
for (i in 1:t) {
M2[i,seq(1,n)+n*(i-1)] <- 1/(n*t)*(t-1)
}
lc3 = inla.make.lincombs(Predictor = M2)
names(lc3) <- paste("temporal.",as.character(seq(1:t)),sep="")
M3 <- matrix(1/(n*t),n*t,n*t)
k=1
for (j in 1:t) {
for (i in 1:n) {
M3[k,seq(0,n*(t-1),by=n)+i] <- 1/(n*t)*(1-n)
M3[k,seq(1,n)+n*(j-1)] <- 1/(n*t)*(1-t)
M3[k,k] <- (n*t-n-t+1)/(n*t)
k=k+1
}
}
lc4 = inla.make.lincombs(Predictor = M3)
names(lc4) <- paste("spatio.temporal.",as.character(seq(1:(n*t))),sep="")
all.lc <- c(lc1,lc2,lc3,lc4)
}else{
all.lc <- NULL
}
if(best.model==1){
formula = Case ~ f(beta1, MeanTemperature, model="iid", values = c(1:15), hyper = param.beta)+
f(beta2, Rainfall, copy="beta1", fixed=T)+
f(beta3, SunshineDuration, copy="beta1", fixed=T) +
f(beta4, Evaporation, copy="beta1", fixed=T)+
f(beta5, RelativeHumidity, copy="beta1", fixed=T)+
f(beta6, MeanTemperatureL1, copy="beta1", fixed=T)+
f(beta7, RainfallL1, copy="beta1", fixed=T)+
f(beta8, SunshineDurationL1, copy="beta1", fixed=T) +
f(beta9, EvaporationL1, copy="beta1", fixed=T)+
f(beta10, RelativeHumidityL1, copy="beta1", fixed=T)+
f(beta11, MeanTemperatureL2, copy="beta1", fixed=T)+
f(beta12, RainfallL2, copy="beta1", fixed=T)+
f(beta13, SunshineDurationL2, copy="beta1", fixed=T) +
f(beta14, EvaporationL2, copy="beta1", fixed=T)+
f(beta15, RelativeHumidityL2, copy="beta1", fixed=T)+
f(IS1,model = "besag", graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = "rw1", cyclic = T, hyper = igprior5))
model.final<- inla(formula, family="poisson", data=data, E=E,
control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001),
control.predictor=list(compute=TRUE, cdf=c(log(1))),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE),
control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = "eb", strategy = "simplified.laplace"))
}else{
formula = Case ~ f(G1,MeanTemperature, model="iid", constr=T,hyper=pcprior)+
f(G2,Rainfall, model="iid", constr=T,hyper=pcprior)+
f(G3,SunshineDuration, model="iid", constr=T,hyper=pcprior)+
f(G4,Evaporation, model="iid", constr=T,hyper=pcprior)+
f(G5,RelativeHumidity, model="iid", constr=T,hyper=pcprior)+
f(G6,MeanTemperatureL1, model="iid", constr=T,hyper=pcprior)+
f(G7,RainfallL1, model="iid", constr=T,hyper=pcprior)+
f(G8,SunshineDurationL1, model="iid", constr=T,hyper=pcprior)+
f(G9,EvaporationL1, model="iid", constr=T,hyper=pcprior)+
f(G10,RelativeHumidityL1, model="iid", constr=T,hyper=pcprior)+
f(G11,MeanTemperatureL2, model="iid", constr=T,hyper=pcprior)+
f(G12,RainfallL2, model="iid", constr=T,hyper=pcprior)+
f(G13,SunshineDurationL2, model="iid", constr=T,hyper=pcprior)+
f(G14,EvaporationL2, model="iid", constr=T,hyper=pcprior)+
f(G15,RelativeHumidityL2, model="iid", constr=T,hyper=pcprior)+
f(beta1, MeanTemperature, model="iid", values = c(1:15), hyper = param.beta)+
f(beta2, Rainfall, copy="beta1", fixed=T)+
f(beta3, SunshineDuration, copy="beta1", fixed=T) +
f(beta4, Evaporation, copy="beta1", fixed=T)+
f(beta5, RelativeHumidity, copy="beta1", fixed=T)+
f(beta6, MeanTemperatureL1, copy="beta1", fixed=T)+
f(beta7, RainfallL1, copy="beta1", fixed=T)+
f(beta8, SunshineDurationL1, copy="beta1", fixed=T) +
f(beta9, EvaporationL1, copy="beta1", fixed=T)+
f(beta10, RelativeHumidityL1, copy="beta1", fixed=T)+
f(beta11, MeanTemperatureL2, copy="beta1", fixed=T)+
f(beta12, RainfallL2, copy="beta1", fixed=T)+
f(beta13, SunshineDurationL2, copy="beta1", fixed=T) +
f(beta14, EvaporationL2, copy="beta1", fixed=T)+
f(beta15, RelativeHumidityL2, copy="beta1", fixed=T)+
f(IS1,model = "besag", graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = "rw1", cyclic = T, hyper = igprior5))
model.final <- inla(formula, family="poisson", data=data, E=E,
control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001),
control.predictor=list(compute=TRUE, cdf=c(log(1))),
control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE),
control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = "eb", strategy = "simplified.laplace"))
}
if(plot.dic==TRUE){
plot(1:max.cluster, dic.list, type="l", xlab="Number of spatial clusters", ylab="DIC",
main=paste("Model",best.model,"of",max.cluster), ylim=y.lim)
points(which.min(dic.list),min(dic.list), pch=19, col="red")
lines(rep(which.min(dic.list),2),c(0,min(dic.list)), lty=2, col="red")
axis(1, at=which.min(dic.list), labels=which.min(dic.list), col.axis="red")
mtext(paste("Minimun DIC value:",round(min(dic.list,na.rm=T),2)), side=3, line=-2, col="red")
}})
results <- list(model.final=model.final, factor.clust=factor.clust, dic.list=dic.list, waic.list=waic.list, cpu.time=time[3], M=M) return(results) }
data<-FullDataLS
head(data)
Define our uniform prior distributions
n <- length(unique(data\(IS1)) # n <- 30 t <- length(unique(data\)IT1)) # t <- 12 t.from <- min(unique(data\(IT1)) # t.from Jan-Dec t.to <- max(unique(data\)IT1)) # t.to <- Dec Rs <- diag(apply(W,2,sum))-W ## Spatial neighborhood matrix
D1 <- diff(diag(t),differences=1) ## RW1 structure matrix Rt <- t(D1)%%D1 Qs=Rs Rt <- t(D1)%%D1
STcluster.modelRan1 <- model.selectionR(cluster.prior=cluster.conf, data, W=W1, Qs=Rs, Qt=Rt, Carto=BANDUNG, max.cluster=360, lincomb=TRUE, plot.dic=TRUE, strategy=“simplified.laplace”)
names(STcluster.modelRan1)
dic<-as.matrix(STcluster.modelRan1\(dic.list)
write.table(dic,"dic.txt")
model.final<-STcluster.modelRan1\)model.final
CLUSTER CONFIGURATION G=16
n<-30 m<-16 #minimum cluster for(i in 1:m){ j <- n*t-i+1 factor.clust <- as.numeric(as.factor(cluster.prior$cluster.store[j,])) K <- length(unique(factor.clust)) aux <- as.numeric(factor(factor.clust,labels=1:K)) }
ID.clust=aux data\(G1<-ID.clust data\)G2<-ID.clust data\(G3<-ID.clust data\)G4<-ID.clust data\(G5<-ID.clust data\)G6<-ID.clust data\(G7<-ID.clust data\)G8<-ID.clust data\(G9<-ID.clust data\)G10<-ID.clust data\(G11<-ID.clust data\)G12<-ID.clust data\(G13<-ID.clust data\)G14<-ID.clust data$G15<-ID.clust
formula16 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)+ f(beta1, MeanTemperature, model=“iid”, values = c(1:15), hyper = param.beta)+ f(beta2, Rainfall, copy=“beta1”, fixed=T)+ f(beta3, SunshineDuration, copy=“beta1”, fixed=T) + f(beta4, Evaporation, copy=“beta1”, fixed=T)+ f(beta5, RelativeHumidity, copy=“beta1”, fixed=T)+ f(beta6, MeanTemperatureL1, copy=“beta1”, fixed=T)+ f(beta7, RainfallL1, copy=“beta1”, fixed=T)+ f(beta8, SunshineDurationL1, copy=“beta1”, fixed=T) + f(beta9, EvaporationL1, copy=“beta1”, fixed=T)+ f(beta10, RelativeHumidityL1, copy=“beta1”, fixed=T)+ f(beta11, MeanTemperatureL2, copy=“beta1”, fixed=T)+ f(beta12, RainfallL2, copy=“beta1”, fixed=T)+ f(beta13, SunshineDurationL2, copy=“beta1”, fixed=T) + f(beta14, EvaporationL2, copy=“beta1”, fixed=T)+ f(beta15, RelativeHumidityL2, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model16 <- inla(formula16, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.model16 <- rbind(model16\(summary.fixed, model16\)summary.random$beta1[,-1]) round(Coeff.model16,4)
#######Stepwise ###Model Selection My.step.
DataL1<-DataL[-c(1:60), ]
my.variable.list
<-c(“MeanTemperature”,“Rainfall”,“SunshineDuration”,“Evaporation”,“RelativeHumidity”,
“MeanTemperatureL1”,“RainfallL1”,“SunshineDurationL1”,“EvaporationL1”,“RelativeHumidityL1”,
“MeanTemperatureL2”,“RainfallL2”,“SunshineDurationL2”,“EvaporationL2”,“RelativeHumidityL2”)
Step<-My.stepwise.glm(Y = “Case”, myoffset=“E”, variable.list = my.variable.list, data = DataL1, myfamily = “poisson”)
MEC2 GROUP
n<-30 m<-16 #minimum cluster for(i in 1:m){ j <- n*t-i+1 factor.clust <- as.numeric(as.factor(cluster.prior$cluster.store[j,])) K <- length(unique(factor.clust)) aux <- as.numeric(factor(factor.clust,labels=1:K)) }
ID.clust=aux data\(G1<-ID.clust data\)G2<-ID.clust data\(G3<-ID.clust data\)G4<-ID.clust data\(G5<-ID.clust data\)G6<-ID.clust data\(G7<-ID.clust data\)G8<-ID.clust data\(G9<-ID.clust data\)G10<-ID.clust data\(G11<-ID.clust data\)G12<-ID.clust data\(G13<-ID.clust data\)G14<-ID.clust data$G15<-ID.clust
Formula.MEC2G = Case ~ f(G1,MeanTemperature, model=“iid”,
constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”,
constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”,
constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”,
constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”,
constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”,
constr=T,hyper=pcprior)+ f(MeanTemperature1, model=“mec”,scale=s,
constr=F, values=MeanTemperature1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(Rainfall1, model=“mec”,scale=s,
constr=F, values=Rainfall1, hyper = list(beta = list(prior = “gaussian”,
param = prior.beta, fixed = FALSE ), prec.u = list( prior = “loggamma”,
param = prior.prec.u, initial = log(prec.u), fixed = FALSE ), prec.x =
list( prior = “loggamma”, param = prior.prec.x, initial = log(prec.x),
fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDuration1, model=“mec”,scale=s,
constr=F, values=SunshineDuration1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(Evaporation1, model=“mec”,scale=s,
constr=F, values=Evaporation1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidity1, model=“mec”,scale=s,
constr=F, values=RelativeHumidity1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(MeanTemperatureL11, model=“mec”,scale=s,
constr=F, values=MeanTemperatureL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RainfallL11, model=“mec”,scale=s,
constr=F, values=RainfallL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDurationL11,
model=“mec”,scale=s, constr=F, values=SunshineDurationL11, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(EvaporationL11, model=“mec”,scale=s,
constr=F, values=EvaporationL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidityL11,
model=“mec”,scale=s, constr=F, values=RelativeHumidityL11, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(MeanTemperatureL21, model=“mec”,scale=s,
constr=F, values=MeanTemperatureL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RainfallL21, model=“mec”,scale=s,
constr=F, values=RainfallL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDurationL21,
model=“mec”,scale=s, constr=F, values=SunshineDurationL21, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(EvaporationL21, model=“mec”,scale=s,
constr=F, values=EvaporationL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidityL21,
model=“mec”,scale=s, constr=F, values=RelativeHumidityL21, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(IS1,model = “besag”, graph = W1, constr
= T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”,
cyclic = T, hyper = igprior5))
Model.MEC2G <- inla(Formula.MEC2G , family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.00001, prec.intercept=0.00001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
summary(Model.MEC2G)
CLUSTER CONFIGURATION G=14
n<-30 m<-14 #minimum cluster for(i in 1:m){ j <- n*t-i+1 factor.clust <- as.numeric(as.factor(cluster.prior$cluster.store[j,])) K <- length(unique(factor.clust)) aux <- as.numeric(factor(factor.clust,labels=1:K)) }
ID.clust=aux data\(G1<-ID.clust data\)G2<-ID.clust data\(G3<-ID.clust data\)G4<-ID.clust data\(G5<-ID.clust data\)G6<-ID.clust data\(G7<-ID.clust data\)G8<-ID.clust data\(G9<-ID.clust data\)G10<-ID.clust data\(G11<-ID.clust data\)G12<-ID.clust data\(G13<-ID.clust data\)G14<-ID.clust data$G15<-ID.clust
formula14 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)+ f(beta1, MeanTemperature, model=“iid”, values = c(1:15), hyper = param.beta)+ f(beta2, Rainfall, copy=“beta1”, fixed=T)+ f(beta3, SunshineDuration, copy=“beta1”, fixed=T) + f(beta4, Evaporation, copy=“beta1”, fixed=T)+ f(beta5, RelativeHumidity, copy=“beta1”, fixed=T)+ f(beta6, MeanTemperatureL1, copy=“beta1”, fixed=T)+ f(beta7, RainfallL1, copy=“beta1”, fixed=T)+ f(beta8, SunshineDurationL1, copy=“beta1”, fixed=T) + f(beta9, EvaporationL1, copy=“beta1”, fixed=T)+ f(beta10, RelativeHumidityL1, copy=“beta1”, fixed=T)+ f(beta11, MeanTemperatureL2, copy=“beta1”, fixed=T)+ f(beta12, RainfallL2, copy=“beta1”, fixed=T)+ f(beta13, SunshineDurationL2, copy=“beta1”, fixed=T) + f(beta14, EvaporationL2, copy=“beta1”, fixed=T)+ f(beta15, RelativeHumidityL2, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14 <- inla(formula14, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.model14 <- rbind(model14\(summary.fixed, model14\)summary.random$beta1[,-1]) round(Coeff.model14,4)
Formula.MEC2G14 = Case ~ f(G1,MeanTemperature, model=“iid”,
constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”,
constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”,
constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”,
constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”,
constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”,
constr=T,hyper=pcprior)+ f(MeanTemperature1, model=“mec”,scale=s,
constr=F, values=MeanTemperature1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(Rainfall1, model=“mec”,scale=s,
constr=F, values=Rainfall1, hyper = list(beta = list(prior = “gaussian”,
param = prior.beta, fixed = FALSE ), prec.u = list( prior = “loggamma”,
param = prior.prec.u, initial = log(prec.u), fixed = FALSE ), prec.x =
list( prior = “loggamma”, param = prior.prec.x, initial = log(prec.x),
fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDuration1, model=“mec”,scale=s,
constr=F, values=SunshineDuration1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(Evaporation1, model=“mec”,scale=s,
constr=F, values=Evaporation1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidity1, model=“mec”,scale=s,
constr=F, values=RelativeHumidity1, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(MeanTemperatureL11, model=“mec”,scale=s,
constr=F, values=MeanTemperatureL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RainfallL11, model=“mec”,scale=s,
constr=F, values=RainfallL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDurationL11,
model=“mec”,scale=s, constr=F, values=SunshineDurationL11, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(EvaporationL11, model=“mec”,scale=s,
constr=F, values=EvaporationL11, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidityL11,
model=“mec”,scale=s, constr=F, values=RelativeHumidityL11, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(MeanTemperatureL21, model=“mec”,scale=s,
constr=F, values=MeanTemperatureL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RainfallL21, model=“mec”,scale=s,
constr=F, values=RainfallL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(SunshineDurationL21,
model=“mec”,scale=s, constr=F, values=SunshineDurationL21, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(EvaporationL21, model=“mec”,scale=s,
constr=F, values=EvaporationL21, hyper = list(beta = list(prior =
“gaussian”, param = prior.beta, fixed = FALSE ), prec.u = list( prior =
“loggamma”, param = prior.prec.u, initial = log(prec.u), fixed = FALSE
), prec.x = list( prior = “loggamma”, param = prior.prec.x, initial =
log(prec.x), fixed = FALSE ), mean.x = list( prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(RelativeHumidityL21,
model=“mec”,scale=s, constr=F, values=RelativeHumidityL21, hyper =
list(beta = list(prior = “gaussian”, param = prior.beta, fixed = FALSE
), prec.u = list( prior = “loggamma”, param = prior.prec.u, initial =
log(prec.u), fixed = FALSE ), prec.x = list( prior = “loggamma”, param =
prior.prec.x, initial = log(prec.x), fixed = FALSE ), mean.x = list(
prior = “gaussian”,
initial = 0, fixed=TRUE ) ))+ f(IS1,model = “besag”, graph = W1, constr
= T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”,
cyclic = T, hyper = igprior5))
Model.MEC2G14 <- inla(Formula.MEC2G14 , family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.00001, prec.intercept=0.00001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
summary(Model.MEC2G14)
formula14a = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a <- inla(formula14a, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
#THe best model
formula14a = Case ~ f(G1,MeanTemperature, model=“iid”,
constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”,
constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”,
constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”,
constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”,
constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”,
constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”,
constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”,
constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T,
hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic
= T, hyper = igprior5))
model14a <- inla(formula14a, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a1 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior) model14a1 <- inla(formula14a1, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
#THe best model
formula14a1 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a1 <- inla(formula14a1, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14b1 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)
model14b1 <- inla(formula14b1, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a2 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a2 <- inla(formula14a2, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a3 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a3 <- inla(formula14a3, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a4 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a4 <- inla(formula14a4, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a5 = Case ~ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a5 <- inla(formula14a5, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a6 = Case ~ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a6 <- inla(formula14a6, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a7 = Case ~ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a7 <- inla(formula14a7, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a8 = Case ~ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a8 <- inla(formula14a8, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14a9 = Case ~ f(G3,SunshineDuration, model=“iid”,
constr=T,hyper=pcprior)+
f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+
f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model =
“besag”, graph = W1, constr = T, hyper = hcprior, group = IT2,
control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14a9 <- inla(formula14a9, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14b = Case ~ f(G1,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G2,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G3,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G4,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(beta1, Rainfall, model=“iid”, values = c(1:4), hyper = param.beta)+ f(beta2, SunshineDuration, copy=“beta1”, fixed=T) + f(beta3, RelativeHumidity, copy=“beta1”, fixed=T)+ f(beta4, EvaporationL2, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14b <- inla(formula14b, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.model14b <- rbind(model14b\(summary.fixed, model14b\)summary.random$beta1[,-1]) round(Coeff.model14b,4)
formula14c = Case ~ f(G1,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G2,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G3,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G4,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14c <- inla(formula14c, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.model14c <- rbind(model14c\(summary.fixed, model14c\)summary.random$beta1[,-1]) round(Coeff.model14c,4)
formula14c = Case ~ f(G1,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G2,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G3,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G4,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14c <- inla(formula14c, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula14 = Case ~ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)+ f(beta2, Rainfall, copy=“beta1”, fixed=T)+ f(beta3, SunshineDuration, copy=“beta1”, fixed=T) + f(beta4, Evaporation, copy=“beta1”, fixed=T)+ f(beta5, RelativeHumidity, copy=“beta1”, fixed=T)+ f(beta7, RainfallL1, copy=“beta1”, fixed=T)+ f(beta8, SunshineDurationL1, copy=“beta1”, fixed=T) + f(beta9, EvaporationL1, copy=“beta1”, fixed=T)+ f(beta10, RelativeHumidityL1, copy=“beta1”, fixed=T)+ f(beta12, RainfallL2, copy=“beta1”, fixed=T)+ f(beta13, SunshineDurationL2, copy=“beta1”, fixed=T) + f(beta14, EvaporationL2, copy=“beta1”, fixed=T)+ f(beta15, RelativeHumidityL2, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model14 <- inla(formula14, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
formula15 = Case ~f(beta1, MeanTemperature, model=“iid”, values = c(1:15), hyper = param.beta)+ f(beta2, Rainfall, copy=“beta1”, fixed=T)+ f(beta3, SunshineDuration, copy=“beta1”, fixed=T) + f(beta4, Evaporation, copy=“beta1”, fixed=T)+ f(beta5, RelativeHumidity, copy=“beta1”, fixed=T)+ f(beta6, MeanTemperatureL1, copy=“beta1”, fixed=T)+ f(beta7, RainfallL1, copy=“beta1”, fixed=T)+ f(beta8, SunshineDurationL1, copy=“beta1”, fixed=T) + f(beta9, EvaporationL1, copy=“beta1”, fixed=T)+ f(beta10, RelativeHumidityL1, copy=“beta1”, fixed=T)+ f(beta11, MeanTemperatureL2, copy=“beta1”, fixed=T)+ f(beta12, RainfallL2, copy=“beta1”, fixed=T)+ f(beta13, SunshineDurationL2, copy=“beta1”, fixed=T) + f(beta14, EvaporationL2, copy=“beta1”, fixed=T)+ f(beta15, RelativeHumidityL2, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5))
model15 <- inla(formula15, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.model15 <- rbind(model15\(summary.fixed, model15\)summary.random$beta1[,-1]) round(Coeff.model15,4)
A <- matrix(1, ncol = 30, nrow = 1) e <- matrix(0, ncol = 1)
formula15 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)+ f(beta1, MeanTemperature, model=“iid”, values = c(1:15), hyper = param.beta)+ f(beta2, Rainfall, copy=“beta1”, fixed=T)+ f(beta3, SunshineDuration, copy=“beta1”, fixed=T) + f(beta4, Evaporation, copy=“beta1”, fixed=T)+ f(beta5, RelativeHumidity, copy=“beta1”, fixed=T)+ f(beta6, MeanTemperatureL1, copy=“beta1”, fixed=T)+ f(beta7, RainfallL1, copy=“beta1”, fixed=T)+ f(beta8, SunshineDurationL1, copy=“beta1”, fixed=T) + f(beta9, EvaporationL1, copy=“beta1”, fixed=T)+ f(beta10, RelativeHumidityL1, copy=“beta1”, fixed=T)+ f(beta11, MeanTemperatureL2, copy=“beta1”, fixed=T)+ f(beta12, RainfallL2, copy=“beta1”, fixed=T)+ f(beta13, SunshineDurationL2, copy=“beta1”, fixed=T) + f(beta14, EvaporationL2, copy=“beta1”, fixed=T)+ f(beta15, RelativeHumidityL2, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5),extraconstr = list(A = A, e = e))
model15 <- inla(formula15, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.model15 <- rbind(model15\(summary.fixed, model15\)summary.random$beta1[,-1]) round(Coeff.model15,5)
#####Aditional extraconstr
t<-12 n<-30
A1 <- kronecker(matrix(1,1,t),diag(n))
A2 <- kronecker(diag(t),matrix(1,1,n))
A.constr <- rbind(A1,A2) dim(A.constr)
e=rep(0,n+t)
A1 <- matrix(1, ncol = 30, nrow = 12)
A <- kronecker(Diagonal(12, 1), Matrix(1, ncol = 30, nrow = 1)) e = rep(0, 12)
k<-12 A <- kronecker(Diagonal(k, 1), Matrix(1, ncol = 30, nrow = 1)) A<-as.matrix(A) e <- matrix(0, ncol = k)
dim(A) dim(e)
A <- matrix(1, ncol = 30, nrow = 12) e <- matrix(0, ncol = 12)
formula15.1 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)+ f(beta1, MeanTemperature, model=“iid”, values = c(1:15), hyper = param.beta)+ f(beta2, Rainfall, copy=“beta1”, fixed=T)+ f(beta3, SunshineDuration, copy=“beta1”, fixed=T) + f(beta4, Evaporation, copy=“beta1”, fixed=T)+ f(beta5, RelativeHumidity, copy=“beta1”, fixed=T)+ f(beta6, MeanTemperatureL1, copy=“beta1”, fixed=T)+ f(beta7, RainfallL1, copy=“beta1”, fixed=T)+ f(beta8, SunshineDurationL1, copy=“beta1”, fixed=T) + f(beta9, EvaporationL1, copy=“beta1”, fixed=T)+ f(beta10, RelativeHumidityL1, copy=“beta1”, fixed=T)+ f(beta11, MeanTemperatureL2, copy=“beta1”, fixed=T)+ f(beta12, RainfallL2, copy=“beta1”, fixed=T)+ f(beta13, SunshineDurationL2, copy=“beta1”, fixed=T) + f(beta14, EvaporationL2, copy=“beta1”, fixed=T)+ f(beta15, RelativeHumidityL2, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5),extraconstr = list(A = A , e = e))
model15.1 <- inla(formula15.1, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.model14 <- rbind(model14\(summary.fixed, model14\)summary.random$beta1[,-1]) round(Coeff.model14,5)
pD15.1<-model15.1\(dic\)dic D15.1<-model15.1\(dic\)mean.deviance
SumLogcpoM15.1<-sum(log(model15.1\(cpo\)cpo))
y<-DATACLUST\(Case E<-DATACLUST\)E
RRH15.1<-model15.1\(summary.fitted.values\)mean*E
R15.1<-sum((RRH15.1-mean(y))2)/sum((y-mean(y))2)
Ra15.1<-1-((360-1)/(360-pD15.1))*(1-R15.1)
interaction<-model14$summary.random
A <- matrix(1, ncol = 30, nrow = 60) e <- matrix(0, ncol = 60)
formula15.2 = Case ~ f(G1,MeanTemperature, model=“iid”, constr=T,hyper=pcprior)+ f(G2,Rainfall, model=“iid”, constr=T,hyper=pcprior)+ f(G3,SunshineDuration, model=“iid”, constr=T,hyper=pcprior)+ f(G4,Evaporation, model=“iid”, constr=T,hyper=pcprior)+ f(G5,RelativeHumidity, model=“iid”, constr=T,hyper=pcprior)+ f(G6,MeanTemperatureL1, model=“iid”, constr=T,hyper=pcprior)+ f(G7,RainfallL1, model=“iid”, constr=T,hyper=pcprior)+ f(G8,SunshineDurationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G9,EvaporationL1, model=“iid”, constr=T,hyper=pcprior)+ f(G10,RelativeHumidityL1, model=“iid”, constr=T,hyper=pcprior)+ f(G11,MeanTemperatureL2, model=“iid”, constr=T,hyper=pcprior)+ f(G12,RainfallL2, model=“iid”, constr=T,hyper=pcprior)+ f(G13,SunshineDurationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G14,EvaporationL2, model=“iid”, constr=T,hyper=pcprior)+ f(G15,RelativeHumidityL2, model=“iid”, constr=T,hyper=pcprior)+ f(beta1, MeanTemperature, model=“iid”, values = c(1:15), hyper = param.beta)+ f(beta2, Rainfall, copy=“beta1”, fixed=T)+ f(beta3, SunshineDuration, copy=“beta1”, fixed=T) + f(beta4, Evaporation, copy=“beta1”, fixed=T)+ f(beta5, RelativeHumidity, copy=“beta1”, fixed=T)+ f(beta6, MeanTemperatureL1, copy=“beta1”, fixed=T)+ f(beta7, RainfallL1, copy=“beta1”, fixed=T)+ f(beta8, SunshineDurationL1, copy=“beta1”, fixed=T) + f(beta9, EvaporationL1, copy=“beta1”, fixed=T)+ f(beta10, RelativeHumidityL1, copy=“beta1”, fixed=T)+ f(beta11, MeanTemperatureL2, copy=“beta1”, fixed=T)+ f(beta12, RainfallL2, copy=“beta1”, fixed=T)+ f(beta13, SunshineDurationL2, copy=“beta1”, fixed=T) + f(beta14, EvaporationL2, copy=“beta1”, fixed=T)+ f(beta15, RelativeHumidityL2, copy=“beta1”, fixed=T)+ f(IS1,model = “besag”, graph = W1, constr = T, hyper = hcprior, group = IT2, control.group = list(model = “rw1”, cyclic = T, hyper = igprior5),extraconstr = list(A = A , e = e))
model15.2 <- inla(formula15.2, family=“poisson”, data=data, E=E, control.fixed=list(mean=0, mean.intercept=0, prec=0.0001, prec.intercept=0.0001), control.predictor=list(compute=TRUE, cdf=c(log(1))), control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE), control.inla = list(tolerance=1e-20, h=1e-8,int.strategy = “eb”, strategy = “simplified.laplace”))
Coeff.model15.2 <- rbind(model15.2\(summary.fixed, model15.2\)summary.random$beta1[,-1]) round(Coeff.model15.2,5)
pD15.2<-model15.2\(dic\)dic D15.2<-model15.2\(dic\)mean.deviance
SumLogcpoM15.2<-sum(log(model15.2\(cpo\)cpo))
y<-DATACLUST\(Case E<-DATACLUST\)E
RRH15.2<-model15.2\(summary.fitted.values\)mean*E
R15.2<-sum((RRH15.2-mean(y))2)/sum((y-mean(y))2)
Ra15.2<-1-((360-1)/(360-pD15.2))*(1-R15.2)